VIETNAM NATIONAL UNIVERSITY, HO CHI MINH CITY

UNIVERSITY OF TECHNOLOGY
FACULTY OF COMPUTER SCIENCE AND ENGINEERING

MATHEMATICAL MODELLING (CO2011)

Assignment

Stochastic Programming and

Applications

Advisor: Dr. Nguyễn Tiến Thịnh

Students: Ngô Ngọc Tuấn Anh - 2210078
Thái Quốc Bảo - 2210266
Võ Thượng Bảo (Leader) - 2210290
Trần Đức Trí Cường - 2210443
Hoàng Anh Hùng - 2211325

Ho Chi Minh City, December 2023

 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Contents
Member list & Workload 3

1 Theoretical basis 4
1.1 Introduction to Stochastic Programming and Optimization . . . . . . . . . . . . 4
1.2 Stochastic linear programming (SLP) . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Wait-and-see approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2 No waiting approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 One-Stage Stochastic linear programming - No Recourse (1-SLP) . . . . . . . . . 5
1.3.1 APPROACH 1. Use Chance constraint and Acceptable risk . . . . . . . . 6
1.3.2 APPROACH 2. For stochastic constraints T (α) ≤ h(α) . . . . . . . . . . 6
1.4 Generic Stochastic Programming with Recourse . . . . . . . . . . . . . . . . . . . 6
1.5 Two-Stage Stochastic linear programming with Recourse . . . . . . . . . . . . . . 7

2 Problem 1: Industry - Manufacturing 9

2.1 Present the problem and Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Algorithm - Coding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.1 Introduction to GAMSPy . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Source code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Solution to the problem with realistic figures . . . . . . . . . . . . . . . . 12

3 Problem 2: Stochastic Linear Program for Evacuation Planning In Disaster

Responses (SLP-EPDR) 14
3.1 Introduction and Definition of the Issue: . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Problem’s solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2.1 Presenting the evacuation problem in 2-SLP form . . . . . . . . . . . . . . 15
3.2.2 Presenting the evacuation problem with the min-cost flow theory . . . . . 16
3.3 Build the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.1 Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.2 Decision variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.3 Constraints of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3.4 Solutions: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3.4.a Problem 1: Min-cost flow problem . . . . . . . . . . . . . 21
3.3.4.b Problem 2: Min-cost flow problem with time factor . . . . 21
3.4 Ways to solve the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4.1.a Minimum-cost Flow Problem (MCFP) . . . . . . . . . . . 21
3.4.1.b Negative Cost Cycle . . . . . . . . . . . . . . . . . . . . . 22
3.4.2 Residual networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.3 Successive Shortest Path Algorithm . . . . . . . . . . . . . . . . . . . . . 23

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 1/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

3.4.3.a Algorithm Overview . . . . . . . . . . . . . . . . . . . . . 23

3.4.3.b Application of Successive Shortest Path Algorithm . . . . 24
3.4.4 The Cycle Canceling Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4.4.a Algorithm Overview . . . . . . . . . . . . . . . . . . . . . 27
3.4.4.b Application of Cycle Canceling Algorithm . . . . . . . . . 29
3.5 Algorithm Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4 Conclusion - Appreciation 33

5 Appendix 33

References 37

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 2/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Member list & Workload

No. Fullname Student ID Problems Percentage of work

Compile documents, write
1 Ngô Ngọc Tuấn Anh 2210078 reports (Theoretical basis 100%
+ problems 1, 2).
Compile documents, write
2 Thái Quốc Bảo 2210266 reports (Theoretical basis 100%
+ problem 1).

3 Võ Thượng Bảo 2210290 Problems 1,2 100%

4 Trần Đức Trí Cường 2210443 Problem 2 100%

5 Hoàng Anh Hùng 2211325 Problem 1 100%

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 3/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

1 Theoretical basis
1.1 Introduction to Stochastic Programming and Optimization
Stochastic programming is a framework for modelling optimization problems that involve un-
certainty. Whereas deterministic optimization problems are formulated with known parameters,
real world problems almost invariably include some unknown parameters. When the parameters
are known only within certain bounds, one approach to tackling such problems is called robust
optimization. Here the goal is to find a solution which is feasible for all such data and optimal
in some sense. Stochastic programming models are similar in style but take advantage of the
fact that probability distributions governing the data are known or can be estimated. The goal
here is to find some policy that is feasible for all (or almost all) the possible data instances and
maximizes the expectation of some function of the decisions and the random variables.
Briefly, Stochastic programming can be viewed as mathematical programming with random
parameters (e.g., random variables, the variables whose possible values depend on the outcomes
of a chance phenomenon). The existence of random variables in the model being considered will
open up countless ways to pose different problems, reflecting various aspects of the issue being
addressed.

1.2 Stochastic linear programming (SLP)

A Stochastic Linear Programming (SLP) is

M inimize Z = g(x) = f (x) = cT .x, s.t. Ax = b and T x ≥ h

With x = (x1 , x2 , . . . , xn ) are called decision variables, certain matrix A and vector B (for deter-
ministic constraints), and with random parameters T, h in T x ≥ h define chance or probabilistic
constraints.
Oftentimes, decision variables x is dependent on random vector ω, which partially represents
for data uncertainty of the problem. How do we solve this problem, if w is really random variable?

1.2.1 Wait-and-see approach

Suppose it is possible to decide about the decision variables x after the observation of the
random vector ω. After having observed the random vector ω, we can solve the problem with
the decision variables vector x being clearly determined.

1.2.2 No waiting approach

In this approach, we need to decide on x before knowing the values of ω. There are 2 main
ways to do it: (a) Guess at uncertainty, (b) Probabilistic Constraints.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 4/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

(a) Guess at uncertainty: We will guess few reasonable values for ω. Three main ways
are:

• Unbiased: Choose mean values for each random parameters ω.

• Pessimistic: Choose worst case values for ω.

• Optimistic: Choose best case values for ω.

(b) Probabilistic Constraints: Random parameters T , h in T x ≥ h (presented above)

define chance or probabilistic constraints.

1.3 One-Stage Stochastic linear programming - No Recourse (1-SLP)

We now start with One-Stage SLP, shortly means stochastic LP with one-stage, no recourse
actions.

Definition 1 Consider the following program LP (α) that is parameterized by the random vector

α:

Pn
Minimize z = g (x) = f (x) = cT x = j=1 cj xj

s.t Ax = b, (certain constraints)

and T x ≥ h, (stochastic constraints)

with assumptions that:

1. Matrix T = T (α) and (vector) h = h(α) express uncertainty via stochastic constraints.

T (α)x ≥ h(α) ⇔ α1 x1 + . . . + xn xn ≥ h (α)

2. Values (T, h) not known: they are unknown before an instance of model occurs, h(α) de-
pends only on random aj .

3. Uncertainty is expressed by probability distribution of random parameters (αj ) = α so

deterministic LP is the degenerate case of Stochastic LP when αj are constant.

Fundamental assumption. We know a (joint) probability distribution of data. Hence the first
approach gives Probabilistic (Chance) Constraint LP.
The Scenario Analysis. Not perfect, but useful, is the second approach. The scenario ap-
proach assumes that there are a finite number of decisions that nature can make (outcomes of
randomness). Each of these possible decisions is called a scenario.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 5/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

1.3.1 APPROACH 1. Use Chance constraint and Acceptable risk

We can replace T x ≥ h by probabilistic constraints P[T x ≥ h] ≥ p2 for some prescribed

reliability level p ∈ (0.5, 1) (to be determined by problem owner).
Risk then is taken care of explicitly, if define an:

acceptable risk rx := P[N ot(T x ≥ h) = P (T x ≤ h) ≤ 1 − p

The chance constraint T x ≤ h implies that the acceptable risk rx is less than a specified
maximal 1 − p ∈ (0, 1).

Definition 2 Stochastic LP or 1 − SLP with Probabilistic Constraints is defined by a ran-

dom coefficients α = (α1 , α2 , ..., αn) in chance constraints and a linear objective f (x):
n
X
SP : min Z, Z = f (x) = cT x = cj x j , cj ∈ R (1)
x
j=1

Ax = b (x = (x , x , . . . , x ) makes decision variables)
1 2 n
s.t
P = [T x = α)c ≤ h ≤ 1 − p (0 < p < 1)

Note: we use parameter vectorα = [α1 , α2 , ...] in general, and denote ω = [ω1 , ω2 , ..., ωS ]
specifically for states s called scenarios.

1.3.2 APPROACH 2. For stochastic constraints T (α) ≤ h(α)

Use Scenario analysis of (α) ≤ h(α).

For every scenario (T s ; hs ) , s = 1, ..., S, solve:

Minimize{f (x) = cT · x; S.t · Ax = b, T s x ≤ hs }

This kind of program targets a specific linear objective while accounting for a probability
function associated with various scenarios. Hence, we find an overall solution by looking at the
scenario solutions xs (s = 1, ..., S).
Advantage: each scenario problem is an LP.
Disadvantage: discrete distribution →
− mixed-integer LP model.

1.4 Generic Stochastic Programming with Recourse

The two-stage stochastic programm, extended from Definition 2 has the form:

2 − SP : minx g(x) with g(x) = f (x) + Eω [v(x, ω)]

where x = (x1 , x2 , . . . , xn ) is the first stage decision variables, f (x) can be linear or not, a part
of the grand objective function g(x).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 6/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

*The mean Q(x) := Eω [v(x, ω)] of a function v : Rn × RS → R upon influences of scenarios

ω. Q(x) is the optimal value of a certain second-stage problem:

min q y| subject to T · x + W · y = h (2)

y∈Rp

Vectors α = α(ω) and y = y(ω) are named correction, tuning or recourse decision variables,
only known after the experiment e.
Briefly we minimize total expected costs g(x) = f (x) + Q(x) while satisfying W · y(ω) =
h(ω) − T (ω) · x.
Here W is called m × p recourse matrix, and we begin with simple case of m = 1, q is the
unit recourse cost vector, having the same dimension as y, and y = y(ω) ∈ Rp .

1.5 Two-Stage Stochastic linear programming with Recourse

Two-stage stochastic programming is one mostly used kind of stochastic programming prob-
lems where the decision-maker has to make decisions in two stages (two different times) for a
fixed phenomenon with uncertainty. The first stage decision plays a vital role here since this
decision has to anticipate based on some random parameters obtained from past experience or
some sort of survey. And the second stage decision is based on the experiment result or some an-
ticipation. Two-stage stochastic programming with recourse discuss the phenomena with simple
second-stage decisions which could be fixed under some circumstances. In making the first stage
decisions, second stage forecast with some sort of possibility (mathematical expectation) has to
be taken into account.
One way to think about stochastic problems is to require the decision maker to make a
decision now and then to minimize the expected costs of the consequences of that decision. This
paradigm is called the recourse model. The simplest form of the recourse model has two stages:
a decision is made in the first stage, then the realization of the uncertain parameters is revealed
at the start of the second stage and recourse actions can be taken given this new information. A
recourse decision can then be made in the second stage that compensates for any bad effects that
might have been experienced as a result of the first-stage decision. Recourse Program means the
program needs some new decisions or recourse actions after the uncertainty is no longer concealed.

Let x ∈ Rn and y ∈ Rn be two variables and let the set of all realizations of the unknown

data be given by ω, Ω = {ω1 , ..., ωn } ⊆ Rr , where r is the number of the random variables

representing the uncertain parameters. Then the stochastic program is given by

Minx z = cT x + Eω [Q(x, ω)]

s.t Ax = b

x≥0

where Q(x, ω) is the optimal value of the second-stage problem.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 7/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Q(x, ω) = Miny qωT y(ω)

s.t Tω + Wω y(ω) = hω

yω ≥ 0, ∀ω ∈ Ω

The first three lines define the first-stage problem and the last three lines define the second-
stage [Link] the first stage:

• x is the decision variable.

• cT represents the cost coefficients of the objective function.

• Eω [Q(x, ω)] denotes the expected value of the optimal solution of the second stage problem.

• A denotes the coefficients.

• b the right-hand side of the first stage constraints.

In the second stage:

• y is the decision variable.

• T represents the transition matrix.

• W the recourse matrix (cost of recourse).

• h the right-hand side of the second stage constraints.

Note: that all parameters and the decision variable of the second stage are dependent on
the specific realization of the stochastic data ω. The objective variable z is also a random vari-
able, since it is a function of ω. As a random variable cannot be optimized, stochastic solvers
automatically optimize the expected value of the objective variable z.
In the first stage, a decision has to be made before the realization of the uncertain data is
clear. The optimal solution of the first stage is fixed and only then it will become known which
values the uncertain parameters will take. Given the fixed solution of the first stage and the
new data, recourse action can be taken in the second stage and the optimal solution determined.
Each possible realization of the uncertain data is represented by ωs ∈ Ω and is called a scenario.
The objective is to find a feasible solution x that minimizes the total cost, namely the sum of
the first-stage costs and the expected second-stage costs.

The formulation above looks a lot messier than the deterministic LP formulation that we

discuss elsewhere. However, we can express this problem in a deterministic form by introducing

a different second-period y variable for each scenario. This formulation is called the deterministic

equivalent:

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 8/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Pn
Min cT x + i=1 pi qiT yi

s.t Ax =b

Ti x + Wi yi = hi , i = 1, ..., n

x ≥0

yi ≥ 0, i = 1, ..., n

where n is the number of scenarios and pi is the probability of the scenario i ’s occurrence.

For our n-scenario problem, we have:

Min cT x + p1 q1T y1 + p2 q2T y2 + . . . + pn qnT yn

s.t Ax =b

T1 x + W1 y1 = h1

T2 x + W2 y2 = h2
.. ..
. .

Tn x + Wn yn = hn

x ≥0

yi ≥ 0, i = 1, ..., n

Notice that the non-anticipativity constraint is met. There is only one first-period decision,
x , whereas there are n second-period decisions, one for each scenario. he first-period decision
cannot anticipate one scenario over another and must be feasible for each scenario. That is, Ax
and Ti x + Wi yi for i = 1, .., n. Because we solve for all the decisions, x and yi simultaneously, we
are choosing x to be (in some sense) optimal over all the scenarios.
By solving the above formula, we obtain an optimal solution x of the first-stage problem and
optimal solutions yi of the second-stage problem for each scenario ωi , i = 1, ..., n. Given x, each
yi gives an optimal second-stage decision corresponding to a realization ω = ωi of the respective
scenario.

2 Problem 1: Industry - Manufacturing

2.1 Present the problem and Solutions
PROBLEM. Produce n products satisfying production ≤ demand. Use the 2-SLPWR model
given in Equation 4 when n = 8 products, the number of scenarios S = 2 with density ps = 1/2,
the number of parts to be ordered before production m = 5, we randomly simulate data vector
b, l, q, s and matrix A of size n × m. We also assume that the random demand vector ω = D =
(D − 1, D2 , . . . , Dn ) where each ωi with density pi follows the binomial distribution Bin(10, 1/2).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 9/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

REQUEST. Build up the nummerical models of Equation 4 with simulated data. Find the
optimal solution x, y ∈ Rm , and z ∈ Rn by suitable soft (as GAMSPy).
Consider a situation where a manufacturer produces n products. There are in total m different
parts (or sub-assemblies) which have to be ordered from third-party suppliers. A unit of product
i requires aij units of part j, where i = 1, .., n and j = 1, .., m. Of course, aij may be zero for
some combinations of i and j. The demand for the products is modeled as a random vector D =
(D1 , ..., Dn ). Before the demand is known, the manufacturer may pre-order the parts from outside
suppliers at a cost of cj per unit of part j. After the demand D is observed, the manufacturer
may decide which portion of the demand is to be satisfied, so that the available numbers of parts
are not exceeded. d. It costs additionally li to satisfy a unit of demand for product i, and the
unit selling price of this product is qi . The parts not used are assessed salvage values sj < cj .
The unsatisfied demand is lost.
Suppose the numbers of parts ordered are equal to xj , j = 1, ..., m. After the demand D
becomes known, we need to determine how much of each product to make. Let us denote the
numbers of units produced by zi , i = 1, ..., n, and the numbers of parts left in inventory by
yj , j = 1, ..., m. For an observed value (a realization) d = (d1 , ..., dn ) of the random demand
vector D, we can find the best production plan by solving the following linear programming
problem:
The second-stage problem:
For an observed value d = (d1 , d2 , . . . , dn ) of the above random demand vector D, we can
find the best production plan by solving the following stochastic linear program (SLP) with:

• Decision variables z = (z1 , z2 , . . . , zn ) - the number of units produced.

• Decision variables y = (y1 , y2 , . . . , ym ) - the number of parts left in inventory.

• The additional cost to meet a unit of product i, li .

• The selling price of a product i, qi .

• The resale price of the residual part j, (sj ≤ bj ),

Therefore, with the above variables and data, the profit will be equal to the profit margin plus
Pn Pm
the resale amount: i=1 (li − qi )zi + j=1 sj yj . In order to create an only problem between two
stages, we will switch from max to min: max = −min, which means, to maximize profits, the
company will strive to achieve the smallest loss cost.
n
X m
X
LSP : min Z = (li − qi )zi − sj yj (3)
z,y
i=1 j=1

Pn
• i=1 (li − qi )zi represents the production cost.
Pm
• j=1 represents the number of residual parts for resale.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 10/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

With real data qi ≥ li , the first sum will be negative, which leads to the result of the first
equation will be negative, thus, we have some expected profit.
Let xj , j = 1, . . . , m be the number of parts need to be ordered before production, the variables
in the equation will have to satisfy these constraints:

 y ≥ 0, j = 1, . . . , m
 j


Pn
yj = xj − i=1 aij xi , j = 1, . . . , m


0 ≤ x ≤ d , i = 1, . . . , n

i i

In which, yj is the number of residual parts (which is non-negative) and is calculated by the
number of ordered parts minus the parts used for production by demand, zi is the number of
products, which is non-negative and not exceeds demand di .
The whole second-stage model can be represents:

T T
minz,y Z = c z − s y with c = (ci := li − qi ) are cost coefficients



y = x − AT z with A = [aij ] is a m x n matrix (4)


0 ≤ z ≤ d,

y≥0

c, s and A are the matrices of production cost, the cost of reserves for production and the amount
of raw materials needed for each product. The size of each matrix depends on the input data.
The first-stage problem
The 2-SLPWR model is based on a general rule: production ≥ demand. Approaching based
on distribution, we set Q(x) := E[Z(z, y)] = Eω[x, ω] representing the optimal value of problem.
The notation b = (b1 , b2 , . . . , bm ) is constructed according to the pre-order cost bj on a unit
of part j (before knowing the demand). The quantities xj are determined from the following
optimization problem:

min g(x, y, z) = bT x + Q(x) = bT x + E[Z(z)], (5)

Pn
where Q(x) = Eω [Z] = i=1 pi ci zi is taken according to the probability distribution of ω = D.
The first part of the objective function represents the pre-order cost and variable x. Con-
versely, the second part represents the expected cost of the optimal production plan, calculated
by the updated order quantity z, which has used the random demand D = d with their density.
In summary, at this stage, we solve the problem after adding the probability factor ω to the
parts of the objective function. The decision variables include the vectors x, y ∈ Rm and z ∈ Rn .
After observing the demand D, the manufacturer can decide which part of the demand will be
met so that the number of available parts does not exceed. The additional cost is li to meet a
unit demand for product i, and the selling price of this product is qi . After knowing the demand
D, we determine the quantity of each product to be produced. The unused parts are converted
into reserve value sj , for vector s = s1 , s2 , . . . , sm .

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 11/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

2.2 Algorithm - Coding

2.2.1 Introduction to GAMSPy

General algebraic modeling system (GAMS) is a high-level modeling system for mathematical
optimization. It is designed for modeling and solving different kinds of problems—linear, non-
linear, as well as mixed-integer optimization problems. The main characteristics of this program
are also its availability for use on various computer platforms, while models are portable from
one platform to another. Therefore, it is very popular for scientific research. This program allows
users to implement different algorithms and to combine different solvers.
GAMSPy, on the other hand, is the combination of GAMS and the well-known programming
language Python. In this case, Python acts as a communication layer, allow us to easily build
complex mathematical models.
GAMS/GAMSPy model consists of the following main elements:

Element Explanation
Sets Used to define the indices in models
Data Can be defined in different types: parameters, tables, or scalar
Variables Should be defined before solving the model
Equations Define the relationship between the data and variables
Model The mutual name for a set of equations, which also include objective function,

Table 1: Elements in GAMS/GAMSPy

2.2.2 Source code

In the Appendix section, Problem 1’s Source code shows group’s source code for Problem
1. The source code and results can also be found at this Google Colab link.

2.2.3 Solution to the problem with realistic figures

The tables below show the figures product selling prices (q), production costs (l), ordering
costs (b), cost of each part left (s), the number of parts needed to make a product (a) and the
demand of products in each scenario (d):

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 12/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

After running our code, we get these results:

• The value of the objective function:

• Number of ordered parts:

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 13/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

• Number of parts left in inventory:

• Number of productions:

3 Problem 2: Stochastic Linear Program for Evacuation

Planning In Disaster Responses (SLP-EPDR)
3.1 Introduction and Definition of the Issue:
Extreme natural disasters, such as earthquakes, hurricanes, tsunamis,... and unnatural ones,
such as terrorist attacks, political issues, internal war, or the most recent one, pandemic, had and
will affect the global with little to none warning beforehand and leave behind considerable damage
and casualties. Researchers and scientists have worked to find a way to limit the damage as much
as possible. Therefore, emergency response was invented, with the main goal is to provide shelter
and assistance to affected people as soon as we can. The proposed solution involves modeling a
travel network using nodes and links. The goal is to efficiently move people from hazardous zones
to safe areas in the least amount of time using a capacity-cost network G(V, A, C, U, D). Here,
V represents the set of nodes, A represents the set of links with random travel times, C(i, j)
denotes the travel time on link (i, j) ∈ A, denoted by cij , U (i, j) represents the capacity of link
(i, j) denoted by uij and D(i) represents the flow at node i ∈ V , denoted by di .
The evacuation process is divided into two stages. In the first stage, it is assumed that the
affected person cannot obtain information about the extent of the disaster. This part of the path
is affected by real-time monitoring devices. Therefore, the problem is solved in two stages - one
is to predict and plan before the disaster occurs, and the other is to adapt the plan when disaster
information and random situations occur over time.
Estimating the impact and damage caused by disasters is inherently difficult due to the near
impossibility of predicting the full extent of their repercussions. To address this issue, the Two-
Stage Stochastic Linear Programming (2-SLP) model is introduced for evacuating the affected
population. The model’s primary objective is to efficiently relocate people from hazardous areas
to safe zones, minimizing the evacuation time.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 14/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

3.2 Problem’s solution

3.2.1 Presenting the evacuation problem in 2-SLP form

The evacuation process is described by assumptions and is divided into 2 stages depending
on the time people get access to the correct information. The general purpose is to build an
optimized evacuation plan in the first stage by reviewing the plans to cope with the disaster in
stage two.
The problem aimed at resolution involves modeling a travel network with nodes representing
locations and links indicating paths. The evacuation process is outlined through assumptions and
divided into two stages. In the first stage, it is assumed that affected individuals lack access to
information about the disaster’s extent and the path of damage. After reaching a time threshold
T, they can obtain accurate road network information through various real-time monitoring
devices. Therefore, addressing the problem involves predicting a plan in the first stage when a
disaster occurs, and in the second stage, adapting the plan based on known disaster information
to accommodate random situations occurring at the specified time.
The evacuation problem under the 2-SLP form is simply illustrated by the network shown
below:

Figure 1: Example of 2-SLP evacuation planning problem

Assume nodes 1 and 8 correspond to the disaster-affected area and the safe zone, respectively.
Four vehicles (a, b, c, d) are designated for evacuating the population to the safe area.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 15/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

In the initial stages of the event, four evacuation vehicles adhere to a predetermined plan:
Vehicles a and b follow routes 1 → 2 → 3 → 4 → 6 → 8, while vehicles c and d traverse routes
1 → 3 → 5 → 6 → 7 → 8.
Following the event (beyond the time point T ∼), once the actual situation is known, the
plan is adjusted accordingly. Assuming two random scenarios to illustrate the planning:
Evacuation plan in Scenario 1:
Vehicles a and b move along routes 1 → 2 → 3 → 4 → 7 → 8.
Vehicles c and d move along routes 1 → 3 → 5 → 6 → 8.
Evacuation plan in Scenario 2:
Vehicles a and b move along routes 1 → 2 → 3 → 4 → 6 → 8.
Vehicle c moves along routes 1 → 3 → 5 → 6 → 8.
Vehicle d moves along routes 1 → 3 → 5 → 6 → 7 → 8.
Both scenarios showcase adjustments based on the actual conditions observed post-event,
underscoring the necessity for flexibility in the evacuation strategy.

3.2.2 Presenting the evacuation problem with the min-cost flow theory

We can construct a detailed evacuation plan for people from dangerous areas to safe areas
using the min-cost flow problem. The goal is to evacuate people from dangerous areas to safe
areas with the minimum evacuation time in the system G(V, A, U, C, D), where:

• V represents the set of nodes.

• A represents the set of links.

• C(i, j) represents the travel time of the link (i, j) ∈ A, denoted as cij .

• U (i, j) represents the travel time of the link (i, j) ∈ A, denoted as uij .

• D(i) is the flow at node i ∈ V , denoted as di .

Assume that the travel time on the link (i, j) ∈ A is the same if the number of affected people on
the road does not exceed the maximum capacity (travel time does not change according to the
number of people). The sources in the network represent the dangerous areas, the sinks represent
the safe areas, and the other nodes are the intersections of the network. The links in the network
represent the movement of the affected people.
This model differs from the usual min-cost flow problem: The 2-SLP model has recourse,
representing the situation where the first and second stages occur at different times in the same
evacuation network. In the first stage, the travel time and capacity of the link are only predicted
by probability, people must be evacuated from the sources to other nodes before knowing the
travel time and capacity of the link in the second stage. In the second stage, the evacuation plan
will be determined based on the specific travel time and capacity.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 16/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Note that the evacuation plan in the first stage may not be feasible for the actual situation,
we can solve it by allowing shorter evacuation time in the second stage. In this case, the result
function will include punish costs in the first stage and the expected value of recourse costs in
the second stage.

3.3 Build the model

3.3.1 Symbols

For the purpose of building the model, Table 2 summarizes the symbols used in the mathe-
matical formulas:

Symbols Definitions
V Set of nodes
A Set of links
i, j Indexes of node i, j ∈ V
(i, j) Indexes of directional link (i, j) ∈ A
s Index of scenario
S Total of the scenarios that can happen
v Value of source node
T̃ Time mark
T Total of time
uij Flow of link (i, j)
usij (t) Flow of link (i, j) in the scenario s at time t
csij (t) Moving time of link (i, j) in the scenario s at time t
µs Probability in the scenario s

Table 2: Indexes and variables used in mathematical formulas

3.3.2 Decision variables

These two types of decision variables presented in Table 2 will be used to build the 2-SLP
evacuation model:

Symbols Definitions
xij Flow on link (i, j)
s
yij (t) Flow on link (i, j) in the scenario s at time t

Table 3: Decision variables used in mathematical formulas

In the first stage, a binary decision variable xij is used to represent the flow on link (i, j).
s
The second decision variable yij (t) is defined to represent the flow on link (i, j) in scenario s at
time t.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 17/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

3.3.3 Constraints of the problem

*The first stage The problem is described as follows: The purpose is to evacuate people
from dangerous places to safe places. In reality, there are more than one sources and sinks in the
evacuation network. Therefore, the physical network with multiple sources should be converted
to an equivalent network with a single supersource, even to multiple sinks. It means a supersource
k is added to the network, and meanwhile the dummy arcs (k, i) should be added. Where K is
the set of source nodes, and let the capacity on the dummy arc be the value of the supply at node
P
i, i.e., uki = di , i ∈ K. Hence the supply value at the supersource can be set as d − k = i∈K di .
In the first stage, a feasible evacuation should be determined from supersource to supersink. The
flow on each link should satisfy the flow balance constraint below:
X X
xij − xji = di (1)
(i,j)∈A (j,i)∈A

where di is a parameter with the following definition:


 v if i = s
 i


di = −vi if i = t



0 otherwise

Meanwhile, the flow on each link must also satisfy the capacity constraint:

0 ≤ xij ≤ uij , ∀(i, j) ∈ A (2)

On a separate note, the flow balance constraint may generate a path with loops and sub-tours
if there are potential loops in the evacuation network. To eliminate loops on the generated
physical evacuation path, the link penalty pij , (i, j) ∈ A is particularly introduced. With this
consideration, the penalty function can be defined as:
X
f (X) = pij xij (3)
(i,j)∈A

where X = {xij |(i, j) ∈ A}.

Because we need to strive for achieving the smallest value of the objective cost, we must
minimize the value of the penalty function. Therefore, the goal of stage one will be minimizing
f (x).
*The second stage
Because the arc travel time and capacity will vary with the departure time, when adding a
supersource the travel times on the dummy arcs are cki = di (t), t ∈ {0, 1, . . . , T }. Since the
time of arriving at the sink for each flow may be different from each other, we first add a copy
j for each original sink j, j ∈ D and then a supersink l is added to the network. For each sink,

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 18/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

we add the arc (j, j ′ ) to the network with infinite capacity, i.e., uij = ∞, and the travel time for
each dummy arc is csij (t), ∀t ∈ {0, 1, . . . , T }.
To ensure high accuracy in predicting pre-disaster planning to minimize loss of life and prop-
erty when evacuating people, the problem must ensure an important time requirement. The
affected people will be evacuated along the a priori evacuation plan in the first stage; in other
words, the evacuation plan in different scenarios before the time threshold is the same as the a
priori path. So the problem conditions are:
X
yij (t) = xij , (i, j) ∈ A, s = 1, 2, . . . , S (4)
i≤T

Therefore, the purpose of the second stage boils down to finding the minimum total time it takes
for affected people to evacuate from the dangerous area to a safe place in each scenario s:
X
Q(Y, s) = min csij (t)yij
s
(t) (5)
(i,j)∈A

s.t.
X X
s s ′
yij (t) − yji (t ) = dsi (t), ∀i ∈ V, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S (6)
(i,j)∈A (j,i)∈A
s
0 ≤ yij ≤ usij ∀(i, j) ∈ A, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S (7)
X
s
yij (t) = lij , ∀(i, j) ∈ A, s = 1, 2, . . . , S (8)
t≤T

With:
(5) is the total evacuation time of all traffic flow in the most optimistic way.
(6), (7) are the flow balance condition and traffic volume condition respectively.
(8) is a condition to ensure that the evacuation plan in each scenario in stage two before time
T is consistent with the plan in the initial phase.
From the above two stages, we can come up with an adaptive plan model to minimize the
penalty for the previous evacuation plan and the expected minimum time of the adaptive evacu-
ation plan according to each scenario with a probability µs , s = 1, 2, . . . , S corresponds to each

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 19/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

scenario as follows:
  P 
min (i,j)∈A pij xij + Ss=1 µs . (i,j)∈As csij (t)yij s
P P

 (t)






s.t.

P P
 (i,j)∈A xij − (j,i)∈A xij = di , ∀i ∈ V



0 ≤ xij ≤ uij , ∀i ∈ V (9)


s s ′
(t ) = dsi (t), ∀i ∈ V, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S
P P
 (i ,j ′ )∈A yij


 t t s
(t) − (jt′ ,it )∈As yij

s
≤ usij , ∀(i, j) ∈ A, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S




0 ≤ yij

y s (t) = x ∈ A, s = 1, 2, . . . , S (∗)

P
t≤T̃ ij ij

(*) is essentially an integer programming model which contains two types of decision variables,
s
i.e., X := {xij }(i,j)∈A and Y := {yij }(i,j)∈A{0,1,...,T },s=1,2,...,S . In this model, the coupling con-
straint is a complex constraint, which leads to the model cannot be solved in polynomial time.
s
Therefore, we use the Lagrangian multiplier αi,j (t), (i, j) ∈ A, s = 1, 2, . . . , St ≤ T for the
coupling constraint, and then this constraint can be relaxed into the objective function in the
following form:  
S
X X X
s  s
αij yij (t) − xij 
s=1 (i,j)∈A t≤T̃

After the relaxation of the formula, the model can be formulated as follows:
 PS  P  P P 
P s s S P s s


 min (i,j)∈A pij x ij + s=1 µ s . (i,j)∈As c ij (t).yij (t) + s=1 (i,j)∈A αij t≤T̃ y ij (t) − x ij






 s.t.

P
 P
(i,j)∈A xij − (i,j)∈A xij = di , ∀i ∈ V



 0 ≤ xij ≤ uij , ∀(i, j) ∈ A

s s ′
(t ) = dsi (t), ∀i ∈ V, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S
 P P


 yit jt′ ∈As yij (t) − (j ′ ,it )∈As yij
 t

0 ≤ y s ≤ us , ∀(i, j) ∈ A, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S

ij ij
(11)
It is worthwhile to note that the variables X and Y can be separated from each other in above
relaxed model:
   
X S
X X S
X X X
min pij xij + µs csij (t).yij
s
(t) + s 
αij s
yij (t) − xij 
(i,j)∈A s=1 (i,j)∈As s=1 (i,j)∈A t≤T̃

 
X S
X X X S
X X X S
X X
= pij xij +  µs .csij (t).yij
s
(t)+ s
αij s
.yij (t)− s
αij xij
(i,j)∈A s=1 (i,j)∈A t∈{0,1,2,...,T } s=1 (i,j)∈A t≤T̃ s=1 (i,j)∈A

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 20/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

!  
X S
X S
X X X X
s
= pij − αij (t) xij +  µs csij (t) + s
αij s
(t) yij (t)
(i,j)∈A s=1 s=1 (i,j)∈A t∈{0,1,2,...,T } t≤T̃

3.3.4 Solutions:

By combining similar terms, the relaxed model is decomposed into two subproblems as follows:

3.3.4.a Problem 1: Min-cost flow problem

Obviously, the first sub-problem can be regarded as a min-cost flow problem, and its form is
given as follows: 
P PS s


 min SP 1(α) = (i,j)∈A (pij − s=1 αij (t))xij



s.t.
P P (12)
(i,j)∈A xij − (i,j)∈A xji = di , ∀i ∈ V





0 ≤ x ≤ u , ∀(i, j) ∈ A

ij ij

PS P s
The objective function of subproblem 1 can be defined as gij = pij − s=1 t≤T̃ αij to represent
the generalized cost of each link. Therefore, the subproblem 1 can be solved by the successive
shortest path algorithm. For convenience, the optimal objective value of subproblem 1 is abbre-
viated as ZSP 1 (α)∗ .

3.3.4.b Problem 2: Min-cost flow problem with time factor

The second subproblem of the relaxed model is associated with the decision variable Y , and
its optimal objective value of the problem is abbreviated as ZSP 2 (α)∗ , shown as follows:
 PS P P 
µs .csij (t) + t≤T̃ αsij (t) ysij (t)
P


 min SP 2(a) = s=1 (i,j)∈A t=0,1,2,...,T



 s.t.
(13)
s
(t) − (j ′ ,it )∈As ysij (t′ ) = dsi (t), ∀i ∈ V, t ∈ {0, 1, . . . , T }, s = 1, 2, . . . , S
P P
 (i ,j ′ )∈A yij


 t t s t

P ij

y
t≤T̃ s (t) = x ij , (i, j) ∈ A, s = 1, 2, . . . , S

3.4 Ways to solve the problem

3.4.1 Definitions

3.4.1.a Minimum-cost Flow Problem (MCFP)

The minimum-cost flow problem (MCFP) is an optimization and decision problem to find
the cheapest possible way of sending a certain amount of flow through a flow network.
The problem consists of a directed graph G(V, E).

1. Each node carries a value called bi . If bi > 0, it means the node has a supply of bi , and if
bi < 0, it means the node has a demand of bi .

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 21/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

2. Each edge on the graph will have a capacity uij ≥ 0 (the maximum amount of goods that
can pass through that node) and cij ≥ 0 (the cost of one unit passing through it).

3. Suppose that an edge (i, j) in E carries xij units of flow.

4. The problem requires finding the minimum flow, i.e., it needs to minimize the following
function:
X
xij .cij
(i,j)∈E

Also, it needs to satisfy the following requirement: 0 ≤ xij ≤ uij :

X X
xij − xki = bi
j:(i,j)∈E k:(k,i)∈E

This means that the sum of the outflows from node i minus the sum of the inflows to node
i will equal the supply or demand of node i.

3.4.1.b Negative Cost Cycle

Negative Cost Cycle is a cycle in a network such that the total cost of the edges in the cycle
is negative.

3.4.2 Residual networks

When solving a streaming problem on a network, after each augmentation step, you update
the network to reflect changes in the stream. A residual network is a representation of the current
network after each update step.
Consider the following concept: Let G be a network, and x be a feasible solution to the
minimum cost flow problem. The residual capacity of an edge (i, j) is defined as rij = uij − xij ,
signifying the additional flow units that can be sent from vertex i to vertex j. It is also possible
to negate the existing flow xij on the arc by sending xij units of flow from j to i along the
arc (i, j). Notably, sending a unit of flow from i to j along the arc (i, j) increases the objective
function by cij , while sending a unit of flow from j to i on the same arc decreases the flow cost
by cij . Building upon these concepts, the residual network concerning the given flow x is defined
as follows. Consider a transportation network G = (V, E). A feasible solution x results in a new
(residual) transportation network denoted as Gx = (V, Ex ), where Ex represents a set of residual
edges corresponding to the feasible solution x. So, what is Ex ? It involves replacing each arc
(i, j) in E with two arcs: (i, j) and (j, i). The arc (i, j) has a cost of cij and a (residual) capacity
of rij = uij − xij , while the arc (j, i) has a cost of −cij and a (residual) capacity of rji = xij .
The set Ex is then constructed from the new edges with a positive residual capacity. Refer to
the figure below to ensure a clear understanding of the construction of the residual network.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 22/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

(b) The residual network with respect to the found

(a) A feasible solution feasible solution

You can notice immediately that such a definition of the residual network has some technical
difficulties. Let’s sum them up:

• If G contains both the edges (i, j) and (j, i), the residual network may contain four edges
between i and j (two parallel arcs from i to j and two contrary). To avoid this situation
we have two options. First, transform the original network to one in which the network
contains either edge (i, j) or edge (j, i), but not both, by splitting the vertexes i and j.
Second, represent our network by the adjacency list, which is handling parallel arcs. We
could even use two adjacency matrices if it were more convenient.

• Let’s imagine now that we have a lot of parallel edges from i to j with different costs.
Unfortunately, we can’t merge them by summarizing their capacities, as we could do while
we were finding the maximum flow. So, we need to keep each of the parallel edges in our
data structure separate.

3.4.3 Successive Shortest Path Algorithm

3.4.3.a Algorithm Overview

This algorithm can be generalized by three steps:

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 23/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

• Step 1: Select x as a feasible flow between any origin-destination pair. This flow has
the minimum cost among the set of feasible flows with the same flow value |x|.

• Step 2: Terminate the algorithm if |x| = v or if there is no shortest path in the

residual network (V, A(x), C(x), U (x), D). Otherwise, find the shortest path P with
the maximum flow value using a label-correcting algorithm and proceed to Step 3.
The functions A(x), C(x), U (x) are defined as follows:

A(x) = {(i, j) | (i, j) ∈ A, xij < uij } ∪ {(j, i) | (j, i) ∈ A, xij > 0}


c ,
ij (i, j) ∈ A, xij < uij
C(x) =
−c , (j, i) ∈ A, x > 0
ij ji


u , (i, j) ∈ A, x < u
ij ij ij
U (x) =
x , (j, i) ∈ A, x > 0
ji ji

• Step 3: Increase the flow on each edge of P up to the capacity limit of the edge with
the smallest limit on the path, ensuring that |x| ≤ v. If the flow value |x| does not
exceed v, return to Step 1.

3.4.3.b Application of Successive Shortest Path Algorithm

Let us define:

In which: xij represents flow unit, uij represents capacity, and cij represents the cost of path
(i, j).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 24/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

First we add 2 additional nodes s, t, which are the start and end node, and then calculate
the minimum cost by adding 1 flow.

We have the yellow path representing the path with minimum cost, which is 0 + 6 + 3 + 0 = 9.
Increase the flow of each edge up to the capacity limit of the edge with the smallest limit on
the path (min(10, 10, 5) = 5 flows).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 25/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

The current flow is 5, so the cost is: 5.0 + 5.6 + 5.3 + 5.0 = 45.
We proceed to find the Residual Network.

Using Bellman-Ford algorithm to check for negative cycles, we couldn’t find any of them.
Continue to find the path with minimum cost, we have the yellow path, which has the cost:
0 + 6 + 6 + 4 + 0 = 16.
After increasing the flow on each edge up to the capacity limit, we have the following Residual
Network:

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 26/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

At this time, there is no path from s to t so we end the algorithm here. Rebuild the trans-
portation model:

The optimal flow has the cost of: 10.6 + 5.6 + 5.3 + 5.4 = 125.

3.4.4 The Cycle Canceling Algorithm

3.4.4.a Algorithm Overview

This section describes the negative cycle optimal conditions and, as a consequence, cycle-
canceling algorithm.
We are starting with this important theorem:
Theorem 1 (Solution Existence): Let G be a transportation network. Suppose that G
contains no uncapacitated negative cost cycle and there exists a feasible solution of the minimum
cost flow problem. Then the optimal solution exists.
Proof. One can see that the minimum cost flow problem is a special case of the linear pro-
gramming problem. The latter is well known to have an optimal solution if it has a feasible
solution and its objective function is bounded. Evidently, if G doesn’t contain an uncapacitated
negative cycle then the objective function of the minimum cost flow problem is bounded from
below — therefore, the assertion of the theorem follows forthwith.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 27/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Theorem 2 (Negative Cycle Optimality Conditions). Let x∗ be a feasible solution of

a minimum cost flow problem. Then x∗ is an optimal solution if and only if the residual network
Gx∗ contains no negative cost (directed) cycle.
Proof:

• If the solution x∗ is optimal, this means there is no way to improve it further.

• If the surplus network Gx∗ contains a negative cost cycle, you can “go around” this cycle
to reduce the cost of the flow. Therefore, x∗ cannot be optimal. Conversely, if the surplus
network Gx∗ does not contain a negative cost cycle, it means there is no way to “go around”
and reduce the cost of the flow. Therefore, x∗ is optimal. In summary, this relates to how
the cost of the flow changes when we change the solution, and if there is no negative cost
cycle, then the current solution is optimal.

In summary, this relates to how the cost of the flow changes when we change the solution,
and if there is no negative cost cycle, then the current solution is optimal.
This theorem gives the cycle-canceling algorithm for solving the minimum cost flow problem.
First, we use any maximum flow algorithm to establish a feasible flow in the network. Then
the algorithm attempts to improve the objective function by finding negative cost cycles in the
residual network and augmenting the flow on these cycles. Let us specify a program in pseudocode
like it is done:

1. Establish a feasible flow x in the network

2. while (Gx contains a negative cycle) do

3. Identify a negative cycle W

4. δ ← minrij : (i, j) ∈ W

5. Augment δ units of flow along the cycle W

6. Update Gx

How many iterations does the algorithm perform? First, note that all the data is integral. After
line 1 of the program we have an integral feasible solution x. It implies the integrality of Gx .
In each iteration of the cycle in line 2 the algorithm finds the minimum residual capacity in the
found negative cycle. In the first iteration δ will be an integer. Therefore, the modified residual
capacities will be integers, too. And in all subsequent iterations the residual capacities will be
integers again. This reasoning implies:
Theorem 3 (Integrality Property): If all edge capacities and supplies/demands on ver-
texes are integers, then the minimum cost flow problem always has an integer solution. The
cycle-canceling algorithm works in cases when the minimum cost flow problem has an optimal
solution and all the data is integral and we don’t need any other assumptions.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 28/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Now let us denote the maximum capacity of an arc by U and its maximum absolute value of
cost by C. Suppose that m denotes the number of edges in G and n denotes the number of ver-
texes. For a minimum cost flow problem, the absolute value of the objective function is bounded
by mCU . Any cycle canceling decreases the objective function by a strictly positive amount.
Since we are assuming that all data is integral, the algorithm terminates within O(mCU ) itera-
tions. One can use O(nm) algorithm for identifying a negative cycle (for instance, Bellman-Ford’s
algorithm or label correcting algorithm), and obtain complexity O(nm2 CU ) of the algorithm.
In conclusion, The Cycle Canceling algorithm is a well-established method used to solve the
minimum-cost flow problem. The primary objective is to optimize the flow in a network by
iteratively identifying and eliminating negative cycles, adjusting the flow values to improve the
overall cost.

3.4.4.b Application of Cycle Canceling Algorithm

In contrast to the example of Successive Shortest Path Algorithm, which has 1 supply and 1
demand, at this example, we have 2 supplies and 2 demands.
Let us define:

In which: xij represents flow unit, uij represents capacity, and cij represents the cost of path
(i, j).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 29/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

A feasible solution has the cost of: 6.8 + 6.5 + 4.6 + 5.4 = 122. We proceed to build the
following Residual Network:

We have the cycle 1 → 2 → 3 → 4 → 1 is a negative cycle, which has the cost of: 6 + 6 +
(−5) + (−8) = −1 and the capacity of it is: min(6, 1, 6, 6) = 1.
Proceed to augmentation along the cycle (1 cycle).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 30/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

The cycle 3 → 5 → 4 → 3 is a negative cycle that has the cost of 3 − 4 − 6 = −7 and the
capacity is min(5, 5, 5) = 5. Continue to augmentation along the cycle (5 cycles).

Cycle 1 → 3 → 4 → 2 → 1 is a negative cycle that has the cost of 6 + 6 − 5 − 8 = −1 and the

capacity of it is min(5, 5, 5, 5) = 5. Continue to augmentation along the cycle (5 cycles).

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 31/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

Using the Bellman-Ford algorithm to check, we find that there is no negative cycles left, we
can end the algorithm here. Rebuild the transportation network.

The optimal flow has the cost of: 6.6 + 5.6 + 5.3 = 81.

3.5 Algorithm Efficiency

• Successive Shortest Path Algorithm:

– Time Complexity: O(V E 2 log(U )), where V represents the number of vertices, E is
the number of edges, and U is the maximum capacity in the network.
– Pros: Easy to comprehend, simple to implement, effective for solving most straight-
forward graph problems.
– Cons: Not particularly efficient for complex graphs, relies heavily on path selection to
increase flow, potentially yielding inaccurate results.

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 32/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

• Cycle Canceling Path Algorithm:

– Time Complexity: O(V E 2 U ), with V as the number of vertices, E as the number of

edges, and U as the maximum capacity in the network.
– Pros: Particularly effective when dealing with integer capacities.
– Cons: The time complexity might be elevated in certain scenarios compared to alter-
native algorithms.

In conclustion: Depending on specific cases, we may choose different algorithms to solve

problems. In this context, the Successive Shortest Path algorithm remains a stable choice in all
cases, especially for the problem we are exploring, which is finding the widest evacuation route.

4 Conclusion - Appreciation
For the Question 1 of this year’s Assignment, the optimal cost is a negative number, meaning
we are making a profit. Calculating the profit from each product component gives us a production
plan that yields the best profit. Therefore, the 2-SLPWR problem is very useful in various fields
of life, or is this case, both Problem 1 and the Emergency Evacuation Problem - Problem 2.
Regarding Question 2, based on the complexity of each algorithm, we can see that the
Successive Shortest Path algorithm has the highest efficiency among the two. Thus, for the
EPDR problem, the Successive Shortest Path algorithm is the most suitable one.
In our submission, we have attached: (1) this report, (2) a python file that contain our solution
for Problem 1, (3) a C++ file containing source code for 2 Min-cost flow algorithms included in
this report to further illustrate the algorithms, and (4) a weekly progress log.
We can’t complete this report without the support of the teachers in charge of this course. We
want to express our sincere appreciation to the teachers, specifically to Mr. Thịnh, for resolving
our questions in class.

5 Appendix
Below is the group’s source code for Problem 1.
1 from gamspy import Container , Set , Parameter , Variable , Equation , Model , Sum ,
Sense
2 import numpy as np
3 import random
4 import pandas as pd
5

6 n = 8
7 m = 5
8 S = 4
9 p_s = 1/ S

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 33/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

10 model = Container ()
11

12 Products = [ f " Product { i +1} " for i in range ( n ) ]

13 Suppliers = [ f " Supplier { j +1} " for j in range ( m ) ]
14 Scenarios = [ f ’ Scenario { i +1} ’ for i in range ( S ) ]
15

16 i = Set ( container = model , name = " i " , records = Products )

17 j = Set ( container = model , name = " j " , records = Suppliers )
18 k = Set ( container = model , name = " k " , records = Scenarios )
19

20 dataframe_d = pd . DataFrame (
21 [[ sce , pro , np . random . binomial (10 ,0.5 , size =1) [0]] for pro in Products for sce
in Scenarios ] , columns =[ " Scenario " , " Product " , " Demand " ]
22 ) . set_index ([ " Scenario " , " Product " ])
23 d = Parameter ( container = model , name = " d " , domain =[ k , i ] , records = dataframe_d .
reset_index () )
24

25 dataframe_A = pd . DataFrame (
26 [[ pro , sup , random . randint (0 ,10) ] for sup in Suppliers for pro in Products ] ,
columns =[ " Product " , " Supplier " , " Units " ]
27 ) . set_index ([ " Product " , " Supplier " ])
28 A = Parameter ( container = model , name = " A " , domain =[ i , j ] , records = dataframe_A .
reset_index () )
29

30 dataframe_b = pd . DataFrame (
31 [[ sup , random . randint (1000 ,10000) ] for sup in Suppliers ] , columns =[ " Supplier " ,
" value " ]
32 ) . set_index ( " Supplier " )
33 b = Parameter ( container = model , name = " b " , domain =j , records = dataframe_b . reset_index
() )
34

35 dataframe_s = pd . DataFrame (
36 [[ sup , random . randint (1000 , dataframe_b . loc [ sup , ’ value ’ ]) ] for sup in
Suppliers ] ,
37 columns =[ " Supplier " , " Cost " ]
38 ) . set_index ( " Supplier " )
39 s = Parameter ( container = model , name = " s " , domain =j , records = dataframe_s . reset_index
() )
40

41 dataframe_l = pd . DataFrame (
42 [[ pro , random . randint (100 _000 , 1 _000_000 ) ] for pro in Products ] ,
43 columns =[ " Product " , " value " ]
44 ) . set_index ( " Product " )
45

46 l = Parameter ( container = model , name = " l " , domain =i , records = dataframe_l . reset_index
() )
47

48 dataframe_q = pd . DataFrame (
49 [[ pro , random . randint (4 _000_000 , 10 _000_000 ) ] for pro in Products ] ,

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 34/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

50 columns =[ " product " , " value " ]

51 ) . set_index ( " product " )
52 q = Parameter ( container = model , name = " q " , domain =i , records = dataframe_q . reset_index
() )
53

54 cost = Parameter ( container = model , name = " cost " , domain = i )

55 cost [ i ] = l [ i ] - q [ i ]
56

57 x = Variable (
58 container = model ,
59 name = " x " ,
60 domain =j ,
61 type = " Positive " ,
62 )
63

64 y = Variable (
65 container = model ,
66 name = " y " ,
67 domain =[ k , j ] ,
68 type = " Positive " ,
69 )
70

71 z = Variable (
72 container = model ,
73 name = " z " ,
74 domain =[ k , i ] ,
75 type = " Positive " ,
76 )
77

78 supply = Equation (
79 container = model ,
80 name = " supply " ,
81 domain =[ k , j ]
82 )
83

84 supply [k , j ] = y [k , j ] == x [ j ] - Sum (i , A [i , j ] * z [k , i ])
85

86 demand = Equation (
87 container = model ,
88 name = " demand " ,
89 domain =[ k , i ]
90 )
91

92 demand [k , i ] = z [k , i ] <= d [k , i ]
93

94 obj = Sum (j , b [ j ] * x [ j ]) + p_s * ( Sum (( k , i ) , cost [ i ] * z [k , i ]) - Sum (( k , j ) , s [ j ]

* y [k , j ]) )
95

96 prob = Model (

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 35/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

97 container = model ,
98 name = " prob " ,
99 equations = model . getEquations () ,
100 problem = " MIP " ,
101 sense = Sense . MIN ,
102 objective = obj ,
103 )
104

105 prob . solve ()

106

107 print ( prob . ob j ec ti v e_ va l ue )

108

109 print ( x . records )

110

111 print ( y . records )

112

113 print ( z . records )

Listing 1: Problem 1’s Source code

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 36/37
 University of Technology, Ho Chi Minh City
Faculty of Computer Science and Engineering

References
[1] [Link], “A two-stage stochastic programming framework for evacuation planning in disaster
responses”, Computers & Industrial Engineering, vol. 145, p. 106458, 2020.

[2] Alexander Shapiro and Andy Philpott, “A Tutorial on Stochastic Programming ”, 2007.

[3] Hashnayne Ahmed, “Formulation of Two-Stage Stochastic Programming with Fixed Re-
course”.

[4] [Link]

[5] [Link]

Assignment for Mathematical Modelling - Academic year 2023 - 2024 Page 37/37