Bilevel optimization constraints

GAMS
1

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'/;
2

This code does not run due to some formatting issues. Please paste multiline code into triple backticks (```). The backticks are also inserted automatically when clicking the preformatted text (</>) icon. Thanks!

3

Hi Fred,

Thanks for reaching out, I am able to run it however, it is showing some infeasible solutions. I tried to relax some of the constraints but still I am unable to get feasible results. Do you think I use a model other than EMP? Thanks!

4

I am sure you are. But when you copy/paste the code into this forum without enclosing it properly between code tags, the formatting breaks. When other users copy it from the forum, they will get lots of compilation errors. One could manually fix this but if you want to get help, I suggest you make it as easy as possible for others to run your code.

5

Oh ok, I just edited the script. Thank you for highlighting this. That is my first time using the forum.

6

Now you put the falsely formatted code in code tags which still doesn’t help. I suggest to just copy from your gms file and paste.

7

Hopefully it should work now!

1 Like
8

The empinfo specification has multiple issues:

  • By default, the JAMS solver expects the empinfo.dat to be in the scratch directory. This is not the case with the current statement. Using the magic variable %emp.info% achieves that purpose.
  • At the lower level, there is an issue with the objective variable: it is not scalar. Does your model have multiple followers, possible in a Nash equilibrium problem? If yes, each lower level agent must be specified. Currently it is not possible to do that in a static fashion, one needs loop and put statements.
  • While performing this, I noticed that e2 was only defined on a subset of its domain. This needs to be mirrored in the empinfo file, otherwise JAMS will (rightfully) complain about a non-existent equation. To alleviate the burden of keeping the conditional in sync, I would recommend using a dynamic set to restrict the definition of e2 and reuse it in the code defining the empinfo file

Substituting the $onEcho/$OffEchopart with the following results in a model that runs without error. No guarantee on the correctness of the model.

file emp / '%emp.info%' /;
put emp /'bilevel p'/;
loop(i,
   put 'min' objin(i)
   loop(t, put v(i,t) w(i,t) bk(i,t) brdf(i,t);)
   loop(t, put defint(i,t) e3(i,t) e4(i,t) e15(i,t);)
   loop(t$(not t.first), put e2(i,t);)
   put defin(i) /;
)
putclose emp;
1 Like
9

Hi Olivier Huber,

Thank you for reaching out and for supporting me in this, following your comments I have managed to reach to a feasible model, this means that I am on the right track hopefully. I am sharing below the script that is feasible now. I appreciate if you can get the chance to check it again. Thanks again!

* Cleaned and Structured GAMS Bilevel Model with EMP Annotations *

$onempty
option iterlim = 500000;
option limrow = 0, limcol = 0, solprint = on, sysout = on;

*--------------------*
*       SETS        *
*--------------------*
sets
    t /t1, t2, t3, t4/
    subset(t) /t1, t2/
    i /i1, i2/;

alias (t,tt);

set tdef(i,t) /i1.t2*t4, i2.t2*t4/;

*--------------------*
*    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 5/, 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 5/, y /2/, y_par /1/, tax /0/,
    subsidy_f /10/, subsidy_m /10/,
    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.2/,
    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);

* Objective equations *
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);

* Demand and production *
e1(t)..           d(t) =e= dbar(t) - alpha*alpha_par*(p(t)*(1+tax)) + beta*beta_par*x(t);
e2(i,t)$tdef(i,t).. 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);

* Emissions and investments *
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);

* Capital Accumulation *
e6(t)$((not t.first)).. xhm(t) =e= xhm(t-1) + betah*bm(t);
e20(t)$t.first..        xhm(t) =e= betah*bm(t);

* Inventory logic *
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);

* R&D logic *
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));

* Aggregation and feasibility *
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);

* Budget constraints *
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;

* Market mechanics *
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);

* Tech investment lags *
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/;

file emp / '%emp.info%' /;
put emp /'bilevel p'/;
loop(i,
   put 'min' objin(i);
   loop(t, put v(i,t) w(i,t) bk(i,t) brdf(i,t););
   loop(t, put defint(i,t) e3(i,t) e4(i,t) e15(i,t););
   loop(t$(tdef(i,t)), put e2(i,t););
   put defin(i) /;
)
putclose emp;

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'/;
10

After a cursory look, I see no obvious error. For the overall correctness of the model, an expertise in the application area is required.

1 Like
11

Thank you Olivier, do you know an expert that can help me by any chance?

12

As mentioned before, expertise in the application area is required to validate a model. I think it will be hard to find someone who just does it. If you want to talk with others about your model, a good start would be to use meaningful symbol names and/or explanatory text and/or comments to increased the chance that others understand what this is all about. A mathematical description of the model (in a pdf) goes a long way towards this goal as well. Then, it usually makes sense to ask specific questions and not just “can you confirm that my model is correct?”.

If you have a valid license, you might want to follow the advice and then either ask support@gams.com for specific help, or, if you consider this being a bigger project, reach out to consulting@gams.com.