The Linear Multi-commodity min-cost flow (MMCF) problems can be modeled in a number of ways depending on how the term “commodity” is defined.
There are three major options: a commodity may originate at a subset of nodes in the network and be destined for another subset of nodes, or it may originate from a single node and be destined for a subset of nodes, or it may originate from a single node and be destined for a single node.
Let us consider a model of the first case.
In the MMCF problem, a directed graph $G(N,A)$, where $N$ a set of nodes is and $A$ is a set of Arcs is given. $A$ set of $K$ commodities has to be routed on $G$ at minimal total cost while satisfying the usual flow-conservation constraints at the nodes.
Parameters
$c^k_{ij}$ : cost of assigning commodity $k$ to arc $ij$
$d_{ij}$ : capacity of arc $ij$
$b^k_{i}$ : demand/supply of commodity $k$ at node $i$
Decision variable
$x^k_{ij}$ : amount of flow of commodity $k$ to arc $ij$
The arc-formulation according to Cappanera and Frangioni is given as:
$$
\begin{array}{ll}
\mbox{minimize} & \sum_{k\in \mathcal{K}}\sum_{{ij}\in \mathcal{A}} c^k_{ij} x^k_{ij} \quad\quad\quad(1)\\
\mbox{subject to} & \sum_{{ij}\in \mathcal{A}} x^k_{ij} - \sum_{{ji}\in \mathcal{A}} x^k_{ji} = b^k_{i}\quad \forall{{i}\in \mathcal{N}}, \forall{{k}\in \mathcal{k}} \quad(2)\\
& 0\leq x^k_{ij} \leq d^k_{ij} \quad \forall{{ij}\in \mathcal{A}}, \forall{{k}\in \mathcal{k}} \quad(3)\\
& \sum_{{k}\in \mathcal{K}} x^k_{ij} \leq d_{ij} \quad \forall{{ij}\in \mathcal{A}}, \quad(4)\\
\end{array}
$$
Equation (1) defines the objective function and equation (2) is the flow-conservation constraint. Equations (3) and (4) are respectively the individual and mutual capacity constraints.
The goal here is to solve instances from MNETGEN (a popular MMCF benchmark instance set) using the formulation above. Please check the link for full details on the instance generation and file formats.
In the following code, we introduce the various varaibles that hold our instance data
#import packages
import cvxpy as cp
import numpy as np
problemName = 'instance/64-4-1'
nodeFile = problemName+'.nod'
arcFile = problemName+'.arc'
mutFile = problemName+'.mut'
supFile = problemName+'.sup'
instanceInfo = [] # [commodity, nodes, arcs, cap_arcs]
numNodes = 0
numCommodity = 0
numArcs = 0
numCapacitatedArcs = 0
bigM = 99999999
The above code assumes the instance files are stored in a folder named instances. Solving a single instance requires 4 files with details of the file formats at MNETGEN.
In the next code, we read the each of the 4 files.
#Read the node file
nodeFileData = open(nodeFile, "r")
for eachLine in nodeFileData:
instanceInfo.append(int(eachLine))
numCommodity = instanceInfo[0]
numNodes = instanceInfo[1]
numArcs = instanceInfo[2]
nodeFileData.close()
#Read mutual capacity pointer file
mutualCapacityPointer = np.loadtxt(mutFile)
mcpFileData.close()
#Read arc file
uniqueArcs = {}
arcCapacity = {}
mutualCapacity = {}
cost = {}
temp_arc = [0]*numCommodity
temp_cost = [0]*numCommodity
arcFileData = open(arcFile, "r")
for eachLine in arcFileData:
data = eachLine.split(" ")
uniqueArcs[int(data[0])] = [int(data[1]),int(data[2])]
mutualCapacity[int(data[0])] = int(data[6])
inner_index = int(data[3])-1
outer_index = int(data[0])-1
if temp_arc[inner_index]==0:
temp_arc[inner_index] = [0]*numArcs
temp_arc[inner_index][outer_index]=int(data[5])
else:
temp_arc[inner_index][outer_index]=int(data[5])
if temp_cost[inner_index]==0:
temp_cost[inner_index] = [bigM]*numArcs
temp_cost[inner_index][outer_index] = float(data[4])
else:
temp_cost[inner_index][outer_index] = float(data[4])
for i in range(numCommodity):
arcCapacity[i+1] = temp_arc[i]
cost[i+1] = temp_cost[i]
arcFileData.close()
#Read sup file
supplyDemand = {}
temp_supDem = [0]*numCommodity
supFileData = open(supFile, "r")
for eachLine in supFileData:
data = list(map(int, eachLine.split(" ")))
inner_index = data[1]-1
outer_index = data[0]-1
if temp_supDem[inner_index]==0:
temp_supDem[inner_index] = [0]*numNodes
temp_supDem[inner_index][outer_index]= data[2]
else:
temp_supDem[inner_index][outer_index]= data[2]
for i in range(numCommodity):
supplyDemand[i+1] = temp_supDem[i]
supFileData.close()
We skip the details of the above code to read the problem data since there are other ways for it to be accomplished.
In the final section, we introduce the model and solve the MMCF problem.
#---------------------------Solve the math model------------------------------
#
C = [(a,k) for a in uniqueArcs.keys() for k in range(numCommodity)] #decision variable indexer
X_key = {}
my_ind = 0
for (a,k) in C:
X_key[a,k] = my_ind
my_ind+=1
X = cp.Variable(shape=len(C), nonneg=True, name="X")
allConstraints = []
#constraint 1
for i in range(numNodes):
for k in range(numCommodity):
sumIn = 0
sumOut = 0
for akey,value in uniqueArcs.items():
if value[0]==(i+1):
sumIn +=X[X_key[akey,k]]
elif value[1]==(i+1):
sumOut +=X[X_key[akey,k]]
const1=sumIn-sumOut==supplyDemand[k+1][i]
allConstraints.append(const1)
#constraint 2
for (a,k) in C:
maxVal = bigM#sys.float_info.max
if arcCapacity[k+1][a-1]==-1:
const2=X[X_key[a,k]]<=maxVal
allConstraints.append(const2)
else:
const2=X[X_key[a,k]]<=arcCapacity[k+1][a-1]
allConstraints.append(const2)
#constraint 3
for a in uniqueArcs.keys():
if mutualCapacity[a]!=0:
totalVars = [X[X_key[a,k]] for k in range(numCommodity)]
lhs = cp.sum(totalVars)
const3=lhs<=mutualCapacityPointer[mutualCapacity[a]-1][1]
allConstraints.append(const3)
actualCost = np.asarray([cost[k+1][a-1] for (a,k) in C])
objExpr = cp.sum(actualCost@X) #objective function
objFunc = cp.Minimize(objExpr)
prob = cp.Problem(objFunc, allConstraints)
prob.solve()
print("MMCF objective function value is {}".format(prob.value))
MMCF objective function value is 290806.3029314588
It is evident from the math model code that the $x^k_{ij}$ was first deflated to $x^k_{a}$ then to $x$. To emphasize, the change to 1 dimension is to demonstrate a way to get pass the 2 dimension limitation in CVXPY as noted in the manual.To achieve this, we build an indexer for $x$ and then store the mapping of the 2-d format in a dictionary. This way the code is still readable as in the math model and it shows how already coded math models in other platforms can be ported to CVXPY.