$eolcom // set n 'node' /r1*r3,lp1*lp41/; set res(n) 'resource' /r1*r3/; set lp(n) 'load point' /lp1*lp41/; set zn 'zone' /zs1*zs3,z1*z9/; set znL(zn) /z1*z9/; set lps(lp) /lp1*lp10,lp12*lp38,lp40,lp41/; alias (n,n1,n2) ; set line(n1,n2) 'directed graph' / r1.lp1 lp1.lp2 lp2.lp3 lp2.lp4 lp3.lp15 lp4.lp5 lp5.lp6 lp5.lp7 lp7.lp8 lp7.lp9 lp8.lp22 lp9.lp10 lp10.lp11 lp10.lp12 lp11.lp25 lp12.lp13 lp13.lp29 r2.lp14 lp14.lp15 lp14.lp16 lp14.lp17 lp16.lp31 lp17.lp18 lp17.lp19 lp18.lp33 lp19.lp20 lp20.lp21 lp20.lp22 lp20.lp23 lp21.lp36 lp23.lp24 lp24.lp25 lp24.lp26 lp24.lp27 lp26.lp39 lp27.lp28 lp27.lp29 lp28.lp41 r3.lp30 lp30.lp31 lp30.lp32 lp32.lp33 lp32.lp34 lp34.lp35 lp35.lp36 lp35.lp37 lp37.lp38 lp38.lp39 lp38.lp40 lp40.lp41 /; set lines(n1,n2) 'directed lines with switches' / r1.lp1 lp3.lp15 lp4.lp5 lp7.lp9 lp8.lp22 lp11.lp25 lp13.lp29 r2.lp14 lp16.lp31 lp17.lp19 lp18.lp33 lp20.lp23 lp21.lp36 lp26.lp39 lp28.lp41 r3.lp30 lp32.lp34 lp35.lp37 /; parameter ls(n1,n2) 'lines without switches =1; with switches =2' / r1.lp1 2 lp1.lp2 1 lp2.lp3 1 lp2.lp4 1 lp3.lp15 2 lp4.lp5 2 lp5.lp6 1 lp5.lp7 1 lp7.lp8 1 lp7.lp9 2 lp8.lp22 2 lp9.lp10 1 lp10.lp11 1 lp10.lp12 1 lp11.lp25 2 lp12.lp13 1 lp13.lp29 2 r2.lp14 2 lp14.lp15 1 lp14.lp16 1 lp14.lp17 1 lp16.lp31 2 lp17.lp18 1 lp17.lp19 2 lp18.lp33 2 lp19.lp20 1 lp20.lp21 1 lp20.lp22 1 lp20.lp23 2 lp21.lp36 2 lp23.lp24 1 lp24.lp25 1 lp24.lp26 1 lp24.lp27 1 lp26.lp39 2 lp27.lp28 1 lp27.lp29 1 lp28.lp41 2 r3.lp30 2 lp30.lp31 1 lp30.lp32 1 lp32.lp33 1 lp32.lp34 2 lp34.lp35 1 lp35.lp36 1 lp35.lp37 2 lp37.lp38 1 lp38.lp39 1 lp38.lp40 1 lp40.lp41 1 /; parameter Yini(n1,n2) 'Initial status of lines with switches' / r1.lp1 1 lp3.lp15 0 lp4.lp5 1 lp7.lp9 1 lp8.lp22 0 lp11.lp25 0 lp13.lp29 0 r2.lp14 1 lp16.lp31 0 lp17.lp19 1 lp18.lp33 0 lp20.lp23 1 lp21.lp36 0 lp26.lp39 0 lp28.lp41 0 r3.lp30 1 lp32.lp34 1 lp35.lp37 1 /; parameter z2lp(zn,n) 'relation between zone and load point' / zs1.r1 1 zs2.r2 1 zs3.r3 1 z1.lp1 1 z1.lp2 1 z1.lp3 1 z1.lp4 1 z2.lp5 1 z2.lp6 1 z2.lp7 1 z2.lp8 1 z3.lp9 1 z3.lp10 1 z3.lp11 1 z3.lp12 1 z3.lp13 1 z4.lp14 1 z4.lp15 1 z4.lp16 1 z4.lp17 1 z4.lp18 1 z5.lp19 1 z5.lp20 1 z5.lp21 1 z5.lp22 1 z6.lp23 1 z6.lp24 1 z6.lp25 1 z6.lp26 1 z6.lp27 1 z6.lp28 1 z6.lp29 1 z7.lp30 1 z7.lp31 1 z7.lp32 1 z7.lp33 1 z8.lp34 1 z8.lp35 1 z8.lp36 1 z9.lp37 1 z9.lp38 1 z9.lp39 1 z9.lp40 1 z9.lp41 1 / ; parameter cg(n),cz(znL),ci(n1,n2),cib(n1,n2),cb(n1,n2),cr(n1); cg(n)=0.1; cg(res)=0.2; cz(znL)=50; ci(n1,n2)=0.1; cib(n1,n2)=0.1; cb(n1,n2)=5; cr(n1)=10; parameter P_d(n1),Q_d(n1); P_d(n1)=0; P_d('lp2')=0.1; P_d('lp5')=0.1; P_d('lp8')=0.1; P_d('lp10')=0.1; P_d('lp12')=0.1; P_d('lp14')=0.1; P_d('lp20')=0.1; P_d('lp24')=0.1; P_d('lp27')=0.1; P_d('lp30')=0.1; P_d('lp35')=0.2; P_d('lp40')=0.2; P_d('lp41')=0.2; Q_d(n1)=0; Q_d('lp2')=0.02; Q_d('lp5')=0.02; Q_d('lp8')=0.02; Q_d('lp10')=0.02; Q_d('lp12')=0.02; Q_d('lp14')=0.02; Q_d('lp20')=0.02; Q_d('lp24')=0.02; Q_d('lp27')=0.02; Q_d('lp30')=0.02; Q_d('lp35')=0.04; Q_d('lp40')=0.04; Q_d('lp41')=0.04; scalar R,X,M,V_snorm,V_smin,V_smax,I_smax,P_lmax,Q_lmax; scalar Base,Vn,In,Rn,Xn 'base values of voltage,current,resistance,reactance'; Base=1; Vn=7.967; In=Base/Vn; Rn=sqr(Vn)/Base; Xn=sqr(Vn)/Base; P_d(lp)=P_d(lp)/Base; Q_d(lp)=Q_d(lp)/Base; V_snorm=sqr(1); V_smin=sqr(0.95); I_smax=2; P_lmax=2; Q_lmax=2; R=0.9/Rn; X=0.9/Xn; M=1; Variable P_l(n1,n2),Q_l(n1,n2),b(n1,n2); Positive variable V_s(n1),I_s(n1,n2),Q_g(n),sw(n1,n2),rlp(n1); Free variable Z_obj,P_g(n); Binary variable X_z(zn),Y_l(n1,n2),Belta(n1,n2); Equation obj 'objective is to minimize energy cost and the number of switching action, maximizing energized zones', ActivePowerBal 'active power balance for each node', VarPowerBal 'reactive power baalnce for each node', Loadshed 'load shedding for each node', VolDropLine 'voltage drop for lines without switches', VolDropSw 'voltage drop for lines with switches', CalVoltage 'calculate voltage for each node', VolLimLo 'lower limit of voltage', VolLimUp 'upper limit of voltage', CurrentLimL 'lower limit of current for lines without switches', CurrentLimU 'upper limit of current for lines without switches', CurrentLimLs 'lower limit of current for lines with switches', CurrentLimUs 'upper limit of current for lines with switches', ActivePowLimL 'lower limit of active power for lines without switches', ActivePowLimU 'upper limit of active power for lines without switches', ActivePowLimLs 'lower limit of active power for lines with switches', ActivePowLimUs 'upper limit of active power for lines with switches', VarPowLimL 'lower limit of reactive power for lines without switches' VarPowLimU 'uppper limit of reactive power for lines without switches' VarPowLimLs 'lower limit of reactive power for lines with switches' VarPowLimUs 'uppper limit of reactive power for lines with switches' SwAction1 'calculate switching action', SwAction2 'calculate switching action', RadialNode 'radial constraint for node', RadialLine 'radial constraint for line', auxillary1, auxillary2; obj .. Z_obj=e=sum(n,cg(n)*P_g(n))-sum(znL,cz(znL)*X_z(znL))+sum((n1,n2)$(ls(n1,n2)=1),ci(n1,n2)*I_s(n1,n2)*R)+sum((n1,n2)$(ls(n1,n2)=2),cib(n1,n2)*I_s(n1,n2)) +sum(lines(n1,n2),cb(n1,n2)*sw(n1,n2))+sum(n1,cr(n1)*rlp(n1)*P_d(n1)); ActivePowerBal(n1,zn)$(z2lp(zn,n1)=1) .. P_g(n1)+sum(n2$(ls(n2,n1)>0),P_l(n2,n1))-sum(n2$(ls(n1,n2)=1),P_l(n1,n2)+I_s(n1,n2)*R)-sum(n2$(ls(n1,n2)=2),P_l(n1,n2))=e=P_d(n1)*(X_z(zn)-rlp(n1)); VarPowerBal(n1,zn)$(z2lp(zn,n1)=1) .. Q_g(n1)+sum(n2$(ls(n2,n1)>0),Q_l(n2,n1))-sum(n2$(ls(n1,n2)=1),Q_l(n1,n2)+I_s(n1,n2)*X)-sum(n2$(ls(n1,n2)=2),Q_l(n1,n2))=e=Q_d(n1)*(X_z(zn)-rlp(n1)); Loadshed(lp,zn)$(z2lp(zn,lp)=1) .. rlp(lp)=l=X_z(zn); VolDropLine(n1,n2)$(ls(n1,n2)=1) .. V_s(n1)-V_s(n2)=e=2*(P_l(n1,n2)*R+Q_l(n1,n2)*X)+I_s(n1,n2)*(sqr(R)+sqr(X)); VolDropSw(n1,n2)$(ls(n1,n2)=2) .. V_s(n1)-V_s(n2)=e=b(n1,n2); CalVoltage(n1,n2)$(ls(n1,n2)>0).. V_s(n2)*I_s(n1,n2)=g=sqr(P_l(n1,n2))+sqr(Q_l(n1,n2)); VolLimLo(n1,zn)$(z2lp(zn,n1)=1) .. V_s(n1)=g=V_smin*X_z(zn); VolLimUp(n1,zn)$(z2lp(zn,n1)=1) .. V_s(n1)=l=V_snorm*X_z(zn); CurrentLimL(n1,n2)$(ls(n1,n2)=1) .. I_s(n1,n2)=g=-I_smax; CurrentLimU(n1,n2)$(ls(n1,n2)=1) .. I_s(n1,n2)=l=I_smax; CurrentLimLs(n1,n2)$(ls(n1,n2)=2) .. I_s(n1,n2)=g=-I_smax*Y_l(n1,n2); CurrentLimUs(n1,n2)$(ls(n1,n2)=2) .. I_s(n1,n2)=l=I_smax*Y_l(n1,n2); ActivePowLimL(n1,n2)$(ls(n1,n2)=1) .. P_l(n1,n2)=g=-P_lmax; ActivePowLimU(n1,n2)$(ls(n1,n2)=1) .. P_l(n1,n2)=l=P_lmax; ActivePowLimLs(n1,n2)$(ls(n1,n2)=2) .. P_l(n1,n2)=g=-P_lmax*Y_l(n1,n2); ActivePowLimUs(n1,n2)$(ls(n1,n2)=2) .. P_l(n1,n2)=l=P_lmax*Y_l(n1,n2); VarPowLimL(n1,n2)$(ls(n1,n2)=1) .. Q_l(n1,n2)=g=-Q_lmax; VarPowLimU(n1,n2)$(ls(n1,n2)=1) .. Q_l(n1,n2)=l=Q_lmax; VarPowLimLs(n1,n2)$(ls(n1,n2)=2) .. Q_l(n1,n2)=g=-Q_lmax*Y_l(n1,n2); VarPowLimUs(n1,n2)$(ls(n1,n2)=2) .. Q_l(n1,n2)=l=Q_lmax*Y_l(n1,n2); SwAction1(n1,n2)$(ls(n1,n2)=2) .. sw(n1,n2)=g=Y_l(n1,n2)-Yini(n1,n2); SwAction2(n1,n2)$(ls(n1,n2)=2) .. sw(n1,n2)=g=Yini(n1,n2)-Y_l(n1,n2); RadialNode(znL) .. sum(lines(n1,n2)$(z2lp(znL,n1)=1),Belta(n1,n2))+sum(lines(n2,n1)$(z2lp(znL,n1)=1),Belta(n1,n2))=e=X_z(znL); RadialLine(n1,n2)$(ls(n1,n2)=2) .. Belta(n1,n2)+Belta(n2,n1)=e=Y_l(n1,n2); auxillary1(n1,n2)$(ls(n1,n2)=2) .. b(n1,n2)=l=100*(1-Y_l(n1,n2)); auxillary2(n1,n2)$(ls(n1,n2)=2) .. b(n1,n2)=g=-100*(1-Y_l(n1,n2)) option limrow = 15; Model MasterProblem / obj, ActivePowerBal, VarPowerBal, Loadshed, VolDropLine, VolDropSw, CalVoltage, VolLimLo, VolLimUp, CurrentLimL, CurrentLimU, CurrentLimLs, CurrentLimUs, $ontext ActivePowLimL, ActivePowLimU, ActivePowLimLs, ActivePowLimUs, VarPowLimL, VarPowLimU, VarPowLimLs, VarPowLimUs, $offtext SwAction1, SwAction2, RadialNode, RadialLine, auxillary1, auxillary2 / ; option miqcp=cplex; V_s.fx(res)=V_snorm; Belta.fx('r1','lp1')=0; Belta.fx('r2','lp14')=0; Belta.fx('r3','lp30')=0; Q_g.fx(lp)=0; P_g.fx(lps)=0; P_g.up('lp39')=0.00; P_g.up('lp11')=0.10; option decimals=8; X_z.fx('z7')=0; I_s.up('lp21','lp36')=0.10; Option optca=0; Option optcr=0; solve MasterProblem using miqcp minimizing Z_obj; Display P_l.l,Q_l.l,V_s.l,I_s.l,P_g.l,Q_g.l,sw.l,rlp.l,X_z.l,Y_l.l,Belta.l,b.l;