* Unit Committment: Social Welfare Model * Written by H.Shavandi, M. Pirnia and J.D. Fuller * This code represents the numerical analysis parts of the paper * "Near equilibrium electricity prices under non-convexities" * Written by H.Shavandi, M. Pirnia and J.D. Fuller * Please, use proper referencing if using the code $eolcom # $inlinecom { } $offsymxref set tt time periods /0*36/; set t(tt) /1*24/; Set tone(t)/1/; Set i generators /GEN1*GEN12 /; Set ioff(i) generators that should be off for some hours starting hour 1 due to minimum down time /GEN5, GEN6/; Set iwoff(i) generators with no force to be off at starting hour of day /GEN1*GEN4, GEN7*GEN12/; Set n nodes/n1*n24/; Set dn(n) demand nodeds /n1*n10, n13*n16, n18*n20 /; Set tn(n) transfer nodes /n11, n12, n17, n21*n24/; Set l transmission lines /l1*l34/; alias (t,k); alias (t,kk); Parameter toff(ioff)/GEN5 1, GEN6 6/; Parameter Lambda(i,t) marginal cost of generator i for operation at each time period; Parameter LambdaRES(i,t) marginal cost of generator i being in reserve status at each time period; Parameter CSU(i,t) startup cost of generator i at period t; Parameter CSD(i,t) shutdown cost of generator i at period t; Parameter CRP(i,t) ramping cost of generator i at period t ; Parameter RU(i) maximum ramping up for generator i ; Parameter RD(i) maximum ramping down for generator i ; Parameter SU(i) maximum startup rate for generator i ; Parameter SD(i) maximum shutdown rate for generator i; Parameter Gmax(i) maximum power generation for generator i ; Parameter Gmin(i) minimum power generation for generator i ; Parameter Dmin(i) minimum down time for generator i; Parameter Umin(i) minimum up time for generator i ; Parameter IZ(i) initial z var for generator i; Parameter IG(i) initial g var for generator i ; Parameter Ares(i) The subset of reserve generators; Parameter A(n,t) benefit function parameter ; Parameter B(n,t) benefit function parameter ; Parameter Qmax(n,t) maximum demand of node n at time period t ; Parameter Qmin(n,t) minimum demand of node n at time period t ; Parameter Pmin(n,t) minimum demand of node n at time period t; Parameter Pmax(n,t) minimum demand of node n at time period t ; Parameter RES(t) reserve amount at time t ; Parameter Agen(i,n) binary value for maping each generator to each bus node ; Parameter Agrid(l,n) binary value that maps and node to each line ; Parameter Fmax(l) maximum possible flow of line l; Parameter Bsuscep(l) Susceptance of line l ; $onecho > IEEE24Bus9.rsp par=Lambda rng=Lambda!a1 rdim=1 cdim=1 par=LambdaRES rng=LambdaRES! rdim=1 cdim=1 par=CSU rng=CSU!a1 cdim=1 rdim=1 par=CSD rng=CSD!a1 rdim=1 cdim=1 par=CRP rng=RC!a1 rdim=1 cdim=1 par=RU rng=RU!a1 rdim=1 cdim=0 par=RD rng=RD!a1 rdim=1 cdim=0 par=SU rng=SU!a1 rdim=1 cdim=0 par=SD rng=SD!a1 rdim=1 cdim=0 par=Gmax rng=Gmax!a1 rdim=1 cdim=0 par=Gmin rng=Gmin!a1 rdim=1 cdim=0 par=Dmin rng=DT!a1 rdim=1 cdim=0 par=Umin rng=UP!a1 rdim=1 cdim=0 par=IZ rng=IZ!a1 rdim=1 cdim=0 par=IG rng=IG!a1 rdim=1 cdim=0 par=Ares rng=Ares!a1 rdim=1 cdim=0 par=A rng=A!a1 rdim=1 cdim=1 par=B rng=B!a1 rdim=1 cdim=1 par=Qmax rng=Qmax!a1 rdim=1 cdim=1 par=Qmin rng=Qmin!a1 rdim=1 cdim=1 par=Pmin rng=Pmin!a1 rdim=1 cdim=1 par=Pmax rng=Pmax!a1 rdim=1 cdim=1 par=RES rng=RES!a1 rdim=1 cdim=0 par=Agen rng=Agen!a1 rdim=1 cdim=1 par=Agrid rng=Agrid!a1 rdim=1 cdim=1 par=Fmax rng=Fmax!a1 rdim=1 cdim=0 par=Bsuscep rng=BL!a1 rdim=1 cdim=0 $offecho $call 'gdxxrw input=IEEE24Bus9.xlsx o=IEEE24Bus9.gdx log=IEEE24Bus9.log@IEEE24Bus9.rsp'; $gdxin 'IEEE24Bus9.gdx' $load Lambda $load LambdaRES $load CSU $load CSD $load CRP $load RU $load RD $load SU $load SD $load Gmax $load Gmin $load Dmin $load Umin $load IZ $load IG $load Ares $load A $load B $load Qmax $load Qmin $load Pmin $load Pmax $load RES $load Agen $load Agrid $load Fmax $load Bsuscep $gdxin ; scalar m/0.5/; scalar mi/100/; Scalar pi/0.8/; scalar phi/0.05/; scalar h/1/; Fmax(l)=0.85*Fmax(l); LambdaRES(i,t)=0.4*Lambda(i,t); Variables objectiveFun objective function f(l,t) Power flow on line l at period t teta(n,t) State variable of node n at time period t ; Positive variables g(i,t) generation output of generator i at period t gr(i,t) Reserve capacity of generator i at period t r(i,tt) ramping amount of generator i at period t q(dn,t) demand of node n at period t v(i,t) startup status of generator i at time period t w(i,t) shutdown status of generator i at time period t Binary variables z(i,t) operation status of generator i at time period t ; Equations objfun Objective function ramp1(i,t) First constraint of ramping for generator i at period t (2b) ramp2(i,t) Second constraint of ramping for generator i at period t (2c) ramp3(i,t) Third constraint of ramping for generator i at period t (2d) ramp4(i,t) Forth constraint of ramping for generator i at period t (2e) ramp1tone(i,tone) First constraint of ramping for generator i at period t (2b)**separated for t=1 for coding purpose because of initial status** ramp2tone(i,tone) Second constraint of ramping for generator i at period t (2c)**separated for t=1 for coding purpose because of initial status** ramp3tone(i,tone) Third constraint of ramping for generator i at period t (2d)**separated for t=1 for coding purpose because of initial status** ramp4tone(i,tone) Forth constraint of ramping for generator i at period t (2e)**separated for t=1 for coding purpose because of initial status** MINgen(i,t) Minimum generation level for generator i at period t (2f) MAXgen(i,t) Maximum generation level for generator i at period t (2g) downtime(i,t) Minimum down time constraint for generator i at period t (2h) uptime(i,t) Minimum up time constraint for generator i at period t (2i) operation(i,t) Operation status for generator i at period t (2j) mktclr1(dn,t) Market clearing constraint for demand nodes at period t (4c) mktclr2(tn,t) Market clearing constraint for transfer nodes at period t (4c) Reserve(t) Reserve constraint at period t (4b) MinDemand(dn,t) Minimum demand constraint (1d) operation1(i,t) The extra Generator's operational constraint (2q) operationtone(i,tone) Operation status for generator i at period 1 (2j) operation1tone(i,tone) The extra Generator's operational constraint at period 1 (2q) transGRID(t) Constraint of Transmission network operator (3g) transMAX(l,t) Constraint of Transmission network operator (3b) transMIN(l,t) Constraint of Transmission network operator (3c) transPHYC2(l,t) Constraint of Transmission network operator (3d) transPHYC3(n,t) Constraint of Transmission network operator (3e) transPHYC4(n,t) Constraint of Transmission network operator (3e) transPHYC5(n,t) Constraint of Transmission network operator (3f) iniOFF(ioff,t) Generators that are down at initial status (2l) ; objfun.. objectiveFun =e= sum(i, sum(t, Lambda(i,t)*g(i,t)+LambdaRES(i,t)*Ares(i)*gr(i,t)+ CSU(i,t)*v(i,t)+CSD(i,t)*w(i,t)+CRP(i,t)*r(i,t)))-sum(dn, sum(t, A(dn,t)*q(dn,t)-m*B(dn,t)*power(q(dn,t),2))); ramp1(i,t)$(ord(t) ge 2).. -g(i,t)+g(i,t-1)+r(i,t)=g=0; ramp2(i,t)$(ord(t) ge 2).. -g(i,t-1)+g(i,t)+r(i,t)=g=0; ramp3(i,t)$(ord(t) ge 2).. g(i,t)-g(i,t-1)-RU(i)*z(i,t)-(SU(i)-RU(i))*v(i,t)=l=0; ramp4(i,t)$(ord(t) ge 2).. g(i,t-1)-g(i,t)-RD(i)*z(i,t-1)-(SD(i)-RD(i))*w(i,t)=l=0; ramp1tone(i,tone).. g(i,tone)-IG(i)-r(i,tone)=l=0; ramp2tone(i,tone).. IG(i)-g(i,tone)-r(i,tone)=l=0; ramp3tone(i,tone).. g(i,tone)-IG(i)-RU(i)*z(i,tone)-(SU(i)-RU(i))*v(i,tone)=l=0; ramp4tone(i,tone).. IG(i)-g(i,tone)-RD(i)*IZ(i)-(SD(i)-RD(i))*w(i,tone)=l=0; MINgen(i,t).. Gmin(i)*z(i,t)-g(i,t)=l=0; MAXgen(i,t).. g(i,t)+Ares(i)*gr(i,t)-Gmax(i)*z(i,t)=l=0; downtime(i,t)$(Dmin(i)>0).. Dmin(i)*w(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk) <= ord(t)+ Dmin(i)-1)), 1-z(i,kk))=l=0; uptime(i,t)$(Umin(i)>0).. Umin(i)*v(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk)<= ord(t)+ Umin(i)-1)), z(i,kk))=l=0; operation(i,t)$(ord(t)>=2).. z(i,t-1)-z(i,t)+v(i,t)-w(i,t)=e=0; operation1(i,t)$(ord(t)>=2).. z(i,t-1)+v(i,t)=l=1; operationtone(i,tone).. IZ(i)-z(i,tone)+v(i,tone)-w(i,tone)=e=0; operation1tone(i,tone).. IZ(i)+v(i,tone)=l=1; mktclr1(dn,t).. Sum(i, Agen(i,dn)*g(i,t))-q(dn,t)-sum(l, Agrid(l,dn)*f(l,t))=e=0; mktclr2(tn,t).. Sum(i, Agen(i,tn)*g(i,t))-sum(l, Agrid(l,tn)*f(l,t))=e=0; Reserve(t).. Sum(i, Ares(i)*gr(i,t))-sum(dn, phi*q(dn,t))=e=0; MinDemand(dn,t).. q(dn,t)-h*Qmin(dn,t)=g=0; iniOFF(ioff,t)$(ord(t)<=toff(ioff)).. z(ioff,t)=e=0; transGRID(t).. sum(n, sum(l, Agrid(l,n)*f(l,t)))=e=0; transMAX(l,t).. f(l,t)-Fmax(l)=l=0; transMIN(l,t).. -f(l,t)-Fmax(l)=l=0; transPHYC2(l,t).. mi*Bsuscep(l)*sum(n, Agrid(l,n)*teta(n,t))-f(l,t)=e=0; transPHYC3(n,t).. -teta(n,t)-pi=l=0; transPHYC4(n,t).. teta(n,t)-pi=l=0; transPHYC5(n,t)$(ord(n)=1).. teta(n,t)=e=0; model UC_MC /all/; option MIQCP=CPLEX; option optcr=0.00; Solve UC_MC using miqcp minimizing objectiveFun ; display v.l, z.l , w.l , g.l ,q.l; *solving the continous SW model by fixing the Z to find the electricity prices Parameter zV(i,t); zV(i,t)=z.l(i,t); Variables objectiveFunn objective function ff(l,t) power flow on line l at period t tetaa(n,t) State variable of node n at time period t ; Positive variables gg(i,tt) generation output of generator i at period t grr(i,t) reserve capacity of generator i at period t rr(i,tt) ramping amount of generator i at period t qq(dn,t) demand of node n at period t vv(i,tt) startup status of generator i at time period t ww(i,tt) shutdown status of generator i at time period t ; Equations objfunn the objective function rampp1(i,t) the first constraint of ramping for generator i at period t rampp2(i,t) the second constraint of ramping for generator i at period t rampp3(i,t) the third constraint of ramping for generator i at period t rampp4(i,t) the forth constraint of ramping for generator i at period t ramp1tonee(i,tone) the first constraint of ramping for generator i at period 1 ramp2tonee(i,tone) the second constraint of ramping for generator i at period 1 ramp3tonee(i,tone) the third constraint of ramping for generator i at period 1 ramp4tonee(i,tone) the forth constraint of ramping for generator i at period 1 MINgenn(i,t) minimum generation level for generator i at period t MAXgenn(i,t) maximum generation level for generator i at period t downtimee(i,t) minimum down time constraint for generator i at period t uptimee(i,t) minimum up time constraint for generator i at period t operationn(i,t) operation status for generator i at period t mktclrr1(dn,t) market clearing constraint at period t for demand nodes mktclrr2(tn,t) market clearing constraint at period t for transfer nodes Reservee(t) reserve constraint at period t MinDemandd(dn,t) minimum demand constraint operationn1(i,t) Extra operational constraint operationtonee(i,tone) operation status for generator i at period 1 operation1tonee(i,tone) Extra operational constraint at period 1 transGRIDD(t) Constraint of Transmission network operator (3g) transMAXX(l,t) Constraint of Transmission network operator (3b) transMINN(l,t) Constraint of Transmission network operator (3c) transPHYCC2(l,t) Constraint of Transmission network operator (3d) transPHYCC3(n,t) Constraint of Transmission network operator (3e) transPHYCC4(n,t) Constraint of Transmission network operator (3e) transPHYCC5(n,t) Constraint of Transmission network operator (3f) ; objfunn.. objectiveFunn =e= sum(i, sum(t, Lambda(i,t)*gg(i,t)+LambdaRES(i,t)*Ares(i)*grr(i,t)+ CSU(i,t)*vv(i,t)+CSD(i,t)*ww(i,t)+CRP(i,t)*rr(i,t)))-sum(dn, sum(t, A(dn,t)*qq(dn,t)-m*B(dn,t)*power(qq(dn,t),2))); rampp1(i,t)$(ord(t) ge 2).. gg(i,t)-gg(i,t-1)-rr(i,t)=l=0; rampp2(i,t)$(ord(t) ge 2).. gg(i,t-1)-gg(i,t)-rr(i,t)=l=0; rampp3(i,t)$(ord(t) ge 2).. gg(i,t)-gg(i,t-1)-RU(i)*zV(i,t)-(SU(i)-RU(i))*vv(i,t)=l=0; rampp4(i,t)$(ord(t) ge 2).. gg(i,t-1)-gg(i,t)-RD(i)*zV(i,t-1)-(SD(i)-RD(i))*ww(i,t)=l=0; ramp1tonee(i,tone).. gg(i,tone)-IG(i)-rr(i,tone)=l=0; ramp2tonee(i,tone).. IG(i)-gg(i,tone)-rr(i,tone)=l=0; ramp3tonee(i,tone).. gg(i,tone)-IG(i)-RU(i)*zV(i,tone)-(SU(i)-RU(i))*vv(i,tone)=l=0; ramp4tonee(i,tone).. IG(i)-gg(i,tone)-RD(i)*IZ(i)-(SD(i)-RD(i))*ww(i,tone)=l=0; MINgenn(i,t).. Gmin(i)*zV(i,t)-gg(i,t)=l=0; MAXgenn(i,t).. gg(i,t)+Ares(i)*grr(i,t)-Gmax(i)*zV(i,t)=l=0; downtimee(i,t)$(Dmin(i)>0).. Dmin(i)*ww(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk) <= ord(t)+ Dmin(i)-1)), 1-zV(i,kk))=l=0; uptimee(i,t)$(Umin(i)>0).. Umin(i)*vv(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk)<= ord(t)+ Umin(i)-1)), zV(i,kk))=l=0; operationn(i,t)$(ord(t)>=2).. zV(i,t-1)-zV(i,t)+vv(i,t)-ww(i,t)=e=0; operationn1(i,t)$(ord(t)>=2).. zV(i,t-1)+vv(i,t)=l=1; operationtonee(i,tone).. IZ(i)-zV(i,tone)+vv(i,tone)-ww(i,tone)=e=0; operation1tonee(i,tone).. IZ(i)+vv(i,tone)=l=1; mktclrr1(dn,t).. Sum(i, Agen(i,dn)*gg(i,t))-qq(dn,t)-sum(l, Agrid(l,dn)*ff(l,t))=e=0; mktclrr2(tn,t).. Sum(i, Agen(i,tn)*gg(i,t))-sum(l, Agrid(l,tn)*ff(l,t))=e=0; Reservee(t).. Sum(i, Ares(i)*grr(i,t))-sum(dn, phi*qq(dn,t))=e=0; MinDemandd(dn,t).. qq(dn,t)-h*Qmin(dn,t)=g=0; transGRIDD(t).. sum(n, sum(l, Agrid(l,n)*ff(l,t)))=e=0; transMAXX(l,t).. ff(l,t)-Fmax(l)=l=0; transMINN(l,t).. -ff(l,t)-Fmax(l)=l=0; transPHYCC2(l,t).. mi*Bsuscep(l)*sum(n, Agrid(l,n)*tetaa(n,t))-ff(l,t)=e=0; transPHYCC3(n,t).. -tetaa(n,t)-pi=l=0; transPHYCC4(n,t).. tetaa(n,t)-pi=l=0; transPHYCC5(n,t)$(ord(n)=1).. tetaa(n,t)=e=0; model UC_SW2 / objfunn rampp1 rampp2 rampp3 rampp4 ramp1tonee ramp2tonee ramp3tonee ramp4tonee MINgenn MAXgenn downtimee uptimee operationn operationn1 operationtonee operation1tonee mktclrr1 mktclrr2 Reservee MinDemandd transGRIDD transMAXX transMINN transPHYCC2 transPHYCC3 transPHYCC4 transPHYCC5 /; option QCP=CPLEX; option optcr=0.00; Solve UC_SW2 using qcp minimizing objectiveFunn ; Parameter PriceDual(tn,t); Parameter ResPrice(t); PriceDual(tn,t)=mktclrr2.m(tn,t); ResPrice(t)=Reservee.m(t); display PriceDual, ResPrice; Parameters GenCost(i) PriceLinear(n,t) DemComp(dn,t) to check that demand bus minQ complementarities are zero GenPaymentLinear(i) ProfitGenLinear(i) TotalProfitGenLin Tdemand(dn,t) NodePaymentLinear(n) TotalNodePaymentLin TransmissionPaymentLinear ; GenCost(i) = Sum(t, Lambda(i,t)*g.l(i,t)+LambdaRES(i,t)*Ares(i)*gr.l(i,t)+CSU(i,t)*v.l(i,t)+CSD(i,t)*w.l(i,t)+CRP(i,t)*r.l(i,t)); PriceLinear(dn,t)=mktclrr1.m(dn,t); DemComp(dn,t)=(PriceLinear(dn,t)+phi*ResPrice(t)-A(dn,t)+B(dn,t)*qq.l(dn,t))*(qq.l(dn,t)-h*Qmin(dn,t)); PriceLinear(tn,t)=PriceDual(tn,t); GenPaymentLinear(i)=Sum(t, Sum(n, Agen(i,n)*PriceLinear(n,t)*g.l(i,t)))+sum(t, Ares(i)*gr.l(i,t)*ResPrice(t)); ProfitGenLinear(i)=GenPaymentLinear(i)-GenCost(i); TotalProfitGenLin=sum(i, ProfitGenLinear(i)); NodePaymentLinear(dn)=sum(t, (PriceLinear(dn,t)+ResPrice(t)*phi)*q.l(dn,t)); TotalNodePaymentLin=sum(dn, NodePaymentLinear(dn)); TransmissionPaymentLinear=sum(t, sum(l, sum(n, -PriceLinear(n,t)*Agrid(l,n)*f.l(l,t)))); Display DemComp, PriceLinear, ProfitGenLinear, TotalProfitGenLin, TotalNodePaymentLin; Execute_Unload 'SW-Tq-Mar13.gdx',v, z, w, g, gr, q, f, GenCost,PriceDual, PriceLinear, ResPrice, GenPaymentLinear, ProfitGenLinear, NodePaymentLinear, TransmissionPaymentLinear ; Execute 'GDXXRW.EXE SW-Tq-Mar13.gdx O=SW-Tq-Mar13.xls Var=v rng=Startup!A1 rdim=1 cdim=1 Var=z rng=GenStatus!A1 rdim=1 cdim=1 Var=w rng=Shutdown!A1 rdim=1 cdim=1 Var=g rng=Generation!A1 rdim=1 cdim=1 Var=gr rng=ResCapacity!A1 rdim=1 cdim=1 Var=q rng=qDemand!A1 rdim=1 cdim=1 Var=f rng=NetworkFlow!A1 rdim=1 cdim=1 Par=GenCost rng=GenFinantial!A1 rdim=0 cdim=1 Par=PriceDual rng=PriceDual!A1 rdim=1 cdim=1 Par=PriceLinear rng=PriceLinear!A1 rdim=1 cdim=1 Par=ResPrice rng=ResPrice!A1 rdim=0 cdim=1 Par=GenPaymentLinear rng=GenFinantial!A8 rdim=0 cdim=1 Par=ProfitGenLinear rng=GenFinantial!A16 rdim=0 cdim=1 Par=NodePaymentLinear rng=NodePayment!A4 rdim=0 cdim=1 Par=TransmissionPaymentLinear rng=TransPayment!A4'; * Solving the individual generation units problem individualy by considering prices come from SW Parameter PriceGEN(i,t); PriceGEN(i,t)=sum(n, Agen(i,n)*priceLinear(n,t)); Variables objGENFunnn objective function ; Positive variables ggg(i,tt) generation output of generator i at period t rrr(i,tt) ramping amount of generator i at period t grrr(i,t) reserve capacity of generator i at period t vvv(i,tt) startup status of generator i at time period t www(i,tt) shutdown status of generator i at time period t Binary variables zzz(i,tt) operation status of generator i at time period t ; Equations objfunGENN the objective function ramppp1(i,t) the first constraint of ramping for generator i at period t ramppp2(i,t) the second constraint of ramping for generator i at period t ramppp3(i,t) the third constraint of ramping for generator i at period t ramppp4(i,t) the forth constraint of ramping for generator i at period t rampp1tonee(i,tone) the first constraint of ramping for generator i at period 1 rampp2tonee(i,tone) the second constraint of ramping for generator i at period 1 rampp3tonee(i,tone) the third constraint of ramping for generator i at period 1 rampp4tonee(i,tone) the forth constraint of ramping for generator i at period 1 MINgennn(i,t) minimum generation level for generator i at period t MAXgennn(i,t) maximum generation level for generator i at period t downttt(i,t) minimum down time constraint for generator i at period t upttt(i,t) minimum up time constraint for generator i at period t operrr(i,t) operation status for generator i at period t opertoneee(i,tone) operation status for generator i at period 1 oper1toneee(i,tone) extra operational constraint at period 1 operrr1(i,t) extra operational constraint initialOFFFF(ioff,t) initial OFF status ; objfunGENN.. objGENFunnn =e= sum(i, sum(t, Lambda(i,t)*ggg(i,t)+LambdaRES(i,t)*Ares(i)*grrr(i,t)+ CSU(i,t)*vvv(i,t)+CSD(i,t)*www(i,t)+CRP(i,t)*rrr(i,t)-PriceGEN(i,t)*ggg(i,t)-ResPrice(t)*Ares(i)*grrr(i,t))); ramppp1(i,t)$(ord(t) ge 2).. ggg(i,t)-ggg(i,t-1)-rrr(i,t)=l=0; ramppp2(i,t)$(ord(t) ge 2).. ggg(i,t-1)-ggg(i,t)-rrr(i,t)=l=0; ramppp3(i,t)$(ord(t) ge 2).. ggg(i,t)-ggg(i,t-1)-RU(i)*zzz(i,t)-(SU(i)-RU(i))*vvv(i,t)=l=0; ramppp4(i,t)$(ord(t) ge 2).. ggg(i,t-1)-ggg(i,t)-RD(i)*zzz(i,t-1)-(SD(i)-RD(i))*www(i,t)=l=0; rampp1tonee(i,tone).. ggg(i,tone)-IG(i)-rrr(i,tone)=l=0; rampp2tonee(i,tone).. IG(i)-ggg(i,tone)-rrr(i,tone)=l=0; rampp3tonee(i,tone).. ggg(i,tone)-IG(i)-RU(i)*zzz(i,tone)-(SU(i)-RU(i))*vvv(i,tone)=l=0; rampp4tonee(i,tone).. IG(i)-ggg(i,tone)-RD(i)*IZ(i)-(SD(i)-RD(i))*www(i,tone)=l=0; MINgennn(i,t).. Gmin(i)*zzz(i,t)-ggg(i,t)=l=0; MAXgennn(i,t).. ggg(i,t)+Ares(i)*grrr(i,t)-Gmax(i)*zzz(i,t)=l=0; downttt(i,t)$(Dmin(i)>0).. Dmin(i)*www(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk) <= ord(t)+ Dmin(i)-1)), 1-zzz(i,kk))=l=0; upttt(i,t)$(Umin(i)>0).. Umin(i)*vvv(i,t)-sum(kk$((ord(kk)>=ord(t)) and (ord(kk)<= ord(t)+ Umin(i)-1)), zzz(i,kk))=l=0; operrr(i,t)$(ord(t)>=2).. zzz(i,t-1)-zzz(i,t)+vvv(i,t)-www(i,t)=e=0; operrr1(i,t)$(ord(t)>=2).. zzz(i,t-1)+vvv(i,t)=l=1; opertoneee(i,tone).. IZ(i)-zzz(i,tone)+vvv(i,tone)-www(i,tone)=e=0; oper1toneee(i,tone).. IZ(i)+vvv(i,tone)=l=1; initialOFFFF(ioff,t)$(ord(t)<=toff(ioff)).. zzz(ioff,t)=e=0; model GEN /objfunGENN ramppp1 ramppp2 ramppp3 ramppp4 rampp1tonee rampp2tonee rampp3tonee rampp4tonee MINgennn MAXgennn downttt upttt operrr operrr1 opertoneee oper1toneee initialOFFFF /; option MIP=CPLEX; Solve GEN using mip minimizing objGENFunnn ; display vvv.l, zzz.l, www.l, ggg.l, grrr.l; Parameters GenCostResPrice(i) GenProfitResPrice(i) DisEquilGen(i) TotDisEquilGen; GenCostResPrice(i)=Sum(t, Lambda(i,t)*ggg.l(i,t)+LambdaRES(i,t)*Ares(i)*grrr.l(i,t) +CSU(i,t)*vvv.l(i,t)+CSD(i,t)*www.l(i,t)+CRP(i,t)*rrr.l(i,t)); GenProfitResPrice(i)= sum(t, priceGEN(i,t)*ggg.l(i,t)+ResPrice(t)*grrr.l(i,t))-GenCostResPrice(i); DisEquilGen(i)=GenProfitResPrice(i)-ProfitGenLinear(i); TotDisEquilGen=sum(i,DisEquilGen(i)); display DisEquilGen ; Execute_Unload 'SW-GEN-Mar13.gdx',vvv, zzz, www, ggg, GenCostResPrice, GenProfitResPrice, DisEquilGen; Execute 'GDXXRW.EXE SW-GEN-Mar13.gdx O=SW-GEN-Mar13.xls Var=vvv rng=startup!A1 rdim=1 cdim=1 Var=zzz rng=Operation!A1 rdim=1 cdim=1 Var=www rng=Shutdown!A1 rdim=1 cdim=1 Var=ggg rng=generation!A1 rdim=1 cdim=1 Par=GenCostResPrice rng=GenResPrice!A1 rdim=1 cdim=0 Par=GenProfitResPrice rng=GenResPrice!D1 rdim=1 cdim=0 Par=DisEquilGen rng=DisEquilibrium!A1 rdim=1 cdim=0 '; **Soloving the consumers' model with considering the price determined by the MC model** Variable ConsResPriceObj Consumers' objective function; Positive Variable qqq(dn,t) Demand 1uantity; Equations MinQ(dn,t) minimum demand constraint ConsObj objective function value; ConsObj.. ConsResPriceObj=e=sum(dn, sum(t, (PriceLinear(dn,t)+phi*ResPrice(t))*qqq(dn,t)-(A(dn,t)*qqq(dn,t)-m*B(dn,t)*power(qqq(dn,t),2)))); MinQ(dn,t).. qqq(dn,t)-h*Qmin(dn,t)=g=0; Model cons /ConsObj, MinQ /; OPtion QCP=CPLEX; Solve cons Using QCP minimizing ConsResPriceObj; Parameter UtilityIndividual(dn) UtilitySW(dn) DisEquil(dn) TotDisEquilCons; UtilityIndividual(dn)=sum(t, -(PriceLinear(dn,t)+phi*ResPrice(t))*qqq.l(dn,t)+(A(dn,t)*qqq.l(dn,t)-m*B(dn,t)*power(qqq.l(dn,t),2))); UtilitySW(dn)=sum(t, -(PriceLinear(dn,t)+phi*ResPrice(t))*q.l(dn,t)+(A(dn,t)*q.l(dn,t)-m*B(dn,t)*power(q.l(dn,t),2))); DisEquil(dn)=UtilityIndividual(dn)-UtilitySW(dn); TotDisEquilCons = sum(dn, DisEquil(dn)); display DisEquil ; Execute_Unload 'SW-Cons-Mar13.gdx',qq, UtilitySW, UtilityIndividual, DisEquil; Execute 'GDXXRW.EXE SW-Cons-Mar13.gdx O=SW-Cons-Mar13.xls Var=qq rng=ResponsiveDemand!A1 rdim=1 cdim=1 Par=UtilitySW rng=UtilitySW!A1 rdim=0 cdim=1 Par=UtilityIndividual rng=UtilityIndividual!A1 rdim=0 cdim=1 Par=DisEquil rng=DisEquil!A1 rdim=0 cdim=1'; * Solving the individual Transmission problem given the prices come from SW Variables tranPay ff(l,t) tetaa(n,t); Equations tranOBJ tranGRID(t) tranMAX(l,t) tranMIN(l,t) tranPHYC2(l,t) tranPHYC3(n,t) tranPHYC4(n,t) tranPHYC5(n,t) ; tranOBJ.. tranPay=e=sum(t, sum(l, sum(n, PriceLinear(n,t)*Agrid(l,n)*ff(l,t)))); tranGRID(t).. sum(n, sum(l, Agrid(l,n)*ff(l,t)))=e=0; tranMAX(l,t).. ff(l,t)-Fmax(l)=l=0; tranMIN(l,t).. -ff(l,t)-Fmax(l)=l=0; tranPHYC2(l,t).. mi*Bsuscep(l)*sum(n, Agrid(l,n)*tetaa(n,t))-ff(l,t)=e=0; tranPHYC3(n,t).. -tetaa(n,t)-pi=l=0; tranPHYC4(n,t).. tetaa(n,t)-pi=l=0; tranPHYC5(n,t)$(ord(n)=1).. tetaa(n,t)=e=0; Model trans /tranOBJ, tranGRID, tranMAX, tranMIN, tranPHYC2, tranPHYC3, tranPHYC4, tranPHYC5/; OPtion LP=CPLEX; Solve trans Using LP minimizing tranPay; Parameter DisEquilTrans transResPricePay; transResPricePay=-tranPay.l; DisEquilTrans = transResPricePay - TransmissionPaymentLinear; display DisEquilTrans,TotDisEquilCons,TotDisEquilGen; Execute_Unload 'SW-Trans-Mar13.gdx',ff, transResPricePay; Execute 'GDXXRW.EXE SW-Trans-Mar13.gdx O=SW-Trans-Mar13.xls Var=ff rng=FlowResPrice!A1 rdim=1 cdim=1 Par=transResPricePay rng=TransPayment!A1 rdim=0 cdim=0';