Hello Dear GAMS World Forum Members,
I wrote a stochastic programming code with 2 random variables. There is no problem in here and now problem version. However, there is a huge problem in Expected Value (EV) problem. I tried different solutions to solve the issue, but I am not capable of solving. Any help is appreciated. Thanks in advance.
Nazmi
The error is listed below:
The following MIP errors were detected in model trial:
**** 784 Symbol s_mh in dimension three in scenario dictionary has no data
Code is given below:
solver options
option threads=32;
option iterlim=2e9;
option reslim=2e9;
option MIP=cplex;
option optcr=1e-4;
sets
i starthubs /125/
s scenarios /1*10/
;
alias (i,j,k,l);
parameter prob(s),c(i,k),ct(i,k),ort(i,j)
,h(i,j,s),fc(k),c1(k,l),o(i),d(i),ctm(i,k,s);
prob(s)=1/card(s);
parameter mh(i,j),cm(i,k);
$call gdxxrw hfa.xls par=ort rng=A1:Z26 Cdim=1 Rdim=1
$call gdxxrw flcd.xls par=c rng=A1:Z26 Cdim=1 Rdim=1
$gdxin hfa.gdx
$load ort
$gdxin
$gdxin flcd.gdx
$load c
$gdxin
scalar cc; scalar poisson ‘poisson variate’; scalar vw;
cc=5477;
loop(s,
cc=cc+ord(s)*1250;
execseed=round(cc,0);
loop((i,j),
poisson = -1; vw = 0; repeat vw = vw + [-log(uniform(0.001,1))]; poisson = poisson + 1;
until (vw >= ort(i,j)); h(i,j,s) = poisson;
); );
parameter ccp,ortc(i,k),sprime(i,k),sspma,o(i),d(i);
sspma=1;
loop((i,k)$(not sameas(i,k)),ortc(i,j) = log(sqr(c(i,k))/sqrt(sqr(c(i,k))+sqr(c(i,k)))); sprime(i,k) = sqrt(log(sqr((c(i,k)sspma)/c(i,k))+1));
);
ccp=5477;
execseed=round(ccp,0);
loop(s,
ctm(i,k,s) = exp(normal(0,1)(sprime(i,k))+ortc(i,k))/25000;
ctm(i,i,s) = 0;
ccp=ccp+ord(s)*1250;
);
o(i)=sum(j,ort(i,j));
d(i)=sum(j,ort(j,i));
mh(i,j)=sum(s,h(i,j,s)*prob(s));
cm(i,k)=sum(s,ctm(i,k,s));
fc(k)=10000;
*finding maximum link in the system
parameter delta,A(i,k);
delta=smax((i,k),c(i,k));
creating A matrix
loop((i,k),if(c(i,k)<(0.5delta),A(i,k)=1;else A(i,k)=0));
*A(i,k)=1;
scalar starttime; starttime = jnow;
positive variables
y(i,k,l)
z(i,k);
free variables
r sub problem objective
q master objective
;
binary variable
x(k);
equations
oeq
eneq1(i,k)
coneq(i,k)
coneq2(k)
eneq2(i,k)
eneq2a(k)
eneq3(i)
bieq(i,k);
oeq… q=e=sum(k,fc(k)x(k))
+0.8sum((i,k,l),cm(k,l)*y(i,k,l))
+sum((i,k),cm(i,k)z(i,k)(sum(j,mh(i,j))+sum(j,mh(j,i))));
eneq1(i,k)… z(i,k)=l=A(i,k)*x(k);
coneq(i,k)$(not sameas(i,k))… sum(j,mh(i,j))*z(i,k)=e=sum((j),mh(i,j)*z(j,k))
+sum((l),y(i,k,l))
-sum((l),y(i,l,k));
coneq2(k)… sum(j,mh(k,j))*x(k)=e=sum((j),mh(k,j)*z(j,k))
+sum((l),y(k,k,l))
-sum((l),y(k,l,k));
eneq2(i,k)(not sameas(i,k)).. sum((l)((not sameas(l,k))),y(i,k,l))=l=sum(j,mh(i,j))*z(i,k);
eneq2a(k)… sum((l)$((not sameas(l,k))),y(k,k,l))=l=sum(j,mh(k,j))*x(k);
eneq3(i)… sum(k$(not sameas(i,k)),z(i,k))+x(i)=g=1;
bieq(i,k)$(not sameas(i,k))… z(i,k)=l=1;
model trial /all/;
parameter
s_mh(s,i,j),
s_cm(s,i,k),
s_z(s,i,k) ,
s_y(s,i,k,l),
rep solution metric report ,
srep(s,*) scenario probability / #s.prob 0 /;
file emp / ‘%emp.info%’ /
put emp ‘* problem %gams.i%’
/‘ExpectedValue q EV_q’
/‘stage 2 z y mh cm eneq1 coneq coneq2 eneq2 eneq2a eneq3 bieq’
put / 'jrandvar ’
loop((i,j), put mh.tn(i,j) );
loop(s,
put /prob(s); loop((i,j), put h(i,j,s); ); );
put / 'jrandvar ’
loop((i,k), put cm.tn(i,k) );
loop(s,
put /prob(s); loop((i,k), put ctm(i,k,s); ); );
putclose emp;
option limrow=0, limcol=0;
Set dictSP / s. scenario.‘’
mh. randvar. s_mh
cm. randvar. s_cm/;
scalar starttime; starttime = jnow;
- Solve stochastic recourse problem or here-and-now
solve trial using emp min q scenario dictSP;
*Calculate time
scalar elapsed; elapsed = (jnow - starttime)243600;
rep(‘here-and-now’) = q.l;
parameter hubop; hubop=sum(k,x.l(k));
*Display solution of RP
display x.l,elapsed;
-
Solve expected value (EV) problem
solve trial using mip min q; -
Compute expected result of using the EV solution by evaluating first stage decision of EV solution under all scenarios:
x.fx(k) = x.l(k);
display mh;
Parameter
s_cost(s) cost by scenario;
Set dictMS / s .scenario.‘’
mh .param. s_mh
cm .param. s_cm
q .level. s_cost /;
solve trial using mip min q scenario dictMS;
rep(‘Expected Value of EV solution’) = sum(s, s_cost(s)*srep(s,‘prob’));
option rep:2:0:1 display rep;