the equation name is hydro_output and variable ph
Code: Select all
* this code for 6-bus system considering the line constraints with one hydro unit and one thermal unit
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 /
Hy number of hydro units /hy1/ ;
scalar zeta /14/ ;
scalar Q /1/;
Alias (t,h);
Alias (bus,node);
table datahy (Hy,*) hydropower generation coef.
h0 alpha eff
hy1 0.857 0.0001315 8.091;
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.0004 13.51 176.95 0 0 55 55 105 108 4 4 1 4 0 1
g2 2 100 10 0.001 32.63 129.97 0 0 50 50 106 112 2 3 0 0 3 1
g3 6 20 10 0.005 17.70 137.41 0 0 20 20 43 45 1 1 0 0 1 1;
Table charac(Hy,*)
lwmin lwmax Vini Vfin Qmin Qmax Pmin Pmax
hy1 1 10 1000 1000 1 10 10 150;
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 500
1.4 0.258 500
2.3 0.037 500
2.4 0.197 500
3.6 0.018 500
4.5 0.037 500
5.6 0.14 500 ;
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
3.g3
/;
set Bushy(bus,Hy) hydro of each bus
/
4.hy1/;
* -----------------------------------------------------
* 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),lw(Hy,t) , Spill(Hy,t),head(hy,t),ph(hy,t);
Binary Variable u(i,t), y(i,t), z(i,t);
* inequality constraints power limits for thermal and hydro, volume limits and water discharge limits
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");
lw.LO(Hy,t) = charac(Hy,'lwmin');
lw.UP(Hy,t) = charac(Hy,'lwmax');
ph.lo(Hy,t) = charac(Hy,'Pmin');
ph.up(Hy,t) = charac(Hy,'Pmax');
*Q.LO(Hy,t) = charac(Hy,'Qmin');
*Q.UP(Hy,t) = charac(Hy,'Qmax');
Spill.LO(Hy,t)= 0;
Equation
Uptime1, Uptime2, Uptime3, Dntime1, Dntime2, Dntime3, Ramp1, Ramp2,
startc, shtdnc, genconst1, Genconst2, Genconst3, Genconst4,const1,balance,waterhead,hydro_output;
*,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
* the hydro water head equation modeled as shahidehpour paper H = H0 + lw(t)
waterhead(Hy,t,bus)$Bushy(bus,Hy) .. head(Hy,t)=e= datahy(Hy,'h0')+lw(hy,t);
*waterhead(t).. head(t)=e= sum(hy,datahy(Hy,'h0')+lw(hy,t));
*+lw(hy,t)
* the hydro output power = Head * eff * Q
hydro_output(Hy,t,bus)$Bushy(bus,Hy).. Ph(hy,t)=e=Q*datahy(Hy,'eff')*head(hy,t);
*Q*datahy(Hy,'eff')*head(hy,t)
* ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
const1(bus,node,t)$(conex(bus,node))..
Pij(bus,node,t) =e= branch(bus,node,'bij')*(delta(bus,t) - delta(node,t));
balance(bus,t).. sum(i$GB(bus,i), p(i,t))+sum(Hy$Bushy(bus,Hy),ph(Hy,t))- WD(t,'d')*BusData(bus,'pd') =e= sum(node$conex(node,bus), Pij(bus,node,t));
*+sum(Hy$Bushy(bus,Hy),ph(Hy,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')*gendata(i,'fuel') + sum(k, data(k,i,'s')*pk(i,t,k)));
*zeta*sum((Hy,t),Ph(Hy,t)) hydro generation cost
Model UCLP / Uptime1, Uptime2, Uptime3, Dntime1, Dntime2, Dntime3, Ramp1, Ramp2,
startc, shtdnc, genconst1, Genconst2, Genconst3, Genconst4,const1,balance, waterhead,hydro_output/;
option optCr = 0.0;
*solve UCLP minimizing costThermal using RMINLP;
file emp / '%emp.info%' /;
put emp;
$onput
randvar Q discrete 0.25 1
0.25 2
0.25 4
0.25 6
* 0.1151555 0
* 0.3698302 10
* 0.1023296 20
* 0.935213 3
* 0.64016 10
* 0.426233 5
* 0.274741 6
* 0.350147 7
chance hydro_output 0.2
$offput
putclose emp;
Set scen scenarios / s1*s24 /;
Parameter
s_Q(scen)
s_head(scen,hy,t)
;
Set dict / scen .scenario.''
Q .randvar .s_Q
head .level .s_head
/;
solve UCLP min costThermal use emp scenario dict;
display s_Q,s_head;
