fm_gams.gms

用于电力系统的一个很好的分析软件

kg.up =  1;
lambdac.lo = 0;

* -------------------------------------------------------------------------
* objective function
cost.. obj =e= -lambdac;
* -------------------------------------------------------------------------

Peq(B).. (lambdac+kg)*(sum(G,Ps_idx(B,G)*Ps(G)) + X(B,'Pg0')) -
	 lambdac*(sum(C,Pd_idx(B,C)*Pd(C)) + X(B,'Pl0')) -
	 V(B)*sum(BB,V(BB)*(Gh(B,BB)*cos(a(B)-a(BB)) +
	 Bh(B,BB)*sin(a(B)-a(BB)))) =e= 0;

Qeq(B).. Qg(B) - lambdac*(X(B,'Ql0') +
	 sum(C,Pd_idx(B,C)*(D(C,'tgphi')*Pd(C)))) -
	 V(B)*sum(BB,V(BB)*(Gh(B,BB)*sin(a(B)-a(BB)) -
	 Bh(B,BB)*cos(a(B)-a(BB)))) =e= 0;

Thetaref(B)$busref(B).. a(B) =e= 0;

* -------------------------------------------------------------------------
* flow limit equations
* -------------------------------------------------------------------------

$label floweq

$if %flow% == 0 $goto solvestat
$if %flow% == 1 $goto iflows
$if %flow% == 2 $goto pflows
$if %flow% == 3 $goto sflows

$label pflows
Pijeq(L).. Pij(L) =e= N(L,'g0')*sum(B,Li(L,B)*V(B)*V(B)) -
	   N(L,'b')*sum(B,Li(L,B)*V(B)*sin(a(B)))*
	   sum(B,Lj(L,B)*V(B)*cos(a(B))) +
	   N(L,'b')*sum(B,Li(L,B)*V(B)*cos(a(B)))*
	   sum(B,Lj(L,B)*V(B)*sin(a(B))) -
	   N(L,'g')*sum(B,Li(L,B)*V(B)*cos(a(B)))*
	   sum(B,Lj(L,B)*V(B)*cos(a(B))) -
	   N(L,'g')*sum(B,Li(L,B)*V(B)*sin(a(B)))*
	   sum(B,Lj(L,B)*V(B)*sin(a(B)));

Pjieq(L).. Pji(L) =e= N(L,'g0')*sum(B,Lj(L,B)*V(B)*V(B)) -
	   N(L,'b')*sum(B,Lj(L,B)*V(B)*sin(a(B)))*
	   sum(B,Li(L,B)*V(B)*cos(a(B))) +
	   N(L,'b')*sum(B,Lj(L,B)*V(B)*cos(a(B)))*
	   sum(B,Li(L,B)*V(B)*sin(a(B))) -
	   N(L,'g')*sum(B,Lj(L,B)*V(B)*cos(a(B)))*
	   sum(B,Li(L,B)*V(B)*cos(a(B))) -
	   N(L,'g')*sum(B,Lj(L,B)*V(B)*sin(a(B)))*
	   sum(B,Li(L,B)*V(B)*sin(a(B)));

$if %control% == 3 $goto solvestat 
$if %control% == 5 $goto solvestat 

N(L,'b')$(ord(L) eq lambda('line')) = -1E-6;
N(L,'g')$(ord(L) eq lambda('line')) = 0;
N(L,'b0')$(ord(L) eq lambda('line')) = 0;
N(L,'g0')$(ord(L) eq lambda('line')) = 0;

Pijceq(L).. Pijc(L) =e= N(L,'g0')*sum(B,Li(L,B)*Vc(B)*Vc(B)) -
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))*
	    sum(B,Lj(L,B)*Vc(B)*cos(ac(B))) +
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))*
	    sum(B,Lj(L,B)*Vc(B)*sin(ac(B))) -
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))*
	    sum(B,Lj(L,B)*Vc(B)*cos(ac(B))) -
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))*
	    sum(B,Lj(L,B)*Vc(B)*sin(ac(B)));

Pjiceq(L).. Pjic(L) =e= N(L,'g0')*sum(B,Lj(L,B)*Vc(B)*Vc(B)) -
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))*
	    sum(B,Li(L,B)*Vc(B)*cos(ac(B))) +
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))*
	    sum(B,Li(L,B)*Vc(B)*sin(ac(B))) -
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))*
	    sum(B,Li(L,B)*Vc(B)*cos(ac(B))) -
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))*
	    sum(B,Li(L,B)*Vc(B)*sin(ac(B)));

$goto solvestat

$label iflows

Pijeq(L).. Pij(L) =e= sqrt(
	   sqr(N(L,'g')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	       N(L,'g')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	       N(L,'b')*sum(B,Li(L,B)*V(B)*sin(a(B)))+
	       N(L,'b')*sum(B,Lj(L,B)*V(B)*sin(a(B)))-
	       N(L,'b0')*sum(B,Li(L,B)*V(B)*sin(a(B))))+
	   sqr(N(L,'b')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	       N(L,'b')*sum(B,Lj(L,B)*V(B)*cos(a(B)))+
	       N(L,'g')*sum(B,Li(L,B)*V(B)*sin(a(B)))-
	       N(L,'g')*sum(B,Lj(L,B)*V(B)*sin(a(B)))+
	       N(L,'b0')*sum(B,Li(L,B)*V(B)*cos(a(B)))));

Pjieq(L).. Pji(L) =e= sqrt(
	   sqr(N(L,'g')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	       N(L,'g')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	       N(L,'b')*sum(B,Lj(L,B)*V(B)*sin(a(B)))+
	       N(L,'b')*sum(B,Li(L,B)*V(B)*sin(a(B)))-
	       N(L,'b0')*sum(B,Lj(L,B)*V(B)*sin(a(B))))+
	   sqr(N(L,'b')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	       N(L,'b')*sum(B,Li(L,B)*V(B)*cos(a(B)))+
	       N(L,'g')*sum(B,Lj(L,B)*V(B)*sin(a(B)))-
	       N(L,'g')*sum(B,Li(L,B)*V(B)*sin(a(B)))+
	       N(L,'b0')*sum(B,Lj(L,B)*V(B)*cos(a(B)))));

$if %control% == 3 $goto solvestat 
$if %control% == 5 $goto solvestat 

N(L,'b')$(ord(L) eq lambda('line')) = -1E-6;
N(L,'g')$(ord(L) eq lambda('line')) = 0;
N(L,'b0')$(ord(L) eq lambda('line')) = 0;
N(L,'g0')$(ord(L) eq lambda('line')) = 0;

Pijceq(L).. Pijc(L) =e= sqrt(
	    sqr(N(L,'g')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'b0')*sum(B,Li(L,B)*Vc(B)*sin(ac(B))))+
	    sqr(N(L,'b')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))+
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b0')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))));

Pjiceq(L).. Pjic(L) =e= sqrt(
	    sqr(N(L,'g')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'b0')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B))))+
	    sqr(N(L,'b')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))+
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b0')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))));

$goto solvestat

$label sflows

Pijeq(L).. Pij(L) =e= sum(B,Li(L,B)*V(B))*sqrt(
	   sqr(N(L,'g')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	   N(L,'g')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	   N(L,'b')*sum(B,Li(L,B)*V(B)*sin(a(B)))+
	   N(L,'b')*sum(B,Lj(L,B)*V(B)*sin(a(B)))-
	   N(L,'b0')*sum(B,Li(L,B)*V(B)*sin(a(B))))+
	   sqr(N(L,'b')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	   N(L,'b')*sum(B,Lj(L,B)*V(B)*cos(a(B)))+
	   N(L,'g')*sum(B,Li(L,B)*V(B)*sin(a(B)))-
	   N(L,'g')*sum(B,Lj(L,B)*V(B)*sin(a(B)))+
	   N(L,'b0')*sum(B,Li(L,B)*V(B)*cos(a(B)))));

Pjieq(L).. Pji(L) =e= sum(B,Lj(L,B)*V(B))*sqrt(
	   sqr(N(L,'g')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	   N(L,'g')*sum(B,Li(L,B)*V(B)*cos(a(B)))-
	   N(L,'b')*sum(B,Lj(L,B)*V(B)*sin(a(B)))+
	   N(L,'b')*sum(B,Li(L,B)*V(B)*sin(a(B)))-
	   N(L,'b0')*sum(B,Lj(L,B)*V(B)*sin(a(B))))+
	   sqr(N(L,'b')*sum(B,Lj(L,B)*V(B)*cos(a(B)))-
	   N(L,'b')*sum(B,Li(L,B)*V(B)*cos(a(B)))+
	   N(L,'g')*sum(B,Lj(L,B)*V(B)*sin(a(B)))-
	   N(L,'g')*sum(B,Li(L,B)*V(B)*sin(a(B)))+
	   N(L,'b0')*sum(B,Lj(L,B)*V(B)*cos(a(B)))));

$if %control% == 3 $goto solvestat 
$if %control% == 5 $goto solvestat 

N(L,'b')$(ord(L) eq lambda('line')) = -1E-6;
N(L,'g')$(ord(L) eq lambda('line')) = 0;
N(L,'b0')$(ord(L) eq lambda('line')) = 0;
N(L,'g0')$(ord(L) eq lambda('line')) = 0;

Pijceq(L).. Pijc(L) =e= sum(B,Li(L,B)*Vc(B))*sqrt(
	    sqr(N(L,'g')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'b0')*sum(B,Li(L,B)*Vc(B)*sin(ac(B))))+
	    sqr(N(L,'b')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))+
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b0')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))));

Pjiceq(L).. Pjic(L) =e= sum(B,Lj(L,B)*Vc(B))*sqrt(
	    sqr(N(L,'g')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'b0')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B))))+
	    sqr(N(L,'b')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))-
	    N(L,'b')*sum(B,Li(L,B)*Vc(B)*cos(ac(B)))+
	    N(L,'g')*sum(B,Lj(L,B)*Vc(B)*sin(ac(B)))-
	    N(L,'g')*sum(B,Li(L,B)*Vc(B)*sin(ac(B)))+
	    N(L,'b0')*sum(B,Lj(L,B)*Vc(B)*cos(ac(B)))));

$goto solvestat

* -------------------------------------------------------------------------
* Solve Market Problem
* -------------------------------------------------------------------------

$label solvestat
model market /all/;
scalar mstatus  /0/;
scalar sstatus  /0/;

option nlp = minos;
option iterlim = 50000;

$if %control% == 1 $goto linearmodel
$if %control% == 2 $goto linearmodel
$if %control% == 3 $goto nonlinearmodel
$if %control% == 4 $goto nonlinearmodel
$if %control% == 5 $goto nonlinearmodel
$if %control% == 6 $goto nonlinearmodel

$label linearmodel

* option lp = cplex;
solve market using lp minimizing obj;

if ((market.modelstat eq 1), mstatus = 1;);
if ((market.modelstat eq 2), mstatus = 2;);
if ((market.modelstat eq 3), mstatus = 3;);
if ((market.modelstat eq 4), mstatus = 4;);
if ((market.modelstat eq 5), mstatus = 5;);
if ((market.modelstat eq 6), mstatus = 6;);
if ((market.modelstat eq 7), mstatus = 7;);
if ((market.modelstat eq 8), mstatus = 8;);
if ((market.modelstat eq 9), mstatus = 9;);
if ((market.modelstat eq 10), mstatus = 10;);
if ((market.modelstat eq 11), mstatus = 11;);
if ((market.modelstat eq 12), mstatus = 12;);
if ((market.modelstat eq 13), mstatus = 13;);

if ((market.solvestat eq 1), sstatus = 1;);
if ((market.solvestat eq 2), sstatus = 2;);
if ((market.solvestat eq 3), sstatus = 3;);
if ((market.solvestat eq 4), sstatus = 4;);
if ((market.solvestat eq 5), sstatus = 5;);
if ((market.solvestat eq 6), sstatus = 6;);
if ((market.solvestat eq 7), sstatus = 7;);
if ((market.solvestat eq 8), sstatus = 8;);
if ((market.solvestat eq 9), sstatus = 9;);
if ((market.solvestat eq 10), sstatus = 10;);
if ((market.solvestat eq 11), sstatus = 11;);
if ((market.solvestat eq 12), sstatus = 12;);
if ((market.solvestat eq 13), sstatus = 13;);

$goto psatoutput

$label nonlinearmodel

solve market using nlp minimizing obj;

if ((market.modelstat eq 1), mstatus = 1;);
if ((market.modelstat eq 2), mstatus = 2;);
if ((market.modelstat eq 3), mstatus = 3;);
if ((market.modelstat eq 4), mstatus = 4;);
if ((market.modelstat eq 5), mstatus = 5;);
if ((market.modelstat eq 6), mstatus = 6;);
if ((market.modelstat eq 7), mstatus = 7;);
if ((market.modelstat eq 8), mstatus = 8;);
if ((market.modelstat eq 9), mstatus = 9;);
if ((market.modelstat eq 10), mstatus = 10;);
if ((market.modelstat eq 11), mstatus = 11;);
if ((market.modelstat eq 12), mstatus = 12;);
if ((market.modelstat eq 13), mstatus = 13;);

if ((market.solvestat eq 1), sstatus = 1;);
if ((market.solvestat eq 2), sstatus = 2;);
if ((market.solvestat eq 3), sstatus = 3;);
if ((market.solvestat eq 4), sstatus = 4;);
if ((market.solvestat eq 5), sstatus = 5;);
if ((market.solvestat eq 6), sstatus = 6;);
if ((market.solvestat eq 7), sstatus = 7;);
if ((market.solvestat eq 8), sstatus = 8;);
if ((market.solvestat eq 9), sstatus = 9;);
if ((market.solvestat eq 10), sstatus = 10;);
if ((market.solvestat eq 11), sstatus = 11;);
if ((market.solvestat eq 12), sstatus = 12;);
if ((market.solvestat eq 13), sstatus = 13;);

$if not %control% == 4 $goto psatoutput

lambdac.up = lambdac.l;
lambdac.lo = lambdac.l-1.0E-5;
lambda('omega') = 0;
option nlp = minos;
solve market using nlp minimizing obj;

$label psatoutput

$libinclude psatout Ps.l G
$libinclude psatout Pd.l C

$if %control% == 1 $libinclude psatout Pbalance.m
$if %control% == 1 $goto end_output

$if %control% == 2 $libinclude psatout Peq.m B
$if %control% == 2 $goto end_output

$libinclude psatout V.l B
$libinclude psatout a.l B
$libinclude psatout Qg.l B
$libinclude psatout Peq.m B
$libinclude psatout Pij.l L
$libinclude psatout Pji.l L

$if %control% == 5 $goto no_dual

$libinclude psatout V.m B
$libinclude psatout Pij.m L
$libinclude psatout Pji.m L

$label no_dual

$if %control% == 3 $goto end_output

$libinclude psatout lambdac.l
$libinclude psatout kg.l

$if %control% == 5 $goto end_output

$libinclude psatout Vc.l B
$libinclude psatout ac.l B
$libinclude psatout Qgc.l B
$libinclude psatout Pijc.l L
$libinclude psatout Pjic.l L

$libinclude psatout Pceq.m B

$if %control% == 4 $goto end_output

$libinclude psatout lambdac.m

$label end_output

$libinclude psatout mstatus
$libinclude psatout sstatus