# bertsekas.mod
# 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.
set Nodes circular; # nodes in the network
set Types; # types of arcs
param demand {Nodes}; # flow demand, i -> i+3
param c {Nodes,Types}; # delay coefficients for links in the network
param gamma; # gamma is a coupling coefficient. The delay on some of
# the links is increased gamma * delay on merging links.
var x {Nodes}; # outside flow
var z {Nodes}; # inside flow
var y {Nodes}; # dual vars
var fx {i in Nodes} = 1 + x[i] + x[i]^2;
var fz {i in Nodes} = 1 + z[i] + z[i]^2;
s.t. delay_o {i in Nodes}: 0 <= x[i]
complements
c[i,"oon"]*fx[i]
+ gamma * (1 + x[next(i,Nodes,3)] + x[next(i,Nodes,4)]
+ (x[next(i,Nodes,3)]+x[next(i,Nodes,4)])^2)
+ 10 * c[i,"ohy"]*(1 + x[i] + x[next(i,Nodes,3)]
+ x[next(i,Nodes,4)]
+ (x[i] + x[next(i,Nodes,3)]
+ x[next(i,Nodes,4)])^2)
+ 2 * gamma * fx[next(i,Nodes,3)]
+ c[next(i),"oby"] * (1 + x[i] + x[next(i,Nodes,4)]
+ (x[i] + x[next(i,Nodes,4)])^2)
+ 10 * c[next(i),"ohy"]
* (1 + x[i] + x[next(i)] + x[next(i,Nodes,4)]
+ (x[i] + x[next(i)] + x[next(i,Nodes,4)])^2)
+ 2 * gamma * fx[next(i,Nodes,4)]
+ c[next(i,Nodes,2),"oby"]
* (1 + x[i] + x[next(i)] + (x[i]+x[next(i)])^2)
+ 10 * c[next(i,Nodes,2),"ohy"]
* (1 + x[i] + x[next(i)] + x[next(i,Nodes,2)]
+ (x[i] + x[next(i)] + x[next(i,Nodes,2)])^2)
+ 2 * gamma * fx[i]
+ c[next(i,Nodes,3),"ooff"] * fx[i]
>= y[i];
s.t. delay_i {i in Nodes}: 0 <= z[i]
complements
c[i,"ion"]*fz[i]
+ gamma * fz[next(i)]
+ 10 * c[i,"ihy"]*(1 + z[i] + z[next(i)] + (z[i] + z[next(i)])^2)
+ 2 * gamma * fz[next(i)]
+ c[prev(i),"iby"]*fz[i]
+ 10 * c[prev(i),"ihy"] * (1 + z[i] + z[next(i,Nodes,4)]
+ (z[i]+z[next(i,Nodes,4)])^2)
+ 2 * gamma * fz[i]
+ c[prev(i,Nodes,2),"ioff"] * fz[i]
>= y[i];
s.t. dem_cons {i in Nodes}: x[i] + z[i] == demand[i];