*==>freebert.gms
$offsymxref offsymlist offuellist offuelxref
$TITLE Bertsekas problem as a 15-variable MCP
$ontext
This problem is based on a traffic assigment problem, found in:
D.P. Bertsekas and E.M. Gafni, "Projection Methods for Variational
Inequalities with Application to the Traffic Assignment Problem",
Mathematical Programming Study 17(1982), 139-159.
Just so you know, x(i) is the flow (from node(i) to node(i++3)) on the
(longer) outside loop, z(i) the flow on the inside loop, and y the
dual variables.
The source problem here is a VI, where the constraints consist of
non-negative bounds and x(I) + z(I) = demand(I). Without much
fiddling, we can write the problem as a MCP.
We can write the problem in (at least) 2 different ways. In the
first, y is free, and the demand constraints are satisfied exactly.
In the second, we change the demand constraint to an inequality, where
extra flow is allowed. This doesn't occur in solutions, and we can
write the problem as an NCP.
$offtext
sets
run / run1 * run7 /,
I nodes in the network / 1 * 5 /,
types types of arcs / ihy, ion, ioff, iby, ohy, oon, ooff, oby / ;
parameters
includerun (run) /
run1 01
run2 01
run3 01
run4 01
run5 01
run6 01
run7 01
/,
demand(I),
c (I, types);
scalar gamma;
$include 'bert_gaf.dat'
* demand and gamma data included in bertsekas and gafni files
parameters
ipty(I,run),
iptx(I,run),
iptz(I,run);
table iptx(I,run)
run1 run2 run3 run4 run5 run6
1
2
3
4
5 ;
table iptz(I,run)
run1 run2 run3 run4 run5 run6
1
2
3
4
5 ;
iptz (I,"run1") = demand (I) - iptx (I,"run1");
iptz (I,"run2") = demand (I) - iptx (I,"run2");
iptz (I,"run3") = demand (I) - iptx (I,"run3");
ipty (I,"run2") = 4000;
ipty (I,"run5") = 4000;
ipty (I,"run3") = 1;
ipty (I,"run6") = 1;
iptx (I,"run7") = 0.5;
ipty (I,"run7") = 0.5;
iptz (I,"run7") = 0.5;
positive variable
x(I) outside flow,
z(I) inside flow;
free variable
y(I) dual vars;
equations
delay_o(I),
delay_i(I),
dem_cons(I);
delay_o(I) ..
c(I,"oon")*(1 + x(I) + power(x(I),2))
+ gamma * (1 + x(I++3) + x(I++4) + power(x(I++3)+x(I++4),2))
+ 10 * c(I,"ohy")*(1 + x(I) + x(I++3) + x(I++4)
+ power(x(I) + x(I++3) + x(I++4),2))
+ 2 * gamma * (1 + x(I++3) + power(x(I++3),2))
+ c(I++1,"oby") * (1 + x(I) + x(I++4) + power(x(I) + x(I++4),2))
+ 10 * c(I++1,"ohy") * (1 + x(I) + x(I++1) + x(I++4)
+ power(x(I) + x(I++1) + x(I++4), 2))
+ 2 * gamma * (1 + x(I++4) + power(x(I++4),2))
+ c(I++2,"oby") * (1 + x(I) + x(I++1) + power(x(I)+x(I++1),2))
+ 10 * c(I++2,"ohy") * (1 + x(I) + x(I++1) + x(I++2)
+ power(x(I) + x(I++1) + x(I++2), 2))
+ 2 * gamma * (1 + x(I) + power(x(I),2))
+ c(I++3,"ooff") * (1 + x(I) + power(x(I),2))
=g= y(I);
delay_i(I) ..
c(I,"ion")*(1 + z(I) + power(z(I),2))
+ gamma * (1 + z(I++1) + power(z(I++1),2))
+ 10 * c(I,"ihy")*(1 + z(I) + z(I++1) + power(z(I) + z(I++1),2))
+ 2 * gamma * (1 + z(I++1) + power(z(I++1),2))
+ c(I--1,"iby")*(1 + z(I) + power(z(I),2))
+ 10 * c(I--1,"ihy") * (1 + z(I) + z(I++4) + power(z(I)+z(I++4), 2))
+ 2 * gamma * (1 + z(I) + power(z(I),2))
+ c(I--2,"ioff") * (1 + z(I) + power(z(I),2))
=g= y(I);
dem_cons (I) .. x(I) + z(I) =e= demand(I);
MODEL freebert / delay_o.x, delay_i.z, dem_cons.y /;
option limrow = 0;
option limcol = 0;
option iterlim = 1000;
option reslim = 120;
loop (run,
if (includerun(run),
x.l(I) = iptx(I,run);
z.l(I) = iptz(I,run);
y.l(I) = ipty(I,run);
solve freebert using mcp;
));