stochastic programming

anan · May 15, 2019, 2:31pm

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

* 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_Fiand · May 16, 2019, 7:13am

Hi,

Your stage 2 equation is an equality.

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

balance(t).. sum(i, p(i,t)) =g= dataLP(t,'load')-pv;

I hope this helps!

Fred

anan · May 19, 2019, 3:56pm

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_Fiand · May 20, 2019, 10:50am

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

randvar pv beta 4 7

gives you values between 0 and 1.

I hope this helps!

Fred

anan · May 21, 2019, 1:57am

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 · May 28, 2019, 6:50am

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?

$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_Fiand · May 28, 2019, 1:37pm

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 · May 30, 2019, 4:36pm

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

$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