Considering the following network traffic models and assuming that there are 1000 drivers, how can we calculate the social optimum for the model, by assigning the number of drivers that uses each possible route.
I know the problem can be solved by linear programming using Pulp. However, my code is broken.
from pulp import *
prob = LpProblem("Network",LpMinimize)
a = LpVariable('A',lowBound=0,upBound=None,cat=LpInteger)
b = LpVariable('B',lowBound=0,upBound=None,cat=LpInteger)
c = LpVariable('C',lowBound=0,upBound=None,cat=LpInteger)
# objective function
prob += a * (((a+b)/100) + 20 + 5) + b * (((a+b)/100) + 10 + ((b+c)/100)) + c * (10 + 20 + ((b+c)/100))
#constraints
prob += a + b + c == 1000
prob.solve()
print("Status: ", LpStatus[prob.status])
for v in prob.variables():
print(v.name,'=',v.varValue)
Error:
TypeError: Non-constant expressions cannot be multiplied
Solution:
ABCF = 201
ABEF = 798
ADEF = 1
One option for removing the non-linear terms would be to replace the integer decision variables on the non-linear costed arcs with binary decision variables for each possible number of drivers along that arc (will scale badly - but no issue for this small problem).
Example code given below. Note that the solution this returns is different to the one you give in your question (I get ABCF=250
, ABEF=750
), I make this solution out to have objective value of 29375 VS the one you've quoted having cost of a little over 29399 - unless I've misunderstood?...
from pulp import *
prob = LpProblem("Network",LpMinimize)
nodes = ['A', 'B', 'C', 'D', 'E', 'F']
arc_costs = [[('A', 'D'), 10],
[('B', 'C'), 20],
[('B', 'E'), 10],
[('D', 'E'), 20],
[('C', 'F'), 5]]
ab_cost = LpVariable('ab_cost', lowBound=0, upBound=None, cat=LpContinuous)
ef_cost = LpVariable('ef_cost', lowBound=0, upBound=None, cat=LpContinuous)
# Binary variable == 1 iff number of drivers along arc = index value
ab_flow = LpVariable.dicts('ab_flow', range(1001), cat=LpBinary)
ef_flow = LpVariable.dicts('ef_flow', range(1001), cat=LpBinary)
# Variables to contain the selected Number of drivers
ab_flow_val = LpVariable('ab_flow_val', lowBound=0, upBound=None, cat=LpContinuous)
ef_flow_val = LpVariable('ef_flow_val', lowBound=0, upBound=None, cat=LpContinuous)
arc_flow_val = LpVariable.dicts('arc_flow_val', [i[0] for i in arc_costs],
lowBound=0, upBound=None, cat=LpInteger)
# objective function
prob += lpSum([arc_flow_val[i]*j for i, j in arc_costs] + ab_cost + ef_cost)
#constraints
# costs for the non-linear costed arcs:
prob += ab_cost == lpSum([ab_flow[x]*((x**2)/100) for x in range(1001)])
prob += ef_cost == lpSum([ef_flow[x]*((x**2)/100) for x in range(1001)])
# only one of binary variables can be true for each of the non-linear arcs
prob += lpSum([ab_flow[i] for i in range(1001)]) == 1
prob += lpSum([ef_flow[i] for i in range(1001)]) == 1
# set flow values from the binary variables:
prob += ab_flow_val == lpSum([ab_flow[i]*i for i in range(1001)])
prob += ef_flow_val == lpSum([ef_flow[i]*i for i in range(1001)])
# 1000 must leave and 1000 must arrive
prob += ab_flow_val + arc_flow_val[('A', 'D')] == 1000
prob += arc_flow_val[('C', 'F')] + ef_flow_val == 1000
# Non terminal nodes must have flow balance
for n in nodes[1:-1]:
if n == 'B':
prob += ab_flow_val == arc_flow_val[('B', 'C')] + arc_flow_val[('B', 'E')]
elif n == 'E':
prob += arc_flow_val[('B', 'E')] + arc_flow_val[('D', 'E')] == ef_flow_val
else:
# flow into node == flow out of node
prob += lpSum([arc_flow_val[i] for i in arc_flow_val.keys() if i[0] == n]) == lpSum([arc_flow_val[j] for j in arc_flow_val.keys() if j[1] == n])
prob.solve()
print("Status: ", LpStatus[prob.status])
for v in prob.variables():
if ('_val' in v.name) or ('cost' in v.name):
print(v.name,'=',v.varValue)