set
j       /j1*j21/
k      /DN,PM,LR,MR,HR/
ttt/t0,t1*t5/
t(ttt)     /t1*t5/
v    /v1*v5/
kk(k)
kkk(k)/pm/
kkkk(k)/pm,lr/
kkkkk(k)/pm,lr,mr/
kkkkkk(k)/pm,lr,mr,hr/


alias(t,tt);
scalar
gr1/0.5306/
gr2 /1.3855/
iriideal/2.2/
*c1/0.0881/
*c2/0.2995/
w1  /0.5/
w2  /0.25/
w3  /0.25/
parameters
l(j)
s0(j)
delta(j)
gc(k)
gd(k)
y(kkkkkk)
AADT(j)
b(v,t)
delta(j)
AADTT(j)
maxz3
minz1
minz2


l(j)/j1 8,j2 2,j3 6,j4 12,j5 6,j6 4,j7 2,j8 3,j9 5,j10 2,j11 2,j12 5,j13 6,j14 4,j15 4,j16 2,j17 10,j18 6,j19 9,j20 4,j21 3.1/
s0(j)/
j1        80.1
j2        64.4
j3        50.1
j4        82.2
j5        61.4
j6        60.8
j7        86.3
j8        36.2
j9        76.9
j10       46.6
j11       81.2
j12       69.2
j13       70.4
j14       82.1
j15       75.9
j16       78.6
j17       84.8
j18       76.5
j19       92.4
j20       90.9
j21       94.9/
iri0(j)/
j1        3.4
j2        3.3
j3        3.9
j4        3.8
j5        3.8
j6        3.1
j7        3.2
j8        3.1
j9        5.0
j10       4.2
j11       2.4
j12       2.2
j13       2.6
j14       2.7
j15       2.5
j16       2.7
j17       2.3
j18       3.0
j19       1.9
j20       2.2
j21       2.4 /


gc(k)/DN 0,PM 5715,LR 17145,MR 28575,HR 57150/
gd(k)/DN 0,PM 0.317,LR 0.951,MR 1.585,HR 3.170/
c(k)/DN 0,PM 93,LR 1984,MR 4030,HR 8835/
y(kkkkkk)/PM 1, LR 1,MR 1,HR 1/
uc(k)/DN 0,PM 0,LR 1,MR 2,HR 3.5/ ;
delta(j)=uniform(0.5,1) ;
AADT(j) =2000;
AADTT(j)=200;

display delta;

table b(v,t)
           t1          t2          t3          t4          t5
v1        0.09        0.09        0.09        0.09        0.09
v2        0.65        0.65        0.65        0.65        0.65
v3        0.04        0.04        0.04        0.04        0.04
v4        0.04        0.02        0.02        0.02        0.02
v5        0.5         0.5         0.5         0.5          0.5    ;


table I(j,k)
        DN        PM        LR        MR        HR
j1      0         3         15        25        40
j2      0         3         15        25        40
j3      0         3         15        25        40
j4      0         3         15        25        40
j5      0         3         15        25        40
j6      0         3         15        25        40
j7      0         3         15        25        40
j8      0         3         15        25        40
j9      0         3         15        25        40
j10     0         3         15        25        40
j11     0         3         15        25        40
j12     0         3         15        25        40
j13     0         3         15        25        40
j14     0         3         15        25        40
j15     0         3         15        25        40
j16     0         3         15        25        40
j17     0         3         15        25        40
j18     0         3         15        25        40
j19     0         3         15        25        40
j20     0         3         15        25        40
j21     0         3         15        25        40 ;

TABLE IR(j,k)
          DN       PM         LR         MR         HR
j1        0        0.5        1.5        2.5        4
j2        0        0.5        1.5        2.5        4
j3        0        0.5        1.5        2.5        4
j4        0        0.5        1.5        2.5        4
j5        0        0.5        1.5        2.5        4
j6        0        0.5        1.5        2.5        4
j7        0        0.5        1.5        2.5        4
j8        0        0.5        1.5        2.5        4
j9        0        0.5        1.5        2.5        4
j10       0        0.5        1.5        2.5        4
j11       0        0.5        1.5        2.5        4
j12       0        0.5        1.5        2.5        4
j13       0        0.5        1.5        2.5        4
j14       0        0.5        1.5        2.5        4
j15       0        0.5        1.5        2.5        4
j16       0        0.5        1.5        2.5        4
j17       0        0.5        1.5        2.5        4
j18       0        0.5        1.5        2.5        4
j19       0        0.5        1.5        2.5        4
j20       0        0.5        1.5        2.5        4
j21       0        0.5        1.5        2.5        4;


scalar
h/21/;

parameter f(ttt),hh(ttt),l1(j,k,t);
f(ttt)=20000;
hh(ttt)=0.2*h;
l1(j,k,t)=0.0000001;



display hh;
binary variable
x(j,k,t),yy(ttt);

positive variable
s(j,ttt),k2(ttt)
iri(j,t)
oc(j)
dc(j,k)
EMC
ETD
ERR
mins(t)
;

*oc.fx(j)=uniformint(10,20);
s.lo(j,t)=1;
s.up(j,t)=100;
iri.lo(j,t)=0;
*iri.up(j,t)=1000000;

iri.m(j,t)=uniformint(0,100);
s.l(j,'t0')=s0(j);
k2.fx(ttt)=uniformint(1,19999);



free variable
z1
z2
z3
z4;





equation
obj1
obj2
obj3
eq22
eq33
eq44
eq66
eq77
eq1
eq2
eq3
eq4
eq5
eq6
eq7
eq8
eq9
eq10
eq11
eq12
eq13
newobjective
;

obj1  .. z1=e=EMC+ETD+ERR;

eq22  .. EMC=e=sum((j,k,t),gc(k)*x(j,k,t))  ;

eq33  .. ETD=e=sum((j,k,t),gd(k)*AADT(j))*sum((j,k,t),x(j,k,t));

eq44  .. ERR=e=sum((j,t),((gr1*(AADT(j)-AADTT(j)))+(gr2*AADTT(j)))*
(iriideal-iri(j,t)));

obj2     .. z2=e=sum((j,k,t),c(k)*l1(j,k,t))*sum((j,k,t),x(j,k,t));

*eq66(j,t)  .. oc(j)=e=(c1*(AADT(j)-AADTT(j))+c2*(AADTT(j))) ;

eq77(j,k).. Dc(j,k)=e=AADT(j)*UC(k);

obj3  .. z3=e=sum((j,k,t),L(j)*I(j,k)*x(j,k,t));

eq1(j,t)  .. s(j,t) =e=sum((k,tt)$(ord(tt) <= ord(t)),(x(j,k,tt)*I(j,k)*power(delta(j),(ord(t)-ord(tt)))))
+(s0(j)*(power(delta(j),ord(t))));

eq2(j,t)  .. iri(j,t) =e=sum((k,tt)$(ord(tt) <= ord(t)),(x(j,k,tt)*Ir(j,k)*power(delta(j),(ord(t)-ord(tt)))))
+(iri0(j)*(power(delta(j),ord(t))));


*eq2(j,t,k,ttt).. iri(j,t)=l=ir(j,k)**power(delta(j),(ord(t)-ord(ttt)))+iri(j,'t1')*(power(delta(j),(ord(t))));

eq3(v,t)  .. sum((j,kk,k),c(kk)*l1(j,k,t))=l=b(v,t);

eq4(kkkkkk,j)  .. sum(t,x(j,kkkkkk,t))=l=Y(kkkkkk);

eq5(j,t)  .. sum(k,x(j,k,t))=e=1;

*eq6(ttt).. f(ttt)=g=f(ttt+1);
eq7(ttt).. hh(ttt)*yy(ttt)=g=1;
*eq8(ttt).. f(ttt)=g=k2(ttt);
eq9(j)$(ord(j)>=50).. s0('j21')=e=94.9;

eq10(t).. mins(t)=g=sum(j,s(j,t))/21;
eq11(t)$(ord(t) < 5).. mins(t)=g=50;
eq12(t)$(ord(t) eq 5).. mins(t)=g=70;

eq13(t).. sum((j,k),x(j,k,t))=g=0.2*card(j);

newobjective  ..  z4=e=w1*(Z1/minz1)+w2*(z2/minz2)+w3*(-z3/maxz3);
loop((j,t),

if(s.l(j,t-1)<=100,
kk(k)=kkk(k);
elseif(s.l(j,t-1)>=90 and s.l(j,t)<=100),
kk(k)=kkkk(k)

elseif(s.l(j,t-1)>=70 and s.l(j,t)<90),
kk(k)=kkkk(k)

elseif(s.l(j,t-1)>=50 and s.l(j,t)<=70),
kk(k)=kkkkk(k)

elseif(s.l(j,t-1)>=35 and s.l(j,t)<=50),
kk(k)=kkkkk(k)

elseif(s.l(j,t-1)<=35),
kk(k)=kkkkkk(k)

));


model model1
/
obj1
eq22
eq33
eq44
*eq66
eq77
obj2
obj3
eq1
eq2
eq3
eq4
eq5
*eq6
eq7
*eq8
eq9
eq10
eq11
eq12
eq13
/;



option optcr=0,optca=0;
solve model1 us minlp min z1;
display
z1.l
x.l
*oc.l
Dc.l
s.l
iri.l
k2.l;




model model2
/
obj1
eq22
eq33
eq44
*eq66
eq77
obj2
obj3
eq1
eq2
eq3
eq4
eq5
*eq6
eq7
*eq8
eq9
eq10
eq11
eq12
eq13
/;

solve model2 us minlp min z2;
display
z2.l
x.l
*oc.l
Dc.l
s.l
iri.l
k2.l;


model model3
/
obj1
eq22
eq33
eq44
*eq66
eq77
obj2
obj3
eq1
eq2
eq3
eq4
eq5
*eq6
eq7
*eq8
eq9
eq10
eq11
eq12
eq13
/;
solve model3 us minlp max z3;
display
z3.l
x.l
*oc.l
Dc.l
s.l
iri.l
k2.l;




minz2=z2.l  ;
minz1=z1.l  ;
maxz3=z3.l;


model model4
/
obj1
eq22
eq33
eq44
*eq66
eq77
obj2
obj3
eq1
eq2
eq3
eq4
eq5
*eq6
eq7
*eq8
eq9
eq10
eq11
eq12
eq13
newobjective
/;

solve model4 us minlp max z4;
display
z4.l
x.l
*oc.l
Dc.l
s.l
iri.l
k2.l
mins.l;