как моделировать IEEE 6-шинный системный блок? - PullRequest
Новая
       51

как моделировать IEEE 6-шинный системный блок?

1 голос
/ 02 мая 2019

Я пытаюсь смоделировать систему с 6 шинами IEEE, чтобы минимизировать стоимость, код работает хорошо, но у меня есть проблема в моделировании уравнения начальных затрат, которое можно найти в уравнениях № 12 и 13 https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1664974

проблема в этой строке const5 (gen, t) .. suc (Gen, t) = e = GD (Gen, 'Scost') * (сумма (t, u (Gen, t)))); 

* 6 bus system with 24 hours demand
set bus number of buses      /1*6/
    Gen   number of generators /g1*g3/
    t   number of hours      /t1*t24/
    k                        /1*100/
    nl       number of segment /1*3/
    slack(bus) / 1     /
    j start up cost intervals /j1*j8/
    n             /1*24/;
alias(t,tt);

Scalar
   Sbase /   100 /
   VOLL  / 10000 /
   VOLW  /    50 /;

Alias (bus,node);

Table GD(Gen,*) 'generating units characteristics'
        Pmax Pmin a      b     c     Scost  Dcost fuelc MU    MD   RU   RD
   g1   220  100  0.05   10    100   100    10    1     4     4    55   55
   g2   100  10   0.001  40.66 162   200    20    1     2     3    50   50
   g3   20   10   0.006  22.06 171   0      0     1                20   20 ;


Set GB(bus,Gen) 'connectivity index of each generating unit to each bus'
/
1.g1
2.g2
6.g3/;

Table branch(bus,node,*) line information
       x       limit
1.2    0.170   200
1.4    0.258   100
2.3    0.037   100
2.4    0.197   100
3.6    0.018   100
4.5    0.037   100
5.6    0.140   100;

Table WD(t,*) hourly load demand
        d
   t1   175
   t2   165
   t3   158
   t4   154
   t5   155
   t6   160
   t7   168
   t8   177
   t9   186
   t10  206
   t11  228
   t12  236
   t13  242
   t14  243
   t15  248
   t16  255
   t17  256
   t18  246
   t19  245
   t20  237
   t21  236
   t22  232
   t23  208
   t24  195;

Table BusData(bus,*) percent of hourly load of each bus
    pd
3   0.20
4   0.40
5   0.40 ;

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')$branch(bus,node,'Limit')       = 1/branch(bus,node,'x');
branch(bus,node,'z')$branch(bus,node,'Limit')         = branch(bus,node,'x');
branch(node,bus,'z')                                  = branch(bus,node,'z');


binary variable u(Gen,t) status of unit i in period t;
Variable OF, Pg(Gen,t),pk(Gen,t,k);
positive variable suc(Gen,t);
Parameter conex(bus,node);
conex(bus,node)$(branch(bus,node,'limit') and branch(node,bus,'limit')) = 1;
conex(bus,node)$(conex(node,bus)) = 1;
Parameter data (k, Gen ,*) ;
data(k,Gen,'DP')=(GD(Gen ,"Pmax")-GD( Gen ,"Pmin")) / card ( k ) ;
data (k , Gen , 'Pini')= ( ord (k)-1)*data( k , Gen , 'DP' ) + GD( Gen , "Pmin" ) ;
data ( k , Gen , 'Pfin' ) = data( k , Gen , 'Pini' ) + data( k , Gen , 'DP' ) ;
data (k, Gen , 'Cini')=GD (Gen ,"a")*( data (k , Gen , 'Pini' ))*( data (k , Gen , 'Pini' ))
+GD (Gen ,"b")*data (k, Gen , 'Pini')+GD (Gen ,"c") ;
data (k, Gen , 'Cfin')=GD (Gen ,"a")*(data (k , Gen , 'Pfin' ))*(data(k , Gen , 'Pfin' ))
+GD (Gen ,"b")*data (k, Gen , 'Pfin')+GD (Gen ,"c") ;
data (k , Gen , 's') =( data (k , Gen , 'Cfin')-data( k , Gen , 'Cini' ) ) / data ( k , Gen , 'DP' ) ;

Pk.up ( Gen , t , k ) = data ( k , Gen , 'DP' ) ;
Pk.lo ( Gen , t , k ) =0;

Equation const1, const2, const3,const4,const5;

const1.. OF =e=sum ((bus,Gen,t)$GB(bus,Gen) , GD (Gen , 'a')*GD( Gen ,"Pmin")*GD( Gen ,"Pmin")
+GD (Gen , 'b')*GD ( Gen ,"Pmin") +GD ( Gen , 'c' )
+ sum (k, data (k, Gen , 's')*pk ( Gen , t , k) ) ) ;
*const1.. OF =e= sum((bus,Gen,t)$GB(bus,Gen), Pg(Gen,t)*Pg(Gen,t)*GD(Gen,'a')*Sbase+Pg(Gen,t)*GD(Gen,'b')*Sbase+GD(Gen,'c')) ;
const2(Gen,t).. pg(Gen,t+1) - pg(Gen,t) =l= GD(Gen,'RU')/Sbase;

const3(Gen,t).. pg(Gen,t-1) - pg(Gen,t) =l= GD(Gen,'RD')/Sbase;
const4(bus,t).. sum(Gen,pg(Gen,t)) =g= WD(t,'d')*BusData(bus,'pd')/Sbase;
const5(gen,t).. suc(Gen,t)=e= GD (Gen , 'Scost')*(sum(t,u(Gen,t)));


Pg.lo(Gen,t) = GD(Gen,'Pmin')/Sbase;
Pg.up(Gen,t) = GD(Gen,'Pmax')/Sbase;


Model loadflow /all/;


Option LP = Cplex;
solve loadflow us mip min OF;
...