* This is idential to MC3 except the Monte Carlo was changed to count events per year. * Note that this system has a well defined peak each day which makes the LOLE ~= LOLEV. * * 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 LOLEV=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 LOLEV ALOLP MWHEUE' WRITE(*,*) ' MONTECARLO SEQUENTIAL SOLUTION' WRITE(*,*) ' YRS LOLH LOLEV ALOLP MWHEUE' * BEGIN THE LOOP SIMULATING EACH YEAR WITH SEQUENTIAL HOURS WITHIN THAT YEAR 6 YRS=YRS+1.D0 NEV=0 ! COUNT THE NUMBER OF EVENTS PER YEAR IEVPH=0 ! 0 if previous hour was not an event or 1 if it was an event DO 8 J=1,365 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 IEVPH=0 ! NO LOSS OF LOAD AND ANY PREVIOUS EVENTS HAVE ENDED 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(IEVPH.EQ.0) THEN NEV=NEV+1 ! THIS IS THE FIRST HOUR OF AN EVENT, ADD ONE TO THE NUMBER OF EVENTS COUNTING IEVPH=1 ! A LOSS OF LOAD EVENT HAS STARTED, SO TURN ON THE FLAG SAYING IT IS IN PROGRESS ENDIF ENDIF 8 CONTINUE LOLEV=LOLEV+NEV ! SUM UP THE EVENTS FOR EACH YEAR AND LOLEV STORES THE RUNNING SUM IF(NEV.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,LOLEV/YRS,1.-ALLR/YRS,EUE/YRS WRITE(*,'(F8.0,3F11.6,F10.3)') & YRS,LOLH/YRS,LOLEV/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 COPT SOLUTION COPT = capacity outage probability table LOLH LOLE ALLR 0.225305 0.099956 0.095577 MWHEUE MWHTOTAL puEUEppm 33.39505 28573440. 1.16874 RESERVE MARGIN = 13.57% TIME = 0.047 SEC Frequency and Duration Monte Carlo Solution YRS LOLH LOLEV ALOLP MWHEUE 1000. 0.196000 0.097000 0.021000 20.065 TIME = 17.922 SEC 10000. 0.232100 0.104200 0.022300 34.014 TIME = 179.391 SEC 100000. 0.234020 0.104210 0.022790 35.379 TIME = 1792.844 SEC 1000000. 0.225614 0.100805 0.022271 33.696 TIME = 17962.500 SEC <- 5 hrs 380,000 times slower than direct! The F&D MC LOLEV is identical to the F&D MC LOLE. In LOLE we count each occurrence per day once. In LOLEV we count each occurrence without worrying about whether the event aligns with the day or not. This test shows that they are identical for a daily summer peaking demand profile and lower demands for the rest of the year. LOLEV would be higher than LOLE if there were two peak demands during the peak demand days. LOLEV might be lower than LOLE if loss of load extends from one day to the next day. LOLEV might depart from LOLE if repair times to the generators is less than the duration of load shedding events. 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