Large scale column generation for multi commodity flow
problems: Application to transport optimization
Benedikt Lienkampa,∗ , Maximilian Schiffera,ba TUM School of Management, Technical University of Munich, Munich, Germany

benedikt.lienkamp@tum.de, schiffer@tum.deb Munich Data Science Institute, Technical University of Munich, Munich, Germany∗ Corresponding author

Extended abstract submitted for presentation at the 11th Triennial Symposium on
Transportation Analysis conference (TRISTAN XI)
June 19-25, 2022, Mauritius IslandApril 4, 2022
Keywords: multi commodity flow, column generation, time expanded networks1INTRODUCTION

According to the United Nations, the global urban population will increase by 900 million from2020 to 2030, reaching a total of 5.2 billion (UN, 2018). This population increase in combination
with already high congestion in cities all around the world calls for novel transportation concepts
to keep up with future urban travel demand. The rise of autonomous vehicles and shared intermodal transportation modes offer an opportunity to mitigate the rising congestion by smarterroute choice. To enable smart utilization of these novel transportation modes it is imperative
to be able to optimize route planning in such a large system. In this context, we will focus on
optimizing passenger flow in a capacitated network with fixed vehicle routes.
A number of approaches exist to optimize route planning in large systems. The concept ofpartially time-expanded networks has been successfully applied to service network design problems (Boland et al. , 2017) and allows for iterative refinement to obtain a reduced model size.
Even though this concept can be used to solve service network design problems optimally, it is
not directly applicable to solve large instances. Simulation-based models capture transportation
systems with high accuracy, but are generally not amenable to efficient optimization. Networkflow models are amenable to efficient optimization and capture transportation systems with highaccuracy. Accordingly, they have for example been used to study the interaction of Autonomous
Mobility-on-Demand (AMoD) with the public transportation system (PTS) (Salazar et al. ,2020) and the control of AMoD systems in congested road networks (Rossi et al. , 2018) on amesoscopic level. To the best of our knowledge, so far no computational efficient algorithm existsthat can be applied to optimize such network flow models for manifold passenger movements invery large systems.
Against this background we develop an algorithmic framework which uses a time expanded
network as well as a spatial expansion to model multiple transportation modes. We define aminimum cost multi commodity flow problem (MCMCFP) to solve the routing of a set of travel
requests optimally. The novelty in our approach lies in presenting a scalable column generation
approach which solves the defined MCMCFP optimally. This column generation leverages acyclic pricing approach in combination with a filtering mechanism to a improve computation
time. Our preliminary results show that this algorithm is able to solve multi commodity flow
problems with 3.137 trips in 23 seconds and with 31.371 trips in 9 minutes optimally.
TRISTAN XI Symposium

Original abstract submittal

B. Lienkamp and M. Schiffer

22Methodology

We seek to find optimal routes for a set of passengers in a capacitated transportation systemwith fixed vehicle routes, e.g., bus and subway lines. An instance of our problem contains theschedule of the fixed vehicle routes in our system and a set of passengers. For each route driven
by a vehicle, the vehicle schedule contains information about the stops of the route, arriving
times at the stops and capacity of the vehicle performing the route. For each passenger, trip
information about the origin/destination coordinates and the departure time is known. We wantto find paths for all passengers, such that the sum over all passenger travel times is minimized
and all vehicle capacities constraints are kept.
We model such a transportation system as a multilayered digraph G = (V, A) with a set oftemporal vertices V and a set of temporal arcs A ⊆ V × V. Let R be the set of all routes
and S be the set of stops in our transportation network. The graph G contains a route layer
Gr = (Vr , Ar ) for each
S route r ∈ R and a waiting layer Gs = (Vs , As ) for each stop s ∈ S. Theroute layers GR = r∈R Gr contain temporal vertices andSarcs which represent the routes driven
in our transportation network. The waiting layers GS = s∈S Gs contain temporal vertices andarcs which represent passengers waiting at a stop. Additionally we add a set of transit arcs
AT and walking arcs AW . Transit arcs can be used to switch between different routes in thetransportation network by moving from a route layer Gr , to a waiting layer Gs or from a waiting
layer Gs0 to a route layer Gr0 , with r, r0 ∈ R and s, s0 ∈ S. Walking arcs allow passengers to walk
between stops. Here, temporal vertices of waiting layers Gs and Gs0 of different stops s, s0 ∈ S
are connected. The capacity of the route arcs corresponds to the capacity of the vehicle operating
the route. All other arcs are unbounded. Collecting all definitions, it holds that V = VR ∪ VSand A = AR ∪ AS ∪ AT ∪ AW . This graph represents the temporal and spatial expansion of thetransportation network. Figure 1 shows a simplified example of the network graph’s structure.
To formulate the passenger flow problem as a MCMCFP, we add a source and a sink vertex foreach passenger trip to the graph G. Hereby, we connect the source (S + ) and sink vertices (S − )to vertices of the waiting layers GS . This extends the vertex set of G to V = VR ∪ VS ∪ VS + ∪ VS −and the arc set to A = AR ∪ AS ∪ AT ∪ AW ∪ AS + ∪ AS − . Each passenger travel demand is
now represented as an individual flow going from a source vertex vs+ ∈ VS + to a sink vertex
vs− ∈ VS − with flow demand one.
Let cij be the travel time and uij the capacity of arc (i, j) ∈ A, dpi the vertex demand of vertex
i ∈ V and passenger p ∈ P and xpij the decision variable, whether passenger p ∈ P uses arc(i, j) ∈ A. The vertex demand is defined as
if vertex i ∈ V is the origin of passenger p ∈ P,
1p
di = −1 if vertex i ∈ V is the destination of passenger p ∈ P,0otherwise

Figure 1 – Scematic illustration of the multilayered network structure
TRISTAN XI Symposium

Original abstract submittal

B. Lienkamp and M. Schiffer

3

We can now formulate the continuous relaxation of finding the minimal sum of all passenger trip
travel times as
X Xminxcij xpij
(1a)s.t.

p∈P (i,j)∈AXX
xpji = dpi ,xpij −j∈N − (i)j∈N + (i)X pxij ≤ uij ,p∈P
xpij ≥ 0,i ∈ V, p ∈ P(1b)(LP 1)(i, j) ∈ A(1c)(i, j) ∈ A, p ∈ P (1d)Here, the objective function models the sum over all passenger travel times. The constraint(1b) ensures flow conservation in the network with N + (i)/N − (i) being the outgoing/ingoingneighbourhood of i ∈ V. The remaining constraints enforce the capacity constraints of allvehicles and ensure positive passenger flows.
Since this LP 1 is not scalable for larger sets of passengers because of the exponential growth ofthe model formulation, we reformulate it to apply a column generation approach.
minλ cT Y λs.t.
Yλ≤uΛλ = 1|P|λ ≥ 0|P|

(2a)
(2b)
(2c)
(2d)

(LP 2)In this LP 2 all possible flows for the passengers are stored as columns in the matrix Y ∈ R|A|×n ,where n is the number of all possible flows. The vector u ∈ R|A| is the capacity vector of allarcs a ∈ A and the matrix Λ ∈ R|P|×n is the incidence matrix between the flows of Y and thepassengers p ∈ P. This means that Λ(i, j) = 1, if the flow of the j-th column of Y corresponds topassenger i and 0 otherwise. The decision variable λ ∈ Rn is used to select a convex combinationof all possible flows for each passenger. A convex formulation of valid passenger flows of passengerp ∈ P yields a new valid flow for passenger p.
We add dummy variables with high objective costs for each passenger, initiate the column generation algorithm without knowing any columns of matrix Y and iteratively add new columns
to the linear program by solving the pricing problem LP 3 for each passenger p ∈ P.
minX ps.t.

(c − w∗ )T X p − αp∗
Xp ≤ uBX p = dp
Xp ≥ 0(3a)
(3b)
(3c)
(3d)

(LP 3)

Here, w∗ is the vector of dual variables for the capacity constraints (2b), and α∗ is the vectorof dual variables for the convexity constraint (2c). The matrix B ∈ R|V|×|A| is the vertex-edge
incidence matrix of graph G. This pricing problem can be solved as a shortest path problem if
we neglect constraint (3b), which is already enforced in constraint (2b) of the restricted masterproblem. We solve these shortest path problems via the A∗ algorithm.
To reduce the number of pricing problems called in each iteration, we add a filtering mechanism
to the algorithm, which decides whether the pricing problem of passenger p ∈ P should be solved.
Here, we find all trips, which use or could use saturated arcs and resolve their pricing problems.
If no new paths are found, we recalculate all pricing problems to ensure optimality.
To ensure integral solutions in LP 2, we propose a branch-and-price (B&P) approach. Here, ourB&P algorithm solves LP 2 at each vertex of the branch-and-bound (B&B) tree. In the vast
majority of tested instances, integral solutions where found in the root node of the B&B tree
TRISTAN XI Symposium

Original abstract submittal

B. Lienkamp and M. Schiffer

4or solving the last updated RMP of our algorithm with λ fixed as an integral variable returned
solutions with an optimality gap of under 0.05%. In the case of no integral solution in the root
node, we use the just described approach to find an upper bound and perform a depth first search
approach on our B&B tree.

3

Preliminary results

We test our algorithm on a case study for the city of Munich, where the network graph consists of1008 stops and 362 routes, which results in a network graph G = (V, A) with |V| = 145.119 and|A| = 510.224. The total trip demand set consists of 31.371 travel requests. The passenger datais generated by the travel demand modelling tool MITO (Moeckel et al. , 2019) through a MonteCarlo sampling followed by a nested mode choice model. Figure 2(a) shows the computation time
for subsets of all travel requests ranging from 10-90% with and without the filter. All subset sizes
have been run on 10 different instances. Our algorithm is able to solve multi commodity flow
problems with 2.370 trips in 23 seconds and with 21.326 trips in 10 minutes. Additionally, one
can see the stable computation time on the different subsets, which shows the generalizability of
our algorithm. Figure 2(b) analogously shows the effectiveness of our filter with respect to thenumber of solved pricing problems. Solving these instances with a simple LP (see LP 1) without
a column generation approach is not possible, since even for the smallest instance we need todefine nearly 500 million constraints, such that a standard desktop computer with 16GB of RAM
runs out of memory, even for initializing the model. Our code was written in Python 3.9 using
the callable library of Gurobi 9.1 to solve the RMP.

(a) Distribution of computation time(b) Number of pricing problems solvedFigure 2 – Computational results
Acknowledgements: This work was supported by the German Research Foundation (DFG)
under grant no. 458548755.

References
Boland, N., Hewitt, M., Marshall, L., & Savelsbergh, M. 2017. The continuous-time service network
design problem. Operations Research, 65(5), 1303–1321.
Moeckel, R., Kuehnel, N., Llorca, C., Moreno, A. T., & Rayaprolu, H. 2019. Microscopic Travel demandmodeling: using the agility of agent-based modeling without the complexity of activity-based models.
In: Annual Meeting of the Transportation Research Board.
Rossi, F., Zhang, R., Hindy, Y., & Pavone, M. 2018. Routing autonomous vehicles in congested transportation networks: Structural properties and coordination algorithms. Autonomous Robots, 42(7),1427–1442.
Salazar, M., Lanzetti, N., Rossi, N., Schiffer, M., & Pavone, M. 2020. Intermodal Autonomous Mobilityon-Demand. IEEE Transactions on Intelligent Transportation Systems, 21(9), 3946–3960.
UN, DESA. 2018. Revision of world urbanization prospects. New York: United Nations Department ofEconomic and Social Affairs.

TRISTAN XI Symposium

Original abstract submittal