Hi I am having troubles with the below script where the solution is infeasible and the model is always giving me problems is it possible to support me? Please find it below:
* Cleaned and Structured GAMS Bilevel Model with EMP Annotations *
$onempty
option iterlim = 500000;
*--------------------*
* SETS *
*--------------------*
sets
t /t1, t2, t3, t4/
subset(t) /t1, t2/
i /i1, i2/;
*--------------------*
* VARIABLES *
*--------------------*
variables
p(t), x(t), w(i,t), v(i,t), n(t), bk(i,t), brdf(i,t), bm(t), brd(t),
objin(i), objint(i,t), objout, objoutt(t), q(t), xk(i,t), xrdf(i,t),
er(i,t), em(t), inv(t), xhm(t), xrm(t), xrd(t), d(t);
positive variables p, x, w, v, n, bk, brdf, bm, brd, q, xk, xrdf, er, em, inv, xhm, xrm, xrd, d;
* Lower bounds to prevent division by zero *
q.lo(t) = 1e-6;
v.lo(i,t) = 1e-6;
*--------------------*
* PARAMETERS *
*--------------------*
parameters
alpha /0.5/, alpha_par /1/, beta /0.3/, beta_par /0.3/,
dbar(t) /t1*t4 1/, delta /0.6/, g /0/, g_par /0/,
l(t) /t1*t4 0/, l_par /0/, h /0.11/, h_par /1/,
sigma(t) /t1*t4 0.2/, phi(t) /t1*t4 0.15/,
eta(t) /t1*t4 0.05/, eta_par /1/, omega(t) /t1*t4 0.07/,
psi /0.2/, betah /0.1/, betard /0.1/,
gamma(i) /i1 0.6, i2 0.6/, gamma_par /1/,
gammafrd(i) /i1*i2 0.2/, gammam /0.5/, gammard /0.2/,
m(t) /t1*t4 0.5/, y /0.6/, y_par /1/, tax /0/,
subsidy_f /1/, subsidy_m /1/,
emm(t) /t1*t4 0.5/,
tbar(t) /t1*t4 0.6/, bbarm(t) /t1*t4 0.5/,
embar(t) /t1*t4 0.3/,
zbar(i,t) /i1.t1 0.2, i2.t1 0.2, i1.t2 0.2, i2.t2 0.2, i1.t3 0.2, i2.t3 0.2, i1.t4 0.2, i2.t4 0.2/,
s(i,t) /i1.t1 0.4, i2.t1 0.2, i1.t2 0.4, i2.t2 0.2, i1.t3 0.4, i2.t3 0.2, i1.t4 0.4, i2.t4 0.2/,
c(i,t) /i1.t1 0.3, i2.t1 0.4, i1.t2 0.3, i2.t2 0.4, i1.t3 0.3, i2.t3 0.4, i1.t4 0.3, i2.t4 0.4/,
e(i,t) /i1.t1 0.8, i2.t1 0.6, i1.t2 0.6, i2.t2 0.6, i1.t3 0.8, i2.t3 0.6, i1.t4 0.8, i2.t4 0/,
erbar(i,t) /i1.t1 0.5, i2.t1 0.5, i1.t2 0.5, i2.t2 0.5, i1.t3 0.5, i2.t3 0.5, i1.t4 0.5, i2.t4 0.5/,
bkbar(i,t) /i1.t1 0.75, i2.t1 0.5, i1.t2 0.75, i2.t2 0.5, i1.t3 0.75, i2.t3 0.5, i1.t4 0.75, i2.t4 0.5/,
betak(i) /i1 0.15, i2 0.1/, betakrd(i) /i1*i2 0.1/,
xo(i,t) /i1.t1 0.1, i2.t1 0.1, i1.t2 0.1, i2.t2 0.1, i1.t3 0.1, i2.t3 0.1, i1.t4 0.1, i2.t4 0.1/;
*--------------------------*
* EQUATION DEFINITIONS *
*--------------------------*
Equations
defout, defoutt(t), defin(i), defint(i,t),
e1(t), e2(i,t), e3(i,t), e4(i,t), e5(t),
e6(t), e7(t), e8(t), e9(t), e10(t),
e11(t), e12(t), e13(t), e14(t), e15(i,t),
e16(t), e17(t), e18(t), e19(i,t), e20(t),
e21(t), e22(t), e23(i,t), e24(i,t);
* Equations (define the logic of the model) *
defout.. objout =e= sum(t, objoutt(t));
defoutt(t).. objoutt(t) =e= p(t)*d(t) - (sigma(t)*q(t) + (1-gammam)*bm(t) + (1-gammard)*brd(t) + (l(t)+l_par)*em(t)*q(t));
defin(i).. objin(i) =e= sum(t, objint(i,t));
defint(i,t).. objint(i,t) =e= w(i,t)*v(i,t) - (1-gamma(i)*gamma_par)*bk(i,t) - (1-gammafrd(i))*brdf(i,t) - c(i,t)*v(i,t) - (g+g_par)*er(i,t)*v(i,t);
e1(t).. d(t) =e= dbar(t) - alpha*alpha_par*(p(t)*(1+tax)) + beta*beta_par*x(t);
e2(i,t)$(not t.first).. xk(i,t) =e= xk(i, t-1) + betak(i)*bk(i,t);
e19(i,t)$(t.first).. xk(i,t) =e= betak(i)*bk(i,t);
e3(i,t).. er(i,t)*v(i,t) =e= e(i,t)*v(i,t) - xk(i,t) - xo(i,t)*v(i,t);
e4(i,t).. (g*g_par*e(i,t) + c(i,t))*v(i,t) + (1-gamma(i)*gamma_par)*bk(i,t) + (1-gammafrd(i))*brdf(i,t) =l= zbar(i,t)*subsidy_f;
e5(t).. em(t)*q(t) =e= emm(t)*q(t) - xhm(t) - xrm(t) - xrd(t);
e6(t)$(not t.first).. xhm(t) =e= xhm(t-1) + betah*bm(t);
e20(t)$(t.first).. xhm(t) =e= betah*bm(t);
e16(t)$(t.first).. n(t) =e= eta(t)*eta_par*q(t);
e7(t)$(not t.first).. n(t) =e= eta(t)*eta_par*q(t) + omega(t-1)*inv(t-1);
e8(t)$(t.first).. xrm(t) =e= psi*(n(t)-n(t));
e18(t)$(not t.first).. xrm(t) =e= psi*(n(t-1)-n(t));
e9(t).. q(t) =e= sum(i, v(i,t));
e10(t).. (1-eta(t)*eta_par)*q(t) =g= d(t);
e11(t)$(t.first).. inv(t) =e= (1-eta(t)*eta_par)*q(t) - d(t);
e17(t)$(not t.first).. inv(t) =e= (1-omega(t-1))*inv(t-1) + (1-eta(t)*eta_par)*q(t) - d(t);
e12(t).. sum(i, w(i,t)*v(i,t)) + q(t)*(l(t)+l_par)*em(t) + (1-gammam)*bm(t) + (1-gammard)*brd(t) + sigma(t)*q(t) =l= m(t)*subsidy_m;
e13(t).. sum(i, gamma(i)*gamma_par*bk(i,t)) + gammam*bm(t) + gammard*brd(t) + sum(i, gammafrd(i)*brdf(i,t)) - q(t)*(l(t)+l_par)*em(t) - sum(i, g*g_par*e(i,t)*v(i,t)) =l= y*y_par;
e14(t).. x(t) =e= (xrm(t)+xhm(t)+xrd(t))/q(t) + sum(i, xk(i,t)/v(i,t) + xo(i,t));
e15(i,t).. w(i,t) =g= delta*p(t);
e21(t)$(not subset(t)).. xrd(t) =e= xrd(t-1) + betard*brd(t-2);
e22(t)$(subset(t)).. xrd(t) =e= 0;
e23(i,t)$(not subset(t)).. xrdf(i,t) =e= betakrd(i)*brdf(i,t-2);
e24(i,t)$(subset(t)).. xrdf(i,t) =e= 0;
*--------------------*
* MODEL & EMP *
*--------------------*
model dairy /all/;
$onecho > empinfo.dat
bilevel p
min objin(i) v(i,t) w(i,t) bk(i,t) brdf(i,t) defint(i,t) e2(i,t) e3(i,t) e4(i,t) e15(i,t) defin(i)
$offecho
solve dairy using emp maximizing objout;
*--------------------*
* OUTPUT RESULTS *
*--------------------*
file results /results.txt/;
put results;
put 'period' @15,'price_f1' @30,'price_f2' @45,'price_m' @60,'qprod_t' @75,'tot_emissP' @90,
'profit_f1' @105,'profit_f2' @120,'profit_m' @135,'emis_redf1' @150,'emis_red_f2' @165,
'emis_m' @180,'emis_tk' @195,'emis_org' @210,'emis_R&D' @225,'cost_emiss' @240,
'tot_emissF1' @255,'tot_emissF2' @270,'saction' @285,'invent_cost' @300,'subsidy_f' @315,
'subsidy_m' @330,'gov_budget' @345,'proc_rent1' @360,'proc_rent2' @375,'gov_techinn' @390,
'cost_emisret' @405,'waste' @420,'elas_price' @435,'elas_co2' @450,'qpord_f1' @465,
'qprod_f2' @480,'q_dem' @495,'tax'/;