pythonpulpgame-theory

How to solve a Network Traffic Problem using Pulp?


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. Network Model

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


Solution

  • 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)