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?
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;
);