Degree Project in Technology

First cycle, 15 credits

Linear Programming Algorithms for

Multi-commodity Flow Problems
Exploring the performance of linear programming algorithms for
finding optimal solutions to multi-commodity flow problems

ISAAC ENQUIST & PHILIP SJÖGREN

Stockholm, Sweden 2022

Linear Programming Algorithms for Multi-commodity Flow
Problems
Isaac Rosenberg Enquist Philip Sjögren
isaacre@kth.se phisjo@kth.se
June 6, 2022

i
ii
Sammanfattning
Ett flödesproblem består av att flytta ett flertal varor från deras respektive källor till
deras sänkor genom ett nätverk där varje kant har olika kostnader och begränsningar.
Denna avhandling utforskar olika linjärprogrammeringsalgoritmer och deras prestanda
när det gäller att hitta en optimal lösning för ett givet flödesproblem. Genom att testa
en mängd olika begräsningar och strukturer för nätverket utreders vilken algoritm som
är mest gångbar för givna nätverk och problemformuleringar. Vi implementerar även
vår egna lösning och jämför dennas prestanda gentemot toppmoderna implementa-
tioner algoritmer såsom Simplex och Interior-Point. Resultaten visar att av de algorit-
mer som testades fanns det inte någon klass som är optimal för att lösa flödesproblem,
utan deras prestanda beror på hur problemet och nätverket är strukturerat.

Nyckelord:
Flödesproblem, Linjärprogrammering, Linjäroptimering, Simplexmetoden, Grafteori,
Optimeringslära, Transportproblem, Minkostnadsflödesproblem

iii
iv
Abstract
A multi-commodity flow problem consists of moving several commodities from their
respective sources to their sinks through a network where each edge has different costs
and capacity constraints. This paper explores different linear programming algorithms
and their performance regarding finding optimal solutions for multi-commodity flow
problems. By testing several different network constraints, it is examined which algo-
rithms are most suitable for specific network and problem structures. Furthermore, an
implementation of our own multi-commodity solver is created and its performance is
compared against state-of-the-art linear programming solvers. The results show that
for the methods we tested it is difficult to discern which class of linear programming
methods are optimal solvers for multi-commodity flow problems and that their perfor-
mance depends on how the network and commodities are structured.

Keywords:
Multi-commodity Flow Problem, Linear Programming, Linear Optimization, Simplex,
Interior-point Method, Graph Theory, Optimization Theory, Transportation Problems

v
vi
Preface
This thesis explores common algorithms to solve multi-commodity flow problems and
examines their efficiency. The authors Isaac Enquist and Philip Sjögren are students
at the Degree Programme in Engineering Physics at the Royal Institute of Technology.
The study was performed in spring 2022 at the institution of Optimization and Sys-
tems Theory within the course Degree Project in Engineering Physics (SA114X). The
authors wish to convey their sincere appreciation to the supervisor Isabel Haasler for
continuous guidance and assistance during process of writing this paper.

vii
viii
Contents
1 Introduction 1

2 Background 2
2.1 Basic graph theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 Commodity flow problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.3 Multi-commodity flow problems . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Linear programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Basic and non-basic variables . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.6 Simplex method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.7 Finding an initial basic solution . . . . . . . . . . . . . . . . . . . . . . . . 9
2.8 Degeneracy and implemented algorithm . . . . . . . . . . . . . . . . . . . 9
2.9 Interior point methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Method 13
3.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Parameter choices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 The algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

4 Results and Analysis 17

4.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.1 Network size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.2 Commodities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2.1 Network size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2.2 Commodities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5 Discussion and Conclusion 24

5.1 implemented Simplex method . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.2 Further Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

ix
x
1 Introduction
A multitude of modern world mechanisms may be interpreted as flows through net-
works, and are reliant on efficient methods for optimizing that flow in relation to some
arbitrary resource. There are numerous example where multi-commodity network flows
are used used as models. Global freight shipping, routing vehicles to minimize traffic,
street planing, messages in telecommunications and forest management are all exam-
ples of areas where implementations of commodity-flow problems are used.1
The purpose of this thesis is to develop and compare efficient algorithms used for solv-
ing multi-commodity network flow problems. A multi-commodity flow problem con-
sists of moving two or more commodities between their respective sources and sinks
through a network where each edge has different costs and capacity constraints. Linear
programming or optimization is a field of methods which attempt to find the optimal
solution to problems whose requirements consist of linear constraints and costs. This
paper aims to examine different types of linear programming methods and explore
their performance when it comes to finding optimal solutions to multi-commodity flow
problems.

1
Cynthia Barnhart, Niranjan Krishnan, and Pamela H. Vance. “Multicommodity flow problems”.
In: Encyclopedia of Optimization. Ed. by Christodoulos A. Floudas and Panos M. Pardalos. Boston,
MA: Springer US, 2009, pp. 2354–2362.

1
2 Background
This section provides a background to commodity flow problems and an introduction
to solutions to them with different linear programming techniques.

2.1 Basic graph theory

In graph theory a graph G(V, E) is defined as a set of vertices V linked together by a
set of edges E where each edge is associated with two vertices. A common synonym for
vertex is node, which will be the preferred term for the remainder of this paper. If the
number of nodes is denoted I then each node in V may be associated with an index
i = 1, . . . , I and a specific edge in E that connects nodes i and j may be denoted (i, j).
The paper will exclusively concern directed graphs which is a specific type of graph
where each edge (i, j) ∈ E only allows movement or flow from node i to node j. This
means that (i, j) ̸= (j, i). The graph will furthermore be assumed to have a dual-arc
formulation, meaning that each existing edge between two nodes has a corresponding
edge in the opposite direction. In notation this may be written as:

∃ (i, j) ∈ E ⇐⇒ ∃ (j, i) ∈ E, i, j ∈ V (1)

From this basic formulation of graph theory a simple directed graph may be illustrated
in the figure below:
A simple directed graph
(2, 3)/(3, 2)
3 2

(1, 3)/(3, 1) (1, 2)/(2, 1)

1

Figure 1: A simple dual-arc directed graph G(V, E) where N = 3. An edge between node i
and j is denoted (i, j).

If for some graph G(V, E) each edge (i, j) ∈ E is associated with a corresponding cost
cij then the graph may be referred to as a network. If for this network V one or more
nodes are assigned as the sources (starting node) and the sinks (end node) for some
flow of size V then the graph may be referred to as a flow network. Flow networks will
be expanded upon in the following subsection. The aim for the remaining sections of
the background is to provide methods with which the cost of a flow through such a
flow network may be minimized.

2.2 Commodity flow problems

A commodity flow problem consists of a flow network in which one or more commodi-
ties should move from one or more starting nodes (sources) to one or more end nodes
(sinks) without breaching any constraints. A constraint could for example be a capac-
ity constraint on an edge. Furthermore each path is associated with a cost which needs

2
to be minimized by taking the cheapest path through the graph. A simple example of
a graph consisting of N = 3 nodes with costs and constraints is illustrated below:
Constraints and costs for the network
-10
(10/0.4)
3 2

(10/0.1) (5/0.2)
1
+10

Figure 2: Constraints and costs for a simple network (constraint/cost).

In Figure 2 a positive and negative total flow, 10 in the figure, is used to indicate the
nodes that acts as the source and the sink of the flow respectively. In addition to this
the associated cost and constraint to the edges, here given by a upper bound, is pre-
sented next to the edge in the figure as (constraint/cost). For this example this means
that both edges in a given dual-arc share a common constraint and cost. Now consider
that the aim is to move 10 units of a commodity from node one to node two. This is a
simple example of a single-commodity flow problem. Two solutions are shown below:
Solutions and costs for different flows
Cost of 5 Cost of 3.5
-10 -10
10 5
3 2 3 2

10 5 5
1 1
+10 +10

Figure 3: Two different solutions to the single-commodity flow problem.

It is clear that both solutions satisfies the conditions of moving 10 units of commod-
ity from node one to node two. If an additional goal is to minimize the cost of moving
commodity, feasible and optimal solutions are important concepts. Since both flows
fulfill the constraints they are both feasible, but the flow with a cost of 3.5 is the opti-
mal solution of the problem. The previous single-commodity flow problem can also be
defined mathematically as a minimization problem:

3
 minimize z = 0.2x12 + 0.2x21 + 0.1x13 + 0.1x31 + 0.4x23 + 0.4x32
x
x12 ≤ 5
x21 ≤ 5
x13 ≤ 10
x31 ≤ 10
x23 ≤ 10 (2)
Subject to x32 ≤ 10
x12 −x12 +x13 −x31 = 10
−x12 +x12 +x23 −x32 = −10
−x13 +x31 −x23 +x32 = 0
x12 , x21 , x13 , x31 , x23 , x32 ≥ 0

where z is the cost and xij is the amount of commodity flowing between edge (i, j). It
is possible to use vector notation to describe a general single-commodity flow problem.
Firstly a cost vector, denoted c is defined. The cost vector for the previous problem
seen in Figure 2 becomes:

cT = 0.2 0.2 0.1 0.1 0.4 0.4

 
(3)

Secondly denote a upper bound constraint vector U , which holds all constraints associ-
ated with inequalities in equation 2:

U T = 5 5 10 10 10 10
 
(4)

Also denote a constraint vector D for the required flow for all nodes, which holds all
constraints associated with equities in equation 2:

DT = 10 −10 0
 
(5)

Thirdly a constraint vector b , which holds all constraints, may be formed by the afore-
mentioned U and D. For the previous problem illustrated in Figure 2 the vector be-
comes:

 T
T U  
b = = 5 5 10 10 10 10 10 −10 0 (6)
D

Finally one has to consider the source and sink of the system. One way to handle this
is to ensure that the sum of all flows from the sources and sinks is equal to the amount
of commodity being created or removed from the system. Furthermore the sum of all
flows for all other nodes should be zero since no commodity should be created or re-
moved from the system at those points. For the problem in Figure 2:

4
 
I
X 10 if i = 1

(xij − xji |(i, j) ∈ E) = −10 if i = 2 i = 1, . . . , I (7)

j=1 
0 otherwise

where I is the number of nodes in the graph, thereby I = 3 for this problem. The
general formulation of single-commodity flow problem therefore becomes:

minimize z = cT x
x

XI V if node i is a source

Subject to (xij − xji |(i, j) ∈ E) = −V if node i is a sink ∀i (8)

j=1 
0 otherwise
U ≥x≥0

where V is the total amount of commodity, z is the cost, c is the cost function, b is the
constraint vector and as previous I is the number of nodes in the network.

2.3 Multi-commodity flow problems

The previous result can easily be generalized to a multi-commodity flow problem by
adding an index k to describe each commodity. This gives:

k
XI V if node i is a source in k

k k
(xij − xji |(i, j) ∈ E) = −V k if node is a sink in k ∀ i, k (9)

j=1 
0 otherwise

where V k is the amount of commodity k being created or removed from the system,
xkij is the amount of commodity k travelling on the edge (i, j) and I is still the number
of nodes in the network. The complete problem formulation for the multi-commodity
flow problem, which this report intends to find and investigate effective algorithms for,
then becomes:

minimize z = cT x
x

k
XI V if node i is a source in k

Subject to (xkij − xkji |(i, j) ∈ E) = −V k if node i is a sink in k ∀ i, k (10)

j=1 
0 otherwise
U ≥x≥0

An example of a multi-commodity flow problem could be presented in graph similar to

that of a simpler single-commodity flow problem, seen in Figure 2. An illustration for
the same network as in the previous section now with two commodities is shown below
in Figure 3. The figure includes both the formulation of the problem and the optimal
solution shown side-by-side:

5
 Problem and solution for a multi-commodity flow problem
Problem Solution
-8 5 -8
(15/15)
3 2 3 2
2
+5 +5
2 3

(10/10) (3/3) 5 3
1 1
+8 -5 +8 -5

Figure 4: Problem formulation with constraints and cost the the left and resulting solution
to the right to the multi-commodity flow problem. The different commodities are assigned the
colors red and blue.

In Figure 4 both the formulation of a multi-commodity flow problem and its optimal
solution is presented. The problem is formulated on the left in the figure where the
constraint and cost of each edge is presented as before such that (constraint/cost) is
used as a label for each edge. For ease of calculation the cost of each edge is set equal
to its capacity. The problem presented consists of two commodities with each com-
modity assigned a different color, red or blue, with labels placed next to their related
source and sink node. A source and sink is designated with a positive or negative sign
respectively. In the solution presented to the right the different paths that each com-
modity is forced to take is presented with directed edges given the color associated
with their respective commodity.

2.4 Linear programming

Linear programming is a subset of optimization theory surrounding optimization prob-
lems which have strictly linear constraints. A general linear program can be written in
standard form as:

minimize z = cT x
x
Ax = b (11)
Subject to
x≥0

The reason that this describes a general linear programming problem despite not in-
cluding any inequalities for the constraint vector b is because one can use slack vari-
ables to transform the problem into the desired form, like bellow:
( (
x1 ≤ 3 x1 + s 1 = 3
x1 , x 2 ≥ 0 becomes x1 , x 2 , s 1 , s 2 ≥ 0 (12)
x2 ≥ 4 x2 − s2 = 4

By adding slack variables to all inequality constraints in the multi-commodity flow

problem it is clear that the problem becomes standard form, as shown in equation 11.
An important consequence of transforming the problem into standard form is that this

6
form is well suited for optimization algorithms like the Simplex method which is de-
scribed further in Section 2.6. Another important detail is the fact that one can apply
the fundamental theorem of linear programming.
Theorem 1 (Fundamental theorem of linear programming). If a linear program has a
bounded optimal solution then there exists an extreme point of the feasible region which
is optimal.2
Since the multi-commodity flow problem is bounded by the the maximum capacity
of the paths and the fact that all commodities are greater or equal to zero, the the-
orem applies if a solution exists. It follows that to find the optimal solution only the
extreme points (corners) of the n-dimensional feasible region need be investigated.

2.5 Basic and non-basic variables

When the matrix A is singular simple linear algebra entails that there exists only one
solution where x = A−1 b. Most linear programming problems, including multi-commodity
optimal flow problems, are not this simple. In most cases the system will be under-
determined and the matrix A will be of size m × n where n > m, since there are less
constraints than variables. In commodity problems this means that several different
routes can accomplish the same flow, for example several routes can be driven to reach
the same goal. In other words, the task is to find the optimal routes for all flows from
a large number of solutions.
One approach is to split the matrix A into two different matrices by choosing m lin-
early independent columns so that:

Ax = b → Aβ xβ + Aν xν = b (13)

where Aβ is square and invertible with dimensions n × n. The indices of the columns of
A used in Aβ are called the basic variables, these are denoted β = (β1 , ..., βm ). Like-
wise, indices of the columns of A used in Aν are called non-basic variables and are
denoted ν = (ν1 , ..., νl ) where l = n − m. Since Aβ is invertible and has the same
dimension as the constraint vector b a solution to the equation is given by xβ = A−1 β b
and xν = 0. Such a solution is called a basic feasible solution.
By finding the correct combination of basic variables which correspond to the optimal
solution, it is possible to solve the problem. This is the basic idea of the Simplex algo-
rithm, which is further explained in the next section.

2.6 Simplex method

The Simplex method utilizes the fact that an optimal solution to a bounded linear
problem lays on the corner points by wandering between different corners/extreme
points in the n-dimensional feasible region. Wandering between these extreme points
is equivalent changing basic variables β. Although there are many steps most con-
sist of elementary matrix operations. Moreover these operations are effective from a
2
Igor Griva. Linear and nonlinear optimization. 2nd ed.. Philadelphia: Society for Industrial and
Applied Mathematics, 2009.

7
computational standpoint. To make the algorithm easier to follow some notation is in-
troduced, the notation and flowchart are taken from a course book Optimizations by
Amol Sasane and Krister Svanberg:3

Notation Element of Definition

b̄ Rm Aβ b̄ = b
a¯j = (j = 1, ...., n) Rm Aβ a¯j = aj
y Rm ATβ y = cβ
rν Rn rν = cν − ATν y
z̄ R z̄ = cT b̄ = y T Aβ b̄ = y T b

Table 1: Definitions of the designations used in the Simplex algorithm.

If there are basic variables β which make up basic feasible solutions, it is possible to
solve the optimization problem by using the Simplex algorithm illustrated in the flowchart
below.
The Simplex algorithm
Initial basic
feasible solution
β = (β1 , . . . , βm )
ν = (ν1 , . . . , νl )

Calculate b, y, rν :
Aβ b = b
ATβ y = cβ
rν = cν − ATν y

yes Optimal solution:

rν ≥ 0 xβ = b Stop
xν = 0
no
Choose q s.t.
rνq = min {rνi : 1 = 1, . . . , l} ;
Calculate aνq : Aβ aνq = aνq

n o
bi
tmax = min ai,νq : ai,νq > 0 ;
no yes Solution does
bp aνq ≤ 0 Stop
Choose p s.t. tmax = ap,ν ; not exist
q
Interchange νq and βp

Figure 5: A description of the Simplex algorithm.

3
Amol Sasane and Krister Svanberg. Optimization. Stockholm: Department of Mathematics, Royal
Institute of Technology.

8
2.7 Finding an initial basic solution
To perform the Simplex algorithm the first step is to find an initial basic feasible solu-
tion. One way is to consider the following problem:

X
minimize yi
i
 
  x (14)
A Im =b
Subject to y
y≥0 x≥0

This problem has a basic feasible solution where x = 0 and y = b. Since y ≥ 0 the
optimal value for the sum is 0, which is fulfilled when y = 0. The key is that when
y = 0, the x values correspond to a basic feasible solution for the original problem. By
executing the Simplex algorithm on the problem stated in equation 15, for which a ba-
sic feasible solution is assured, it is possible to find a feasible solution for the original
problem. This is called solving the phase I problem. It is important to note that if the
objective function does not equate to 0 then the original problem lacks a solution.

2.8 Degeneracy and implemented algorithm

b
In the case of several indices p satisfying the condition tmax = ap,νp , see flowchart
q
in Figure 5, an arbitrary decision has to be made as to which of the multiple vari-
ables corresponding to p should leave the basis β. This does not pose a problem if
tmax > 0 since the cost function improves. But in the case where the transformed
constraint vector b includes zeroes, there is a possibility to get one or more p so that
b
tmax = ap,νp = 0. Such a solution is called a degenerate feasible solution. In the case of
q
several such zeroes, this can lead to cycling, where the algorithm wanders between the
same corners in n-dimensional solution space in an infinite cycle without improving.
This will in turn lead to the algorithm not converging.
Cycling is a big problem when working with multi-commodity flow problems since a
large proportion of the constraint vector consists of zeroes, namely all nodes which are
not sinks of sources in some commodity. There are a plethora of algorithms to solve
these issues, they are collectively named pivot rules because they decide which column
should leave and enter the basis β. In the Simplex algorithm presented Figure 5 the
pivot is performed were ever a choice of q and p has to be made. Many of these pivot
rules are complex and an implemented algorithm may use multiple pivot rules for dif-
ferent cases in order to find a solution in as few iterations as possible.
A simple and reasonably effective pivot rule is Bland’s rule:
1. Choose the highest-numbered (i.e., rightmost) nonbasic column with a negative
(reduced) cost.
2. Now among the rows, choose the one with the lowest ratio between the (trans-
formed) right hand side and the coefficient in the pivot tableau where the coeffi-
cient is greater than zero. If the minimum ratio is shared by several rows, choose
the row with the highest-numbered column (basic variable) in it.

9
It can be shown that Simplex with Bland’s rule converges to the optimal solution from
all feasible solutions. As a result of its simplicity Bland’s Rule has an inherit probabil-
ity for needing a considerable amount of iterations to converge.4
Expressing the Bland’s pivot rule in mathematical notation is also possible. This anti-
cycling Simplex algorithm necessitates consistent indexing over the columns of the ma-
trix A, while retaining the previous independent indexing over β and ν. Below the in-
dex k is introduced to assist in a given arbitrary order of columns.

Given a matrix A ∈ Rm×n

Introduce k = 1, . . . , n (15)
 T
Where A1k . . . Amk = k : th column of A

The pivot steps may be described using the new index k with kβ , kν denoting the cor-
responding matrix index for the basic and non-basic variables. Bland’s rule with the
addition of index k can be written as follows:
1. Choose q s.t.
kνq = max {kνi : rνi < 0} (16)

2. Choose p s.t.  
bi
kβp = max kβi : tmax = (17)
ai,νq

3. Interchange νq and βp
With this new pivot rule the previous Simplex algorithm found in Figure 5 may be
rewritten and the implemented algorithm follows in the figure below:

4
Robert G. Bland. “New Finite Pivoting Rules for the Simplex Method”. In: Mathematics of Op-
erations Research 2.2 (1977), pp. 103–107.

10
 The implemented Simplex algorithm
Initial basic
feasible solution
β = (β1 , . . . , βm )
ν = (ν1 , . . . , νl )

Calculate b, y, rν :
Aβ b = b
ATβ y = cβ
rν = cν − ATν y

yes Optimal solution:

rν ≥ 0 xβ = b Stop
xν = 0
no
Choose q s.t.
kνq = max {kνi : rνi < 0} ;
Calculate aνq : Aβ aνq = aνq
n o
bi
tmax = min ai,ν : ai,νq > 0
q
yes
Choose
n p s.t. o no aνq ≤ 0 Solution does
Stop
bi not exist
kβp = max kβi : tmax = ai,ν ;
q
Interchange νq and βp

Figure 6: A description of the Simplex algorithm with Bland’s rule implemented

The changes made to the algorithm are retained to steps in the algorithm where piv-
oting is performed, which can be seen as the two steps where a choice for the pivots q
and p have to be made.

2.9 Interior point methods

Interior point methods are another widely used class of solvers for convex optimization
problems, which includes all linear programming problems. Contrary to the Simplex
method, which moves between edges the Interior-point method traverses inside the fea-
sible region.5 The main idea is to transform the problem, including the constraints,
into a continuous function which should be minimized. A common method to achieve
this is to use barrier functions and then using standard search methods to locate the
global minimum.

5
Stephen G. Nash and Ariela Sofer. Linear and nonlinear programming. New York : McGraw-Hill,
1996.

11
As an example, see the following problem:

minimize f (⃗x)
x
(18)
Subject to xi ≤ b

Next, a barrier function g(x) is defined that approaches infinity at the boundary, as an
example a logarithmic barrier function is used:

n
X
g(⃗x) = log(b − xi ) (19)
i=1

When µ is a small positive parameter it is possible to rephrase the original problem as:

minimize f (⃗x) − µg(⃗x) (20)

x

Although the problem is not precisely equivalent, as µ approaches zero the problems
become more and more equivalent. Since µ is arbitrarily small the solution will lie
incredibly close to the boundary while not surpassing it because of the divergent na-
ture of the logarithm. Hereon it is possible to find the global minima with more stan-
dard methods such as gradient search and newtons methods. One caveat with using
interior-point methods is the fact that the problems need to be convex, since all lin-
ear problems (and thereby our multi-commodity formulation) are convex this is not an
issue.

12
3 Method
3.1 Problem definition
The comparisons between the different linear programming algorithms were made on
a Manhattan style grid. In such a network all nodes that do not lie on the network’s
sides are connected to four other nodes, see Figure 7. This choice was made to mirror
a standard city layout since a common application of multi-commodity flow problems
is the delivery of different goods. Furthermore the network is simple to create and gen-
eral in the sense that one can create a different network with less connections by set-
ting an upper bound of a path to zero, so that no commodities can utilize that path.
Manhattan grid

1 2 3

4 5 6

7 8 9

Figure 7: A Manhattan-grid of size N = 3.

A grid size of N , which is defined as the amount of nodes on each side of the network,
implies that there are N 2 nodes. The amount of edges in the network equals 4N (N −
1). This can be understood by realizing that each node aside from the nodes on the
sides are connected to four other nodes. Therefore one can take 4N 2 and subtract 4N
because of the four sides on the grid. All the edges are counted twice this way, which
is reasonable since they only connect to two nodes. This gives a total of 4N 2 − 4N =
4N (N − 1) edges.
The multi-commodity flow problem can be divided into two parts. Firstly making sure
each path between two nodes does not exceed its upper bound and secondly making
sure that each commodity is moved from the correct sources to the correct sinks and
that the remaining nodes have a total inflow of 0.
The first problem is easily modelled as the sum of all commodity flows for a specific
edge with an added slack variable which should be equal to the upper bound. This is
done for all edges and results in an identity matrix for each commodity plus an iden-
tity matrix for the slack variable being equal to the upper bound. The size of the iden-
tity matrices is 4N (N − 1) × 4N (N − 1) which equals the amount of paths.
The second problem ensures that all nodes have the correct amount of inflow. The
sum of all flows originating from each node should be zero unless they are a source or

13
sink. In that case the sum of flows should equal the predefined amount being sourced
or sunk. This results in a block matrix in which each block are node constraints for a
specific commodity. A node block is a directed incidence matrix and needs to abide by
the following constraints:

i
XI V if node i is a source

(xij − xji |(i, j) ∈ E) = −V i if node i is a sink ∀ i (21)

j=1 
0 otherwise

This can be represented as a matrix equation where the sum of all flows for a com-
modity and a given node is equal to the sum of all in and outgoing flows for the node
in question. Denoting each respective node-constraint matrix as noden , the complete
matrix denoted A for k commodities on which linear optimization is performed be-
comes:  
I1 I2 I3 ... Ik Islack
node1 0 0 ... 0 0 
 
 0 node2 0 ... 0 0 
A= 0 (22)
 
 0 node3 . . . 0 0  
 . .. .. .. .. .. 
 .. . . . . . 
0 0 0 . . . nodek 0
where each identity block I is of size 4N (N − 1) × 4N (N − 1) and each block of node
constraints noden is a directed incidence matrix of size N 2 × 4N (N − 1). A single com-
modity case of such a matrix A without slack, with a single directed incidence matrix
node1 , can be found in equation 2. The Dimension of the resulting matrix A becomes

N 2 K + 4N (N − 1) × 4N (N − 1)(K + 1)


m×n = (23)

The constraints vector denoted b is made up of upper bounds which correspond with
the identity blocks. This ensures that no upper bounds are violated. The upper bounds
are followed by the sum of flows for all nodes and their respective commodities. This
ensures that all sources, sinks and nodes have the correct in and outgoing flow. In
other words, the constraint vector b is defined by the aforementioned upper bound
constraint vector U and k different node flow constraints vectors Dk , each associated
with a different commodity k. Both U and D are described further in Section 2.2.
The resulting constraint vector b for a multi-commodity flow problem is written in the
equation below:
bT = U T D1T D2T . . . DkT
 
(24)

where the length of the vector b will be m = N 2 K + 4N (N − 1). The final vector used
is the cost vector c. In this formulation each arc has the same cost regardless of which
commodity is travelling on it. This may be formulated as:

clij = cm
ij ∀ l, m (25)

where l and m are commodities and xij signifies an arc between node i and node j.
Each element of the cost vector signifies a cost for a specific arc and commodity. The

14
final 4N (N − 1) variables corresponds with the slack variables of the problems, there-
fore their associated cost will be 0. If the corresponding cost for the edges in a single
commodity is designated ck then the cost vector c for the multi-commodity case be-
comes:

cT = ck ck . . .
 
ck 0 (26)

where each block of ck is identical but corresponds to a different commodity. The vec-
tor c will have the length n = 4N (N − 1)(K + 1). The objective becomes finding a
vector x which fullfills Ax = b and minimizes z = cT x. Equivalent to the earlier seen:

minimize z = cT x
x
Ax = b (27)
Subject to
x≥0

Solving this problem gives an optimal solution to a multi-commodity flow-problem in a

Manhattan-grid.

3.2 Parameter choices

Between the upper bounds, cost vector, amount of commodities and size of sources,
sinks and the network itself, it is evident that the amount of configurations are end-
less, even for small networks. Therefore some quantities will remain constant to reduce
the dimensions of the problem and scope of this research. Other quantities may be
randomly generated or iterated over for testing. This subsection examines how test-
ing is performed and which quantities are placed in which category and the reasoning
behind these choices. Note that the choices are implemented in all results unless other-
wise stated.
Since some quantities are randomized over, and the computer’s internal clock is too in-
accurate to give good results for small problems a choice was made to run several tri-
als of problems. The resulting time to complete the problems is divided by the amount
of trials to get an average time per trial for the current constraints. A number of tri-
als amounting to one hundred was chosen since it is big enough to reduce the variance
created by the randomness while being small enough to avoid enormous calculation
times for larger problems.
The cost vector c defines the costs of each path through the network. The size of the
entries in c, relative to each other, affects which path will be the optimal solution. To
encourage the commodity to not simply take the shortest path, each path received a
random cost cij ∈ [1, 10], so that the most cost effective route was not the shortest one.
The reason that the cost vector is randomized is to force the algorithms to solve differ-
ent problems of the same size for each trial. For example, when iterating over another
quantity the specific cost vector should not be what affects the result but rather the
quantity in question.
The size of the network N increases the amount of elements in A with a factor N 4 .
The size will be held constant since different size networks can affect the results, and

15
therefore testing quantities is done over a constant size. The size used depends on the
quantity tested and is displayed in the respective results and denoted N . Furthermore,
an analysis is performed to examine how the size of the grid affects the computation
time.
The upper bounds of the paths are set to a constant value of 15. While it would be
interesting to explore how different upper bounds affect the performance of the algo-
rithms, there are a couple of reasons why this is difficult to implement. For example,
if the upper bounds are randomized, there would be a chance of creating systems that
do not have a solution (if the bounds are too low or a system where bounds do not
matter if they are set too high). Instead, a upper bound of 15 was set and the com-
modities distribute correspondingly. For example, it is easy to create a source with 16
units of commodity to force the flow to take at least two paths.
The amount of commodities has an effect on the size of the matrix. In the first case,
the amount of commodities will be held constant since, different amount of commodi-
ties can affect the results. In the second case, the amount of commodities used differs
and is displayed clearly in the results denoted K. An in depth analysis is performed to
examine how the amount of commodities affect the computation time.
The multi-commodity flow problem is modeled with one source for each commodity a
general, and arbitrary rule, the number of sinks is chosen to be D
5 which is three for
D = 15. All these sources and sinks are randomly distributed to nodes in the Manhat-
tan grid. A summary of all quantities and their uses is presented below:

Quantity Denoted Generation Value Examined

Cost of edges - Randomized Between [1, 10] -
Network size N Constant Defined in plot Yes
Upper bound of edges - Constant Set to 15 -
Amount of commodities K Constant Defined in plot Yes
Total flow per commodity D Constant Defined in plot -
Amount of sources - Constant Set to 1 -
Amount of sinks - Constant Set to 3 -
Position of sources - Randomized - -
Position of sinks - Randomized - -

Table 2: Definition of quantities and how they are used in results.

3.3 The algorithms

The algorithms used for solving the problem are the standard dual Simplex and Interior-
point solvers included in MATLAB’s linear optimization package.6 Furthermore, our
own implementation of revised primal Simplex is tested. All solvers represent a line in
the result plots. Where ”Simplex” denotes MATLAB’s dual Simplex solver, ”Interior-
point” denotes MATLAB’s Interior-point solver based on the Predictor-corrector method
and ”implemented” denotes our primal Simplex implementation.

6
The MathWorks Inc. linprog Solve linear programming problems. url: https://se.mathworks.
com/help/optim/ug/linprog.html (visited on 05/19/2022).

16
4 Results and Analysis
This section will present the performance of the different linear programming methods,
whilst also providing some analysis on the obtained results.

4.1 Results
In order to measure the performance of the different algorithms relatively unbiased,
their average computation time is calculated over a number of trials. This will focus
the results on the general trend whilst enabling estimation of the error. In addition to
this, the number of internal step-iterations needed for convergence is also calculated.

4.1.1 Network size

The computation times for solving multi-commodity flow problems are averaged over a
number of trials I over an increase in the size of the networks Manhattan grid N . This
increase in network size is performed for a constant number of commodities K = 2.

Figure 8: Average computation time plotted against size of the Manhattan grid. K denotes
the number of commodities, D the total flow for each commodity and I the number of trials.
The error bars signify standard deviation.

As can be seen in Figure 8, the computation time of the implemented algorithm is or-
ders of magnitude larger than that of the built-in versions for larger network sizes. A
smaller network size is therefore needed in order to differentiate between the built-in
methods, which is shown in Figure 9 below.

17
Figure 9: Average computation time plotted against size of the Manhattan grid. K denotes
the number of commodities, D the total flow for each commodity and I the number of trials.
The error bars signify standard deviation.

In Figure 9 it may be observed that the computation time for the implemented algo-
rithms is faster, or on par, with the built-in algorithms for smaller network sizes. The
poor performance of the implemented algorithm prohibits investigation of larger net-
works, which motivates an isolated comparison between the two built-in algorithms,
shown in Figure 10 below.

Figure 10: Average computation time plotted against size of the Manhattan grid. K denotes
the number of commodities, D the total flow for each commodity and I the number of trials.
The error bars signify standard deviation.

In Figure 10 it is shown that the Interior-point algorithm outperforms the Simplex

algorithm for large networks. An isolated observation for medium sized networks is
shown in Figure 11.

18
Figure 11: Average computation time plotted against size of the Manhattan grid. K denotes
the number of commodities, D the total flow for each commodity and I the number of trials.
The error bars signify standard deviation.

Figure 11 indicates that the Interior-point algorithm is more efficient for smaller net-
work sizes in comparison to the Simplex method. Lastly a comparison over the the
three algorithms is to be made where the number of iterations is compared.

Figure 12: Number of iterations plotted against size of the Manhattan grid. K denotes the
number of commodities, D the total flow for each commodity and I the number of trials. The
error bars signify standard deviation.

In Figure 12 the implemented and Simplex algorithms perform comparably and the
Interior-point method is near constant over increasing network size.

19
4.1.2 Commodities
The computation times for solving multi-commodity flow problems are averaged over a
number of trials I over an increase in the number of commodities K in the network.
This increase in complexity of the network is performed over a constant size of its
Manhattan grid of N = 4. A size of N = 4 is motivated by feasible computation time
for the implemented method.

Figure 13: Average computation time plotted against the number of commodities in the
multi-commodity flow problem. N denotes the size of the grid, D the total from for each com-
modity and I the number of trials. The error bars signify standard deviation.

In Figure 13 it can be observed that the computation time of the implemented algo-
rithm once again is orders of magnitude larger than that of the built-in versions. It
may be observed that the computation time for the implemented method increases
drastically even for a relatively low number of commodities. In Figure 14 a small num-
ber of commodities is shown in order to differentiate the built-in methods.

20
Figure 14: Average computation time plotted against the number of commodities in the
multi-commodity flow problem. N denotes the size of the grid, D the total from for each com-
modity and I the number of trials. The error bars signify standard deviation.

It can be observed in Figure 14 that the implemented method has a computation time
in the same order of magnitude for very few commodities. Furthermore the Interior-
point method is more efficient for N = 4. When only observing the built-in methods,
the size of the Manhattan grid may be increased to N = 10, which is shown with a
larger range of commodities in Figure 15 below.

Figure 15: Average computation time plotted against the number of commodities in the
multi-commodity flow problem. N denotes the size of the grid, D the total from for each com-
modity and I the number of trials. The error bars signify standard deviation.

In Figure 15 the computation time for the Simplex algorithm is significantly lower
than that of the Interior-point method. The corresponding results for a smaller span
of commodities is shown is Figure 16 below.

21
Figure 16: Average computation time plotted against the number of commodities in the
multi-commodity flow problem. N denotes the size of the grid, D the total from for each com-
modity and I the number of trials. The error bars signify standard deviation.

The plots in Figure 16 shows that the built-in Interior-point algorithm is more efficient
than the Simplex algorithm for fewer commodities. Lastly a comparison over the the
three algorithms is to be made where the number of iterations is compared.

Figure 17: Number of iterations plotted against the number of commodities in the multi-
commodity flow problem. N denotes the size of the grid, D the total from for each commodity
and I the number of trials. The error bars signify standard deviation.

In Figure 17, similarly to Figure 12, the implemented and Simplex algorithms per-
forms comparably and the Interior-point method is near constant over increasing net-
work size.

22
4.2 Analysis
4.2.1 Network size
It may be noted that in Figure 8 and Figure 9 that the computation time of the im-
plemented method increases exponentially for increasing network sizes N . The imple-
mented method is meanwhile equally as or more efficient for the smallest networks.
An increase in the computation time for the built-in methods can also be observed
in Figure 10 and Figure 11 with network sizes of N ≥ 10. It is also clear that the
Interior-point algorithm performs better than the Simplex algorithm when increasing
the size of the Manhattan grid N for a constant number of commodities.
In Figure 12 the comparison is made over the number of iterations needed for conver-
gence. It is apparent that the implemented algorithm and the built-in Simplex algo-
rithm are comparable. Meanwhile, the Interior-point method has very few iterations,
since it functions in a different way.

4.2.2 Commodities
In Figure 13 and Figure 14 it can be observed that the implemented algorithm neces-
sitates computation times several orders of magnitude higher for higher numbers of
commodities than that of the built-in algorithms. The implemented algorithm is not
comparable with the built-in algorithm for a low number of commodities in the net-
work.
When focusing on the built-in algorithms in Figure 15 and Figure 16, it is observable
that the Simplex algorithm requires lower computation time than the Interior-point
method for networks with many commodities, K ≥ 10. Although, for lower number of
commodities, K ≤ 8, is the Interior-point algorithm is faster than the Simplex algo-
rithm.
In Figure 17 the comparison is made of the number of iterations needed for conver-
gence. The implemented algorithm increases faster than the Simplex algorithm, al-
though they are still comparable in order of magnitude. The Interior-point method is
not distinguishable on the scale, a consequence of it being completely different from
Simplex, and therefore direct comparison is not reasonable.

23
5 Discussion and Conclusion
The main conclusion that can be drawn from the results is that among the algorithms
tested there is no superior algorithm for solving multi-commodity problems. Which
algorithm performs the best is dependent on the parameters of the problem at hand,
and it may prove to be difficult to know which method is most efficient a priori. The
network size and amount of commodities both had a clear effect on the performance
of the solvers, where the interior-point method performed better on bigger network
sizes while the Simplex method performed better when the amount of commodities
increased. The results are reasonable since both algorithms have polynomial-time aver-
age case complexity.7
The reason that one algorithm performs better than another is difficult to understand
due to the complicated nature of the algorithms and the solution-spaces of the prob-
lems. Furthermore, the large amount of quantities and their multivariate relationships
makes it hard to test different hypothesis as to why. Although finding the cause of the
performance difference is difficult and becomes a guessing game, there are some theo-
ries which may be worth mentioning.
A theory as to why the interior-point method performs better when the network size
increases, is the fact that the amount of iterations is constant and the matrix is sparse.
This means that taking many variables into regard at once is easy if they are zero,
since they do not affect things like a gradient. Furthermore, when the size is increased,
the matrix increases with around N 4 , which in practice will be doubled again for the
phase 1 problem. This means that the Simplex method will have many columns to
take into account which could slow down performance. The reason that Simplex per-
forms better with more commodities networks may stem from the fact that the interior-
point method has to consider many variables and for a complex network of commodi-
ties its difficult to compute things like the gradient when there are many nonzero in-
dices.
While solving linear programs in practice, it might be wise to test several methods to
discover which one performs best for the given problem structure. Although caution
must be taken with regard to which problems the choice of algorithm is based on. For
example, if the performance of the multi-commodity problem is tested on a small grid
size, the interior-point method would come out on top. If instead the network size is
increased for the actual simulation it would be considerably slower than using Simplex.
Therefore it might be good idea to test for several network sizes to sense the general
trend of the computation times before deciding which method to use for the linear pro-
gramming problem at hand. An important conclusion is to choose algorithms carefully
when working with big linear problems.
The choice of comparing the algorithms is based on the relatively arbitrary metric of
computation time which was motivated by the lack of other metrics. An other alter-
native could have been the number of iterations needed for convergence. This met-
ric would however fair ever worse than computation time which can be observed in
Figure 12 and Figure 17. The reason for this is that different methods requires differ-
ent number of iterations for convergence and each iteration may demand a different
7
Nash and Sofer, Linear and nonlinear programming.

24
amount of computation time.

5.1 implemented Simplex method

The naive implementation of Simplex did not stand the test and performed signifi-
cantly worse than the other algorithms. This most likely stems from the fact that state
of the art Simplex methods such as MATLAB use smart transformations and tricks
to reduce the computation time. For example, although sparse matrices were used in
the implementation, there probably is a lot of time to be saved when further consider-
ing the results in Figure 8 and Figure 13. An indication that MATLAB probably uses
some sort of pre-processing is shown by Figure 9, in which, for a small network (N=2)
and few commodities, the implemented Simplex performs better. This is probably be-
cause the implemented version is not pre-processed and therefore is faster for small
problems. Furthermore, there are different types of Simplex methods, some of which
utilize duality to increase performance. The amount of iterations was also higher, but
the difference was not as extreme as for the computation time. The reason for this is
that there are quicker and more effective pivot rules than Bland’s rule, which may re-
duce the amount of iterations significantly. As a result of the papers scope, such pivot
rules were not explored.8

5.2 Further Research

On improving the computation time further, there are several things that need con-
sideration. In the case that finding a basic feasible solution is what is slowing down
the Simplex method for large amounts of commodities, it might decrease computation
time if one used interior point to solve the phase I problem. A combination of different
optimization methods may even be faster than the individual methods by themselves.
Furthermore, it would be interesting to look into other types of Simplex and interior
point methods, which might be more suitable for the problem in question. Looking
into other classes of linear program solvers, such as column generation, may also give
interesting results.
When it comes to understanding which method to use for which multi-commodity
problem, there is also much research left to be done to discover and map all multivari-
ate relationships. It would be interesting to look at more of the quantities in detail,
and also change the type of network as to add more complexity, which would more
closely mirror real life. For example, it would be interesting to see how the algorithms
performed for a real road network. Discovering good rules of thumb for which algo-
rithm to use could save computation time when working with this type of problem in
the future.

8
Nash and Sofer, Linear and nonlinear programming.

25
References
[1] Cynthia Barnhart, Niranjan Krishnan, and Pamela H. Vance. “Multicommodity
flow problems”. In: Encyclopedia of Optimization. Ed. by Christodoulos A. Floudas
and Panos M. Pardalos. Boston, MA: Springer US, 2009, pp. 2354–2362.
[2] Igor Griva. Linear and nonlinear optimization. 2nd ed.. Philadelphia: Society for
Industrial and Applied Mathematics, 2009.
[3] Amol Sasane and Krister Svanberg. Optimization. Stockholm: Department of
Mathematics, Royal Institute of Technology.
[4] Robert G. Bland. “New Finite Pivoting Rules for the Simplex Method”. In: Math-
ematics of Operations Research 2.2 (1977), pp. 103–107.
[5] Stephen G. Nash and Ariela Sofer. Linear and nonlinear programming. New York
: McGraw-Hill, 1996.
[6] The MathWorks Inc. linprog Solve linear programming problems. url: https :
//se.mathworks.com/help/optim/ug/linprog.html (visited on 05/19/2022).

26
TRITA –

www.kth.se