Hi; I think this book might come in handy:

Power system Optimization Modeling in GAMS by Alireza Soroudi.

Some of the chapters are available to download, anyway you can find the book in Amazon for $127.20 (Kindle version).

Don't know if the following code will be useful but you igth find exaactly what you are looking for in the book.

Best regards.

GAMS CODE

$ONTEXT

DC-OPF GAMS code for three-bus network.

$OFFTEXT

Sets

bus /1*3/

slack(bus) /3/

Gen /g1*g3/;

Scalars

Sbase /100/;

Alias (bus,node);

Table GenData(Gen ,*) Generating units characteristics

b pmin pmax

g1 10 0 65

g2 11 0 100;

Set

GBconect(bus, Gen) connectivity index of each generating unit to each bus

/ 1 . g1

3 . g2/;

Table BusData(bus ,*) Demands of each bus in MW

Pd

2 100;

Set conex Bus connectivity matrix

/

1 . 2

2 . 3

1 . 3/;

conex(bus,node)$(conex(node,bus))=1;

Table branch (bus, node,*) Network technical characteristics

x Limit

1 . 2 0.2 100

2 . 3 0.25 100

1 . 3 0.4 100;

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');

Variables

OF

Pij(bus, node)

Pg(Gen)

delta(bus);

Equations const1 , const2 , const3;

const1(bus, node)$conex(bus, node).. Pij(bus, node) =E=

branch(bus, node, 'bij')*(delta(bus)-

delta(node));

const2(bus).. +sum(Gen$GBconect(bus, Gen), Pg(Gen))-BusData(bus ,'pd')/Sbase =E=

+sum(node$conex(node, bus), Pij(bus, node));

const3.. OF =G= sum(Gen, Pg(Gen)*GenData(Gen, 'b')*Sbase);

Model loadflow /const1, const2, const3/;

Pg.lo(Gen)=GenData(Gen, 'Pmin')/Sbase;

Pg.up(Gen)=GenData(Gen, 'Pmax')/Sbase;

delta.up(bus)=pi;

delta.lo(bus)=-pi;

delta.fx(slack)=0;

Pij.up(bus, node)$((conex(bus, node)))=1*branch(bus, node, 'Limit')/Sbase;

Pij.lo(bus, node)$((conex(bus, node)))=-1*branch(bus, node, 'Limit')/Sbase;

Solve loadflow minimizing OF using lp;

Parameter report(bus,*), Congestioncost;

report(bus, 'Gen(MW)')= sum(Gen$GBconect(bus, Gen), Pg.l(Gen))*sbase;

report(bus, 'Angle')= delta.l(bus);

report(bus, 'load (MW)')= BusData(bus, 'pd');

report(bus, 'LMP($/MWh)')=const2.m(bus)/sbase;

Congestioncost= sum((bus, node), Pij.l(bus, node)*(-const2.m(bus)+

const2.m(node)))/2;

Display report, Pij.l, Congestioncost;