* Eugene G Preston is not responsible for any errors MC4i may contain. This program is a variation off MC4 * in that the convolution of the down states shifts to the right. Good precision is maintained with real*4 numbers. * These changes allow extremely large networks such as the Eastern grid to be easily modeled with small arrays. * Standard deviation percentage has been added to the MC section. MC4i update was posted 12/30/2015. updated 02/23/2016. * * BEGIN THE DIRECT SOLUTION - based on DCPk.for convolution method PARAMETER MF=3000 ! MAX SIZE OF ARRAY F(X), INCREASE AS NEEDED, MAYBE MF=80000 FOR THE EASTERN GRID CHARACTER*100 A/' '/ ! USED TO READ IN A STRING OF NEW NUMBERS OR TEST TO SEE IF BLANK REAL F(-1:MF)/1.,0.,MF*0./ ! F(X) CUMULATIVE DISTRIBUTION FOR GENERATORS FOR THE TOTAL SYSTEM REAL FO(-1:MF)/1.,0.,MF*0./ ! FO(X) CUMULATIVE DISTRIBUTION FOR GENERATORS FOR THE OUTSIDE AREA 0 REAL F1(-1:MF)/1.,0.,MF*0./ ! F1(X) CUMULATIVE DISTRIBUTION FOR GENERATORS FOR THE INSIDE AREA 1 REAL FT(-1:MF) ! TEMPORARY F1 IS RECREATED EACH HOUR AS F1+ATCS AS AN OUTSIDE GENERATOR FEEDING AREA 1 REAL*8 LOLE,ALLR,ALOLP,EUE,ENGY,LOLH,D8 ! AREA 1 STATS FOR THE DIRECT SOLUTION SUMMED IN DOUBLE PRECISION REAL FOR(10,2)/.15,.10,.09,.07,.06,.05,.04,.04,.04,.04,10*.0/ ! no derated states, replace this with real data INTEGER PMAX(10,2)/5,10,13,22,43,57,101,198,248,276,10*0/ ! no derated states, replace this with real data ! note that PMAX( ,2) is the DERATING from PMAX( ,1), not the rating at a lower level than PMAX( ,1) INTEGER PMX/0/,PMX1/0/,PMXO/0/,PMXT ! PMX PMX1 PMXO ARE TOTAL GENERATION MW FOR THE SYSTEM, AREA 1, AND AREA 0 USED BY EACH DISTRIBUTION INTEGER PMN,PMN1,PMNO ! PMX MW CORRESPONDS TO F(I=1)=LOLP@PMX F1(I=1)=LOLP@PMX1 AND FO(I=1)=LOLP@PMXO REAL HL(24,365),HL1(24,365),DAILY(24) ! HL1 ARE AREA 1 HOURLY LOADS AND HL IS THE TOTAL SYSTEM, REPLACE THIS WITH ACTUAL SYSTEM HOURLY LOADS &/.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/ * ADD THE ATCS OF THE TIE LINES INTO AREA 1, DEFAULT INITIALLY TO 1 ATC LIMIT, BUT THIS CAN BE CHANGED CONVERSATIONALLY AFTER THE FIRST PASS. REAL ATC(2,12)/300.,.99,200.,.01,-300.,1.,18*0./ ! The ATC is +300 MW import power with probability .99, 200 MW probability = .01 ! (2,11) is the probability of the 0 MW +0 MW import state and (2,12) is the probability of the -0 MW export export ! Replace this ATC data with actual study ATC results in a real two area study. REAL ATCP(12) ! This stores the hourly product of ATC probabilty times availability of generation for import/export REAL ATCX(10) ! This is the ATCP MW that is convolved. Too little sending power may reduce this MW below the specified ATC ! note that ATCP(11)&(12) are slack import&export probabilities for 0 MW transfer in both directions * 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, replace this with real data if you have it REAL RR(110),FR(110),GMTTR,AMTTR/0./ ! PROBABILITY OF HOURLY RATE OF REPAIR AND HOURLY PROBABILITY OF FAILURE FOR MONTE CARLO INTEGER ISTATE(110) ! 0=GEN IS OFF 1=GEN IS ON at Pmax( ,1) 2=Gen is on at derated power Pmax( ,2) INTEGER YRS ! COUNT THE NUMBER OF YEARS SIMULATED REAL*8 ATCH(12)/12*0.D0/ ! COUNT THE NUMBER OF ATC HOURS IN MONTE CARLO FOR EACH LINE STATE REAL*8 ATCC(12)/12*0.D0/ ! SAME AS ATCH BUT ONLY COUNTED WHEN IN CONJUNCTION WITH A LOSS OF LOAD REAL*8 LOLE1,ALLR1,ALOLP1,EUE1,ENGY1,LOLH1! AREA 1 STATS FOR MONTE CARLO SUMMED IN DOUBLE PRECISION REAL*8 V(12),S(12),B(12) ! THESE STORE INFO TO CALCULATE % STANDARD DEVIATIONS INTEGER SD(12) OPEN(3,FILE='OP') * CALCULATE F(X) CUMULATIVE CAPACITY DISTRIBUTION CURVE USING THE DIRECT METHOD 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 PMX=PMX+PMAX(I,1) ! TOTAL INSTALLED GENERATION CAPACITY PMX = PMAX SUM IMX=PMX ! IMX WILL BE SHIFTED DOWN LATER AMTTR=AMTTR+PMAX(I,1)*MTTR(I) IF(IMX.GT.MF) IMX=MF IF(K.EQ.1) THEN PMX1=PMX1+PMAX(I,1)! AREA 1 INSTALLED GENERATION CAPACITY IMX1=PMX1 IF(IMX1.GT.MF) IMX1=MF IF(PMAX(I,1).LE.0.OR.PMAX(I,1).LT.PMAX(I,2)) THEN WRITE(*,*) 'gen',I,' PMAX is too small' GOTO 41 ENDIF ENDIF IF(K.LE.10) THEN PMXO=PMXO+PMAX(I,1)! OUTSIDE AREA 1 INSTALLED GENERATION CAPACITY IMXO=PMXO IF(IMXO.GT.MF) IMXO=MF ENDIF DO 1 JX=IMX,0,-1 ! this convolves together all 110 three state generators J1=JX-PMAX(I,1) ! shift the fully outage FOR state to the right J1 amount J2=JX-PMAX(I,2) ! shift the derated power state to the right J2 amount IF(J1.LT.-1) J1=-1 ! this effectively captures F(J1=X<0)=1 for Pr of J1 shift IF(J2.LT.-1) J2=-1 ! this effectively captures F(J2=X<0)=1 for Pr of J2 shift D4=1.-FOR(I,1)-FOR(I,2) ! the probability of full power at Pmax FJ1=F(J1) ! this is a precaution if J1=JX FJ2=F(J2) ! this is a precaution if J2=JX F(JX)= D4*F(JX) ! the full power state is not shifted to the right F(JX)=F(JX)+FOR(I,1)*FJ1 ! the full outage state is shifted to the right J1 F(JX)=F(JX)+FOR(I,2)*FJ2 ! the partial power state is shifted to the right J2 IF(K.EQ.1) F1(JX)=F(JX) ! AREA 1 GENERATORS, note in an actual system we convolve the F F1 FO separately IF(K.EQ.10) FO(JX)=F(JX) ! total system excluding area 1 generators (technically should be the last 10, but first 10 is the same set in this example) 1 CONTINUE IMTTR=(24.*AMTTR/PMX)+.5 AMTTR=IMTTR/24. ! AVERAGE (WEIGHTED BY PMAX) MTTR OF THE POWER PLANTS, YOU COULD ENCODE A NEW AMTTR HERE * shrink IMX to the range of significance for F(X) LOLP DO 33 JX=0,MF D4=1.-F(JX) IF(D4.LT..99999996) THEN ! test for significance IMX=JX ! IMX is as far out on the curve as we need to go ELSE DO 32 K=JX,MF ! the LOLP beyond IMX is 0 32 F(K)=0. GOTO 34 ENDIF 33 CONTINUE 34 PMN=PMX-IMX IF(IMX.EQ.MF) THEN ! cumulative distribution is too wide, increase the MF value WRITE(*,'(13Hwarning: 1-F(,I5,3H) =,F12.8)') MF,F(MF) ENDIF * J=IMX+1 ! this is to make sure the print statement doesn't run greater than MF * IF(J.GT.MF) J=MF ! writes MC4iF(X) to compare with MC4F(X) double precision & not inverted F(x) * WRITE(3,*) 'PMX=',PMX,' PMN=',PMX-IMX,' IMX=',IMX * WRITE(3,*) ' I MW 1 - F(X) F(X) = LOLP' * WRITE(3,'(I6,I8,2F12.8)') (J-I,PMX-I,1.D0-F(I),F(I),I=J,-1,-1) * shrink IMX1 to the range of significance for F1(X) LOLP DO 36 JX=0,MF D4=1.-F1(JX) IF(D4.LT..99999996) THEN ! test for significance IMX1=JX ! IMX1 is as far out on the curve as we need to go ELSE DO 35 K=JX,MF ! the LOLP beyond IMX1 is 0 35 F1(K)=0. GOTO 37 ENDIF 36 CONTINUE 37 PMN1=PMX1-IMX1 IF(IMX1.EQ.MF) THEN ! cumulative distribution is too wide, increase the MF value WRITE(*,'(13Hwarning: 1-F(,I5,3H) =,F12.8)') MF,F1(MF) ENDIF * J=IMX1+1 ! this is to make sure the print statement doesn't run greater than MF * IF(J.GT.MF) J=MF * WRITE(3,*) 'PMX1=',PMX1,' PMN1=',PMX1-IMX1,' IMX1=',IMX1 * WRITE(3,*) ' I MW 1 - F1(X) F1(X) = LOLP' * WRITE(3,'(I6,I8,2F12.8)') (J-I,PMX1-I,1.D0-F1(I),F1(I),I=J,-1,-1) * shrink IMXO to the range of significance for FO(X) LOLP DO 39 JX=0,MF D4=1.-FO(JX) IF(D4.LT..99999996) THEN ! test for significance IMXO=JX ! IMXO is as far out on the curve as we need to go ELSE DO 38 K=JX,MF ! the LOLP beyond IMXO is 0 38 FO(K)=0. GOTO 40 ENDIF 39 CONTINUE 40 PMNO=PMXO-IMXO IF(IMXO.EQ.MF) THEN ! cumulative distribution is too wide, increase the MF value WRITE(*,'(13Hwarning: 1-F(,I5,3H) =,F12.8)') MF,FO(MF) ENDIF * J=IMXO+1 ! this is to make sure the print statement doesn't run greater than MF * IF(J.GT.MF) J=MF * WRITE(3,*) 'PMXO=',PMXO,' PMNO=',PMXO-IMXO,' IMXO=',IMXO * WRITE(3,*) ' I MW 1 - FO(X) FO(X) = LOLP' * WRITE(3,'(I6,I8,2F12.8)') (J-I,PMXO-I,1.D0-FO(I),FO(I),I=J,-1,-1) * 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 D5=300. ! MIN PEAK LOAD IN AREA 1 D6=700. ! WIN PEAK LOAD IN AREA 1 D7=846. ! SUM PEAK LOAD IN AREA 1 note area 1 is summer peaking 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 PD=0. ! REMEMBER THE PEAK DEMAND PD1=0. ! REMEMBER THE AREA 1 PEAK DEMAND PDO=0. ! REMEMBER OUTSIDE AREA 1 PEAK DEMAND * SWEEP THROUGH THE YEAR CALCULATING LOADS AND INDICES DO 3 J=1,365 DO 3 I=1,24 HL(I,J)=IFIX(D1*DAILY(I)+.01) ! LIGHT LOAD PERIOD HL1(I,J)=IFIX(D5*DAILY(I)+.01) ! LIGHT LOAD PERIOD AREA 1 IF(J.LE.10) HL(I,J)=IFIX(D2*DAILY(I)+.01) ! 10 WINTER PEAK DEMANDS IF(J.LE.10) HL1(I,J)=IFIX(D6*DAILY(I)+.01) ! 10 WINTER PEAK DEMANDS IN AREA 1 IF(J.GE.201.AND.J.LE.214) HL(I,J)=IFIX(D3*DAILY(I)+.01) ! 14 SUMMER PEAK DEMANDS IF(J.GE.201.AND.J.LE.214) HL1(I,J)=IFIX(D7*DAILY(I)+.01) ! 14 SUMMER PEAK DEMANDS IN AREA 1 HL(I,J)=HL(I,J)+D4 ! ADD THE CONSTANT LOAD IF(HL(I,J).GT.PD) PD=HL(I,J) ! REMEMBER THE PEAK DEMAND IF(HL1(I,J).GT.PD1) PD1=HL1(I,J) ! REMEMBER THE PEAK DEMAND IN AREA 1 IF(HL(I,J)-HL1(I,J).GT.PDO) PDO=HL(I,J)-HL1(I,J) ! REMEMBER THE PEAK DEMAND OUTSIDE AREA 1 IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR JX=IX(HL(I,J),PMX,MF) LOLE=LOLE+F(JX) ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*(1.D0-F(JX)) ! PROBABILITY OF DAILY GEN UP STATES ENDIF IF(PMN.LE.HL(I,J)) THEN ! WHERE LOLP > 0 JX=IX(HL(I,J),PMX,MF) DO 2 K=PMN,HL(I,J) EUE=EUE+F(JX) ! EXPECTED UNSERVED ENERGY - THE AREA UNDER F(X) 2 JX=JX+1 ENDIF JX=IX(HL(I,J),PMX,MF) LOLH=LOLH+F(JX) ! LOLH = RUNNING SUM OF THE HOURLY LOLP 3 ENGY=ENGY+HL(I,J) ALLR=1.D0-ALLR * SHOW HOURLY LOADS FOR THE TOTAL SYSTEM AND AREA 1 WRITE(3,*)' ' WRITE(3,*)' TOTAL SYSTEM HOURLY LOADS:' 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, & HL(I, 1),F(IX(HL(I, 1),PMX,MF)), & HL(I,201),F(IX(HL(I,201),PMX,MF)), & HL(I, 11),F(IX(HL(I, 11),PMX,MF)) WRITE(3,'(/'' D1 MIN, D2 WIN, D3 SUM, D4 ADD ='',4I7/)') & IFIX(D1),IFIX(D2),IFIX(D3),IFIX(D4) WRITE(3,*)' AREA 1 HOURLY LOADS:' WRITE(3,*)' WINTER SUMMER OFF PEAK' WRITE(3,*)'HR MW LOLP MW LOLP MW LOLP' DO 24 I=1,24 24 WRITE(3,'(I3,3(F7.0,F12.8))') I, & HL1(I, 1),F1(IX(HL1(I, 1),PMX1,MF)), & HL1(I,201),F1(IX(HL1(I,201),PMX1,MF)), & HL1(I, 11),F1(IX(HL1(I, 11),PMX1,MF)) WRITE(3,'(/'' D5 MIN, D6 WIN, D7 SUM ='',3I7/)') & IFIX(D5),IFIX(D6),IFIX(D7) * DISPLAY RESULTS ALOLP=LOLE/AMTTR**.4 ! APPROXIMATE ANNUAL LOLP, WHEN MTTR=1DAY ALOLP=LOLE AND ALOLP IS OBSERVED TO DROP OFF AS AVERAGE MTTR IS INCREASED CALL CPU_TIME(T1) WRITE(3,10) ' SINGLE AREA SA DIRECT SOLUTN TIME =',T1-T0, & LOLH,LOLE,ALLR,ALOLP,EUE,ENGY,EUE/ENGY*1.D6 10 FORMAT(/A36,F5.2,' SEC'/ & ' --------------------------------------------'/ & 4X,' LOLH LOLE ALLR ~ALOLP'/4F11.6// & 4X,'MWHEUE MWHTOTAL puEUEppm'/F11.5,F12.0,F10.5) WRITE(3,'(/19H RESERVE MARGIN =,F8.2,2H %)') & (PMX-PD)/PD*100. * SHOW THE DIRECT SOLUTION LOLP FOR WINTER, SUMMER, AND OFF PEAK HOURS 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(*,'(/'' TIME ='',F10.3,'' SEC''/)') T1-T0 T0=T1 WRITE(*,*) 'ENTER NEW D1,D2,D3,D4 total area, D5,D6,D7 area 1' WRITE(*,*) 'OR RETURN TO CONTINUE' READ(*,'(A100)') A IF(A.NE.' ') THEN READ(A,*,ERR=41) D1,D2,D3,D4,D5,D6,D7 GOTO 9 ENDIF * Perform 2 area model calculations - follow the steps below to see the logic 14 LOLE=0.D0 ! LOSS OF LOAD EXPECTATION FOR AREA 1 LOLH=0.D0 ! LOSS OF LOAD HOURS FOR AREA 1 ALLR=1.D0 ! ANNUAL LOAD LOSS RISK FOR AREA 1 EUE=0.D0 ! EXPECTED UNSERVED ENERGY FOR AREA 1 ENGY=0.D0 ! TOTAL LOAD ENERGY FOR AREA 1 WRITE(3,*) ' ' WRITE(3,*) '2 AREA INTERCONNECTION' WRITE(*,*) '2 AREA INTERCONNECTION' WRITE(3,*) '+MW=IMPORT -MW=EXPORT' WRITE(*,*) '+MW=IMPORT -MW=EXPORT' WRITE(3,*) ' ATCMW PROBLTY' WRITE(*,*) ' ATCMW PROBLTY' ATC(2,11)=1. ! initialize the probability of the +0 MW ATC (import) to area 1, ATC(1,11) always equals 0 MW ATC(2,12)=1. ! initialize the probability of the -0 MW ATC (export) from area 1, ATC(1,12) always equals 0 MW DO 13 K=1,10 ! only the first 10 are user defined, 11 and 12 are calculated 0 MW import and export power flow states IF(ATC(1,K).GT.0.) ATC(2,11)=ATC(2,11)-ATC(2,K) IF(ATC(1,K).LT.0.) ATC(2,12)=ATC(2,12)-ATC(2,K) IF(ATC(1,K).EQ.0.) GOTO 13 WRITE(3,'(I8,F11.5)') IFIX(ATC(1,K)),ATC(2,K) WRITE(*,'(I8,F11.5)') IFIX(ATC(1,K)),ATC(2,K) 13 CONTINUE IF(ABS(ATC(2,11)).LT.1.E-6) ATC(2,11)=0. ! THIS IS A ROUNDOFF PROBLEM IF(ABS(ATC(2,12)).LT.1.E-6) ATC(2,12)=0. IF(ATC(2,11).LT.0.) THEN ! THE SUM OF THE ATC IMPORT PROBABILITIES MUST SUM TO 1 AT MOST WRITE(3,*) 'WARNING, ATC IMPORT PROBABILITY ERROR, TRY AGAIN' GOTO 20 ENDIF IF(ATC(2,12).LT.0.) THEN ! THE SUM OF THE ATC EXPORT PROBABILITIES MUST SUM TO 1 AT MOST WRITE(3,*) 'WARNING, ATC EXPORT PROBABILITY ERROR, TRY AGAIN' GOTO 20 ENDIF IF(ATC(1,11).LT.1.) WRITE(3,'(A8,F11.5)') ' +0',ATC(2,11) IF(ATC(1,11).LT.1.) WRITE(*,'(A8,F11.5)') ' +0',ATC(2,11) IF(ATC(1,12).LT.1.) WRITE(3,'(A8,F11.5)') ' -0',ATC(2,12) IF(ATC(1,12).LT.1.) WRITE(*,'(A8,F11.5)') ' -0',ATC(2,12) * BEGIN THE HOURLY ANALYSIS FOR TWO AREAS * CALCULATE AREA 1 FIRST DO 11 J=1,365 ! SWEEP THROUGH EVERY DAY DO 11 I=1,24 ! SWEEP THROUGH EVERY HOUR ATCP(11)=1. ! WE HAVE TO PICK UP THE PROBABILITY OF THE 0 MW ATC IMPORT STATE DO 12 K=1,10 ! GO THROUGH ALL THE ATC STATES ATCX(K)=0. ! THIS IS THE ATC MW THAT WILL BE APPLIED IF ATCS ARE SPECIFIED THAT ARE TOO LARGE ATCP(K)=0. ! ATCP WILL STORE THE PROBABILITY OF THE Kth ATC MW TAKING BOTH TRANSMISSION AND OUTSIDE GENERATION INTO ACCOUNT IF(ATC(1,K).LE.0.) GOTO 12 ! only the import constraints > 0 are used IAS=0 ! we are going to send the area 0 outside power for this hour through the ATC constraits IAMX=0 ! this will be the largest ATC DO 19 IA=1,ATC(1,K) ! WE WILL DO THIS ATC(1,K) TIMES AND THEN DIVIDE THE ATCP BY THAT VALUE TO GET THE AVERAGE IF(IA.GT.PD1) GOTO 19 ! THE SPECIFIED ATC EXCEEDS THE PEAK DEMAND X=HL(I,J)-HL1(I,J)+IA JX=IX(X,PMXO,MF) IF(FO(JX).GT..3) GOTO 19 ! THERE IS NOT ENOUGH GENERATION LEFT IN THE SENDING AREA IAS=IAS+IA ! SUM OF IA'S ATCP(K)=ATCP(K)+(1.D0-FO(JX))*IA ATCX(K)=IA IF(IA.GT.IAMX) IAMX=IA ! CAPTURE THE LARGEST ATC FOR THIS HOUR AS IAMX (INTEGER ATC MAX) 19 CONTINUE ATCP(K)=ATC(2,K)*ATCP(K)/IAS ! PRODUCT OF THE TRANSMISSION PROBABILITY TIMES THE WEIGHTED OUTSIDE GENERATION AVAILABILITY ATCP(11)=ATCP(11)-ATCP(K) ! THIS IS THE PROBABILITY OF THE 0 MW ATC STATE, I.E. IT GOES ALONG WITH ATC(1,11)=0 MW 12 CONTINUE PMXT=IAMX+PMX1 ! NOW CONVOLVE THE NEWLY CREATED ATC STATES INTO F1 IMXT=IMX1+IAMX IF(IMXT.GT.MF) IMXT=MF FT(-1)=1. DO 15 JX=IMXT,0,-1 ! WE WILL CREATE FT ONE LAYER AT A TIME. BEGIN WITH THE NO ATC MW STATE. J1=JX-IAMX ! UNLIKE MC4, THIS NO ATC STATE SHIFTS THE F1 TO THE RIGHT THE MAX AMOUNT, IAMX IF(J1.LT.-1) J1=-1 15 FT(JX)=F1(J1)*ATCP(11) DO 16 K=1,10 ! PICK UP THE REST OF THE ATC +MW IMPORT OF POWER STATES IF(ATCX(K).GT.0.) THEN ! ONLY POSITIVE ATCS ARE USED. THEIR SUM OF PR SHOULD BE 1. DO 18 JX=IMXT,0,-1 J1=(JX-IAMX)+ATCX(K)+.001 ! THE SHIFT TO THE RIGHT IS IAMX-ATCX, THE DIFFERENCE IN THE MAX ATC AND THE ATC, I.E. THE MW IF THE ATC FAILS. IF(J1.LT.-1) J1=-1 18 FT(JX)=FT(JX) + F1(J1)*ATCP(K) ! ATCS ARE TREATED AS A MULTI STATE GENERATOR, EACH MW LEVEL BEING SEPARATELY SHIFTED AND ADDED ENDIF 16 CONTINUE IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR FOR AREA 1 D8=FT(IX(HL1(I,J),PMXT,MF))! DAILY PEAK LOLP LOLE=LOLE+D8 ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*(1.D0-D8) ! PROBABILITY OF DAILY GEN UP STATES IN AREA 1 PLUS THE ATCS TIE LINE ENDIF IF(PMN1.LE.HL1(I,J)) THEN ! THE LOLP IS ZERO BELOW PMN1 DO 17 K=PMN1,HL1(I,J) X=K 17 EUE=EUE+FT(IX(X,PMXT,MF)) ! EXPECTED UNSERVED ENERGY - THE AREA 1 BELOW FT(X) TO THE RIGHT OF THE LOAD X ENDIF LOLH=LOLH+FT(IX(HL1(I,J),PMXT,MF)) ! LOLH = RUNNING SUM OF THE HOURLY LOLP ENGY=ENGY+HL1(I,J) 11 CONTINUE ALLR=1.D0-ALLR * DISPLAY RESULTS ALOLP=LOLE/AMTTR**.4 CALL CPU_TIME(T1) WRITE(3,*) WRITE(3,10) ' AREA 1 A1 DIRECT SOLUTION TIME = ',T1-T0, & LOLH,LOLE,ALLR,ALOLP,EUE,ENGY,EUE/ENGY*1.D6 ENGYA1=ENGY ! WILL USE THIS LATER IN THE MC PRINTOUT SECTION WRITE(3,'(/19H RESERVE MARGIN =,F8.2,2H %)') & (PMX1-PD1)/PD1*100. T0=T1 WRITE(*,*) ' ' WRITE(*,*) ' AREA 1 DIRECT SOLUTION INDICES' WRITE(*,*) ' LOLH LOLE puEUEppm' WRITE(*,'(2F11.6,F11.3)') LOLH,LOLE,EUE/ENGY*1.D6 * CALCULATE THE OUTSIDE OF AREA 1 SECONDLY LOLE=0.D0 ! LOSS OF LOAD EXPECTATION FOR OUTSIDE AREA 1 LOLH=0.D0 ! LOSS OF LOAD HOURS FOR OUTSIDE AREA 1 ALLR=1.D0 ! ANNUAL LOAD LOSS RISK FOR OUTSIDE AREA 1 EUE=0.D0 ! EXPECTED UNSERVED ENERGY FOR OUTSIDE AREA 1 ENGY=0.D0 ! TOTAL LOAD ENERGY FOR OUTSIDE AREA 1 DO 21 J=1,365 ! SWEEP THROUGH EVERY DAY DO 21 I=1,24 ! SWEEP THROUGH EVERY HOUR ATCP(12)=1. ! WE HAVE TO PICK UP THE PROBABILITY OF THE 0 MW ATC EXPORT STATE DO 22 K=1,10 ! GO THROUGH ALL THE ATC STATES ATCX(K)=0. ! THIS IS THE ATC MW THAT WILL BE APPLIED IF ATCS ARE SPECIFIED THAT ARE TOO LARGE ATCP(K)=0. ! ATCP WILL STORE THE PROBABILITY OF THE Kth ATC MW TAKING BOTH TRANSMISSION AND OUTSIDE GENERATION INTO ACCOUNT IF(ATC(1,K).GE.0.) GOTO 22 ! only the export constraints < 0 are used IAS=0 IAMX=0 ! this will be the largest ATC DO 29 IA=-1,ATC(1,K),-1 ! DO THE ATC(1,K) TIMES AND THEN DIVIDE BY THE ATCP BY THAT VALUE TO GET THE AVERAGE IF(-IA.GT.PDO) GOTO 29 ! THE SPECIFIED ATC EXCEEDS THE PEAK DEMAND X=HL1(I,J)-IA ! KEEP IN MIND THAT THIS ATC IS NEGATIVE, SO A '-' MEANS ADD THE ATC TO THE LOAD JX=IX(X,PMX1,MF) IF(F1(JX).GT..3) GOTO 29 ! THERE IS NOT ENOUGH GENERATION LEFT IN THE SENDING AREA IAS=IAS-IA ATCP(K)=ATCP(K)-(1.D0-F1(JX))*IA ATCX(K)=-IA IF(-IA.GT.IAMX) IAMX=-IA ! CAPTURE THE LARGEST ATC FOR THIS HOUR AS IAMX (INTEGER ATC MAX) 29 CONTINUE ATCP(K)=ATC(2,K)*ATCP(K)/IAS ! PRODUCT OF THE TRANSMISSION PROBABILITY TIMES THE WEIGHTED OUTSIDE GENERATION AVAILABILITY ATCP(12)=ATCP(12)-ATCP(K) ! THIS IS THE PROBABILITY OF THE 0 MW ATC STATE, I.E. IT GOES ALONG WITH ATC(1,12)=0 MW 22 CONTINUE PMXT=IAMX+PMXO ! NOW CONVOLVE THE NEWLY CREATED ATC STATES INTO FO IMXT=IMXO+IAMX IF(IMXT.GT.MF) IMXT=MF FT(-1)=1. DO 25 JX=IMXT,0,-1 ! WE WILL CREATE FT ONE LAYER AT A TIME. BEGIN WITH THE NO ATC MW STATE. J1=JX-IAMX ! UNLIKE MC4, THIS NO ATC STATE SHIFTS THE F1 TO THE RIGHT THE MAX AMOUNT, IAMX IF(J1.LT.-1) J1=-1 25 FT(JX)=FO(J1)*ATCP(12) DO 26 K=1,10 ! PICK UP THE REST OF THE ATC +MW IMPORT OF POWER STATES IF(ATCX(K).GT.0.) THEN ! ONLY POSITIVE ATCS ARE USED. THEIR SUM OF PR SHOULD BE 1. DO 28 JX=IMXT,0,-1 J1=(JX-IAMX)+ATCX(K)+.001 ! THE SHIFT TO THE RIGHT IS IAMX-ATCX, THE DIFFERENCE IN THE MAX ATC AND THE ATC, I.E. THE MW IF THE ATC FAILS. IF(J1.LT.-1) J1=-1 28 FT(JX)=FT(JX) + FO(J1)*ATCP(K) ! ATCS ARE TREATED AS A MULTI STATE GENERATOR, EACH MW LEVEL BEING SEPARATELY SHIFTED AND ADDED ENDIF 26 CONTINUE JX=HL(I,J)-HL1(I,J)+.001 ! LOAD OF THE OUTSIDE AREA X=JX IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR FOR AREA 1 D8=FT(IX(X,PMXT,MF)) LOLE=LOLE+D8 ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*(1.D0-D8) ! PROBABILITY OF DAILY GEN UP STATES IN AREA 1 PLUS THE ATCS TIE LINE ENDIF IF(PMNO.LE.JX) THEN ! THE LOLP IS ZERO BELOW PMNO DO 27 K=PMNO,JX XX=K 27 EUE=EUE+FT(IX(XX,PMXT,MF))! EXPECTED UNSERVED ENERGY - THE AREA 1 BELOW FT(X) TO THE RIGHT OF THE LOAD X ENDIF LOLH=LOLH+FT(IX(X,PMXT,MF)) ! LOLH = RUNNING SUM OF THE HOURLY LOLP ENGY=ENGY+X 21 CONTINUE ALLR=1.D0-ALLR * DISPLAY RESULTS ALOLP=LOLE/AMTTR**.4 CALL CPU_TIME(T1) WRITE(3,*) WRITE(3,10) ' AREA 0 A0 DIRECT SOLUTION TIME = ',T1-T0, & LOLH,LOLE,ALLR,ALOLP,EUE,ENGY,EUE/ENGY*1.D6 WRITE(3,'(/19H RESERVE MARGIN =,F8.2,2H %)') & (PMXO-PDO)/PDO*100. WRITE(*,*) ' ' WRITE(*,*) ' AREA 0 DIRECT SOLUTION INDICES' WRITE(*,*) ' LOLH LOLE puEUEppm' WRITE(*,'(2F11.6,F11.3)') LOLH,LOLE,EUE/ENGY*1.D6 WRITE(*,'(/'' TIME ='',F10.3,'' SEC''/)') T1-T0 T0=T1 20 WRITE(*,*) 'ENTER NEW ATC MW,PR PAIRS OR RETURN TO CONTINUE' READ(*,'(A100)') A A(61:100)=',0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0' IF(A(1:60).NE.' ') THEN READ(A,*,ERR=41) ((ATC(I,J),I=1,2),J=1,10) GOTO 14 ENDIF * * BEGIN THE (Area 1) SEQUENTIAL MONTECARLO SIMULATION, WILL BE REPORTED IN INTERVALS OF 1000, 10K, 100K, ETC YEARS * In this simulation the ATCs for import power are used and the export ATCs are not used. The first +MW,Pr entry for the ATC * should be the N-0 all lines in service state. The ATC Monte Carlo will assume that there is a transition from this up state * to a lower state and the repair time for that lower state is assumed to be 8 hours to get back to the upper state. There can * be several lower MW states of lower probability. The probabilities of all these ATC states are summed and what is remaining that * is less than 1 is assigned to a 0 MW probability state in array ATC(2,11). The 2 is for probability and the 11 is for import. * Just keep in mind that ATC(1,1) is the MW import power of the N-0 state and its probability is ATC(2,1). You entered this * earlier. All the data entering the Monte Carlo solution below comes from the last run using the direct solution. OPEN(5,STATUS='SCRATCH',FORM='UNFORMATTED',ACCESS='SEQUENTIAL') ! SINGLE AREA SA EVENTS OPEN(6,STATUS='SCRATCH',FORM='UNFORMATTED',ACCESS='SEQUENTIAL') ! AREA 1 EVENTS OPEN(7,STATUS='SCRATCH',FORM='UNFORMATTED',ACCESS='SEQUENTIAL') ! ATC EVENTS WRITE(3,*) ' ' WRITE(3,*) ' ' WRITE(3,*) ' ' WRITE(3,*) 'FD SEQUENTIAL MONTECARLO SOLUTION' WRITE(*,*) 'SETTING UP THE SEQUENTIAL MONTE CARLO SIMULATION' * WRITE(*,*)'ENTER A GLOBAL MTTR (DAYS) OR RETURN FOR MTTR DEFAULTS' GMTTR=0. * READ(*,'(A100)') A * IF(A.NE.' ') READ(A,*,ERR=41) GMTTR M=3 ! SEED FOR RANDOM NUMBER GENERATOR 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 D4=FOR(I,1)+FOR(I,2) ! PROBABILITY OF ALL THE FAILURE STATES R=RAN(M) IF(R.LT.D4) ISTATE(K1+I)=0 ! DETERMINE GENERATORS IN A FAILURE STATE AT START OF THE RUN IF(R.LT.FOR(I,2)) ISTATE(K1+I)=2 ! PARTIAL OUTAGE STATE, NOTE THAT MOST FOR(-,2)=0 SO R WOULD NEVER BE LESS THAN 0 IF(K.EQ.1) THEN IF(GMTTR.GT.0.) MTTR(I)=GMTTR ! USE NEW GLOBAL MTTR MEAN TIME TO REPAIR VALUE MTTR(I)=MTTR(I)*24. ! CONVERT MEAN TIME TO REPAIR FROM DAYS TO HOURS ENDIF RR(K1+I)=1./MTTR(I) ! CALCULATE THE REPAIR RATE TRANSITION PROBABILITY FOR EACH HOUR 5 FR(K1+I)=RR(K1+I)*D4/(1.-D4) ! CALCULATE THE FAILURE RATE FROM THE REPAIR RATE AND FORCED OUTAGE RATE, NOTE THAT PARTIAL FOR IS ALSO INCLUDED RRT=8. ! 8 HOURS TO REPAIR TRANSMISSION IF(GMTTR.GT.0.) AMTTR=GMTTR WRITE(*,'(20H TRANSMISSION MTTR =,F9.3,5H DAYS)') RRT/24. WRITE(*,'(20H GENERATION MTTR =,F9.3,5H DAYS)') AMTTR * WRITE(*,*) 'ENTER A NEW MTTR HOURS VALUE OR RETURN' * READ(*,'(A100)') A * IF(A.NE.' ') READ(A,*,ERR=41) RRT WRITE(3,'(20H TRANSMISSION MTTR =,F9.3,5H DAYS)') RRT/24. WRITE(3,'(20H GENERATION MTTR =,F9.3,5H DAYS)') AMTTR WRITE(3,*) &'---------------------------------------------------------' RRT=1./RRT IST=1 ! ASSUME THE TRANSMISSION IS IN STATE 1 OR THE N-0 STATE INITIALLY ATC21=1.-ATC(2,1) ! THIS IS THE PROBABITY OF NOT BEING IN N-0 AND IS USED FREQUENTLY SO WE SAVE IT HERE FRT=RRT*ATC21/ATC(2,1) ! FAILURE RATE OF TRANSMISSION IS CALCULATED FROM THE RRT AND PROBABILITY OF BEING IN N-0 STATE LOLE=0.D0 ! NUMBER OF LOSS OF LOAD EXPECTATION EVENTS, NO MORE THAN 1 PER DAY WILL BE COUNTED. LOLE1=0.D0 ! NUMBER OF LOSS OF LOAD EXPECTATION EVENTS FOR AREA 1, 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 LOLH1=0.D0 ! LOSS OF LOAD HOURS FOR AREA 1, ADD 1 TO EVERY HOUR THERE IS AN OBSERVED LOSS OF LOAD ALOLP=0.D0 ! ANNUAL LOAD LOSS RISK, WILL COUNT THE NUMBER OF YEARS THERE ARE NO LOSS OF LOAD EVENTS ALOLP1=0.D0 ! ANNUAL LOAD LOSS RISK FOR AREA 1, 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 EUE1=0.D0 ! EXPECTED UNSERVED ENERGY FOR AREA 1, WILL SUM THE DIFFERENCE BETWEEN THE LOAD AND GENERATION AVAILABLE EACH HOUR MULT=1000 ! START OUT WITH 1000 YEARS OF SIMULATION BEFORE REPORTING RESULTS WRITE(*,*) 'NUMBER OF MONTECARLO ITERATIONS?' READ(*,*) ITERMAX IF(ITERMAX.LT.MULT) STOP YRS=0 * BEGIN THE LOOP SIMULATING EACH YEAR WITH SEQUENTIAL HOURS WITHIN THAT YEAR 6 YRS=YRS+1 NED=0 ! COUNT THE NUMBER OF EVENT DAYS JUST FOR THIS ONE YEAR NED1=0 ! COUNT THE NUMBER OF AREA 1 EVENT DAYS JUST FOR THIS ONE YEAR DO 8 J=1,365 NEED=0 ! COUNT THE NUMBER OF EVENTS EACH DAY NEED1=0 ! COUNT THE NUMBER OF EVENTS EACH DAY FOR AREA 1 DO 8 I=1,24 GEN=0. ! THE MW GENERATION ON LINE THIS HOUR GEN1=0. ! THE MW GENERATION ON LINE IN AREA 1 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) K1=K-((K-1)/10)*10 ! K1 IS THE ORIGINAL GENERATOR NUMBER, RANGE IS 1-10 IF(ISTATE(K).EQ.1) THEN IF(R.LE.FR(K)) THEN ! WE HAVE DETERMINED THAT WE ARE EITHER GOING INTO FAILURE OR PARTIAL FAILURE, WHICH? ISTATE(K)=0 ! ASSUME WE FLIP THE GENERATOR TO A FAILURE STATE R=RAN(M) ! BUT DRAW ANOTHER RANDOM NUMBER TO SEE IF WE ARE GOING INTO A PARTIAL FAILURE STATE IF(R.LE.FOR(K1,2)) ISTATE(K)=2 ! INSTEAD THE K1 GENERATOR GOES INTO A PARTIAL POWER STATE IF THIS TEST IS MET ENDIF ELSE ! GETTING TO HERE MEANS WE WERE NOT IN FULL POWER STATE BUT EITHER PARTIAL OR FULL OUTAGE IF(R.LE.RR(K)) ISTATE(K)=1 ! GENERATOR IS REPAIRED FROM EITHER A FULL OUTAGE STATE OR A DERATED STATE ENDIF K1=K-((K-1)/10)*10 ! THIS FINDS THE CORRECT GENERATOR PMAX IF(ISTATE(K).EQ.1) GEN=GEN+PMAX(K1,1)! SUM UP THE POWER OF THE GENS ON LINE FULLY UP FOR THE TOTAL SYSTEM IF(ISTATE(K).EQ.2) GEN=GEN+PMAX(K1,2)! SUM UP THE POWER OF THE GENS ON LINE PARTIAL POWER STATE IF(K.GT.10) GOTO 7 ! THIS WILL TEST ALL 10 GENERATORS IN AREA 1 AS WELL AS TEST THE 10 ATC STATES IF(ISTATE(K).EQ.1) GEN1=GEN1+PMAX(K1,1)! SUM UP THE POWER OF THE GENS ON LINE IF(ISTATE(K).EQ.2) GEN1=GEN1+PMAX(K1,2)! SUM UP THE POWER OF THE GENS ON LINE PARTIAL POWER STATE 7 CONTINUE * NOW GO THROUGH THE TRANSMISSION STATES TO FIND WHICH STATE WE ARE IN R=RAN(M) ! A NEW RANDOM NUMBER FOR THE TRANSMISSION STATE IF(IST.EQ.1) THEN ! WE ARE IN STATE 1 IF(R.GT.FRT) THEN ! STAY IN STATE 1 ATCH(1)=ATCH(1)+1.D0 ! COUNT THE NUMBER OF ATC HOURS IN STATE 1, LATER THIS ARRAY IS USED ONLY WHEN TRANSMISSION WAS NEEDED GOTO 30 ENDIF R=RAN(M) ! WE ARE GETTING READY TO EXIT OUT OF STATE 1 TO ONE OF THE OTHER STATES PRS=0. ! PROBABILITY SUM FROM THE BOTTOM UP DO 23 K2=11,1,-1 IF(ATC(2,K2).LE.0..OR.ATC(1,K2).LT.0.) GOTO 23 ! SKIP THIS IT'S NOT A VALID IMPORT CONSTRAINT PRS=PRS+ATC(2,K2) IF(R.LE.PRS/ATC21) THEN ! TRANSITION TO STATE K2 IST=K2 ! REMEMBER THIS STATE WE ARE IN ATCH(K2)=ATCH(K2)+1.D0 ! KEEP A COUNT OF THE HOURS IN THIS STATE GOTO 30 ENDIF 23 CONTINUE WRITE(*,*)'CANNOT BE HERE AT 1' STOP ELSE ! IS IS GREATER THAN 1 IF(R.GT.RRT) THEN ! STAY IN THIS OTHER LOWER STATE FOR NOW ATCH(IST)=ATCH(IST)+1.D0 GOTO 30 ELSE IST=1 ! TRANSITION BACK TO STATE 1 ATCH(1)=ATCH(1)+1.D0 GOTO 30 ENDIF WRITE(*,*)'CANNOT BE HERE AT 2' STOP ENDIF WRITE(*,*)'CANNOT BE HERE AT 3' STOP 30 IF(GEN+.0001.GE.HL(I,J)) THEN ! THERE IS ENOUGH GENERATION, THE .0001 INSURES THE IF TEST IS ACCURATE IN CASE REAL NUMBERS HAPPEN TO END IN .999 ELSE EUE=EUE+HL(I,J)-GEN ! HL EXCEEDS THE GENERATION MW, ADD THE MW TIMES ONE HOUR, OR UNSERVED ENERGY IN MWH LOLH=LOLH+1.D0 ! 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 * see what generation is available to area 1 through the tie lines GEN=GEN-GEN1 ! MAKE GEN THE OUTSIDE AREA 1 GENERATION GEN=GEN-(HL(I,J)-HL1(I,J)) ! SUBTRACT THE OUTSIDE AREA 1 LOAD FROM GEN TO SEE WHATS LEFT OVER IN OUTSIDE GENERATION IF(GEN.GT.ATC(1,IST)) GEN=ATC(1,IST) ! POWER CANNOT BE GREATER THAN WHAT CAN GET THROUGH THE TIE LINES IF(GEN.GT.0.) GEN1=GEN1+GEN ! ADD THE TIE LINE POWER FLOW TO THE AREA 1 GENERATION AND THEN DO THE LOAD TEST IF(GEN1+.0001.GE.HL1(I,J)) THEN ! REPEAT EVERYTHING FOR AREA 1 ELSE IF(GEN.EQ.ATC(1,IST)) ATCC(IST)=ATCC(IST)+1.D0 ! COUNT WHICH TRANSMISSION STATES ARE CAUSING LOSS OF LOAD AND NOT JUST LACK OF OUTSIDE GENERATION EUE1=EUE1+HL1(I,J)-GEN1 LOLH1=LOLH1+1.D0 ! COUNT THE HOURS THERE IS LOSS OF LOAD IN AREA 1 IF(NEED1.EQ.0) THEN NED1=NED1+1 ! THIS IS THE FIRST OCCURENCE OF LOSS OF LOAD THIS DAY IN AREA 1 NEED1=1 ! DON'T ALLOW ANY MORE TO BE COUNTED THIS DAY IN AREA 1 ENDIF ENDIF 8 CONTINUE LOLE=LOLE+NED ! SUM UP THE EVENTS FOR EACH YEAR AND LOLE STORES THE RUNNING SUM LOLE1=LOLE1+NED1 ! SUM UP THE EVENTS FOR EACH YEAR AND LOLE STORES THE RUNNING SUM FOR AREA 1 IF(NED.EQ.0) ALOLP=ALOLP+1.D0 ! COUNT THE NUMBER OF YEARS THERE WERE NO EVENTS IF(NED1.EQ.0) ALOLP1=ALOLP1+1.D0 ! COUNT THE NUMBER OF YEARS THERE WERE NO EVENTS IN AREA 1 * WRITE ANNUAL RESULTS TO SCRATCH FILE 5 * V1 V2 V3 V4 WRITE(5) LOLH/YRS,LOLE/YRS,1.-ALOLP/YRS,EUE/YRS ! THIS IS OUTPUT FOR THE SINGLE AREA WRITE(6) LOLH1/YRS,LOLE1/YRS,1.-ALOLP1/YRS,EUE1/YRS ! THIS IS OUTPUT FOR AREA 1 WRITE(7) (ATCC(I)/YRS,I=1,12) ! AVERAGE TRANSMISSION EVENTS PER YEAR (USED ONLY FOR THE STD DEVIATION CALCULATION) IF(YRS.LT.MULT.AND.YRS.LT.ITERMAX) GOTO 6 * DISPLAY RESULTS FOR CHOSEN NUMBER OF YEARS 1K 10K 100K ETC REWIND(5) REWIND(6) REWIND(7) * SINGLE AREA STANDARD DEVIATION %'S V(1)=LOLH/YRS V(2)=LOLE/YRS V(3)=1.-ALOLP/YRS V(4)=EUE/YRS DO 42 I=1,4 42 S(I)=0. DO 43 J=1,YRS READ(5) (B(K),K=1,4) DO 43 I=1,4 43 S(I)=S(I)+((B(I)-V(I))/V(I))**2 ! SD formulas corrected 02/23/2016 DO 44 I=1,4 S(I)=S(I)/(YRS-1.) 44 SD(I)=SQRT(S(I))*100.+.5 * MONTECARLO SEQUENTIAL SOLUTION * LOLH SD% LOLE SD% ALOLP SD% MWHEUE SD% * 0.196000 15 0.097000 15 0.021000 15 20.065 15 WRITE(3,*) &' LOLH SD% LOLE SD% ALOLP SD% MWHEUE SD%' WRITE(*,*) &' LOLH SD% LOLE SD% ALOLP SD% MWHEUE SD%' WRITE(3,'(F10.6,I3,2(F11.6,I3),F10.3,I3,4H SA)') & (V(I),SD(I),I=1,4) WRITE(*,'(F10.6,I3,2(F11.6,I3),F10.3,I3,4H SA)') & (V(I),SD(I),I=1,4) * AREA 1 STANDARD DEVIATION %'S V(1)=LOLH1/YRS V(2)=LOLE1/YRS V(3)=1.-ALOLP1/YRS V(4)=EUE1/YRS DO 45 I=1,4 45 S(I)=0. DO 46 J=1,YRS READ(6) (B(K),K=1,4) DO 46 I=1,4 46 S(I)=S(I)+((B(I)-V(I))/V(I))**2 DO 47 I=1,4 S(I)=S(I)/(YRS-1.) 47 SD(I)=SQRT(S(I))*100.+.5 WRITE(3,'(F10.6,I3,2(F11.6,I3),F10.3,I3,4H A1)') & (V(I),SD(I),I=1,4) WRITE(*,'(F10.6,I3,2(F11.6,I3),F10.3,I3,4H A1)') & (V(I),SD(I),I=1,4) WRITE(3,'(41X,F10.3,'' ppm A1'')') V(4)/ENGYA1*1.E6 * TRANSMISSION ATCS STANDARD DEVIATION %'S DO 48 I=1,12 V(I)=ATCC(I)/YRS ! ATCC IS THE TOTAL EVENTS/YRS TO GET THE AVERAGE NR EVENTS PER YEAR BUT THIS IS NOT PRINTED OUT 48 S(I)=0. DO 49 J=1,YRS READ(7) (B(K),K=1,12) ! NOW READ THIS RUNNING AVERAGE FOR ALL YEARS DO 49 I=1,12 IF(ABS(V(I)).LT.1.E-20) GOTO 49 S(I)=S(I)+((B(I)-V(I))/V(I))**2 ! CALCULATE THE DEVIATION OVER TIME COMPARED TO THE PRESENT 49 CONTINUE DO 50 I=1,12 S(I)=S(I)/(YRS-1.) ! STANDARD DEVIATION AVERAGE 50 SD(I)=SQRT(S(I))*100.+.5 *ATC LOLTV LOSS OF LOAD TRANSMISSION EVENTS * # ATCMW PROB LOLTV SD% CALCPR * 1 300. 0.0100 821. 10 0.0100 WRITE(3,*) ' ATCS LOSS OF LOAD TRANSMISSION EVENTS' WRITE(3,*) ' # ATCMW PROB LOLTV SD% CALCPR' WRITE(*,*) ' ' WRITE(*,*) ' ATCS LOSS OF LOAD TRANSMISSION EVENTS' WRITE(*,*) ' # ATCMW PROB LOLTV SD% CALCPR' DO 31 I=1,11 IF(ATCH(I).GT.0.D0) WRITE(3,'(I3,F7.0,F8.4,F8.0,I4,F9.4)') & I,ATC(1,I),ATC(2,I),ATCC(I),SD(I),ATCH(I)/(YRS*8760.) IF(ATCH(I).GT.0.D0) WRITE(*,'(I3,F7.0,F8.4,F8.0,I4,F9.4)') & I,ATC(1,I),ATC(2,I),ATCC(I),SD(I),ATCH(I)/(YRS*8760.) 31 CONTINUE CALL CPU_TIME(T1) WRITE(3,'(/F10.2,4H min,I10,4H yrs)') (T1-T0)/60.,YRS WRITE(*,'(/F10.2,4H min,I10,4H yrs)') (T1-T0)/60.,YRS WRITE(3,*) &'---------------------------------------------------------' MULT=MULT*10 IF(YRS.LT.ITERMAX) GOTO 6 41 WRITE(*,*) ' ' WRITE(*,*) 'FILE OP CONTAINS THE OUTPUT REPORT' WRITE(*,*) 'PRESS ENTER TO END THIS PROGRAM' READ(*,'(A1)') A END * * X = MW load, PMX = sum of all Pmax MW, MF = the max array size of F(x) FUNCTION IX(X,PMX,MF) INTEGER PMX IX=PMX-X IF(IX.LT.-1) IX=-1 IF(IX.GT.MF) IX=MF RETURN END