*Optimal power flow for 9 bus sets i number of bus /1*9/ r Reference bus/4/ g(i) number of generator /1*6/ dg(i) dispatchable DG units/1*4/ sg(i) stochastic DG units/5,6/ ng(i) number of buses with no generator /7*9/ k Hourly time interval/1*24/ ; alias (i,j); Table BData(i,j,*) Branch data R X 1.4 0 0.0576 4.5 0.017 0.092 5.6 0.039 0.17 3.6 0 0.0586 6.7 0.0119 0.1008 7.8 0.0085 0.072 8.2 0 0.0625 8.9 0.032 0.161 9.4 0.01 0.085 ; * Generating G and B BData(i,j,"R")=BData(j,i,"R"); BData(i,j,"X")=BData(j,i,"X"); *Generalized Branch data BData(j,i,"R")$(BData(i,j,"R") gt 0) = BData(i,j,"R"); BData(j,i,"X")$(BData(i,j,"X") gt 0) = BData(i,j,"X"); *construction of admitance matrix: Parameter Z(i,j), GG(i,j), BB(i,j); Z(i,j) = (BData(i,j,"R"))**2 + (BData(i,j,"X"))**2 ; GG(i,j)$(Z(i,j) ne 0.00) = BData(i,j,"R")/Z(i,j) ; BB(i,j)$(Z(i,j) ne 0.00) = -BData(i,j,"X")/Z(i,j); BB(j,i)$(Z(i,j) ne 0.00) = -BData(i,j,"X")/Z(i,j); Parameter Gbus(i,j),Bbus(i,j),Ybus(i,j); Bbus(i,i) = sum(j,BB(i,j)); Gbus(i,i) = sum(j,GG(i,j)); Gbus(i,j)$(ord(i) ne ord(j)) = -GG(i,j); Bbus(i,j)$(ord(i) ne ord(j)) = -BB(i,j); display Gbus,Bbus; parameter Gbus(i,j),Bbus(i,j); Table BusData(i,*) bus data Vmin Vmax PLoad QLoad 1 0.9 1.1 0 0 2 0.9 1.1 0 0 3 0.9 1.1 0 0 4 0.9 1.1 0 0 5 0.9 1.1 0.90 0.30 6 0.9 1.1 0 0 7 0.9 1.1 0 0.35 8 0.9 1.1 1.00 0 9 0.9 1.1 1.25 0 ; Table Ddata(dg,*) Pmin Pmax Qmax Qmin 1 0 4 3 -3 2 0 5 3 -3 3 0 5.5 3 -3 4 0 7 3 -3 Table CData(i,*) Cost_dg SUC_dg SDC_dg 1 37 20 25 2 40 20 25 3 35 50 25 4 45 50 25 ; Table Sdata(sg,*) Characteristics of stochastic units included in VPP portfolio Pmin Pmax Cost_SG 5 0 9 55 6 0 7 65 ; Table NData(i,*) Cost_sg 5 55 6 65 ; *Table CData(i,*) * C2 C1 C0 * 1 0.11 5 450 * 2 0.085 1.2 600 * 3 0.1225 1 335 * 4 0.1356 1.30 440 * ; Table Price_DA (k,*) Characteristics of dispatchable loads and forecasted market prices for 24 hours F_price Cost_DL P_DLmin P_DLmax 1 60 50 0 0.575 2 40 30 0 0.540 3 42.2 35.5 0 0.549 4 43 41 0 0.551 5 50 50 0 0.575 6 60 65 0 0.610 7 64.5 68 0 0.626 8 70 75 0 0.645 9 75 78 0 0.636 10 76 76 0 0.666 11 65 65 0 0.628 12 80.5 85 0 0.682 13 85 87 0 0.698 14 75 80 0 0.663 15 65 70 0 0.628 16 55 65 0 0.593 17 57 65 0 0.600 18 56 80 0 0.596 19 72.5 80 0 0.654 20 80 85 0 0.680 21 87 89 0 0.705 22 72 75 0 0.652 23 52.5 65 0 0.584 24 60 65 0 0.610 ; Parameters Charge_DMP(i,K), alpha4 ; Scalar a1 is the first model estimation parameter for distribution network demand /0.07/; Scalar a2 is the second model estimation parameter for distribution network demand /8/; *Scalar BM is a positive large number exceeding any max feasible value /1000/; *Scalar HD Bilateral hourly tolerance (in *5%) /0.5/; *alpha1=0.95 ; *alpha2=1.05 ; *alpha3=1.00 ; alpha4=1.00 ; *LMP('5',i,k)= alpha1 * price_DA(K,i) ; *LMP('7',i,K)= alpha2 * price_DA(K,i) ; *LMP('9',i,K)= alpha3 * price_DA(K,i) ; Charge_DMP(i,K) = alpha4*price_DA(K,i) ; *Pdemand(i,k)= a1*price_DA(k,i)+a2 ; variables profit cost revenue Pdg(i) active power output for generator g in bus i Psg(i) Stochastic units output Qg(i) reactive power output for generator in bus i Ploss(i) power losses in bus i Pdl(k) Dispatchable load V(i) voltage magnitude at bus i Teta(i) voltage angle at bus i Pld(i) Active power demand in bus i Qld(i) Reactive power demand in bus i Cost_dg(i) Generation cost of dispatchable units SUC_dg(i) Start-up cost of dg units SDC_dg(i) Shut-down cost of dg units Cost_sg(i) Generation cost of stochastic units *CC2(i) *CC1(i) *CC0(i) Binary variables alpha_dg(k) On and off status of dispatchable DG units beta_dg(k) Start-up status of dispatchableDG units gamma_dg(k) Shut-down status of DG units alpha_sg(k) On and off status of stochastic units ; Teta.fx("1")=0; Pdg.lo(dg)=Ddata(dg,"Pmin"); Pdg.up(dg)=Ddata(dg,"Pmax"); Psg.lo(sg)=Sdata(sg,"Pmin"); Psg.up(sg)=Sdata(sg,"Pmax"); PDL.lo(k)=Price_DA (k,"P_DLmin"); PDL.up(k)=Price_DA (k,"P_DLmax"); Qg.up(dg )=1.00*Ddata(dg,'Qmax'); Qg.lo(dg )=1.00*Ddata(dg,'Qmin'); Pdg.fx(ng)=0; Qg.fx(ng)=0; V.lo(i)=BusData(i,"Vmin"); V.up(i)=BusData(i,"Vmax"); Pld.fx(i)=BusData(i,"PLoad"); Qld.fx(i)=BusData(i,"QLoad"); Cost_dg.fx(i)=CData(i,"Cost_dg"); SUC_dg.fx(i)=CData(i,"SUC_dg"); SDC_dg.fx(i)=CData(i,"SDC_dg"); Cost_sg.fx(i)=Ndata(i,"Cost_sg"); *CC2.fx(i)=CData(i,"C2"); *CC1.fx(i)=CData(i,"C1"); *CC0.fx(i)=CData(i,"C0"); equation eq1 eq2 eq3 eq4 eq5 P_bus(i) active power balance constraints Q_bus(i) reactive power balance constraints; eq1.. profit =e= (revenue)-(cost); eq2(k).. Revenue =E= Sum((i),Pld(i) + Ploss(i)*price_DA(k,"F_PRICE")) ; *eq3(i).. cost=e=sum((i),CC2(i)*100*100*Pdg(i)**2 + 100*CC1(i)*Pdg(i) + CC0(i)); eq3(i,k).. Cost =E= Sum((dg),Pdg(i)*Cdata(i,"cost_dg")*alpha_dg(k)) + (Cdata(i,"SUC_dg")*beta_dg(k))+(Cdata(i,"SDC_dg")*gamma_dg(k)) + Sum((sg),Psg(i)*Ndata(i,"Cost_sg")*alpha_sg(k))+(Price_DA(k,"Cost_DL")*PDL(k)) ; eq4(i,k).. Sum((dg),Pdg(i) + PDL(k)) =E= Pld(i) + Ploss(i) ; eq5(k).. PDL(k)=L= Price_DA(k,"P_DLmax"); P_Bus(i) .. (Pdg(i)+psg(i))-(Pld(i))=e= sum(j,V(i )*(V(j ))*(Gbus(i,j)*cos(Teta(i )-Teta(j ))+Bbus(i,j)*sin(Teta(i )-Teta(j )))); Q_Bus(i ) .. Qg(i)-(Qld(i))=e= sum(j,V(i )*(V(j ))*(Gbus(i,j)*sin(Teta(i )-Teta(j ))-Bbus(i,j)*cos(Teta(i )-Teta(j )))); model OPF /all/; solve OPF using minlp maximizing Profit; display profit.l,revenue.l,cost.l,pdg.l,qg.l,pld.l,qld.l,v.l,teta.l,bbus,gbus,P_bus.l,q_bus.l;