* Your use of this program means you agree Eugene G Preston is not responsible for any errors it may contain. * The MC4 program was written by Gene Preston 09/06/2015 as a variation off MC3 by splitting off area 1. * Area 1 has its own hourly loads that are made to be summer peaking whereas the total system is winter peaking. * The generation in area 1 is ten generators of the original set, i.e. area 1 is for set K=1. The outside generators * are convolved together in FO(x) as in MC3 and the inside area 1 generators are convolved together into F1(x) as * in program MC2. The main difference between MC2 and MC3 is the number of generators. If there were no tie lines * we could simply calculate and look up the indices for area 1 as a stand alone system. Adding a tie line allows * access to area 1 of the generation outside area 1 to improve the reliability of area 1. Also area 1 may have * excess generation it desires to export out of area 1. The effect of the tie line on the import and export capability * is calculated. The area 1 reliabilty indices for import are calculated. The ATCs for power into and out of area 1 * are found from load flow studies with outage contingencies which are given probabilities. So in the ATC array below * we enter pairs of MW,Pr power and probability states for power into area 1 on the tie line which is +MW and then * any ATCs for exporting power out of area 1 are entered as -MW,Pr states. Note that entering any +MW states requires * that the sum of probabilities for power into area 1 sum to 1 so if the sum of MW probabilities is less than 1, the * program assmues there is a 0 MW import probability making the sum of power import state's probabilities equal 1. * Likewise if there are -MW states they will test for limits to power out of the area (such as generation exporting * power out of an area) then if the probabilities of those negative flow states are less than 1, then the program * assumes there is a 0 MW state with probabilty so they sum to 1. Here is how we put all this together in the program. * First we run the single area MC3 and report indices. Next we reset the arrays and create two sets of generation * curves, FO(x) for generation outside area 1 and F1(x) for generation inside area 1. These generation curves are * starting points each hour as the load is constantly changing between the two areas. Storing the FO(x) and F1(x) we * do not have to recalculate the generation curves each hour which saves a whole lot of computer time. So we ask: what * power is available from the larger area to help area 1? Well, we simply look up the probability on FO of the load * outside area 1 which is FO(x=HL-HL1+IA) where IA is varied from 1 to the ATC MW level and FO probabilities noted for * each IA level. A weighted probability of generation available is stored in ATCP. It's possible to specify ATCs * that are much too large so ATCX stores the maximum MW of generation available to send through the transmission link * if the link capacity that is specified as input exceeds this level of available generation. The link MW and probability * states are treated as a multi state generator feeding the load. Once this is done, the indices of area 1 load can be * determined by looking at the generation availability for the desired area 1 load. The outside and inside area 1 loads * are constantly changing so the LOLP is recalculated for each hour. The outside and area 1 LOLP process is inverted * and the LOLP and reliability indices outside area 1 are also calculated. The sequential Monte Carlo program performs * the same calculations for serving area 1. The purpose of including MC in this program is to test its ability to * capture transmission constraint events that affect the area 1 reliability indices and how many iterations are needed * to capture the rare transmission events. The question I want to know is whether MC is too slow and imprecise to be * used as an engineering tool when trying to size the import capacity and reliabity of area 1. Note that the Monte Carlo * solution needs to have the first ATC MW,Pr pair be the import N-0 up state in which all transmission lines are in service. * * BEGIN THE DIRECT SOLUTION CHARACTER*100 A/' '/ REAL*8 F(0:10703)/1.D0,10703*0.D0/ ! THIS IS THE ORIGINAL 973 MW GENERATION DUPLICATED 11 TIMES IN MC3 REAL*8 FO(0:9730)/1.D0, 9730*0.D0/ ! SAVE THE ORIGINAL F(X) FOR THE LAST 10 SETS GO GENERATORS OUTSIDE AREA 1 REAL*8 F1(0: 973)/1.D0, 973*0.D0/ ! 973 MW INTERNAL AREA 1 GENERATION AS IN MC2 REAL*8 FT(0:973),LOLE,ALLR,EUE,ENGY,LOLH ! TEMPORARY F1 RECREATED EACH HOUR AS F1+ATCS AS AN OUTSIDE GENERATOR FEEDING AREA 1 INTEGER PMAX(10)/5,10,13,22,43,57,101,198,248,276/ ! TEN GENERATORS PMAX RATINGS INTEGER PMX/0/,PMX1/0/,PMXO/0/ ! PMX PMX1 PMXO ARE TOTAL GENERATION FOR THE SYSTEM, AREA 1, OUTSIDE AREA 1 REAL*8 FOR(10) & /.15D0,.10D0,.09D0,.07D0,.06D0,.05D0,.04D0,.04D0,.04D0,.04D0/ REAL HL(24,365),HL1(24,365),DAILY(24) ! HL1 ARE AREA 1 HOURLY LOADS AND HL IS THE TOTAL SYSTEM &/.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 REAL ATCP(12) ! this stores the hourly product of ATC probabilty times availability of generation for import/export REAL ATCX(10) ! this is the ATC MW that is convolved. too little sending power may reduce this MW below the specified ATC * 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 FOR GENERATION INTEGER ISTATE(110) ! 0, GEN IS OFF 1, GEN IS ON REAL*8 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,EUE1,ENGY1,LOLH1 ! AREA 1 STATS FOR MONTE CARLO 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) ! TOTAL INSTALLED GENERATION CAPACITY PMX = PMAX SUM IF(K.EQ.1) PMX1=PMX1+PMAX(I) ! AREA 1 INSTALLED GENERATION CAPACITY IF(K.LE.10) PMXO=PMXO+PMAX(I) ! OUTSIDE AREA 1 INSTALLED GENERATION CAPACITY 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 F(IX)=(1.-FOR(I))*F(J)+FOR(I)*F(IX) IF(K.EQ.1) F1(IX)=F(IX) ! AREA 1 GENERATORS IF(K.EQ.10) FO(IX)=F(IX) ! total system excluding area 1 generators (technically should be the last 10, but first 10 is the same set in this example) 1 CONTINUE * 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 LOLE=LOLE+(1.-F(HL(I,J))) ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*F(HL(I,J)) ! PROBABILITY OF DAILY GEN UP STATES ENDIF DO 2 K=1,HL(I,J) 2 EUE=EUE+(1.-F(K)) ! EXPECTED UNSERVED ENERGY - THE AREA ABOVE F(X) AND < PR=1 LOLH=LOLH+(1.-F(HL(I,J))) ! 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),1.-F(HL(I, 1)), & HL(I,201),1.-F(HL(I,201)), & HL(I, 11),1.-F(HL(I, 11)) 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),1.-F1(HL1(I, 1)), & HL1(I,201),1.-F1(HL1(I,201)), & HL1(I, 11),1.-F1(HL1(I, 11)) WRITE(3,'(/'' D5 MIN, D6 WIN, D7 SUM ='',3I7/)') & IFIX(D5),IFIX(D6),IFIX(D7) * DISPLAY RESULTS CALL CPU_TIME(T1) WRITE(3,10) ' SINGLE AREA DIRECT SOLUTION TIME =',T1-T0, & LOLH,LOLE,ALLR,EUE,ENGY,EUE/ENGY*1.D6 10 FORMAT(/A35,F5.2,' SEC'/ & ' -------------------------------------------'/ & 4X,' LOLH LOLE ALLR'/3F11.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 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,*) D1,D2,D3,D4,D5,D6,D7 CALL CPU_TIME(T0)! RESET THE TIME CLOCK 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 CALL CPU_TIME(T0) ! RESET THE TIME CLOCK * 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 IAS=0 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 IX=X+.5 IF(IX.LT.0) IX=0 IF(IX.GT.9730) IX=9730 IF(FO(IX).LT..7) GOTO 19 ! THERE IS NOT ENOUGH GENERATION LEFT IN THE SENDING AREA IAS=IAS+IA ! SUM OF IA'S ATCP(K)=ATCP(K)+FO(IX)*IA ATCX(K)=IA 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 FT(0)=1.D0 DO 15 IX=1,973 ! WE WILL CREATE FT ON LAYER AT A TIME. BEGIN WITH THE NO ATC STATE. 15 FT(IX)=F1(IX)*ATCP(11) ! THIS IS SIMPLY SCALING DOWN F1 TO THE ATCP(11) PROBABILITY AND INITIALIZING FT(X) DO 16 K=1,10 ! PICK UP THE REST OF THE ATC +MW IMPORT OF POWER STATES IF(ATCX(K).LE.0.) GOTO 16 DO 18 IX=1,973 ! WE DONT HAVE TO SWEEP FROM RIGHT TO LEFT BECAUSE F1 IS NOT BEING REPLACED I1=IX-ATCX(K)+.001 IF(I1.LT.0) I1=0 18 FT(IX)=FT(IX) + F1(I1)*ATCP(K) ! ATCS ARE TREATED AS A MULTI STATE GENERATOR, EACH MW LEVEL BEING SEPARATELY SHIFTED AND ADDED 16 CONTINUE IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR FOR AREA 1 LOLE=LOLE+(1.-FT(HL1(I,J))) ! LOLE = RUNNING SUM OF THE DAILY LOLP ALLR=ALLR*FT(HL1(I,J)) ! PROBABILITY OF DAILY GEN UP STATES IN AREA 1 PLUS THE ATCS TIE LINE ENDIF DO 17 K=1,HL1(I,J) 17 EUE=EUE+(1.-FT(K)) ! EXPECTED UNSERVED ENERGY - THE AREA 1 ABOVE FT(X) AND < PR=1 LOLH=LOLH+(1.-FT(HL1(I,J))) ! LOLH = RUNNING SUM OF THE HOURLY LOLP ENGY=ENGY+HL1(I,J) 11 CONTINUE ALLR=1.D0-ALLR * DISPLAY RESULTS CALL CPU_TIME(T1) WRITE(3,*) WRITE(3,10) ' AREA 1 DIRECT SOLUTION TIME =',T1-T0, & LOLH,LOLE,ALLR,EUE,ENGY,EUE/ENGY*1.D6 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 IAS=0 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 IX=X+.5 IF(IX.LT.0) IX=0 IF(IX.GT.973) IX=973 ! THIS IS ALL THE GENERATION THERE IS INSIDE AREA 1 IF(F1(IX).LT..7) GOTO 29 ! THERE IS NOT ENOUGH GENERATION LEFT IN THE SENDING AREA IAS=IAS-IA ATCP(K)=ATCP(K)-F1(IX)*IA ATCX(K)=-IA 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 F(0)=1.D0 ! HAVE TO USE THE LARGER ARRAY F(X) THIS TIME DO 25 IX=1,10703 ! WE WILL CREATE F(X) ON LAYER AT A TIME. BEGIN WITH THE NO ATC STATE. F(IX)=0.D0 IF(IX.GT.9730) GOTO 25 F(IX)=FO(IX)*ATCP(12) ! THIS IS SIMPLY SCALING DOWN FO TO THE ATCP(12) PROBABILITY AND INITIALIZING F(X) 25 CONTINUE DO 26 K=1,10 ! PICK UP THE REST OF THE ATC +MW IMPORT OF POWER STATES IF(ATCX(K).LE.0.) GOTO 26 DO 28 IX=1,9730 ! WE DONT HAVE TO SWEEP FROM RIGHT TO LEFT BECAUSE F1 IS NOT BEING REPLACED I1=IX-ATCX(K)+.001 IF(I1.LT.0) I1=0 28 F(IX)=F(IX) + FO(I1)*ATCP(K) ! ATCS ARE TREATED AS A MULTI STATE GENERATOR, EACH MW LEVEL BEING SEPARATELY SHIFTED AND ADDED 26 CONTINUE IF(I.EQ.12) THEN ! THIS IS THE PEAK HOUR FOR OUTSIDE AREA 1 LOLE=LOLE+(1.-F(HL(I,J)-HL1(I,J))) ! LOLE = RUNNING SUM OF THE DAILY LOLP OUTSIDE AREA 1 ALLR=ALLR*F(HL(I,J)-HL1(I,J)) ! PROBABILITY OF DAILY GEN UP STATES OUTSIDE AREA 1 PLUS THE ATCS TIE LINE ENDIF DO 27 K=1,HL(I,J)-HL1(I,J) 27 EUE=EUE+(1.-F(K)) ! EXPECTED UNSERVED ENERGY - THE AREA 1 ABOVE FT(X) AND < PR=1 LOLH=LOLH+(1.-F(HL(I,J)-HL1(I,J))) ! LOLH = RUNNING SUM OF THE HOURLY LOLP ENGY=ENGY+HL(I,J)-HL1(I,J) 21 CONTINUE ALLR=1.D0-ALLR * DISPLAY RESULTS CALL CPU_TIME(T1) WRITE(3,*) WRITE(3,10) ' BIG GRID DIRECT SOLUTION TIME =',T1-T0, & LOLH,LOLE,ALLR,EUE,ENGY,EUE/ENGY*1.D6 WRITE(3,'(/19H RESERVE MARGIN =,F8.2,2H %)') & (PMXO-PDO)/PDO*100. WRITE(*,*) ' ' WRITE(*,*) ' BIG GRID 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,*) ((ATC(I,J),I=1,2),J=1,10) IF(ATC(1,1).LT.0.) THEN WRITE(*,*)'THE FIRST ATC MW MUST NOT BE NEGATIVE' ! THIS IS BECAUSE MC NEEDS THE FIRST VALUE TO BE THE N-0 IMPORT STATE GOTO 20 ENDIF IF(ATC(2,1).LE.0.) THEN WRITE(*,*)'THE FIRST ATC PROBABILITY MUST BE > 0' ! THIS IS TO KEEP THE MC SOLUTION FROM BLOWING UP GOTO 20 ENDIF 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. CALL CPU_TIME(T0) ! reset the clock 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 RRT=.125 ! HARD WIRE IN THE RATE OF REPAIR FOR TRANSMISSION = .125 i.e. MTTR=8 HOURS FROM ALL LOWER ATC STATES BACK TO THE FIRST N-0 STATE 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 ALLR=0.D0 ! ANNUAL LOAD LOSS RISK, WILL COUNT THE NUMBER OF YEARS THERE ARE NO LOSS OF LOAD EVENTS ALLR1=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 M=3 ! SEED FOR RANDOM NUMBER GENERATOR 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.D0 WRITE(3,*) ' ' WRITE(3,*) ' ' WRITE(3,*) ' MONTECARLO SEQUENTIAL SOLUTION' WRITE(3,*) ' ' WRITE(3,*) ' YRS LOLH LOLE ALOLP MWHEUE' WRITE(*,*) ' ' 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 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) 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 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)! SUM UP THE POWER OF THE GENS ON LINE 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) ALLR=ALLR+1.D0 ! COUNT THE NUMBER OF YEARS THERE WERE NO EVENTS IF(NED1.EQ.0) ALLR1=ALLR1+1.D0 ! COUNT THE NUMBER OF YEARS THERE WERE NO EVENTS IN AREA 1 IF(YRS.LT.MULT.AND.YRS.LT.ITERMAX) 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(3,'(A8,3F11.6,F10.3)') ' AREA 1', & LOLH1/YRS,LOLE1/YRS,1.-ALLR1/YRS,EUE1/YRS WRITE(*,'(F8.0,3F11.6,F10.3)') & YRS,LOLH/YRS,LOLE/YRS,1.-ALLR/YRS,EUE/YRS WRITE(*,'(A8,3F11.6,F10.3)') ' AREA 1', & LOLH1/YRS,LOLE1/YRS,1.-ALLR1/YRS,EUE1/YRS MULT=MULT*10 IF(YRS.LT.ITERMAX) GOTO 6 CALL CPU_TIME(T1) WRITE(3,'(/'' TIME ='',F10.2,'' MIN'')') (T1-T0)/60. WRITE(*,'(/'' TIME ='',F10.2,'' MIN'')') (T1-T0)/60. WRITE(3,*) ' ' WRITE(3,*) 'ATCS: LOLTV=LOSS OF LOAD TRANSMISSION EVENTS' DO 31 I=1,11 IF(ATCH(I).GT.0.D0) & WRITE(3,'(8H ATC MW=,F6.0,4X,5HPROB=,F7.4,4X,6HLOLTV=,F8.0,4X, & 5HPROB=,F7.4,11H CALCULATED)') ATC(1,I),ATC(2,I),ATCC(I),ATCH(I)/ & (YRS*8760.) 31 CONTINUE WRITE(*,*) 'FILE OP CONTAINS THE OUTPUT REPORT' WRITE(*,*) 'PRESS ENTER TO END THIS PROGRAM' READ(*,'(A1)') A END