I'm trying to understand what's the proper way to make a recursive-dynamic model in GAMS.
For this I've made a toy example of 2 linear equations (aka cobweb model). Now it starts at random point P(t=1)=40
and move to equilibrium (solution of the system).
It has appeared to be tricky to run a straightforward code with index t
inside the loop. GAMS always complains (!):
Error 50: Loop controlling index appears in model equation(s)
I've been looking in other code snippets (like so), but have not found the complete example to reproduce. Playing with code I've just put another z
and loop with z
. I works ok.
So far so good, but note how t
is still inside the loop. I'm a bit confused. (At least I've no idea about side-effects of such coding.)
Set t /t1*t10/
* HERE IS THE TRICK
z /z1*z10/;
Parameters a /100/, b /2/, c /20/, d /1/;
Variables P(t), Qd(t), Qs(t), obj;
Positive Variables P, Qd, Qs;
Equations Demand(t), Supply(t), Equilibrium(t), Objective;
Demand(t).. Qd(t) =E= a - b * P(t);
Supply(t).. Qs(t) =E= c + d * P(t-1);
Equilibrium(t).. Qd(t) =E= Qs(t);
Objective.. obj =E= 1;
P.l('t1') = 40;
Model cobwebModel /all/;
Loop(z,
Solve cobwebModel using LP maximizing obj;
P.l(t+1) = (a - c + d * P.l(t)) / (b + d);
);
The question is: how to refactor the code above to make it proper GAMS style and ready for scaling?
As @Lutz correctly admitted, my trick was a computational disaster, as the loop is completely redundant here. When I've figured this out, decided to post another snippet. Irony is that GAMS uses t
to kind of loop through the code w/o a loop.
* This is a correct setting. No need to loop.
Set t /1*10/;
Parameters a /100/, b /2/, c /20/, d /1/;
Variables P(t), Qd(t), Qs(t), obj;
Positive Variables P, Qd, Qs;
Equations Demand(t), Supply(t), Equilibrium(t), Objective;
Demand(t).. Qd(t) =E= a - b * P(t);
Supply(t).. Qs(t) =E= c + d * P(t-1);
Equilibrium(t).. Qd(t) =E= Qs(t);
Objective..obj =E= 1;
* Initial point
P.l('1') = 40;
Model cobwebModel /all/;
Solve cobwebModel using LP maximizing obj;
Display P.l, Qd.l, Qs.l;
As a result it return the series of 10 values of variables that approach to equilibrium (from a starting point). Usual oscillation is in place.
I do not understand completely what you want to solve in the end, but if I got you right, you want to solve for just one element of t
at a time? In your code, that would not happen, instead, you solve 10 times for all 10 elements in t
. But you could use a dynamic subset, to do it differently, something like this:
Set t /t1*t10/
tt(t);
Alias (t,ta);
Parameters a /100/, b /2/, c /20/, d /1/;
Variables P(t), Qd(t), Qs(t), obj;
Positive Variables P, Qd, Qs;
Equations Demand(t), Supply(t), Equilibrium(t), Objective;
Demand(tt(t)).. Qd(t) =E= a - b * P(t);
Supply(tt(t)).. Qs(t) =E= c + d * P(t-1);
Equilibrium(tt(t)).. Qd(t) =E= Qs(t);
Objective.. obj =E= 1;
P.l('t1') = 40;
Model cobwebModel /all/;
Loop(ta,
tt(ta)=yes;
Solve cobwebModel using LP maximizing obj;
P.l(ta+1) = (a - c + d * P.l(ta)) / (b + d);
tt(ta)=no;
);