$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;