$Ontext This model has been coded April to ... 2012 by A.J $Offtext Sets i "number of buses" /1*14/ g "state number" /1*5/ Gen "index of generators" /G1/ Map(Gen,i) "associates generators with buses" /G1.1/ year the horizon year for the planning program /1*3/ ; Alias(i,j); Alias(g,gp); *make a random number for each bus to show the growth rate of the consumption in each bus for 30 yers horizon parameter GR(i,year); execseed = 1e8*(frac(jnow)); loop(i, loop(year, GR(i,year) = uniform(0,0.05); ); ); *candidate buses for allocation of DG's in this system Sets Candidate(i) "Candidate buses for DG placement" /2,3,5,7,9,11,14/; * CandidateP(i) "Candidate buses for DG placement" /1/ ; *SET iter(g); Table G_Data(Gen,*) "generators data" Pmin PMax Qmin QMax * (KW) (KW) (KVAr) (KVAr) G1 0.0 inf -inf inf ; *the parametrs of the system is p.u Table Line_Data(i,j,*) "lines data" Y PHI G * (Ohm^-1) (rad) (Ohm^-1) 1.1 20.3590 -1.2704 6.0250 1.2 16.0609 1.8873 -4.9991 1.5 4.3575 1.8085 -1.0259 2.2 31.7342 -1.2661 9.5213 2.3 4.9147 1.8038 -1.1350 2.4 5.3865 1.8892 -1.6860 2.5 5.4654 1.8873 -1.7011 3.3 10.3063 -1.2631 3.1210 3.4 5.4440 1.9442 -1.9860 4.4 40.0583 -1.3052 10.5130 4.5 22.6370 1.8778 -6.8410 4.7 4.8895 1.5708 0 4.9 1.8555 1.5708 0 5.5 36.7993 -1.3078 9.5680 5.6 4.2574 1.5708 0 6.6 18.5471 -1.2081 6.5799 6.11 4.5369 2.0163 -1.9550 6.12 3.5235 2.0187 -1.5260 6.13 6.8445 2.0407 -3.0989 7.7 19.5490 -1.5708 0 7.8 5.6770 1.5708 0 7.9 9.0901 1.5708 0 8.8 5.6770 -1.5708 0 9.9 24.6742 -1.3532 5.3261 9.10 11.0755 1.9308 -3.9020 9.14 3.3471 2.0103 -1.4240 10.10 15.8602 -1.1976 5.7829 10.11 4.7879 1.9745 -1.8809 11.11 9.3227 -1.1467 3.8359 12.12 6.7515 -0.9339 4.0150 12.13 3.3566 2.4062 -2.4890 13.13 12.6122 -1.0084 6.7249 13.14 2.5791 2.0273 -1.1370 14.14 5.9260 -1.1239 2.5610 ; * The set of buses N must be duplicated to refer to different elements in the same constraint Line_Data(i,j,'Y')$(Ord(i) gt Ord(j))=Line_Data(j,i,'Y'); Line_Data(i,j,'PHI')$(Ord(i) gt Ord(j))=Line_Data(j,i,'PHI'); Line_Data(i,j,'G')$(Ord(i) gt Ord(j))=Line_Data(j,i,'G'); *this dataes are for the base year Table Bus_Data(i,*) "loads power in each + voltage margin" VMax Vmin PL QL * (p.u) (p.u) (MW) (MVAr) 1 1.06 0.94 0 0 2 1.06 0.94 21.7 12.7 3 1.06 0.94 94.2 19 4 1.06 0.94 47.8 -3.9 5 1.06 0.94 7.6 1.6 6 1.06 0.94 11.2 7.5 7 1.06 0.94 0 0 8 1.06 0.94 0 0 9 1.06 0.94 29.5 16.6 10 1.06 0.94 9 5.8 11 1.06 0.94 3.5 1.8 12 1.06 0.94 6.1 1.6 13 1.06 0.94 13.5 5.8 14 1.06 0.94 14.9 5 ; *wind turbin capacity which the allocation in integer solution is based on it Scalars Wind_Turbine Kw /1100/ ; Table matrix_C(g,*) "Table 8" wind_power_state load_state 1 1 1 2 1 0.853 3 1 0.774 4 1 0.713 5 1 0.65 ; Table P_C(g,*) "Probability of each state(Table 9)" P_c 1 0.000761 2 0.004338 3 0.008348 4 0.012587 5 0.012587 ; * Constant PI is used to state limits on bus angles Scalar PI /3.1416/; * Optimization variables are declared Variables Z "the objective function value" P_loss(g,year) "total power loss during state g" V(g,i,year) "voltage of bus i in state g" Delta(g,i,year) "pgase of bus i in state g" Igij(g,i,j,year) "branch current from i to j" Pgij(g,i,j,year) P(g,Gen,year) active power output for generator G Q(g,Gen,year) reactive power output for generator G ; Positive variables P_DG(i) "P of DG with kind of n in bus i" Q_DG(i) "Q of DG with kind of n in bus i"; Integer variables N_DG(i) "number of DG kind of n at bus i" ; *base parameters which used in this program because our equations are p.u!!! Scalars Sbase "base power (KVA)" /100000/ Ibase "base current (A)" /6250/ Vbase "base voltage (V)" /16000/ Zbase "base ampedance ohm" /2.56/ CF "Capacity Factor of kind n DG" /0.2209/; *Scalars *Sbase "base power (KVA)" /100000000/ *Ibase "base current (A)" /434.78/ *Vbase "base voltage (KV)" /230000/ *Zbase "base ampedance ohm" /529/ *CF "Capacity Factor of kind n DG" /0.22/; * Limits are stated on variables using previously defined data P.lo(g,Gen,year)=G_Data(Gen,'Pmin'); P.up(g,Gen,year)=G_Data(Gen,'PMax'); Q.lo(g,Gen,year)=G_Data(Gen,'Qmin'); Q.up(g,Gen,year)=G_Data(Gen,'QMax'); V.lo(g,i,year)=Bus_Data(i,'Vmin'); V.up(g,i,year)=Bus_Data(i,'VMax'); Delta.lo(g,i,year)=-PI; Delta.up(g,i,year)=PI; N_DG.up(i)=10; *Igij.lo(g,i,j)=-300/Ibase; *Igij.up(g,i,j)=300/Ibase; * Bus 1 is considered as the reference bus V.fx(g,'1',year)=1.025; Delta.fx(g,'1',year)=0; *N_DG.fx('3')=0.25; V.l(g,i,year)=1; Delta.l(g,i,year)=0; *P.l(g,Gen)=7000; *Q.l(g,Gen)=2; *P_loss.l(g)= 0.5*Sum((i,j),(abs(Line_Data(i,j,'G')*Zbase))*(power(V.l(g,i),2) + power(V.l(g,j),2) * - 2*V.l(g,i)*V.l(g,j)* cos(Delta.l(g,j)-Delta.l(g,i)))) ;; *Igij.l(g,i,j)=(Line_Data(i,j,'Y')*Zbase)*sqrt(abs(power(V.l(g,i),2) + power(V.l(g,j),2) * - 2*V.l(g,i)*V.l(g,j)* cos(Delta.l(g,j)-Delta.l(g,i)))) ; *N_DG.l(i)=1; *P_DG.l(i)= N_DG.l(i)*Wind_Turbine; * * Constraints Declaration Equations Obj "objective function definition" PowerFlow_P(i,g,year) "active power flow equation for bus i in state g" PowerFlow_Q(i,g,year) "reactive power flow equation for bus i in state g" PowerLoss(g,year) "power loss equation" BranchCurrent(g,i,j,year) "branch current equation" BranchPower(g,i,j,year) "branch power equation" BranchCurrentLimit(g,i,j,year) "limitation of lines" Bus_Penetration(i) "Sum of Power of DG in bus i" Sys_Penetration "Penetration of DG in system" Discret_DG_Power(i) "Discret power for each kind of DG" * Tot_Balance_P(g) "balance of total P because we remove slack bus eq." * Tot_Balance_Q(g) "balance of total Q because we remove slack bus eq." ZeroDiscret_DG_Power(i) *vvvvvv(i) *rrr(i) *rr(i) ; * Constraints Definition * sum((g),P_loss(g)) *matrix_C(g,'load_state')* load_state *Objective Function *objective function is annual energy loss of the system Obj.. Z =e=sum((g,year),P_C(g,'P_c')*P_loss(g,year)*8760) ; *Cons. 1-1 P(G1) Rmoved! and don't use power flow eq. for bus 1 PowerFlow_P(i,g,year).. Sum(Gen$Map(Gen,i),P(g,Gen,year)/Sbase)+matrix_C(g,'wind_power_state')*P_DG(i)/Sbase -matrix_C(g,'load_state')*(Bus_Data(i,'PL')/Sbase)*(1+GR(i,year))=e= Sum((j),V(g,i,year)*V(g,j,year)*(Line_Data(i,j,'Y')*Zbase)*cos(Line_Data(i,j,'PHI')+Delta(g,j,year)-Delta(g,i,year))); *Cons. 1-2 Q(G1) Rmoved! and don't use power flow eq. for bus 1 PowerFlow_Q(i,g,year).. Sum(Gen$Map(Gen,i),Q(g,Gen,year)/Sbase)- matrix_C(g,'load_state')*(Bus_Data(i,'QL')/Sbase)*(1+GR(i,year))=e= -Sum((j),V(g,i,year)*V(g,j,year)*(Line_Data(i,j,'Y')*Zbase)*sin(Line_Data(i,j,'PHI')+Delta(g,j,year)-Delta(g,i,year))); *Cons 2 PowerLoss(g,year).. P_loss(g,year) =e= 0.5*Sum((i,j),(abs(Line_Data(i,j,'G')*Zbase))*(power(V(g,i,year),2) + power(V(g,j,year),2) - 2*V(g,i,year)*V(g,j,year)* cos(Delta(g,j,year)-Delta(g,i,year))))*Sbase ; *Cons 3 BranchCurrent(g,i,j,year).. Igij(g,i,j,year) =e= (Line_Data(i,j,'Y')*Zbase)*sqrt(abs(power(V(g,i,year),2) + power(V(g,j,year),2) - 2*V(g,i,year)*V(g,j,year)* cos(Delta(g,j,year)-Delta(g,i,year)))) ; *appendix constraint BranchPower(g,i,j,year).. Pgij(g,i,j,year) =e=V(g,i,year)*V(g,j,year)*Line_Data(i,j,'Y')*Zbase*cos(Line_Data(i,j,'PHI')+Delta(g,j,year)-Delta(g,i,year)) -power(V(g,i,year),2)*Line_Data(i,j,'Y')*Zbase*cos(Line_Data(i,j,'PHI')); *Cons 4-1 ok *V.fx(g,'1')=1.025; *Cons 4-2 ok *Delta.fx(g,'1')=0; *Cons 5 ok *V.lo(g,i)=Bus_Data(i,'Vmin'); *V.up(g,i)=Bus_Data(i,'VMax'); *Cons 6 ok the branch power limit is 40MW BranchCurrentLimit(g,i,j,year).. abs(Pgij(g,i,j,year)) =l= 0.4; *BranchCurrentLimit(g,i,j).. abs(Igij(g,i,j)) =l= 300/Ibase; *Cons 7 ok Discret_DG_Power(i)$Candidate(i).. P_DG(i) =e= N_DG(i)*Wind_Turbine; ZeroDiscret_DG_Power(i)$(not candidate(i)).. P_DG(i) =e= 0; *rrr(i)$(CandidateP(i))..P_DG(i) =e=0; *rr(i)$(CandidateP(i))..N_DG(i)=e=0; *Cons 8 0.2 assumed as max peneteration in each bus !!!!!!! Bus_Penetration(i).. P_DG(i) =l= 10000; *CONS 9 Sys_Penetration.. Sum((i),CF*P_DG(i)) =l= 0.3*Sum((i),Bus_Data(i,'PL')); *Appendix Cons.A1 *Tot_Balance_P(g).. Sum((Gen,i)$Map(Gen,i),P(g,Gen)/Sbase)+Sum((i),P_DG(i)/Sbase) =e= P_loss(g)/Sbase+Sum((i),matrix_C(g,'load_state')*Bus_Data(i,'PL')/Sbase); *Appendix Cons.A2 *Tot_Balance_Q(g).. Sum((Gen,i)$Map(Gen,i),Q(g,Gen)/Sbase) =e= Sum((i),matrix_C(g,'load_state')*Bus_Data(i,'QL')/Sbase); model AJ / all /; AJ.reslim = 12000000; AJ.iterlim = 10000000; OPTION work = 4000; AJ.workfactor=1.5; AJ.workspace = 4096; *AJ.optfile = 1 ; option minlp = couenne ; *option minlp = baron ; *option minlp=sbb *minlp=lindoglobal *option lstcrs=t *option rtmaxv = 1.e8 *Solve AJ using rminlp Minimizing Z ; *loop(gp, *iter(g) = no; *iter(gp) = yes; SOLVE AJ USING minlp MINIMIZING Z; *); execute_unload "continous+limit+integer.gdx" Z.l,P_loss.l,V.l,Delta.l,P.l,Q.l,P_DG.l,N_DG.l,PowerFlow_P.m,PowerFlow_Q.m execute 'gdxxrw.exe continous+limit+integer.gdx var=Z.l rng=Z!a1:b2' execute 'gdxxrw.exe continous+limit+integer.gdx var=P_loss.l rng=P_loss!a1:dq2' execute 'gdxxrw.exe continous+limit+integer.gdx var=V.l rng=V!a1:ap121' execute 'gdxxrw.exe continous+limit+integer.gdx var=Delta.l rng=Delta!a1:ap121' execute 'gdxxrw.exe continous+limit+integer.gdx var=P.l rng=P!a1:b121' execute 'gdxxrw.exe continous+limit+integer.gdx var=Q.l rng=Q!a1:b121' execute 'gdxxrw.exe continous+limit+integer.gdx var=P_DG.l rng=P_DG!a1:ap2' execute 'gdxxrw.exe continous+limit+integer.gdx var=N_DG.l rng=N_DG!a1:ap2' execute 'gdxxrw.exe continous+limit+integer.gdx equ=PowerFlow_P.m rng=PowerFlow_P!a1:dq42' execute 'gdxxrw.exe continous+limit+integer.gdx equ=PowerFlow_Q.m rng=PowerFlow_Q!a1:dq42'