Stochastic Modelling with 2 Random Variables

Problems with modeling
Post Reply
nsener
User
User
Posts: 4
Joined: 7 years ago

Stochastic Modelling with 2 Random Variables

Post by nsener »

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 /1*25/
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.5*delta),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.8*sum((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)*24*3600;

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;
User avatar
bussieck
Moderator
Moderator
Posts: 1033
Joined: 7 years ago

Re: Stochastic Modelling with 2 Random Variables

Post by bussieck »

The error message indicates that you need to fill the parameter s_mh (and probably s_cm) with some values before the some. This is required by the compiler. Usually there is a clear correspondence between the values of the random variable you write and the scenario parameter you capture the values in. You write e.g. ctm to the EMP info file but collect s_mh. This is possible but in this case you need to initialize these parameters before the solve:

s_mh(s,i,j) = 0;
s_cm(s,i,k) = 0;

I like that you post the code but in order to run the model and reproduce the behavior for others with helps if you also post the data.

-Michael
Post Reply