$title Optimal power flow for a three-bus system $onText For more details please refer to Chapter 6 (Gcode6.2), of the following book: Soroudi, Alireza. Power System Optimization Modeling in GAMS. Springer, 2017. -------------------------------------------------------------------------------- Model type: LP -------------------------------------------------------------------------------- $offText Set bus / 1*3 / slack(bus) / 1 / Gen / g1*g2 / t / t1*t24 /; Scalar Sbase / 100 /; Alias (bus,node); Table GenData(Gen,*) 'generating units characteristics' b pmin pmax RU RD g1 20 100 250 50 50 g2 30 100 250 40 40; * ----------------------------------------------------- Set GBconect(bus,Gen) 'connectivity index of each generating unit to each bus' / 1.g1 2.g2 /; Table BusData(bus,*) 'demands of each bus in MW' Pd 2 200 3 60; Table branch(bus,node,*) 'network technical characteristics' x Limit 1.2 0.1 200 1.3 0.2 60 2.3 0.2 80; Table WD(t,*) w d t1 0.0786666666666667 1 t2 0.0866666666666667 0.944122690036197 t3 0.117333333333333 0.91306915602972 t4 0.258666666666667 0.499733282530006 t5 0.361333333333333 0.888874071251667 t6 0.566666666666667 0.8980186702229 t7 0.650666666666667 0.826786054486569 t8 0.566666666666667 0.851743189178891 t9 0.484 0.806039245570585 t10 0.548 0.887007048961707 t11 0.757333333333333 0.839016955610593 t12 0.710666666666667 0.852733854067441 t13 0.870666666666667 0.870642027052772 t14 0.932 0.834254143646409 t15 0.966666666666667 0.816536483139646 t16 1 0.819394170318156 t17 0.869333333333333 0.874071251666984 t18 0.665333333333333 1 t19 0.656 0.983615926843208 t20 0.561333333333333 0.936368832158506 t21 0.565333333333333 0.887597637645266 t22 0.556 0.809297008954087 t23 0.724 0.74585635359116 t24 0.84 0.733473042484283; Set conex 'bus connectivity matrix' / 1.2, 1.3, 2.3 /; conex(bus,node)$(conex(node,bus)) = 1; branch(bus,node,'x')$(branch(bus,node,'x')=0) = branch(node,bus,'x'); branch(bus,node,'Limit')$(branch(bus,node,'Limit')=0) = branch(node,bus,'Limit'); branch(bus,node,'bij')$conex(bus,node) = 1/branch(bus,node,'x'); ***************************************************** Variable OF, Pij(bus,node,t), Pg(Gen,t), delta(bus,t); Equation const1, const2, const3, const4, const5; const1(bus,node,t)$(conex(bus,node)).. Pij(bus,node,t) =e= branch(bus,node,'bij')*(delta(bus,t) - delta(node,t)); * active flow calculation between each pair of connected buses const2(bus,t).. sum(Gen$GBconect(bus,Gen), Pg(Gen,t)) - WD(t,'d')*BusData(bus,'Pd')/Sbase =e= sum(node$conex(node,bus), Pij(bus,node,t)); * nodal active power balance in each bus const3.. OF =e= sum((bus,Gen,t)$GBconect(bus,Gen), Pg(Gen,t)*GenData(Gen,'b')*Sbase); * objective function calculation const4(gen,t).. Pg(gen,t+1) - Pg(gen,t) =l= GenData(gen,'RU')/Sbase; const5(gen,t).. Pg(gen,t-1) - Pg(gen,t) =l= GenData(gen,'RD')/Sbase; Model loadflow / const1, const2, const3 /; Pg.lo(Gen,t) = GenData(Gen,'Pmin')/Sbase; Pg.up(Gen,t) = GenData(Gen,'Pmax')/Sbase; delta.up(bus,t) = pi; delta.lo(bus,t) =-pi; delta.fx(slack,t) = 0; Pij.up(bus,node,t)$((conex(bus,node))) = 1*branch(bus,node,'Limit')/Sbase; Pij.lo(bus,node,t)$((conex(bus,node))) =-1*branch(bus,node,'Limit')/Sbase; solve loadflow minimizing OF using lp; Parameter report(t,bus,*), Congestioncost,lmp(bus,t); report(t,bus,'Gen(MW)') = sum(Gen$GBconect(bus,Gen), Pg.l(Gen,t))*sbase; report(t,bus,'Angle') = delta.l(bus,t); report(t,bus,'load(MW)') = BusData(bus,'pd')*WD(t,'d'); report(t,bus,'LMP($/MWh)') = const2.m(bus,t)/sbase; Congestioncost = sum((bus,node,t)$conex(bus,node), Pij.l(bus,node,t)*(-const2.m(bus,t) + const2.m(node,t)))/2; lmp(bus,t) = report(t,bus,'LMP($/MWh)'); display report, Pij.l, Congestioncost,OF.l;