stochastic programming Topic is solved

Problems with modeling
Post Reply
anan
User
User
Posts: 20
Joined: 5 years ago

stochastic programming

Post by anan »

I was trying to minimize the cost of three generating units and the code is working well but when I added the stochastic part random variable (PV) it became infeasible no solution please anyone can help me

Code: Select all

* 6-bus system UC without line constraints only generating units and load 
Set
   t 'time'               / t1*t24   /
   i 'generators indices' / g1*g3   /
   k 'cost segments'      / sg1*sg20 /
   char                   / ch1*ch2  /;
scalar pv /4/;
Alias (t,h);

Table gendata(i,*) 'generator cost characteristics and limits'
       bus Pmax Pmin a      b      c       CostsD  costst RU   RD   SU   SD   UT   DT   uini U0  S0 fuel
   g1  1   220  100  0.050  10     100     0       100    55   55   105  108  4    4    1    4   0  1
   g2  2   100  10   0.001  40.66  162     0       200    50   50   106  112  2    3    0    0   3  1
   g3  6   20   10   0.006  22.06  171     0       0      20   20   43   45   1    1    0    0   1  1;
* -----------------------------------------------------
* linearizing the quadratic fuel cost FC = a*p*p + b*p + c;
Parameter data(k,i,*);
data(k,i,'DP')       =(gendata(i,"Pmax") - gendata(i,"Pmin"))/card(k);
data(k,i,'Pini')     =(ord(k) - 1)*data(k,i,'DP') + gendata(i,"Pmin");
data(k,i,'Pfin')     = data(k,i,'Pini') + data(k,i,'DP');
data(k,i,'Cini')     = gendata(i,"a")*power(data(k,i,'Pini'),2)
                     + gendata(i,"b")*data(k,i,'Pini') + gendata(i,"c");
data(k,i,'Cfin')     = gendata(i,"a")*power(data(k,i,'Pfin'),2)
                     + gendata(i,"b")*data(k,i,'Pfin') + gendata(i,"c");
data(k,i,'s')        =(data(k,i,'Cfin') - data(k,i,'Cini'))/data(k,i,'DP');
gendata(i,'Mincost') = gendata(i,'a')*power(gendata(i,"Pmin"),2)
                     + gendata(i,'b')*gendata(i,"Pmin") + gendata(i,'c');
* -----------------------------------------------------
* 24 hours load profile

Table dataLP(t,*)
        load
   t1   175.19
   t2   165.15
   t3   158.67
   t4   154.73
   t5   155.06
   t6   160.48
   t7   168.39
   t8   177.6
   t9   186.81
   t10  206.96
   t11  228.61
   t12  236.1
   t13  242.18
   t14  243.6
   t15  248.86
   t16  255.79
   t17  256
   t18  246.74
   t19  245.97
   t20  237.35
   t21  236.31
   t22  232.67
   t23  208.93
   t24  195.6;

* eq. (5.5d) page 122
Parameter unit(i,char);
unit(i,'ch1') = 24;
unit(i,'ch2') =(gendata(i,'UT') - gendata(i,'U0'))*gendata(i,'Uini');
gendata(i,'Lj') = smin(char,unit(i,char));

*eq. (5.6d) page 123
Parameter unit2(i,char);
unit2(i,'ch1')  = 24;
unit2(i,'ch2')  =(gendata(i,'DT') - gendata(i,'S0'))*(1 - gendata(i,'Uini'));
gendata(i,'Fj') = smin(char,unit2(i,char));

Variable costThermal;
Positive Variable pu(i,t), p(i,t), StC(i,t), SDC(i,t), Pk(i,t,k);
Binary Variable u(i,t), y(i,t), z(i,t);

p.up(i,t)    = gendata(i,"Pmax");
p.lo(i,t)    = 0;
Pk.up(i,t,k) = data(k,i,'DP');
Pk.lo(i,t,k) = 0;
p.up(i,t)    = gendata(i,"Pmax");
pu.up(i,h)   = gendata(i,"Pmax");

Equation
   Uptime1, Uptime2, Uptime3, Dntime1, Dntime2, Dntime3, Ramp1, Ramp2,
   startc, shtdnc, genconst1, Genconst2, Genconst3, Genconst4, balance;
* minimum up time constraints
Uptime1(i)$(gendata(i,"Lj")>0)..
   sum(t$(ord(t)<(gendata(i,"Lj")+1)), 1 - U(i,t)) =e= 0;

Uptime2(i)$(gendata(i,"UT")>1)..
   sum(t$(ord(t)>24-gendata(i,"UT")+1), U(i,t) - y(i,t)) =g= 0;

Uptime3(i,t)$(ord(t)>gendata(i,"Lj") and ord(t)<24-gendata(i,"UT")+2 and not(gendata(i,"Lj")>24-gendata(i,"UT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"UT"))), U(i,h)) =g= gendata(i,"UT")*y(i,t);

* minimum down time constraints 
Dntime1(i)$(gendata(i,"Fj")>0)..
   sum(t$(ord(t)<(gendata(i,"Fj")+1)), U(i,t)) =e= 0;

Dntime2(i)$(gendata(i,"DT")>1)..
   sum(t$(ord(t)>24-gendata(i,"DT")+1), 1 - U(i,t) - z(i,t)) =g= 0;

Dntime3(i,t)$(ord(t)>gendata(i,"Fj") and ord(t)<24-gendata(i,"DT")+2 and not(gendata(i,"Fj")>24-gendata(i,"DT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"DT"))), 1-U(i,h)) =g= gendata(i,"DT")*z(i,t);
   
* start-up and shut-down constraints 
startc(i,t).. StC(i,t) =g= gendata(i,"costst")*y(i,t);
shtdnc(i,t).. SDC(i,t) =g= gendata(i,"CostsD")*z(i,t);

* eq. (5.3d) P lower >= Pmin which PL = Pmin + Pk
genconst1(i,h)..
   p(i,h)    =e= u(i,h)*gendata(i,"Pmin") + sum(k, Pk(i,h,k));

* Eq. (5.4a)
Genconst2(i,h)$(ord(h)>0)..
   U(i,h)    =e= U(i,h-1)$(ord(h)>1) + gendata(i,"Uini")$(ord(h)=1) + y(i,h) - z(i,h);
* Eq. (5.2a)
Genconst3(i,t,k)..
   Pk(i,t,k) =l= U(i,t)*data(k,i,'DP');
   

* ramp rate constraints 
Ramp1(i,t)..             p(i,t+1) - p(i,t) =l= gendata(i,'RU');;
Ramp2(i,t)..             p(i,t-1) - p(i,t) =l= gendata(i,'RD');;

* Demand-Generation Balance constraint
balance(t).. sum(i, p(i,t)) =e= dataLP(t,'load')-pv;

* Objective Function
Genconst4..
   costThermal =e= sum((i,t), StC(i,t) + SDC(i,t))
                +  sum((t,i), u(i,t)*gendata(i,'Mincost') + sum(k, data(k,i,'s')*pk(i,t,k)));
                
Model UCLP / all /;
file emp / '%emp.info%' /; put emp '* problem %gams.i%'/;
$onput
randvar pv beta 4 7

stage 2 pv 
stage 2 balance
$offput
putclose emp;

Set scen        Scenarios / s1*s6 /;
Parameter
    s_pv(scen)   pv by scenario;

Set dict / scen .scenario.''
           pv    .randvar .s_pv /;
           
option emp=lindo;
*option optCr = 0.0;
option reslim = 86400;
*option Savepoint=1;
option optcr=0.005;



solve UCLP min costThermal use emp scenario dict;

Display s_pv;
Fred
Posts: 372
Joined: 7 years ago

Re: stochastic programming

Post by Fred »

Hi,

Your stage 2 equation is an equality.

Code: Select all

balance(t).. sum(i, p(i,t)) =e= dataLP(t,'load')-pv;
With p being a stage 1 variable and pv being a stage 2 random variable I don't see how this is supposed to work.
Maybe you want something like

Code: Select all

balance(t).. sum(i, p(i,t)) =g= dataLP(t,'load')-pv;
I hope this helps!

Fred
anan
User
User
Posts: 20
Joined: 5 years ago

Re: stochastic programming

Post by anan »

Dear Fred,
sorry for the late reply,

thank you so much this worked perfectly.
all the load was satisfied by the units
but at time t23 the load is 208.93 and the generated power is 218.7528 which is much higher and affects the cost
logically after adding the PV the cosThermal should decrease not increase is there any other ways to solve this?


thank you,
Anan
Fred
Posts: 372
Joined: 7 years ago

Re: stochastic programming

Post by Fred »

Hi,

Not sure if I understand what you are asking for.

Of course, sum(i, p(i,t)) can be greater than dataLP(t,'load')-pv, even for all realizations of random variable pv for example because you link subsequent time steps via Ramp1 and Ramp2 equations.

Also not sure what is meant with after adding the PV the cosThermal should decrease not increase. costThermal should decrease compared to what? To the deterministic model? Note that you define pv=4 at the beginning but

Code: Select all

randvar pv beta 4 7 
gives you values between 0 and 1.

I hope this helps!

Fred
anan
User
User
Posts: 20
Joined: 5 years ago

Re: stochastic programming

Post by anan »

yes thank you for your help this was helpful

and yes I'm comparing the new result after adding PV with the deterministic solution without stochastic
I think I had a problem with the distribution fitting of the historical dataset giving me the wrong parameters of the Beta PDF

Best Regards,
Anan
anan
User
User
Posts: 20
Joined: 5 years ago

Re: stochastic programming

Post by anan »

My final code looks like this and it is working but I was trying to index PV over time and I couldn't figure out how to do it
could you help me, please?

Code: Select all

$ontext
this code for 6-bus system considering the line constraints when solved as deterministic the cost=110915

after solving the deterministic problem this code has PV uncertainity model considered in Irradiance variability
modeld as two stage stochastic
and considered as Weibull PDF
$offtext


Set
   t 'time'               / t1*t24   /
   i 'generators indices' / g1*g3   /
   k 'cost segments'      / sg1*sg20 /
   char                   / ch1*ch2  /
   l number of lines  /1*7/
   bus        / 1*6   /
   r          /pv1/;
Alias (t,h);
Alias (bus,node);
scalar  eff /0.75/
        irr /0/
        mw   /1000/;
        
Table gendata(i,*) 'generator cost characteristics and limits'
       bus Pmax Pmin a      b      c       CostsD  costst RU   RD   SU   SD   UT   DT   uini U0  S0 fuel
   g1  1   220  100  0.050  10     100     0       100    55   55   105  108  4    4    1    4   0  1
   g2  2   100  10   0.001  40.66  162     0       200    50   50   106  112  2    3    0    0   3  1
   g3  6   20   10   0.006  22.06  171     0       0      20   20   43   45   1    1    0    0   1  1;
   
Table WD(t,*)
        d
   t1   175.19
   t2   165.15
   t3   158.67
   t4   154.73
   t5   155.06
   t6   160.48
   t7   168.39
   t8   177.6
   t9   186.81
   t10  206.96
   t11  228.61
   t12  236.1
   t13  242.18
   t14  243.6
   t15  248.86
   t16  255.79
   t17  256
   t18  246.74
   t19  245.97
   t20  237.35
   t21  236.31
   t22  232.67
   t23  208.93
   t24  195.6;
   
Table branch(bus,node,*)
       x        limit  
1.2    0.17     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.14     100 ;

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

Set GB(bus,i)
/
   1.g1
   2.g2
   6.g3/;

Parameter pvarea(bus) / 5 43380 /;

* -----------------------------------------------------
* linearizing the quadratic fuel cost FC = a*p*p + b*p + c;
Parameter data(k,i,*);
data(k,i,'DP')       =(gendata(i,"Pmax") - gendata(i,"Pmin"))/card(k);
data(k,i,'Pini')     =(ord(k) - 1)*data(k,i,'DP') + gendata(i,"Pmin");
data(k,i,'Pfin')     = data(k,i,'Pini') + data(k,i,'DP');
data(k,i,'Cini')     = gendata(i,"a")*data(k,i,'Pini')*data(k,i,'Pini')
                     + gendata(i,"b")*data(k,i,'Pini') + gendata(i,"c");
data(k,i,'Cfin')     = gendata(i,"a")*data(k,i,'Pfin')*data(k,i,'Pfin')
                     + gendata(i,"b")*data(k,i,'Pfin') + gendata(i,"c");
data(k,i,'s')        =(data(k,i,'Cfin') - data(k,i,'Cini'))/data(k,i,'DP');
gendata(i,'Mincost') = gendata(i,'a')*gendata(i,"Pmin")*gendata(i,"Pmin")
                     + gendata(i,'b')*gendata(i,"Pmin") + gendata(i,'c');
* -----------------------------------------------------
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')       = 100*(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');

Parameter conex(bus,node);
conex(bus,node)$(branch(bus,node,'limit') and branch(node,bus,'limit')) = 1;
conex(bus,node)$(conex(node,bus)) = 1;

* eq. (5.5d) page 122
Parameter unit(i,char);
unit(i,'ch1') = 24;
unit(i,'ch2') =(gendata(i,'UT') - gendata(i,'U0'))*gendata(i,'Uini');
gendata(i,'Lj') = smin(char,unit(i,char));

*eq. (5.6d) page 123
Parameter unit2(i,char);
unit2(i,'ch1')  = 24;
unit2(i,'ch2')  =(gendata(i,'DT') - gendata(i,'S0'))*(1 - gendata(i,'Uini'));
gendata(i,'Fj') = smin(char,unit2(i,char));

Variable costThermal,delta(bus,t),Pij(bus,node,t);
Positive Variable pu(i,t), p(i,t), StC(i,t), SDC(i,t), Pk(i,t,k),pv;
Binary Variable u(i,t), y(i,t), z(i,t);

p.up(i,t)    = gendata(i,"Pmax");
p.lo(i,t)    = 0;
Pk.up(i,t,k) = data(k,i,'DP');
Pk.lo(i,t,k) = 0;
p.up(i,t)    = gendata(i,"Pmax");
pu.up(i,h)   = gendata(i,"Pmax");

Equation
   Uptime1, Uptime2, Uptime3, Dntime1, Dntime2, Dntime3, Ramp1, Ramp2,
   startc, shtdnc, genconst1, Genconst2, Genconst3, Genconst4,const1,const2,pvout;
   
*,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
pvout(bus).. pv(bus) =e= irr*eff*pvarea(bus)/mw;




const1(bus,node,t)$(conex(bus,node))..
   Pij(bus,node,t) =e= branch(bus,node,'bij')*(delta(bus,t) - delta(node,t));

const2(bus,t).. sum(i$GB(bus,i), p(i,t)) - sum(node$conex(node,bus), Pij(bus,node,t))=g= (BusData(bus,'pd')*WD(t,'d'))-(pv(bus));


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







* minimum up time constraints
Uptime1(i)$(gendata(i,"Lj")>0)..
   sum(t$(ord(t)<(gendata(i,"Lj")+1)), 1 - U(i,t)) =e= 0;

Uptime2(i)$(gendata(i,"UT")>1)..
   sum(t$(ord(t)>24-gendata(i,"UT")+1), U(i,t) - y(i,t)) =g= 0;

Uptime3(i,t)$(ord(t)>gendata(i,"Lj") and ord(t)<24-gendata(i,"UT")+2 and not(gendata(i,"Lj")>24-gendata(i,"UT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"UT"))), U(i,h)) =g= gendata(i,"UT")*y(i,t);

* minimum down time constraints 
Dntime1(i)$(gendata(i,"Fj")>0)..
   sum(t$(ord(t)<(gendata(i,"Fj")+1)), U(i,t)) =e= 0;

Dntime2(i)$(gendata(i,"DT")>1)..
   sum(t$(ord(t)>24-gendata(i,"DT")+1), 1 - U(i,t) - z(i,t)) =g= 0;

Dntime3(i,t)$(ord(t)>gendata(i,"Fj") and ord(t)<24-gendata(i,"DT")+2 and not(gendata(i,"Fj")>24-gendata(i,"DT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"DT"))), 1-U(i,h)) =g= gendata(i,"DT")*z(i,t);
   
* start-up and shut-down constraints 
startc(i,t).. StC(i,t) =g= gendata(i,"costst")*y(i,t);
shtdnc(i,t).. SDC(i,t) =g= gendata(i,"CostsD")*z(i,t);

* eq. (5.3d) P lower >= Pmin which PL = Pmin + Pk
genconst1(i,h)..
   p(i,h)    =e= u(i,h)*gendata(i,"Pmin") + sum(k, Pk(i,h,k));

* Eq. (5.4a)
Genconst2(i,h)$(ord(h)>0)..
   U(i,h)    =e= U(i,h-1)$(ord(h)>1) + gendata(i,"Uini")$(ord(h)=1) + y(i,h) - z(i,h);
* Eq. (5.2a)
Genconst3(i,t,k)..
   Pk(i,t,k) =l= U(i,t)*data(k,i,'DP');
   

* ramp rate constraints 
Ramp1(i,t)..             p(i,t+1) - p(i,t) =l= gendata(i,'RU');;
Ramp2(i,t)..             p(i,t-1) - p(i,t) =l= gendata(i,'RD');;


* Objective Function
Genconst4..
   costThermal =e= sum((i,t), StC(i,t) + SDC(i,t))
                +  sum((t,i), u(i,t)*gendata(i,'Mincost') + sum(k, data(k,i,'s')*pk(i,t,k)));
                
Model UCLP / all /;
file emp / '%emp.info%' /; put emp '* problem %gams.i%'/;
$onput
randvar irr weibull 3 0.78

stage 2 irr pv   
stage 2 const2 pvout
$offput
putclose emp;

Set scen        Scenarios / s1*s6 /;
Parameter
    s_irr(scen)   irr by scenarioo;
Set dict / scen .scenario.''
           irr    .randvar .s_irr /;
option emp=lindo;
*option optCr = 0.0;
option reslim = 86400;
*option Savepoint=1;
option optcr=0.005;



solve UCLP min costThermal use emp scenario dict;
Display s_irr;

Parameter report(t,bus,*),solar(bus,*);
report(t,bus,'Gen(MW)')    = 1*sum(i$GB(bus,i), P.l(i,t));
solar(bus,'PV output(MW)')    = pv.l(bus)
$ifI %system.fileSys%==Unix $exit
$call MSAppAvail Excel
$ifThen not errorLevel 1
   execute_unload "case2_3.gdx" report
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls par=report rng=report!A1'
   execute_unload "case2_3.gdx" P.l
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls var=P rng=PGen!A1'
   execute_unload "case2_3.gdx" solar
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls par=solar rng=PV!A1'
   execute_unload "case2_3.gdx" costThermal.l
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls var=costThermal rng=costThermal!A1'
$endIf
Fred
Posts: 372
Joined: 7 years ago

Re: stochastic programming

Post by Fred »

Where exactly are struggling? PV is a variable that is currently indexed with bus. If you want to index it over time as well, add a t to the index PV(bus,t) and make sure to adjust your code at all places where the variable occurs accordingly.

Regards,
Fred
anan
User
User
Posts: 20
Joined: 5 years ago

Re: stochastic programming

Post by anan »

thank you for your help.
I'm confused now becuase i want PV to be in equation Const2 as generation so I put everything in stage 2 to make the equation const 2 equality constraint but there is something wrong i cant find it

Code: Select all

$ontext
this code for 6-bus system considering the line constraints when solved as deterministic the cost=110915

after solving the deterministic problem this code has PV uncertainity model considered in Irradiance variability
modeld as two stage stochastic
and considered as Weibull PDF
$offtext


Set
   t 'time'               / t1*t24   /
   i 'generators indices' / g1*g3   /
   k 'cost segments'      / sg1*sg20 /
   char                   / ch1*ch2  /
   l number of lines  /1*7/
   bus        / 1*6   /
   r          /pv1/;
Alias (t,h);
Alias (bus,node);
scalar  eff /0.75/
        irr /0/
        mw   /1000/;
        
Table gendata(i,*) 'generator cost characteristics and limits'
       bus Pmax Pmin a      b      c       CostsD  costst RU   RD   SU   SD   UT   DT   uini U0  S0 fuel
   g1  1   220  100  0.050  10     100     0       100    55   55   105  108  4    4    1    4   0  1
   g2  2   100  10   0.001  40.66  162     0       200    50   50   106  112  2    3    0    0   3  1
   g3  6   20   10   0.006  22.06  171     0       0      20   20   43   45   1    1    0    0   1  1;
   
Table WD(t,*)
        d       solar
   t1   175.19  0
   t2   165.15  0
   t3   158.67  0
   t4   154.73  0
   t5   155.06  0
   t6   160.48  0
   t7   168.39  1
   t8   177.6   1
   t9   186.81  1
   t10  206.96  1
   t11  228.61  1
   t12  236.1   1
   t13  242.18  1
   t14  243.6   1
   t15  248.86  1
   t16  255.79  1
   t17  256     1
   t18  246.74  0
   t19  245.97  0
   t20  237.35  0
   t21  236.31  0
   t22  232.67  0
   t23  208.93  0
   t24  195.6   0   ;
   
Table branch(bus,node,*)
       x        limit  
1.2    0.17     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.14     100 ;

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

Set GB(bus,i)
/
   1.g1
   2.g2
   6.g3/;

Parameter pvarea(bus) / 1 10000 /;

* -----------------------------------------------------
* linearizing the quadratic fuel cost FC = a*p*p + b*p + c;
Parameter data(k,i,*);
data(k,i,'DP')       =(gendata(i,"Pmax") - gendata(i,"Pmin"))/card(k);
data(k,i,'Pini')     =(ord(k) - 1)*data(k,i,'DP') + gendata(i,"Pmin");
data(k,i,'Pfin')     = data(k,i,'Pini') + data(k,i,'DP');
data(k,i,'Cini')     = gendata(i,"a")*data(k,i,'Pini')*data(k,i,'Pini')
                     + gendata(i,"b")*data(k,i,'Pini') + gendata(i,"c");
data(k,i,'Cfin')     = gendata(i,"a")*data(k,i,'Pfin')*data(k,i,'Pfin')
                     + gendata(i,"b")*data(k,i,'Pfin') + gendata(i,"c");
data(k,i,'s')        =(data(k,i,'Cfin') - data(k,i,'Cini'))/data(k,i,'DP');
gendata(i,'Mincost') = gendata(i,'a')*gendata(i,"Pmin")*gendata(i,"Pmin")
                     + gendata(i,'b')*gendata(i,"Pmin") + gendata(i,'c');
* -----------------------------------------------------
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')       = 100*(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');

Parameter conex(bus,node);
conex(bus,node)$(branch(bus,node,'limit') and branch(node,bus,'limit')) = 1;
conex(bus,node)$(conex(node,bus)) = 1;

* eq. (5.5d) page 122
Parameter unit(i,char);
unit(i,'ch1') = 24;
unit(i,'ch2') =(gendata(i,'UT') - gendata(i,'U0'))*gendata(i,'Uini');
gendata(i,'Lj') = smin(char,unit(i,char));

*eq. (5.6d) page 123
Parameter unit2(i,char);
unit2(i,'ch1')  = 24;
unit2(i,'ch2')  =(gendata(i,'DT') - gendata(i,'S0'))*(1 - gendata(i,'Uini'));
gendata(i,'Fj') = smin(char,unit2(i,char));

Variable costThermal,delta(bus,t),Pij(bus,node,t);
Positive Variable pu(i,t), p(i,t), StC(i,t), SDC(i,t), Pk(i,t,k),pv(bus,t),ir;
Binary Variable u(i,t), y(i,t), z(i,t);

p.up(i,t)    = gendata(i,"Pmax");
p.lo(i,t)    = 0;
Pk.up(i,t,k) = data(k,i,'DP');
Pk.lo(i,t,k) = 0;
p.up(i,t)    = gendata(i,"Pmax");
pu.up(i,h)   = gendata(i,"Pmax");

Equation
   Uptime1, Uptime2, Uptime3, Dntime1, Dntime2, Dntime3, Ramp1, Ramp2,
   startc, shtdnc, genconst1, Genconst2, Genconst3, Genconst4,const1,const2,pvout;
   
*,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,


pvout(bus,t).. pv(bus,t)=e= irr*WD(t,'solar')*pvarea(bus)/mw;


const1(bus,node,t)$(conex(bus,node))..
   Pij(bus,node,t) =e= branch(bus,node,'bij')*(delta(bus,t) - delta(node,t));

const2(bus,t).. sum(i$GB(bus,i), p(i,t)) - sum(node$conex(node,bus), Pij(bus,node,t)) =e= BusData(bus,'pd')*(WD(t,'d'))- pv(bus,t);


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







* minimum up time constraints
Uptime1(i)$(gendata(i,"Lj")>0)..
   sum(t$(ord(t)<(gendata(i,"Lj")+1)), 1 - U(i,t)) =e= 0;

Uptime2(i)$(gendata(i,"UT")>1)..
   sum(t$(ord(t)>24-gendata(i,"UT")+1), U(i,t) - y(i,t)) =g= 0;

Uptime3(i,t)$(ord(t)>gendata(i,"Lj") and ord(t)<24-gendata(i,"UT")+2 and not(gendata(i,"Lj")>24-gendata(i,"UT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"UT"))), U(i,h)) =g= gendata(i,"UT")*y(i,t);

* minimum down time constraints 
Dntime1(i)$(gendata(i,"Fj")>0)..
   sum(t$(ord(t)<(gendata(i,"Fj")+1)), U(i,t)) =e= 0;

Dntime2(i)$(gendata(i,"DT")>1)..
   sum(t$(ord(t)>24-gendata(i,"DT")+1), 1 - U(i,t) - z(i,t)) =g= 0;

Dntime3(i,t)$(ord(t)>gendata(i,"Fj") and ord(t)<24-gendata(i,"DT")+2 and not(gendata(i,"Fj")>24-gendata(i,"DT")))..
   sum(h$((ord(h)>ord(t)-1) and (ord(h)<ord(t)+gendata(i,"DT"))), 1-U(i,h)) =g= gendata(i,"DT")*z(i,t);
   
* start-up and shut-down constraints 
startc(i,t).. StC(i,t) =g= gendata(i,"costst")*y(i,t);
shtdnc(i,t).. SDC(i,t) =g= gendata(i,"CostsD")*z(i,t);

* eq. (5.3d) P lower >= Pmin which PL = Pmin + Pk
genconst1(i,h)..
   p(i,h)    =e= u(i,h)*gendata(i,"Pmin") + sum(k, Pk(i,h,k));

* Eq. (5.4a)
Genconst2(i,h)$(ord(h)>0)..
   U(i,h)    =e= U(i,h-1)$(ord(h)>1) + gendata(i,"Uini")$(ord(h)=1) + y(i,h) - z(i,h);
* Eq. (5.2a)
Genconst3(i,t,k)..
   Pk(i,t,k) =l= U(i,t)*data(k,i,'DP');
   

* ramp rate constraints 
Ramp1(i,t)..             p(i,t+1) - p(i,t) =l= gendata(i,'RU');;
Ramp2(i,t)..             p(i,t-1) - p(i,t) =l= gendata(i,'RD');;


* Objective Function
Genconst4..
   costThermal =e= sum((i,t), StC(i,t) + SDC(i,t))
                +  sum((t,i), u(i,t)*gendata(i,'Mincost') + sum(k, data(k,i,'s')*pk(i,t,k)));
                
Model UCLP / all /;
file emp / '%emp.info%' /; put emp '* problem %gams.i%'/;
$onput
randvar irr weibull 3 0.78


stage 2 irr pv p Pij pk delta StC SDC u y z      
stage 2  pvout const2 const1 Ramp1 Ramp2 genconst1 genconst2 genconst3 Uptime2 Uptime3 Dntime2 Dntime3 startc shtdnc
$offput
putclose emp;

Set scen        Scenarios / s1*s2 /;
Parameter
    s_irr(scen)   irr by scenarioo;
Set dict / scen .scenario.''
           irr    .randvar .s_irr /;
option emp=lindo;
*option optCr = 0.0;
option reslim = 86400;
*option Savepoint=1;
option optcr=0.005;



solve UCLP min costThermal use emp scenario dict;
Display s_irr;
$ontext
Parameter report(t,bus,*),solar(bus,*);
report(t,bus,'Gen(MW)')    = 1*sum(i$GB(bus,i), P.l(i,t));
solar(bus,'PV output(MW)')    = pv.l(bus)
$ifI %system.fileSys%==Unix $exit
$call MSAppAvail Excel
$ifThen not errorLevel 1
   execute_unload "case2_3.gdx" report
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls par=report rng=report!A1'
   execute_unload "case2_3.gdx" P.l
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls var=P rng=PGen!A1'
   execute_unload "case2_3.gdx" solar
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls par=solar rng=PV!A1'
   execute_unload "case2_3.gdx" costThermal.l
   execute 'gdxxrw.exe  case2_3.gdx o=case2_3.xls var=costThermal rng=costThermal!A1'
$endIf
$offtext
Post Reply