Description
Three different models are used to compute traffic equilibria. These are a mixed complementarity formulation and a primal and dual formulation using NLPs.
Large Model of Types : MCP nlp
Category : GAMS Model library
Main file : traffic.gms
$title Traffic Equilibrium Problem (TRAFFIC,SEQ=169)
$onText
Three different models are used to compute traffic equilibria. These are
a mixed complementarity formulation and a primal and dual formulation
using NLPs.
Ferris, M C, Meeraus, A, and Rutherford, T F, Computing Wardropian
Equilibria in a Complementarity Framework. Optimization Methods and
Software 10 (1999), 669-685.
Keywords: mixed complementarity problem, nonlinear programming, traffic equilibria,
transportation, multicommodity flows
$offText
Set
n 'nodes'
a(n,n) 'directed arcs';
Parameter
trip(n,n) 'trip table'
ca (n,n) 'cost coef A'
cb (n,n) 'cost coef B'
ck (n,n) 'cost coef K';
Alias (n,i,j,k);
Variable
t(i,j) 'time to get from node i to node j'
v(i,j) 'time to traverse arc form i to j'
y(i,j,k) 'flow to k along arc i-j'
x(i,j) 'aggregate flow on arc i-j'
objpnlp 'objective for nlp formulation'
objdnlp 'objective for nlp formulation';
Positive Variable y;
Equation
balance(i,j) 'material balance'
vdef(i,j) 'arc travel time definition'
rational(i,j,k) 'cost minimization condition'
xdef(i,j) 'aggregate flow definition'
defpnlp 'defines objective for primal nlp formulation'
defdnlp 'defines objective for dual nlp formulation';
balance(i,k)$(not sameas(i,k)).. sum(a(i,j), y(a,k)) =e= sum(a(j,i), y(a,k)) + trip(i,k);
rational(a(i,j),k)$(not sameas(i,k)).. v(a) + t(j,k) =g= t(i,k);
vdef(a).. v(a) =e= ca(a) + cb(a)*power(x(a)/ck(a),4);
xdef(a).. x(a) =e= sum(k, y(a,k));
defpnlp.. objpnlp =e= sum(a, ca(a)*x(a) + cb(a)*power(x(a)/ck(a),5)*ck(a)/5);
defdnlp.. objdnlp =e= sum((i,k), trip(i,k)*t(i,k))
- sum(a, (4/5)*(ck(a)/cb(a)**(1/4))*(v(a) - ca(a))**1.25);
Model
pnlp 'primal nlp formulation' / defpnlp, balance, xdef /
dnlp 'dual nlp formulation' / defdnlp, rational /
mcp 'mcp formulation' / rational.y, balance.t, xdef, vdef.v /;
$sTitle Sioux Falls Test Data
$eolCom //
Set
n 'node definition' / 1*24 /
param 'arc cost table headers' / a, b, k /;
Alias (i,j,k,n);
Table arc_cost(n,n,param) 'arc cost data'
a b k
1.2 6 .90 25.9002
2.6 5 .75 4.9582
3.12 4 .60 23.4035
4.11 6 .90 4.9088
5.9 5 .75 10.0000
7.8 3 .45 7.8418
8.9 10 1.50 5.0502
9.10 3 .45 13.9158
10.15 6 .90 13.5120
10.17 8 1.20 4.9935
11.14 4 .60 4.8765
13.24 4 .60 5.0913
14.23 4 .60 4.9248
15.22 4 .60 10.3150
16.18 3 .45 19.6799
18.20 4 .60 23.4035
20.21 6 .90 5.0599
21.22 2 .30 5.2299
22.23 4 .60 5.0000
1.3 4 .60 23.4035
3.4 4 .60 17.1105
4.5 2 .30 17.7828
5.6 4 .60 4.9480
6.8 2 .30 4.8986
7.18 2 .30 23.4035
8.16 5 .75 5.0458
10.11 5 .75 10.0000
10.16 5 .75 5.1335
11.12 6 .90 4.9088
12.13 3 .45 25.9002
14.15 5 .75 5.1275
15.19 4 .60 15.6508
16.17 2 .30 5.2299
17.19 2 .30 4.8240
19.20 4 .60 5.0026
20.22 5 .75 5.0757
21.24 3 .45 4.8854
23.24 2 .30 5.0785;
arc_cost(i,j,param) $= arc_cost(j,i,param); // mirror the values
option arc_cost:4; display arc_cost; // check values
Parameter ca(i,j), cb(i,j), ck(i,j);
ca(i,j) = arc_cost(i,j,"a");
cb(i,j) = arc_cost(i,j,"b");
ck(i,j) = arc_cost(i,j,"k");
Table trip(i,j) 'symmetric trip matrix from i to j'
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
1 0 1 1 5 2 3 5 8 5 13 5 2 5 3 5 5 4 1 3 3 1 4 3 1
2 0 1 2 1 4 2 4 2 6 2 1 3 1 1 4 2 0 1 1 0 1 0 0
3 0 2 1 3 1 2 1 3 3 2 1 1 1 2 1 0 0 0 0 1 1 0
4 0 5 4 4 7 7 12 14 6 6 5 5 8 5 1 2 3 2 4 5 2
5 0 2 2 5 8 10 5 2 2 1 2 5 2 0 1 1 1 2 1 0
6 0 4 8 4 8 4 2 2 1 2 9 5 1 2 3 1 2 1 1
7 0 10 6 19 5 7 4 2 5 14 10 2 4 5 2 5 2 1
8 0 8 16 8 6 6 4 6 22 14 3 7 9 4 5 3 2
9 0 28 14 6 6 6 9 14 9 2 4 6 3 7 5 2
10 0 40 20 19 21 40 44 39 7 18 25 12 26 18 8
11 0 14 10 16 14 14 10 1 4 6 4 11 13 6
12 0 13 7 7 7 6 2 3 4 3 7 7 5
13 0 6 7 6 5 1 3 6 6 13 8 8
14 0 13 7 7 1 3 5 4 12 11 4
15 0 12 15 2 8 11 8 26 10 4
16 0 28 5 13 16 6 12 5 3
17 0 6 17 17 6 17 6 3
18 0 3 4 1 3 1 0
19 0 12 4 12 3 1
20 0 12 24 7 4
21 0 18 7 5
22 0 21 11
23 0 7
24 0;
trip(i,j) $= trip(j,i); // mirror the values
trip(i,j) = trip(i,j)*0.11; // get it back to original values
option trip:2; display trip; // has to match the values in the article
a(i,j) = ca(i,j); // identify arcs using flow cost parameter a
display a;
$sTitle Bound Definitions and Solutions
t.fx(i,i) = 0;
v.lo(a) = ca(a);
y.fx(a(i,j),i) = 0;
option domLim = 100000;
pnlp.workFactor = 3;
Parameter rep(i,k,*) 'summary report';
solve mcp using mcp; rep(i,j,'mcp ') = t.l(i,j);
solve pnlp using nlp minimizing objpnlp; rep(i,j,'primal') = balance.m(i,j);
solve dnlp using nlp maximizing objdnlp; rep(i,j,'dual ') = t.l(i,j);
display rep;