python-3.xscipylinear-programming
Translating a mixed-integer programming formulation to Scipy
Problem Summary
To match passengers (the number of passengers) to capacitated vehicles such that the profit is increased. All the vehicles have the same capacity c. It is not important to track which passenger is matched to which vehicle.
Problem Background
I'm trying to find a solution to the following passenger-matching problem:
The network is represented by graph G=(V,E). Graph G is a complete graph, The Edge e_ij is the edge between nodes i and j. V is the set of nodes/stations. p_{ij} is the profit of traveling through an edge (i,j). Let N be the number of vehicles, all the vehicles have the same capacity c.
At each (discreet) time step, the passengers arrive at their origin station i, and need to be transported to their destination station j. d_{ij} is the demand. f_{ij} <= d_{ij} is the passenger flow, i.e., the number of passengers that are traveling from i to j at time t and are successfully matched to a vehicle. The unmatched passers will leave the system. X_i is the total number of available vehicles at station i at time t.
Objective
The objective is to maximize profit.
I would like to solve the above formulation in Scipy and solve it using milp(). Eq. (2) ensures that the number of vehicles leaving node i is not more than the available vehicles at this node. Equations 3 and 4 bounds the flow to the maximum demand and capacity of the vehicles.
I have difficulty in translating the formulation to Scipy milp code. I would appreciate it if anyone could give me some pointers.
What I have done:
The code for equation (1):
f_obj = [p[i] for i in Edge]
x_obj = [0]*len(Edge)
obj = f_obj + v_obj
Integrality:
f_cont = [0 for i in Edge] # continous
x_int = [1]*len(Edge) # integer
integrality = f_cont + x_int
Equation (2):
def constraints(self):
b = []
A = []
const = [0]*len(Edge) # for f_ij
for i in region: # for x_ij
for e in Edge:
if e[0] == i:
const.append(1)
else:
const.append(0)
A.append(const)
b.append(self.accInit[i])
const = [0]*len(Edge) # for f_ij
return A, b
Equation (4):
[(0, demand[e]) for e in Edge]
Update
CPLEX code:
# Flow[e] == f_ij
# vehicleCount[e] == x_ij
maximize(sum(e in Edge) Flow[e]*prices[e]);
subject to
{
forall(e in Edge)
Flow[e] <= demand[e];
forall(i in region)
sum(e in Edge: e.i==i)vehicleCount[e] <= init_car[i];
forall(e in Edge)
Flow[e] <= capacity*vehicleCount[e];
}
Minimum sample data:
Edge = [(0, 1), (0, 5), (0, 6), (1, 0), (1, 3), (1, 5), (1, 7), (1, 10), (1, 13), (2, 0), (2, 1), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (3, 1), (3, 2), (3, 5), (3, 6), (3, 9), (3, 11), (3, 13), (4, 2), (4, 3), (4, 5), (4, 7), (4, 8), (4, 11), (4, 13), (5, 0), (5, 1), (5, 2), (5, 4), (5, 6), (5, 7), (5, 8), (5, 10), (5, 12), (6, 1), (6, 2), (6, 3), (6, 5), (6, 7), (6, 9), (6, 10), (6, 13), (7, 0), (7, 2), (7, 5), (7, 6), (7, 9), (7, 10), (7, 11), (8, 4), (8, 5), (8, 6), (8, 9), (8, 12), (8, 13), (8, 14), (9, 0), (9, 2), (9, 4), (9, 5), (9, 6), (9, 7), (9, 11), (9, 12), (9, 14), (9, 15), (10, 3), (10, 4), (10, 9), (10, 14), (11, 3), (11, 6), (11, 7), (11, 9), (11, 13), (11, 15), (12, 4), (12, 9), (12, 10), (13, 6), (13, 7), (13, 8), (13, 9), (13, 15), (14, 6), (14, 7), (14, 15), (15, 2), (15, 3), (15, 7), (15, 10), (15, 13), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (9, 9)]
# demand[Edge]={(src, dst): demand}
demand = {(0, 1): 1, (0, 5): 1, (0, 6): 1, (1, 0): 1, (1, 3): 1, (1, 5): 2, (1, 7): 1, (1, 10): 1, (1, 13): 1, (2, 0): 2, (2, 1): 1, (2, 3): 3, (2, 4): 1, (2, 5): 1, (2, 6): 5, (2, 7): 2, (2, 8): 2, (2, 9): 3, (2, 10): 2, (2, 11): 2, (2, 12): 1, (3, 1): 1, (3, 2): 1, (3, 5): 1, (3, 6): 1, (3, 9): 2, (3, 11): 1, (3, 13): 1, (4, 2): 1, (4, 3): 1, (4, 5): 5, (4, 7): 4, (4, 8): 1, (4, 11): 3, (4, 13): 1, (5, 0): 1, (5, 1): 2, (5, 2): 1, (5, 4): 3, (5, 6): 3, (5, 7): 2, (5, 8): 1, (5, 10): 1, (5, 12): 1, (6, 1): 2, (6, 2): 1, (6, 3): 1, (6, 5): 1, (6, 7): 5, (6, 9): 4, (6, 10): 1, (6, 13): 1, (7, 0): 1, (7, 2): 1, (7, 5): 1, (7, 6): 3, (7, 9): 4, (7, 10): 1, (7, 11): 1, (8, 4): 2, (8, 5): 2, (8, 6): 1, (8, 9): 2, (8, 12): 2, (8, 13): 2, (8, 14): 1, (9, 0): 2, (9, 2): 2, (9, 4): 3, (9, 5): 5, (9, 6): 1, (9, 7): 1, (9, 11): 1, (9, 12): 3, (9, 14): 1, (9, 15): 2, (10, 3): 1, (10, 4): 1, (10, 9): 1, (10, 14): 1, (11, 3): 1, (11, 6): 6, (11, 7): 4, (11, 9): 1, (11, 13): 1, (11, 15): 1, (12, 4): 1, (12, 9): 1, (12, 10): 1, (13, 6): 1, (13, 7): 1, (13, 8): 1, (13, 9): 1, (13, 15): 1, (14, 6): 1, (14, 7): 2, (14, 15): 1, (15, 2): 1, (15, 3): 1, (15, 7): 1, (15, 10): 1, (15, 13): 2, (1, 1): 2, (2, 2): 3, (3, 3): 4, (4, 4): 1, (5, 5): 1, (6, 6): 2, (7, 7): 4, (9, 9): 1}
# prices[Edge] = {(src, dst): price}
prices = {(0, 1): 7.1, (0, 5): 8.626999999999999, (0, 6): 11.568749999999998, (1, 0): 8.120000000000001, (1, 3): 8.425, (1, 5): 10.23125, (1, 7): 11.500000000000004, (1, 10): 13.633333333333331, (1, 13): 13.999999999999998, (2, 0): 7.558333333333335, (2, 1): 8.733333333333333, (2, 3): 6.610999999999999, (2, 4): 8.899999999999999, (2, 5): 9.883333333333331, (2, 6): 9.159523809523812, (2, 7): 9.978749999999998, (2, 8): 12.54, (2, 9): 11.6625, (2, 10): 13.3, (2, 11): 12.800000000000002, (2, 12): 15.6, (3, 1): 9.136000000000001, (3, 2): 6.163103448275866, (3, 5): 13.050000000000006, (3, 6): 10.143928571428571, (3, 9): 15.787499999999998, (3, 11): 11.0, (3, 13): 20.4, (4, 2): 9.042499999999999, (4, 3): 7.249999999999998, (4, 5): 7.123913043478264, (4, 7): 11.058518518518518, (4, 8): 6.87142857142857, (4, 11): 14.36714285714286, (4, 13): 8.120000000000001, (5, 0): 6.35, (5, 1): 8.975000000000001, (5, 2): 8.329999999999998, (5, 4): 6.177419354838713, (5, 6): 7.2321428571428585, (5, 7): 9.6675, (5, 8): 7.130000000000001, (5, 10): 9.049999999999999, (5, 12): 13.999999999999998, (6, 1): 8.66, (6, 2): 7.873571428571429, (6, 3): 7.374999999999999, (6, 5): 6.726818181818181, (6, 7): 6.7074509803921485, (6, 9): 9.031200000000002, (6, 10): 8.143846153846155, (6, 13): 11.754999999999999, (7, 0): 14.3, (7, 2): 7.551818181818183, (7, 5): 9.354545454545454, (7, 6): 9.2275, (7, 9): 10.295625, (7, 10): 9.450000000000001, (7, 11): 9.237000000000002, (8, 4): 7.270952380952383, (8, 5): 8.730909090909092, (8, 6): 10.31923076923077, (8, 9): 9.78, (8, 12): 8.796153846153846, (8, 13): 9.3, (8, 14): 14.87333333333333, (9, 0): 9.892857142857142, (9, 2): 12.05909090909091, (9, 4): 7.990909090909093, (9, 5): 9.074333333333334, (9, 6): 9.919354838709681, (9, 7): 11.73782608695652, (9, 11): 8.61, (9, 12): 8.033333333333335, (9, 14): 9.086666666666668, (9, 15): 8.366666666666665, (10, 3): 12.275000000000002, (10, 4): 10.214285714285715, (10, 9): 7.441428571428574, (10, 14): 6.711111111111111, (11, 3): 9.25, (11, 6): 8.576451612903226, (11, 7): 8.644333333333334, (11, 9): 9.09857142857143, (11, 13): 9.299999999999999, (11, 15): 5.973571428571427, (12, 4): 9.899999999999999, (12, 9): 10.266666666666667, (12, 10): 11.299999999999997, (13, 6): 12.000000000000002, (13, 7): 12.4, (13, 8): 11.416666666666664, (13, 9): 8.0375, (13, 15): 6.5, (14, 6): 12.323999999999998, (14, 7): 12.277777777777775, (14, 15): 5.373333333333334, (15, 2): 15.249999999999998, (15, 3): 17.999999999999996, (15, 7): 10.375, (15, 10): 9.9625, (15, 13): 9.199999999999998, (1, 1): 15.12111111111111, (2, 2): 5.354666666666667, (3, 3): 6.814285714285714, (4, 4): 5.949999999999999, (5, 5): 6.486842105263157, (6, 6): 8.116086956521741, (7, 7): 8.7425, (9, 9): 6.1}
# init_car[region] = {region: initial_no_cars}
init_car = {0: 3.0, 1: 34.0, 2: 82.0, 3: 57.0, 4: 16.0, 5: 28.0, 6: 39.0, 7: 43.0, 8: 11.0, 9: 25.0, 10: 13.0, 11: 25.0, 12: 94.0, 13: 32.0, 14: 9.0, 15: 15.0}
Transformed the above-data to be used for Reinderien code
Price= [7.1, 8.626999999999999, 11.568749999999998, 8.120000000000001, 8.425, 10.23125, 11.500000000000004, 13.633333333333331, 13.999999999999998, 7.558333333333335, 8.733333333333333, 6.610999999999999, 8.899999999999999, 9.883333333333331, 9.159523809523812, 9.978749999999998, 12.54, 11.6625, 13.3, 12.800000000000002, 15.6, 9.136000000000001, 6.163103448275866, 13.050000000000006, 10.143928571428571, 15.787499999999998, 11.0, 20.4, 9.042499999999999, 7.249999999999998, 7.123913043478264, 11.058518518518518, 6.87142857142857, 14.36714285714286, 8.120000000000001, 6.35, 8.975000000000001, 8.329999999999998, 6.177419354838713, 7.2321428571428585, 9.6675, 7.130000000000001, 9.049999999999999, 13.999999999999998, 8.66, 7.873571428571429, 7.374999999999999, 6.726818181818181, 6.7074509803921485, 9.031200000000002, 8.143846153846155, 11.754999999999999, 14.3, 7.551818181818183, 9.354545454545454, 9.2275, 10.295625, 9.450000000000001, 9.237000000000002, 7.270952380952383, 8.730909090909092, 10.31923076923077, 9.78, 8.796153846153846, 9.3, 14.87333333333333, 9.892857142857142, 12.05909090909091, 7.990909090909093, 9.074333333333334, 9.919354838709681, 11.73782608695652, 8.61, 8.033333333333335, 9.086666666666668, 8.366666666666665, 12.275000000000002, 10.214285714285715, 7.441428571428574, 6.711111111111111, 9.25, 8.576451612903226, 8.644333333333334, 9.09857142857143, 9.299999999999999, 5.973571428571427, 9.899999999999999, 10.266666666666667, 11.299999999999997, 12.000000000000002, 12.4, 11.416666666666664, 8.0375, 6.5, 12.323999999999998, 12.277777777777775, 5.373333333333334, 15.249999999999998, 17.999999999999996, 10.375, 9.9625, 9.199999999999998, 15.12111111111111, 5.354666666666667, 6.814285714285714, 5.949999999999999, 6.486842105263157, 8.116086956521741, 8.7425, 6.1]
Demand= [1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 3, 1, 1, 5, 2, 2, 3, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 5, 4, 1, 3, 1, 1, 2, 1, 3, 3, 2, 1, 1, 1, 2, 1, 1, 1, 5, 4, 1, 1, 1, 1, 1, 3, 4, 1, 1, 2, 2, 1, 2, 2, 2, 1, 2, 2, 3, 5, 1, 1, 1, 3, 1, 2, 1, 1, 1, 1, 1, 6, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 2, 3, 4, 1, 1, 2, 4, 1]
Availability= [3.0, 34.0, 82.0, 57.0, 16.0, 28.0, 39.0, 43.0, 11.0, 25.0, 13.0, 25.0, 94.0, 32.0, 9.0, 15.0]
Solution
I'm going to do some wild guessing, given how much you've left open to interpretation. Let's assume that
- this is a maximisation problem, since the minimisation problem is trivial
- Expression (1) is actually the maximisation objective function, though you failed to write it as such
panddare floating-point vectorsXis an integer vectorcis a floating-point scalar- the graph edges, since you haven't described them at all, do not matter for problem setup
The variable names are not well-chosen and hide what they actually contain. I demonstrate potential replacements.
import numpy as np
from numpy.random._generator import Generator
from scipy.optimize import milp, Bounds, LinearConstraint
import scipy.sparse
from numpy.random import default_rng
# Edges between a source and destination location
edges = (( 0, 1), ( 0, 5), ( 0, 6), ( 1, 0), ( 1, 3), ( 1, 5), ( 1, 7), ( 1, 10), ( 1, 13),
( 2, 0), ( 2, 1), ( 2, 3), ( 2, 4), ( 2, 5), ( 2, 6), ( 2, 7), ( 2, 8), ( 2, 9),
( 2, 10), ( 2, 11), ( 2, 12), ( 3, 1), ( 3, 2), ( 3, 5), ( 3, 6), ( 3, 9), ( 3, 11),
( 3, 13), ( 4, 2), ( 4, 3), ( 4, 5), ( 4, 7), ( 4, 8), ( 4, 11), ( 4, 13), ( 5, 0),
( 5, 1), ( 5, 2), ( 5, 4), ( 5, 6), ( 5, 7), ( 5, 8), ( 5, 10), ( 5, 12), ( 6, 1),
( 6, 2), ( 6, 3), ( 6, 5), ( 6, 7), ( 6, 9), ( 6, 10), ( 6, 13), ( 7, 0), ( 7, 2),
( 7, 5), ( 7, 6), ( 7, 9), ( 7, 10), ( 7, 11), ( 8, 4), ( 8, 5), ( 8, 6), ( 8, 9),
( 8, 12), ( 8, 13), ( 8, 14), ( 9, 0), ( 9, 2), ( 9, 4), ( 9, 5), ( 9, 6), ( 9, 7),
( 9, 11), ( 9, 12), ( 9, 14), ( 9, 15), (10, 3), (10, 4), (10, 9), (10, 14), (11, 3),
(11, 6), (11, 7), (11, 9), (11, 13), (11, 15), (12, 4), (12, 9), (12, 10), (13, 6),
(13, 7), (13, 8), (13, 9), (13, 15), (14, 6), (14, 7), (14, 15), (15, 2), (15, 3),
(15, 7), (15, 10), (15, 13), ( 1, 1), ( 2, 2), ( 3, 3), ( 4, 4), ( 5, 5), ( 6, 6),
( 7, 7), ( 9, 9))
# demand[Edge]={(src, dst): demand} aka. d
demand = {( 0, 1): 1, ( 0, 5): 1, ( 0, 6): 1, ( 1, 0): 1, ( 1, 3): 1, ( 1, 5): 2, ( 1, 7): 1,
( 1, 10): 1, ( 1, 13): 1, ( 2, 0): 2, ( 2, 1): 1, ( 2, 3): 3, ( 2, 4): 1, ( 2, 5): 1,
( 2, 6): 5, ( 2, 7): 2, ( 2, 8): 2, ( 2, 9): 3, ( 2, 10): 2, ( 2, 11): 2, ( 2, 12): 1,
( 3, 1): 1, ( 3, 2): 1, ( 3, 5): 1, ( 3, 6): 1, ( 3, 9): 2, ( 3, 11): 1, ( 3, 13): 1,
( 4, 2): 1, ( 4, 3): 1, ( 4, 5): 5, ( 4, 7): 4, ( 4, 8): 1, ( 4, 11): 3, ( 4, 13): 1,
( 5, 0): 1, ( 5, 1): 2, ( 5, 2): 1, ( 5, 4): 3, ( 5, 6): 3, ( 5, 7): 2, ( 5, 8): 1,
( 5, 10): 1, ( 5, 12): 1, ( 6, 1): 2, ( 6, 2): 1, ( 6, 3): 1, ( 6, 5): 1, ( 6, 7): 5,
( 6, 9): 4, ( 6, 10): 1, ( 6, 13): 1, ( 7, 0): 1, ( 7, 2): 1, ( 7, 5): 1, ( 7, 6): 3,
( 7, 9): 4, ( 7, 10): 1, ( 7, 11): 1, ( 8, 4): 2, ( 8, 5): 2, ( 8, 6): 1, ( 8, 9): 2,
( 8, 12): 2, ( 8, 13): 2, ( 8, 14): 1, ( 9, 0): 2, ( 9, 2): 2, ( 9, 4): 3, ( 9, 5): 5,
( 9, 6): 1, ( 9, 7): 1, ( 9, 11): 1, ( 9, 12): 3, ( 9, 14): 1, ( 9, 15): 2, (10, 3): 1,
(10, 4): 1, (10, 9): 1, (10, 14): 1, (11, 3): 1, (11, 6): 6, (11, 7): 4, (11, 9): 1,
(11, 13): 1, (11, 15): 1, (12, 4): 1, (12, 9): 1, (12, 10): 1, (13, 6): 1, (13, 7): 1,
(13, 8): 1, (13, 9): 1, (13, 15): 1, (14, 6): 1, (14, 7): 2, (14, 15): 1, (15, 2): 1,
(15, 3): 1, (15, 7): 1, (15, 10): 1, (15, 13): 2, ( 1, 1): 2, ( 2, 2): 3, ( 3, 3): 4,
( 4, 4): 1, ( 5, 5): 1, ( 6, 6): 2, ( 7, 7): 4, ( 9, 9): 1}
# prices[Edge] = {(src, dst): price}, aka. p aka. profit
prices = {(0, 1): 7.1, (0, 5): 8.626999999999999, (0, 6): 11.568749999999998,
(1, 0): 8.120000000000001, (1, 3): 8.425, (1, 5): 10.23125, (1, 7): 11.500000000000004,
(1, 10): 13.633333333333331, (1, 13): 13.999999999999998, (2, 0): 7.558333333333335,
(2, 1): 8.733333333333333, (2, 3): 6.610999999999999, (2, 4): 8.899999999999999,
(2, 5): 9.883333333333331, (2, 6): 9.159523809523812, (2, 7): 9.978749999999998,
(2, 8): 12.54, (2, 9): 11.6625, (2, 10): 13.3, (2, 11): 12.800000000000002,
(2, 12): 15.6, (3, 1): 9.136000000000001, (3, 2): 6.163103448275866,
(3, 5): 13.050000000000006, (3, 6): 10.143928571428571, (3, 9): 15.787499999999998,
(3, 11): 11.0, (3, 13): 20.4, (4, 2): 9.042499999999999, (4, 3): 7.249999999999998,
(4, 5): 7.123913043478264, (4, 7): 11.058518518518518, (4, 8): 6.87142857142857,
(4, 11): 14.36714285714286, (4, 13): 8.120000000000001, (5, 0): 6.35,
(5, 1): 8.975000000000001, (5, 2): 8.329999999999998, (5, 4): 6.177419354838713,
(5, 6): 7.2321428571428585, (5, 7): 9.6675, (5, 8): 7.130000000000001,
(5, 10): 9.049999999999999, (5, 12): 13.999999999999998, (6, 1): 8.66,
(6, 2): 7.873571428571429, (6, 3): 7.374999999999999, (6, 5): 6.726818181818181,
(6, 7): 6.7074509803921485, (6, 9): 9.031200000000002, (6, 10): 8.143846153846155,
(6, 13): 11.754999999999999, (7, 0): 14.3, (7, 2): 7.551818181818183,
(7, 5): 9.354545454545454, (7, 6): 9.2275, (7, 9): 10.295625, (7, 10): 9.450000000000001,
(7, 11): 9.237000000000002, (8, 4): 7.270952380952383, (8, 5): 8.730909090909092,
(8, 6): 10.31923076923077, (8, 9): 9.78, (8, 12): 8.796153846153846, (8, 13): 9.3,
(8, 14): 14.87333333333333, (9, 0): 9.892857142857142, (9, 2): 12.05909090909091,
(9, 4): 7.990909090909093, (9, 5): 9.074333333333334, (9, 6): 9.919354838709681,
(9, 7): 11.73782608695652, (9, 11): 8.61, (9, 12): 8.033333333333335,
(9, 14): 9.086666666666668, (9, 15): 8.366666666666665, (10, 3): 12.275000000000002,
(10, 4): 10.214285714285715, (10, 9): 7.441428571428574, (10, 14): 6.711111111111111,
(11, 3): 9.25, (11, 6): 8.576451612903226, (11, 7): 8.644333333333334,
(11, 9): 9.09857142857143, (11, 13): 9.299999999999999, (11, 15): 5.973571428571427,
(12, 4): 9.899999999999999, (12, 9): 10.266666666666667, (12, 10): 11.299999999999997,
(13, 6): 12.000000000000002, (13, 7): 12.4, (13, 8): 11.416666666666664, (13, 9): 8.0375,
(13, 15): 6.5, (14, 6): 12.323999999999998, (14, 7): 12.277777777777775,
(14, 15): 5.373333333333334, (15, 2): 15.249999999999998, (15, 3): 17.999999999999996,
(15, 7): 10.375, (15, 10): 9.9625, (15, 13): 9.199999999999998, (1, 1): 15.12111111111111,
(2, 2): 5.354666666666667, (3, 3): 6.814285714285714, (4, 4): 5.949999999999999,
(5, 5): 6.486842105263157, (6, 6): 8.116086956521741, (7, 7): 8.7425, (9, 9): 6.1}
# it's a mystery
rand: Generator = default_rng(seed=0)
capacity = rand.integers(low=1, high=100)
# init_car[region] = {region: initial_no_cars} aka. X aka. accInit
init_car = {0: 3.0, 1: 34.0, 2: 82.0, 3: 57.0, 4: 16.0, 5: 28.0, 6: 39.0, 7: 43.0,
8: 11.0, 9: 25.0, 10: 13.0, 11: 25.0, 12: 94.0, 13: 32.0, 14: 9.0, 15: 15.0}
N = len(prices) # number of non-zero edges
c = np.zeros(2*N) # f (passenger flow) and x (number of vehicles???)
# (1) f maximized with coefficients of 'p'
c[:N] = [-prices[edge] for edge in edges]
# x not optimized
CONTINUOUS = 0
INTEGER = 1
integrality = np.empty_like(c, dtype=int)
integrality[:N] = CONTINUOUS # f
integrality[N:] = INTEGER # x
upper = np.empty_like(c)
# (4) upper bound of f
upper[:N] = [demand[edge] for edge in edges]
# (2) upper bound of x
upper[N:] = [init_car[source] for source, dest in edges]
eye_N = scipy.sparse.eye(N)
# (3) 0 <= -f + cx
A = scipy.sparse.hstack((-eye_N, capacity*eye_N))
result = milp(
c=c, integrality=integrality,
bounds=Bounds(lb=np.zeros_like(c), ub=upper),
constraints=LinearConstraint(lb=np.zeros(N), A=A),
)
print(result.message)
flow = result.x[:N]
vehicles = result.x[N:].astype(int)