* The MC3 program was written by Gene Preston 08/15/2015 as an expansion of MC2 by duplicating the generator data * eleven times. A clock routine is added. Please distribute and use as you wish. In using this program you * agree that Gene Preston is not responsible for any consequences that may arise from the use of this program. * This program calculates the annual LOLE and ALLR using the direct approach. * The sequential Monte Carlo is also added to calculate the annual LOLP or ALOLP. * A conversational feature is added in which D1 D2 D3 D4 can be adjusted to achieve a target index. * * BEGIN THE DIRECT SOLUTION CHARACTER*50 A/' '/ REAL*8 F(0:10703)/1.D0,10703*0.D0/,LOLE,ALLR,EUE,ENGY,LOLH REAL FOR(10)/.15,.10,.09,.07,.06,.05,.04,.04,.04,.04/ INTEGER PMAX(10)/5,10,13,22,43,57,101,198,248,276/ REAL HLOAD(24,365),DAILY(24) &/.56,.60,.64,.68,.72,.76,.80,.84,.88,.92,.96,1.0, & 1.0,.96,.92,.88,.84,.80,.76,.72,.68,.64,.60,.56/ * ADDITIONAL INFORMATION FOR THE SEQUENTIAL MONTECARLO SOLUTION, 110 INDIVIDUAL GENERATORS ARE MODELED REAL MTTR(10)/5.,5.,5.,7.,7.,14.,21.,28.,28.,42./ ! MEAN TIME TO REPAIR IN DAYS REAL RR(110),FR(110) ! PROBABILITY OF HOURLY REPAIR AND FAILURE RATES INTEGER ISTATE(110) ! 0, GEN IS OFF 1, GEN IS ON REAL*8 YRS ! COUNT THE NUMBER OF YEARS SIMULATED OPEN(3,FILE='OP') * CALCULATE F(X) CUMULATIVE CAPACITY DISTRIBUTION CURVE USING THE DIRECT METHOD T1=0. ! T1 is the end of the first run time, also used to detect the first pass CALL CPU_TIME(T0) ! this is Compaq Fortran at t=0, T1-T0 is the run time in seconds DO 1 K=1,11 ! repeat 11 times DO 1 I=1,10 ! ten original generators DO 1 IX=K*973,1,-1 ! total capacity of all 110 generators, do it in batches of 10 J=IX-PMAX(I) IF(J.LT.0) J=0 1 F(IX)=(1.-FOR(I))*F(J)+FOR(I)*F(IX) * INITIALIZE CONSTANTS AND VARIABLES, D1-D4 SELECTED TO MAKE LOLE = 0.1 d/y D1=2000. ! 341 MIN LD DAYS PEAK DEMANDS D2=8000. ! 10 WINTER DAYS PEAK DEMANDS D3=7000. ! 14 SUMMER DAYS PEAK DEMANDS D4=1424. ! CONSTANT MW ADDER EVERY HOUR 9 LOLE=0.D0 ! LOSS OF LOAD EXPECTATION LOLH=0.D0 ! LOSS OF LOAD HOURS ALLR=1.D0 ! ANNUAL LOAD LOSS RISK EUE=0.D0 ! EXPECTED UNSERVED ENERGY ENGY=0.D0 ! TOTAL LOAD ENERGY PEAKD=0. ! REMEMBER THE PEAK DEMAND * SWEEP THROUGH THE YEAR CALCULATING LOADS AND INDICES DO 3 J=1,365 DO 3 I=1,24 HLOAD(I,J)=IFIX(D1*DAILY(I)+.01) ! LIGHT LOAD PERIOD IF(J.LE.10) HLOAD(I,J)=IFIX(D2*DAILY(I)+.01) ! 10 WINTER PEAK DEMANDS IF(J.GE.201.AND.J.LE.214) HLOAD(I,J)=IFIX(D3*DAILY(I)+.01) ! 14 SUMMER PEAK DEMANDS HLOAD(I,J)=HLOAD(I,J)+D4 ! ADD THE CONSTANT LOAD IF(HLOAD(I,J).GT.PEAKD) PEAKD=HLOAD(I,J) ! REMEMBER THE PEAK DEMAND IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR LOLE=LOLE+(1.-F(HLOAD(I,J))) ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*F(HLOAD(I,J)) ! PROBABILITY OF DAILY GEN UP STATES ENDIF DO 2 K=1,HLOAD(I,J) 2 EUE=EUE+(1.-F(K)) ! EXPECTED UNSERVED ENERGY - THE AREA ABOVE F(X) AND < PR=1 LOLH=LOLH+(1.-F(HLOAD(I,J))) ! LOLH = RUNNING SUM OF THE HOURLY LOLP 3 ENGY=ENGY+HLOAD(I,J) ALLR=1.D0-ALLR * DISPLAY RESULTS WRITE(3,*) ' DIRECT SOLUTION' WRITE(3,*) ' ' WRITE(3,*) ' LOLH LOLE ALLR' WRITE(3,'(3F11.6)') LOLH,LOLE,ALLR WRITE(3,*) ' ' WRITE(3,*) ' MWHEUE MWHTOTAL puEUEppm' WRITE(3,'(F11.5,F12.0,F10.5)') EUE,ENGY,EUE/ENGY*1.D6 WRITE(3,'(/19H RESERVE MARGIN =,F8.2,1H%)') & (10703.-PEAKD)/PEAKD*100. * SHOW THE DIRECT SOLUTION LOLP FOR WINTER, SUMMER, AND OFF PEAK HOURS WRITE(3,*) ' ' WRITE(3,*)' WINTER SUMMER OFF PEAK' WRITE(3,*)'HR MW LOLP MW LOLP MW LOLP' DO 4 I=1,24 4 WRITE(3,'(I3,3(F7.0,F12.8))') I, & HLOAD(I, 1),1.-F(HLOAD(I, 1)), & HLOAD(I,201),1.-F(HLOAD(I,201)), & HLOAD(I, 11),1.-F(HLOAD(I, 11)) WRITE(3,'(/''D1 MIN, D2 WIN, D3 SUM, D4 ADD ='',4I7)') & IFIX(D1),IFIX(D2),IFIX(D3),IFIX(D4) IF(T1.LE.0.) THEN CALL CPU_TIME(T1) WRITE(3,'(''TIME ='',F10.3,'' SEC'')') T1-T0 ENDIF WRITE(*,*) ' LOLH LOLE puEUEppm' WRITE(*,'(2F11.6,F11.3)') LOLH,LOLE,EUE/ENGY*1.D6 WRITE(*,'(''D1 MIN, D2 WIN, D3 SUM, D4 ADD ='',4I10)') & IFIX(D1),IFIX(D2),IFIX(D3),IFIX(D4) WRITE(*,*) 'ENTER NEW D1,D2,D3,D4 OR RETURN TO CONTINUE' READ(*,'(A50)') A IF(A.NE.' ') THEN READ(A,*) D1,D2,D3,D4 GOTO 9 ENDIF * * BEGIN THE SEQUENTIAL MONTECARLO SIMULATION, WILL BE REPORTED IN INTERVALS OF 1000, 10K, 100K, ETC YEARS T0=0. CALL CPU_TIME(T0) DO 5 K=1,11 ! 11 SETS OF GENERATORS K1=(K-1)*10 DO 5 I=1,10 ! 10 GENERATORS PER SET ISTATE(K1+I)=1 ! START OUT WITH ALL GENERATORS TURNED ON IF(K.EQ.1) MTTR(I)=MTTR(I)*24. ! CONVERT MEAN TIME TO REPAIR FROM DAYS TO HOURS RR(K1+I)=1./MTTR(I) ! CALCULATE THE REPAIR RATE TRANSITION PROBABILITY FOR EACH HOUR 5 FR(K1+I)=RR(K1+I)*FOR(I)/(1.-FOR(I)) ! CALCULATE THE FAILURE RATE FROM THE REPAIR RATE AND FORCED OUTAGE RATE LOLE=0.D0 ! NUMBER OF LOSS OF LOAD EXPECTATION EVENTS, NO MORE THAN 1 PER DAY WILL BE COUNTED. * NOTE THAT UNITS OF DAYS PER YEAR MEANS IF A DAY HAD AN OUTAGE IT HAS AN OUTAGE. * HAVING TWO OUTAGES IN THE SAME DAY DOESN'T CHANGE THE FACT THAT IT HAD AN OUTAGE * IN THAT DAY. THIS IS TO MAINTAIN CONSISTENCY WITH THE DEFINITION OF DAYS PER YEAR. LOLH=0.D0 ! LOSS OF LOAD HOURS, ADD 1 TO EVERY HOUR THERE IS AN OBSERVED LOSS OF LOAD ALLR=0.D0 ! ANNUAL LOAD LOSS RISK, WILL COUNT THE NUMBER OF YEARS THERE ARE NO LOSS OF LOAD EVENTS EUE=0.D0 ! EXPECTED UNSERVED ENERGY, WILL SUM THE DIFFERENCE BETWEEN THE LOAD AND GENERATION AVAILABLE EACH HOUR M=3 ! SEED FOR RANDOM NUMBER GENERATOR MULT=1000. ! START OUT WITH 1000 YEARS OF SIMULATION BEFORE REPORTING RESULTS YRS=0.D0 WRITE(3,*) ' ' WRITE(3,*) ' ' WRITE(3,*) ' MONTECARLO SEQUENTIAL SOLUTION' WRITE(3,*) ' ' WRITE(3,*) ' YRS LOLH LOLE ALOLP MWHEUE' WRITE(*,*) ' MONTECARLO SEQUENTIAL SOLUTION' WRITE(*,*) ' YRS LOLH LOLE ALOLP MWHEUE' * BEGIN THE LOOP SIMULATING EACH YEAR WITH SEQUENTIAL HOURS WITHIN THAT YEAR 6 YRS=YRS+1.D0 NED=0 ! COUNT THE NUMBER OF EVENT DAYS JUST FOR THIS ONE YEAR DO 8 J=1,365 NEED=0 ! COUNT THE NUMBER OF EVENTS EACH DAY DO 8 I=1,24 GEN=0. ! THE MW GENERATION ON LINE THIS HOUR DO 7 K=1,110 ! SWEEP THROUGH ALL THE GENERATORS EACH HOUR CHANGING STATES AND SUMMING THE TOTAL GENERATION AVAILABLE R=RAN(M) IF(ISTATE(K).EQ.1) THEN IF(R.LE.FR(K)) ISTATE(K)=0 ! FLIP THE GENERATOR TO A FAILURE STATE ELSE IF(R.LE.RR(K)) ISTATE(K)=1 ! FLIP THE GENERATOR BACK ON LINE, ITS FIXED ENDIF K1=K-((K-1)/10)*10 ! THIS FINDS THE CORRECT GENERATOR PMAX IF(ISTATE(K).EQ.1) GEN=GEN+PMAX(K1)! SUM UP THE POWER OF THE GENS ON LINE 7 CONTINUE IF(GEN+.0001.GE.HLOAD(I,J)) THEN ! THERE IS ENOUGH GENERATION, THE .0001 INSURES THE IF TEST IS ACCURATE IN CASE REAL NUMBERS HAPPEN TO END IN .999 * INSERT SOMETHING HERE FOR THE OK CONDITION IF THERE IS ANYTHING TO CALCULATE OR REMEMBER ELSE EUE=EUE+HLOAD(I,J)-GEN ! HLOAD EXCEEDS THE GENERATION MW, ADD THE MW TIMES ONE HOUR, OR UNSERVED ENERGY IN MWH LOLH=LOLH+1. ! COUNT THE HOURS THERE IS LOSS OF LOAD IF(NEED.EQ.0) THEN NED=NED+1 ! THIS IS THE FIRST OCCURENCE OF LOSS OF LOAD THIS DAY NEED=1 ! DON'T ALLOW ANY MORE TO BE COUNTED THIS DAY ENDIF ENDIF 8 CONTINUE LOLE=LOLE+NED ! SUM UP THE EVENTS FOR EACH YEAR AND LOLE STORES THE RUNNING SUM IF(NED.EQ.0) ALLR=ALLR+1.D0 ! COUNT THE NUMBER OF YEARS THERE WERE NO EVENTS IF(YRS.LT.MULT) GOTO 6 * DISPLAY RESULTS FOR CHOSEN NUMBER OF YEARS 1K 10K 100K ETC WRITE(3,'(F8.0,3F11.6,F10.3)') & YRS,LOLH/YRS,LOLE/YRS,1.-ALLR/YRS,EUE/YRS WRITE(*,'(F8.0,3F11.6,F10.3)') & YRS,LOLH/YRS,LOLE/YRS,1.-ALLR/YRS,EUE/YRS CALL CPU_TIME(T1) WRITE(3,'(''TIME ='',F10.3,'' SEC'')') T1-T0 WRITE(*,'(''TIME ='',F10.3,'' SEC'')') T1-T0 MULT=MULT*10. WRITE(*,*) 'CONTINUE MONTECARLO? 1=YES' READ(*,'(A50)') A IF(A(1:1).EQ.'1') GOTO 6 WRITE(*,*) 'FILE OP CONTAINS THE OUTPUT REPORT' WRITE(*,*) 'PRESS ENTER TO END THIS PROGRAM' READ(*,'(A1)') A END DIRECT SOLUTION LOLH LOLE ALLR 0.225305 0.099956 0.095577 MWHEUE MWHTOTAL puEUEppm 33.39505 28573440. 1.16874 RESERVE MARGIN = 13.57% WINTER SUMMER OFF PEAK HR MW LOLP MW LOLP MW LOLP 1 5904. 0.00000000 5344. 0.00000000 2544. 0.00000000 2 6224. 0.00000000 5624. 0.00000000 2624. 0.00000000 3 6544. 0.00000000 5904. 0.00000000 2704. 0.00000000 4 6864. 0.00000000 6184. 0.00000000 2784. 0.00000000 5 7184. 0.00000000 6464. 0.00000000 2864. 0.00000000 6 7504. 0.00000000 6744. 0.00000000 2944. 0.00000000 7 7824. 0.00000001 7024. 0.00000000 3024. 0.00000000 8 8144. 0.00000035 7304. 0.00000000 3104. 0.00000000 9 8464. 0.00000663 7584. 0.00000000 3184. 0.00000000 10 8784. 0.00009990 7864. 0.00000002 3264. 0.00000000 11 9104. 0.00116223 8144. 0.00000035 3344. 0.00000000 12 9424. 0.00998912 8424. 0.00000465 3424. 0.00000000 13 9424. 0.00998912 8424. 0.00000465 3424. 0.00000000 14 9104. 0.00116223 8144. 0.00000035 3344. 0.00000000 15 8784. 0.00009990 7864. 0.00000002 3264. 0.00000000 16 8464. 0.00000663 7584. 0.00000000 3184. 0.00000000 17 8144. 0.00000035 7304. 0.00000000 3104. 0.00000000 18 7824. 0.00000001 7024. 0.00000000 3024. 0.00000000 19 7504. 0.00000000 6744. 0.00000000 2944. 0.00000000 20 7184. 0.00000000 6464. 0.00000000 2864. 0.00000000 21 6864. 0.00000000 6184. 0.00000000 2784. 0.00000000 22 6544. 0.00000000 5904. 0.00000000 2704. 0.00000000 23 6224. 0.00000000 5624. 0.00000000 2624. 0.00000000 24 5904. 0.00000000 5344. 0.00000000 2544. 0.00000000 D1 MIN, D2 WIN, D3 SUM, D4 ADD = 2000 8000 7000 1424 TIME = 0.047 SEC MONTECARLO SEQUENTIAL SOLUTION YRS LOLH LOLE ALOLP MWHEUE 1000. 0.196000 0.097000 0.021000 20.065 TIME = 17.188 SEC 10000. 0.232100 0.104200 0.022300 34.014 TIME = 171.781 SEC 100000. 0.234020 0.104210 0.022790 35.379 TIME = 1715.734 SEC