NISTIR 7907

Flow Control in Time-Varying, 

Random Supply Chains

Ion Matei

Assane Gueye

John S. Baras

http://dx.doi.org/10.6028/NIST.IR.7907
NISTIR 7907

Flow Control in Time-Varying, 

Random Supply Chains

Ion Matei

Manufacturing Systems Integration Division

Engineering Laboratory

Assane Gueye

Applied and Computational Mathematics Division

Information Technology Laboratory

John S. Baras

University of Maryland

http://dx.doi.org/10.6028/NIST.IR.7907

January 2013

U.S. Department of Commerce 

Rebecca Blank, Acting Secretary

National Institute of Standards and Technology 

Patrick D. Gallagher, Under Secretary of Commerce for Standards and Technology and Director 
1

Flow Control in Time-Varying, Random Supply

Chains

Ion Matei, Assane Gueye and John S. Baras

Abstract

This paper focuses on the logistics aspect of supply chain management. It proposes a randomized flow management

algorithm for a time-varying, random, supply chain network. A constrained stochastic optimization problem that

maximizes the profit function in terms of the long-run, time-average of the flows in the supply chain is formulated.

The algorithm is distributed and based on queueing theory and stochastic Lyapunov analysis concepts. The long-run,

time average of the flows generated by the algorithm can get arbitrarily close to the solution of the aforementioned

optimization problem. In support of the theoretical results, numerical simulations are also presented.

I. Introduction

Among many possible definitions, the supply chain can be defined as a network of interrelated activities of

procurement, production, distribution, vendition, and consumption of one of more products [19]. Manufacturing is

often outsourced around the world, with each component made in locations chosen for their expertise and low costs

[17]. Consequently, today’s supply chains are increasingly complex and rely on critical infrastructures such as roads,

railways, and airports to move goods [16], and therefore they exhibit the co-existence of operational optimization

with operational vulnerability [17]. This was most recently and dramatically demonstrated in the aftermath of

several accidents and natural disasters. For example, a fire in the Phillips Semiconductor plant in Albuquerque,

New Mexico caused its major customer, Ericsson, to lose $400 million in potential revenues. Another example

concerns the impact of Hurricane Katrina. This storm halted 10% - 15% of the total U.S. gasoline production,

raising both domestic and overseas oil prices [4]. More recently, the tragic earthquake of March 13, 2011, off the

northeastern coast of Japan and the devastating tsunami that followed have shattered the nation, with immense

loss of life and property. In addition, it brought uncertainty of the future, not the least of which is the expected

decades-long impact of the nuclear reactors in Fukushima [17].

This material is based in part upon work supported by the NIST-ARRA Measurement Science and Engineering Fellowship Program award

70NANB10H026, through the University of Maryland, and in part upon work supported by the National Institute of Standards and Technology

(NIST) grant award 70NANB11H148.

Ion Matei and Assane Gueye are with the National Institute of Standards and Technology, Gaithersburg, MD 20899 (ion.matei@nist.gov,

assane.gueye@nist.gov).

John S. Baras is with Electrical Engineering Department and the Institute for Systems Research at the University of Maryland, College Park,

MD 20742 (baras@umd.edu).
2

As the world’s economies become increasingly interconnected into a global economy, supply chain networks face

many new types of risk, including natural disasters, political/social instability, cultural/communication inconsistency,

exchange rate fluctuation, and local legislations [2]. These risks forced the supply chains’ stakeholders to go beyond

the operational optimization and to recognize the operational vulnerabilities of the supply chains and to underline

their time-varying and random nature.

This paper focuses on the logistics aspect of the supply chain management. Logistics plans, implements, and

controls efficient and effective product storage and flows (forward/reverse). Logistics starts from the point of origin

to the point of consumption, with the goal of meeting customer requirements [3]. The paper addresses the flow

management in a supply chain that exhibit stochastic behavior in both links and demands, and in addition it

responds the need for decentralized decisions as point out in [5]. A randomized and decentralized algorithm for the

management of the flow of the product in a time-varying, random supply chain aimed at maximizing the profit of

a firm is proposed. Due to the random nature of the supply chain, the profit function is defined to be dependent on

the (long-run, time) averages of the flows, since the flows are random processes. Hence, the optimization problem

becomes stochastic. The approach for solving the optimization problem is as follows. First, the satisfiability of the

supply chain’ constraints is transformed into a stability condition on a set of queues associated with the supply

chain’s components. Second, a Lyapunov drift analysis technique is used to generate an algorithm that ensures

the stability of the queues, and at the same time maximizes the profit function. This approach avoids the need of

a realization of the stochastic parameters, as it is the case in a stochastic approximation approach. At each time

instant, the algorithm produces decisions on the flows that are implementable (that is, take into account the current

the state of the supply chain). More importantly, the resulting long-run, time averages of the flows get arbitrarily

close to the solution of the stochastic optimization problem. In addition, the algorithm does not require knowledge

of the probability distribution of the random process that drives the supply chain and deals with both supply changes

and demand variability. Furthermore, the actions taken by a specific decision maker are based only on a localized

view of the state of the supply chain. This localized view consists of the state of all the links that have at one end

the decision maker. In other words, the algorithm is distributed.

The operational research literature emphasizes importance of flows’ management in supply chains, with approaches varying from linear, non-linear or mixed-integer programming [1],[18] to game theory [11], [12],[13],[19].

The role of a supply chain, the key strategic drivers of its performance, and the analytical methodologies for its

analysis are extensively treated in [5].

The study of a supply chain under a stochastic setup have been addressed in several works in the literature, however

there are some significant differences compared to the approach presented in this paper. In [10], the authors propose

an algorithm for determining the system reliability with respect to the maximum flow of a network achieving a given

demand. Although the network studied by the authors have nodes that can fail randomly, the demand is assumed

deterministic and no cost/profit functions are considered in their analysis. In [6], the goal is to determine how much

of a particular product a plant should produce, given a (possible random) demand and based on maximizing a utility

function. The authors use a simplified model for a supply chain, formed by plants and retailers only, the resulting
3

network topology being a deterministic bipartite graph. The authors propose a heuristic scheme for determining the

assignment policy and focus most of their attention to a particular type of graph, called expander graphs. Expanders

graphs are interesting due to there spectral properties, that is, they do not degrade by increasing the number of

nodes. Compared to this work, although the profit function can be interpreted as a utility function, in the current

paper the graph is arbitrary and stochastic and the proposed algorithm is based on a rigorous, mathematical analysis.

Another example of supply chain analysis under random demands is introduced in [7]. Similar to the current paper,

the authors focus on determining the flow on the supply chain links based on optimizing an objective function, but

the supply chain is assumed deterministic. Another formulation for the analysis of a supply chain under a stochastic

setup is presented in [15], where the authors consider the processing/transportation costs, demands, supplies, and

capacities to be stochastic parameters. The goal is to minimize the expectation of a cost function and the authors

chose a stochastic approximation strategy to solve the optimization problem. This approach consists of using a

realization of the stochastic parameters to approximate the expectation cost and then use deterministic optimization

techniques to solve the resulting problem. The main disadvantage of this approach is that the accuracy of the solution

depends on the number on samples the joint probability distribution of the stochastic parameters must be knows.

In the current paper, the approach for solving the stochastic optimization problem is not based on a approximation

of the expected cost and there is no need for the probability distribution to be known.

The paper is organized as following. Section II introduces the model for the time-varying supply chain network

considered in this paper. Section III introduces the notion of the capacity region of a supply chain and formulates a

constrained stochastic optimization problem, aimed at maximizing the profit function in terms of the long-run timeaverage of the flows. Section IV describes a randomized, dynamic flow control algorithm for solving the stochastic

optimization problem, using queuing theory concepts to model the constraints. Section V presents a performance

analysis of the flow control algorithm, which shows that the solution of the algorithm can get arbitrarily close to

the solution of the optimization problem described in Section III. The paper ends with numerical simulations of

the proposed algorithm (Section VI) and some concluding remarks (Section VIII).

II. Supply chain model

A firm involved in the production, storage and distribution of a homogeneous product is considered. The firm

uses a set of manufacturing facilities, a set of warehouses and serves a set of retail outlets/demand markets.

The supply chain model used in this paper is similar to the one used in [12], with the main difference that

the network is time-varying and random. An example of a supply chain network is given in Figure 1, where

node 1 represents the firm, nodes {2,3,4} represent the set of manufacturing facilities, nodes {(5,5

′

),(6,6

′

)} are the

warehouses and nodes {7,8,9} designate the retail outlets/demand markets.

A supply chain with only one firm is considered. The single-firm scenario is suitable for a dominant-firm model,

where a single firm controls a dominant share of the market [14]. The sets of firms, manufacturers, warehouses and

retailers are denoted by F , M, W and R, respectively. In addition, let N be the set of all nodes in the network

(with a typical node denoted by i), i.e., N = {F ∪M ∪W ∪ R} ∪ {i

′

|i ∈ W}, with cardinality N = |N|. Note that
4

Fig. 1: Example of supply chain network

similarly to [12], a warehouse i is represented by two nodes in the network (by using i

′

as well) in order to clearly

emphasize the flow of the product passing through the warehouse, i.e., through the link (i,i

′

). The set of links of

the supply chain is denoted by L = {(i, j),i , j ∈ N}, where products “flow” from node i to node j for each (i, j) ∈ L

and where the flow of the product in the chain is driven by the demand at the retailers/markets. It is assumed that

links of the form (i,i

′

) are also included in L.

The supply chain operates in slotted time, with slots normalized to integral units so that slot times occur at

times t ∈ {0,1,2,...}. The state of the supply chain at time t is denoted by S (t). The state process S (t) incorporates

the stochastic/nondeterministic behavior of the supply chain, such as possible disruptions in manufacturing and

transportation due to natural disasters, power outages, technical and malfunctions. For example, the transport or

manufacturing capacity can be at full capacity or at zero capacity in case of uncontrollable events. For simplicity,

throughout the rest of the paper, we assume that the links of the supply chain can be either active or inactive, as

described by S (t). This means that a transportation link may become unavailable at some time slot. The following

assumptions about the statistical properties of S (t) are made.

Assumption 2.1: The process S (t) belongs to a finite set S and evolves according to an identically, independently
5

distributed random process, with stationary distribution given by π = (πs), where

πs = lim

t→∞

1

t

Xt−1

τ=0

1{S (τ)=s}, ∀s ∈ S, (1)

with 1{S (τ)=s} being the indicator function that takes value one whenever S (t) = s, and zero otherwise.

The amount of product flowing through the link (i, j) during time slot t is denoted by μi, j(t) . Without loss of

generality it is assumed that the flows are measured in (final) product units; to recover other units (raw materials

for example) the flows are multiplied by the process rate of the economic unit generating the flow. The random

process di(t) for i ∈ R represents the demand at market i. It is reasonable to assume that the quantity of product

flowing between different entities is upper-bounded, and hence the following assumption is made.

Assumption 2.2: The flows μi, j(t) are non-negative for all time-slots t and there exist positive scalars μ

max

i

such

that

X

b

μi,b(t) ≤ μ

max

i

, ∀i ∈ N, ∀t, (2)

where all pairs (i,b) belong to the set L.

The above inequalities limit the total flow of the product leaving any node, which can be thought of as production,

transportation or storage capabilities limitations.

The following definitions introduce the time averages of the product flows in the supply chain.

Definition 2.1: The time average flows of product in the supply chain are given by

μ ̄i, j(t) =

1

t

Xt−1

τ=0

E{μi, j(τ)}, (3)

and the long-run time averages of flow product are given by

μ ̄i, j = lim

t→∞

μ ̄i, j(t), (4)

for all (i, j) ∈ L.

Additionally, the market demands satisfy the following assumption.

Assumption 2.3: The random processes di(t) are independent and identically distributed, with mean given by

d ̄

i = E{di(t)}, ∀i ∈ R. (5)

The aggregate vectors of product flows and market demands are denoted by μ(t) =



μi, j(t)



for (i, j) ∈ L, and

d(t) = (di(t)) for i ∈ R, respectively.

III. Formulation of the stochastic optimization problem

In this section, the optimization problem the firm needs to solve to maximize its profit is presented. The profit

function is defined as the difference between the revenue from selling the product and the cost for producing the

product. Since the supply chain is assumed random, the profit function is defined in terms of the long-run, time

averages of the product flows. The flows of the product must satisfy a set of constraints induced by the supply
6

chain network. These constraints define the capacity region of a supply chain, which tells how much demand the

supply chain can support.

Definition 3.1: The capacity region Λ of a supply chain is the closure of all vector of demands x = (xi) that can

be supported by the supply chain network, considering all possible strategies for choosing the flows of product,

under the limitations introduced by Assumption 2.2.

In the following, a more detailed characterization of the capacity region of a supply chain is given. To that end,

let Ci, j(s) be the set of flows on link (i, j) satisfying Assumption 2.2, when the supply network is in state s, and

under all possible flow control policies. Let C(s) be the set of all link sets, i.e., C(s) =



Ci, j(s)



for (i, j) ∈ L. Let

co{C(s)} denote the convex hull of the set of all possible values of C(s). Recalling that the state of the supply chain

is an i.i.d. random process, the set of the average convex hull of all possible flows on links, given all possible states

can be defined. This average set can be formally written as a family of graphs Γ, given by

Γ ,

X

s∈S

πsco{C(s)}.

A matrix G = (Gi, j) is said to belong to Γ if there exits a randomized flow control policy that depends on the state

of the network, such that

G =

X

s∈S

πsE{μ(t)|S (t) = s},

where E{μ(t)|S (t) = s} is the expected flow matrix under the considered policy, given that the supply chain is in

state s.

The following Theorem inspired by [8] gives a mathematical characterization of the capacity region of the supply

chain.

Theorem 3.1: The capacity region of a supply chain is given by the set Λ of all demand vectors x = (xi) such

that there exits a flow matrix G = (Gi, j) belonging to the closure of Γ, together with flow variables fi, j such that

fi, j ≥ 0, ∀(i, j) ∈ L, fi, j = 0, ∀(i, j) < L, (6)

X

a∈F

fa,i =

X

b∈W

fi,b, ∀i ∈ M, (7)

X

a∈M

fa,i = fi,i

′ , ∀i ∈ W, (8)

fi,i

′ =

X

b∈R

fi

′

,b, ∀i ∈ W, (9)

X

a∈W

fa

′

,i = xi

, ∀i ∈ R, (10)

fi, j ≤ Gi, j, ∀(i, j) ∈ L. (11)

In the particular case where the process S (t) is i.i.d. (which in fact is the assumption throughout this paper), the

next Corollary presents a further characterization of the capacity region, where Cl(A) is used to denote the closure

of the set A.
7

Corollary 3.1 (adaptation of Corollary 3.9, [8]): If Γ is a closed set and if the state process S (t) is i.i.d. from

slot to slot, the demand vector x is within the capacity region Λ if and only if there exists a stationary (randomized)

policy that chooses μ(t) based only on the current topology state S (t), such that

E







X

a∈F

μa,i(t)







= E







X

b∈W

μi,b(t)







,∀ i ∈ M,

E







X

a∈M

μa,i(t)







= E



μi,i

′ (t)

	

,∀ i ∈ W,

E



μi,i

′ (t)

	

= E







X

b∈R

μi

′

,b(t)







,∀ i ∈ W,

E







X

a∈W

μa

′

,i(t)







= xi,∀ i ∈ R,

where the expectation is taken with respect to the random process S (t) and the (potentially) random policy based

on S (t).

Note that if x ∈ Λ, then any  ̃x such that  ̃x ≤ x entrywise, also belongs to Λ. In addition, it can be shown that the

set Λ is convex, closed and bounded and it contains the vector of all zeros, (i.e., 0 ∈ Λ).

The previous Corollary gives the constraints induced by the supply chain network that the flows of product must

satisfy. Next, a stochastic optimization problem is formulated; problem that describes the objective of the firm under

the network constraints introduced above.

The goal of the firm is to maximize its profit, that is the difference between the revenue and the cost functions.

The revenue function of the firm depends on the quantity of products that reach the retailers/markets in the long-run.

The revenue function is denoted by

f(μ ̄) =

X

i∈W, j∈R

fi

′

, j( ̄μi

′

, j),

where (i

′

, j) represent valid warehouse-retailer pairs, i.e., i ∈W, j ∈ R and (i

′

, j) ∈ L. Cost functions associated with

each link (i, j) ∈ L are also considered, and are denoted by gi, j( ̄μi, j). These cost functions depend on the flow of the

product on the links and are generated by activities such as acquiring raw materials, manufacturing, transportation

or warehouse usage. The total cost function is given by

g(μ ̄) =

X

i∈F , j∈M

gi, j( ̄μi, j)+

X

i∈M, j∈W

gi, j( ̄μi, j)+

X

i∈W

gi,i

′ ( ̄μi,i

′ )+

+

X

i∈W, j∈R

gi

′

, j( ̄μi

′

, j).

Assumption 3.1: The functions fi, j are non-negative, continuously differentiable and concave, while the functions

gi, j are non-negative, continuously differentiable and convex.

The profit function h is the difference between the revenue and the cost functions, i.e.,

h(μ ̄) = f(μ ̄)− g(μ ̄).
8

The firm’s objective is to maximize the profit under the flow constraints induced by the (capacity region of the)

supply chain network. Let xi denote the long-run, average flow of product arriving at market (retailer) i, that is,

xi =

X

a∈W

μ ̄a,i, ∀i ∈ R.

The following stochastic optimization problem is considered:

max

μ ̄,x

h(μ ̄) (12)

subject to: x ∈ Λ,

x ≤ d ̄.

The first constraint introduced above ensures that the average product flows arriving at the markets (retailers) are

within the capacity region of the supply chain network, i.e., can be supported by the network. The second inequality

ensures that the long term flow of the product arriving at the markets are no larger than the demands at the markets.

By Corollary 3.1, the above stochastic optimization problem can be equivalently represented as

max

μ ̄

h(μ ̄) (13)

subject to: P

a∈F μ ̄a,i =

P

b∈W μ ̄i,b,∀i ∈ M,

P

a∈M μ ̄a,i = μ ̄i,i

′ ,∀i ∈ W,

μ ̄i

′

,i =

P

b∈R μ ̄i

′

,b,∀i ∈ W,

P

a∈W μ ̄a

′

,i ≤ d ̄i

,∀i ∈ R,

where  ̄μi, j = E{μi, j(t)} for all (i, j) ∈ L, with μi, j(t) being chosen by some stationary, randomized control algorithm,

based only on the current state S (t).

Assumption 3.2 (Interior point): There exist positive scalars ǫ1 and ǫ2 and two stationary randomized flow control

policies based on the current state S (t), corresponding to ǫ1 and ǫ2, respectively, such that

E{μ

ǫ1

1,i

(t)}+ǫ1 = E







X

b

μ

ǫ1

i,b

(t)







,∀ i ∈ M,

X

a

E

n

μ

ǫ1

a,i

(t)

o

+ǫ1 = E

n

μ

ǫ1

i,i

′ (t)

o

,∀ i ∈ W,

E

n

μ

ǫ1

i,i

′ (t)

o

+ǫ1 = E







X

b

μ

ǫ1

i

′

,b

(t)







,∀ i ∈ W,

E







X

a

μ

ǫ1

a

′

,i

(t)







+ǫ1 ≤ d ̄

i

,∀ i ∈ R,

and

E







X

b

μ

ǫ2

i,b

(t)







+ǫ2 = E{μ

ǫ2

1,i

(t)},∀ i ∈ M,

E

n

μ

ǫ2

i,i

′ (t)

o

+ǫ2 =

X

a

E

n

μ

ǫ2

a,i

(t)

o

,∀ i ∈ W,
9

E







X

b

μ

ǫ2

i

′

,b

(t)







+ǫ2 = E

n

μ

ǫ2

i,i

′ (t)

o

,∀ i ∈ W,

E







X

a

μ

ǫ2

a

′

,i

(t)







≤ d ̄

i

,∀ i ∈ R,

The above Assumption basically states that the optimal solution of (13) is not on the boundary of the capacity

region. In particular, ǫ1 can be viewed as an additional flow on one of the links that arrives at a node and is produced

by a source outside the supply chain, while ǫ2 can be viewed as an additional flow leaving a node on one of the

links but that fails to reach the destination node.

From the numerical optimization point of view, it is more advantageous to work with inequality constraints

rather than equality constraints. As a consequence, each equality constraint in (13) is replaced by two inequality

constraints, as shown in the following:

max

μ ̄

h(μ ̄) (14)

subject to: P

a∈F μ ̄a,i ≤

P

b∈W μ ̄i,b,∀i ∈ M,

P

a∈F μ ̄a,i ≥

P

b∈W μ ̄i,b,∀i ∈ M,

P

a∈M μ ̄a,i ≤ μ ̄i,i

′ ,∀i ∈ W,

P

a∈M μ ̄a,i ≥ μ ̄i,i

′ ,∀i ∈ W,

μ ̄i

′

,i ≤

P

b∈R μ ̄i

′

,b,∀i ∈ W,

μ ̄i

′

,i ≥

P

b∈R μ ̄i

′

,b,∀i ∈ W,

P

a∈W μ ̄a

′

,i ≤ d ̄i

,∀i ∈ R.

In the following sections a mathematical approach for solving the optimization problem (14) is introduced. This

approach is based on queueing theory and on drift analysis.

IV. Flow control algorithm

In this section a flow control algorithm which ensures that the long-run, time-average flows in the supply chain

get arbitrarily close to the optimal solution of (13) is presented. The main idea behind the algorithm is to associate

to each of the inequality constraints a (virtual) queue. As shown in what follows, the inequality constraints are

satisfied if the queues associated to them are stable, in some sense that is about to be defined. By taking advantage

of this property, an algorithm that stabilizes the queues and gets arbitrarily close to the optimal solution of (14) is

proposed. The algorithm is derived as a result of a drift analysis approach on the (virtual) queues. This approach

is closely related to the stochastic Lyapunov theory [9] and avoids using a realization of the stochastic parameters

for approximating the objective function.
10

A. Modeling inequality constraints using queues

This subsection shows why the feasibility of the inequality constraints defined in the optimization problem (14)

can be connected to the stability of a set of queues associated to them.

Consider a queue U(t) (Figure 2) with (possibly random) input λ(t) and output μ(t), whose dynamics is given by

U(t+1) = max{U(t)−μ(t),0}+λ(t).

Fig. 2: Queue schematics

Definition 4.1: The queue U(t) is said to be strongly stable if

lim sup

t→∞

1

t

Xt−1

τ=0

E{U(τ)} < ∞.

Let us now assume that there exists λ ̄ and  ̄μ such that

λ ̄ = lim

t→∞

1

t

Xt−1

τ=0

E[λ(τ)],

and

μ ̄ = lim

t→∞

1

t

Xt−1

τ=0

E[μ(τ)].

Proposition 4.1 (Queue stability): A necessary condition for the strong stability of the queue U(t) is

λ ̄ ≤ μ.  ̄

The necessary condition is quite intuitive. Indeed, if λ > ̄ μ ̄, then the expected queue backlog grows to infinity,

leading to instability. Under additional assumptions on the processes λ(t) and μ(t), it can be shown that λ < ̄ μ ̄ is

also a sufficient condition (see [8] for more details).

As previously mentioned, a set of (virtual) queues are associated to the constraints of the optimization problem

(14), whose dynamics are given in the following.
11

In the case of a manufacturing unit, the dynamics of the queue levels are given by

U

1

i

(t+1) = max







U

1

i

(t)−

X

b

μi,b(t),0







+

X

a

μa,i(t),∀i ∈ M, (15)

U

2

i

(t+1) = max







U

2

i

(t)−

X

a

μa,i(t),0







+

X

b

μi,b(t),∀i ∈ M. (16)

The queues associated to the wharehouses evolve in time according to

U

1

i

(t+1) = maxnU

1

i

(t)−μi,i

′(t),0

o

+

X

a

μa,i(t),∀i ∈ W, (17)

U

2

i

(t+1) = max







U

2

i

(t)−

X

a

μa,i(t),0







+μi,i

′ (t),∀i ∈ W. (18)

and

U

1

i

′ (t+1) = max







U

1

i

′(t)−

X

b

μi

′

,b(t),0







+μi,i

′ (t),∀i ∈ W, (19)

U

2

i

′ (t+1) = maxn

U

2

i

′ (t)−μi,i

′(t),0

o

+

X

b

μi

′

,b(t),∀i ∈ W. (20)

The dynamics of the queues corresponding to the retailers is given by

Ui(t+1) = maxnU

1

i

(t)−di(t),0

o

+

X

a

μa,i(t),∀i ∈ R. (21)

Remark 4.1: In the previous expressions, P

b μi,b(t) represents the summation over all active links carrying

products from node i, at time slot t, as per the state of the supply chain state S (t). A similar interpretation can be

given to the term P

a μa,i(t).

From Proposition 4.1 it can be inferred that any flow control algorithm stabilizing the queues produces a solution

that satisfies the flow constraints defined in the optimization problem (13). Therefore, it makes sense to look for

an algorithm that stabilizes the queues defined above and in the same time maximizes the profit function.

B. Algorithm description

This section introduces a randomized flow control algorithm that can get arbitrarily close to the optimal solution of

(13). The algorithm stabilizes the (virtual) queues and therefore ensures that the inequality constraints are satisfied,

but, most importantly, it shows how the economic entities in the supply chain dynamically adapt their flows based

on the changes in the network.

The algorithm consists of actions taken by the entities involved in the economic activities of the firm, at each

time slot t. Let δ be a positive scalar, that affects the performance of the algorithm. For simplicity, the set of firms

F contains only one firm, say node 1 in the network. In the following the flow control algorithm is described.

• Control of the raw material flow: At every time slot, the firm observes the current levels of the manufacturers’

queues, U

1

b

(t) and U

2

b

(t). Then, at each time t it chooses the amount μ1,b of raw material sent to manufacturer

b, where μ1,b is the solution of the following optimization problem:

min

μ1,b

P

b∈M 

δg1,b(μ1,b)+

h

U

1

b

(t)− U

2

b

(t)

i

μ1,b



(22)

subject to: P

b∈M μ1,b ≤ μ

max

1

, μ1,b ≥ 0,∀b. (23)
12

• Control of the flow of product from manufacturers to warehouses: At every time slot, each manufacturer i

observes the current level of its queues U

1

i

(t) and U

2

i

(t) and the current levels of the queues of the warehouse

b to which product is possible to be sent to (as per the state of S (t)), i.e., U

1

b

(t) and U

2

b

(t). The amount

of product sent to each warehouse b at time slot t is given by μi,b, obtained as solution of the following

optimization problem:

min

μi,b

P

b δgi,b(μi,b)−

hU

1

i

(t)− U

1

b

(t)

i

+

h

U

2

b

(t)− U

2

i

(t)

iμi,b (24)

subject to: P

b∈W μi,b ≤ μ

max

i

, μi,b ≥ 0,∀b, (25)

for all i ∈ M, b ∈ W and (i,b) ∈ L which are active at time t, as per the state of the supply chain given by

S (t).

• Control of the flow of product within the warehouses: At every time slot, each warehouse i observes the current

level of its queues U

1

i

(t), U

1

i

′ (t), U

2

i

(t) and U

2

i

′ (t). The amount of product allowed in the warehouse at time

slot t is given by μi,i

′ , obtained as solution of the following optimization problem:

min

μ

δgi,i

′(μ)−

hU

1

i

(t)− U

1

i

′(t)

i

+

h

U

2

i

′ (t)− U

2

i

(t)

iμ (26)

subject to: 0 ≤ μ ≤ μ

max

i

(27)

for all i ∈ W and (i,i

′

) ∈ L which are active at time t, as per the state of the supply chain given by S (t).

• Control of the flow of product from warehouses to retailers: At every time slot, each warehouse i observes

the current level of its queues backlog U

1

i

′(t) and U

2

i

′(t) and the current level of the queue of the retailer b to

which the product is sent to, i.e., U

1

b

(t). The amount of product sent to retailer b at time slot t is given by

μi

′

,b, where μi

′

,b are obtained as solution of the following optimization problem:

min

μi

′

,b

P

b∈R δgi

′

,b(μi

′

,b)−δ fi

′

,b(μi

′

,b)−

hU

1

i

′ (t)− U

1

b

(t)



− U

2

i

′(t)

i

μi

′

,b (28)

subject to: P

b∈R μi

′

,b ≤ μ

max

i

′

, μi

′

,b ≥ 0, ∀b, (29)

for all i ∈ W, b ∈ R and (i

′

,b) ∈ L which are active at time t, as per the state of the supply chain given by

S (t).

Note that the optimization problems (22)-(28) are convex constrained optimization problems, which can be solved

efficiently at each time slot. Also, note that each of the entities involved in the economic activities does not need to

know the entire state of the network, nor the probability distribution of S (t). Indeed, in the case of the manufacturers,

the raw material flow is determined only by the level of the queues’ backlogs and the cost. When a manufacturer

must decide the flow of the product sent to warehouses, it looks at the current valid links, and it makes the decision

based on the cost of utilizing the respective links, and based on the difference between the queues’ levels of the

manufacturer and warehouses. In the case of the amount of product allowed in a warehouse, the decision is based

on the cost of keeping the product in the warehouse and on the difference between the levels of the (virtual) queues.

Finally, the amount of product sent to retailers from a warehouse is based on the current available links, on the

(localized) profit obtained from sending products to a specific retailer and on the difference between the queues’
13

levels of the warehouse and retailers. This limited need of information for implementing the algorithm makes it

advantageous for controlling the flow of product in increasingly complex and globalized supply chains. Another

important observation is that the manufacturers, warehouses and retailers do not need to know the entire state of

the network at a time slot, nor the statistics of the state process S (t). They only need to observe the state of links

which connect them to their neighbors. In addition, the virtual queues U

1

i

(t) can find an analogy in reality. Indeed,

in the case of a manufacturer for example, the queue can be viewed as a deposit for the raw material waiting to

be processed.

V. Derivation of the algorithm and performance analysis

This section shows the considerations behind the development of the algorithm and analyzes its performance. The

algorithm is derived as a result of a tradeoff between maximizing the profit function and maintaining the stability

of the queues introduced above. Stability of the queues ensures that the constraints introduced by the supply chain

are satisfied. By putting more weight on maximizing the profit function, the flows generated by the algorithm get

closer to the optimal solution. However, the backlogs of the queues are increased as well.

A. Derivation of the algorithm

The algorithm is derived as a result of a tradeoff between a drift function and the profit function. The drift is a

measure of the increase in the queues’ backlogs.

Let U(t) =



U

j

i

(t),i ∈ M,U

j

i

(t),U

j

i

′ (t),i ∈ W, j ∈ {1,2},Ui(t),i ∈ R

be the vector of queues. Using the quadratic

Lyapunov function

V(U(t)) ,

1

2

X

j∈{1,2}





X

i∈M

U

j

i

(t)

2 +

X

i∈W



U

j

i

(t)

2 + U

j

i

′(t)

2







+

1

2

X

i∈R

Ui(t)

2

.

the queues’ drift is given by:

∆(U(t)) , E[V(U(t+1))−V(U(t))|U(t)],

The flow control algorithm for the supply chain results from minimizing an upper bound of the following quantity

∆(U(t))−δE {h(μ(t))|U(t)}, (30)

for each time slot t. Note that minimizing the previous expression means a trade-off between the stability of the

queues through the Lyapunov drift ∆(U(t)) and the firm’s profit through the profit function h, where δ is a weighing

factor. In fact, making δ large enough implies focusing on maximizing the profit (and getting arbitrarily close to

the optimal solution), but at a cost in terms of an increased product congestion in the queues.

Let us Y,U, μ,A be three non-negative reals so that

Y ≤ max{U −μ,0}+ A.

It is not difficult to show that the following inequality holds:

Y

2 ≤ U2 +μ2 + A2 −2U(μ− A). (31)
14

Using the previous inequality, an upper-bound for (30) is as follows:

∆(U(t))−δE {h(μ(t))|U(t)} ≤ BN ̄ − E







X

i∈M

U

1

i

(t)





X

b

μi,b(t)−μ1,i(t)





|U(t)







−

−E







X

i∈M

U

2

i

(t)





−

X

b

μi,b(t)+μ1,i(t)





|U(t)







− E







X

i∈W

U

1

i

(t)





μi,i

′ (t)−

X

a

μa,i(t)





|U(t)







−

−E







X

i∈W

U

2

i

(t)



−μi,i

′ (t)+

X

a

μa,i(t)





|U(t)







− E







X

i∈W

U

1

i

′ (t)





X

b

μi,b(t)−μi,i

′(t)





|U(t)







−

−E







X

i∈W

U

2

i

′ (t)





−

X

b

μi,b(t)+μi,i

′(t)





|U(t)







− E







X

i∈R

Ui(t)





di(t)−

X

a

μa

′

,i(t)





|U(t)







−

−δE







X

(i

′

, j)

fi

′

, j(μi

′

, j(t))|U(t)







+δE







X

i∈M

gi(ri(t))|U(t)







+δE







X

(i, j)

gi, j(μi, j(t))|U(t)







+

+δE







X

(i,i

′

)

gi,i

′ (μi,i′ (t))|U(t)







+δE







X

(i

′

, j)

gi

′

, j(μi

′

, j(t))|U(t)







,

where

B ,

1

N ̄

X

i∈N

2



μ

max

i

2

,

and where N ̄ is the number of all queues.

A rearrangement of the sums in the previous inequality further produces

∆(U(t))−δE {h(μ(t))|U(t)} ≤

≤ BN + E







X

i∈R

Ui(t)di(t)|U(t)







+ E







X

i∈M

δg1,i(μ1,i(t))+

h

U

1

i

(t)− U

2

i

(t)

i

μ1,i(t)|U(t)







+

E







X

(i,b),i∈M,b∈W

δgi,b(μi,b(t))−

hU

1

i

(t)− U

1

b

(t)

i

+

h

U

2

b

(t)− U

2

i

(t)

iμi,b(t)|U(t)







+

+E







X

(i∈W

δgi,i

′ (μi,i′ (t))−

hU

1

i

(t)− U

1

i

′(t)

i

+

h

U

2

i

′ (t)− U

2

i

(t)

iμi,i

′ (t)|U(t)







+

+ E







X

(i

′

,b),i∈W,b∈R

δgi

′

,b(μi

′

,b(t))−δ fi

′

,b(μi

′

,b(t))−

hU

1

i

′(t)− Ub(t)



− U

2

i

′(t)

i

μi

′

,b(t)|U(t)







. (32)

From the above inequality, the derivation of the algorithm is evident. Given queue levels U(t), the flow control

algorithm follows from greedily minimizing the right-hand side of the inequality (32), in terms of the control

variables μ(t) over all possible flow options satisfying the constraints introduced in Assumption 2.2.
15

B. Performance analysis

This subsection shows that the dynamic flow control algorithm introduced above gets arbitrarily close to the

optimal solution of (14). The next theorem proves to be useful in the analysis of the algorithm.

Theorem 5.1: Let Assumptions 2.1 through 3.2 hold and assume that there exist positive constants δ, ǫ and B

such that for all timeslots t and all backlog queue levels U(t), the Lyapunov drift satisfies:

∆(U(t))−δE{h(μ(t))|U(t)} ≤ B−ǫ

X

N ̄

i=1

Ui(t)−δh

∗

, (33)

where h

∗

is the optimal cost function of the stochastic optimization problem (13). Then the follwing inequalities

are satisfied

lim sup

t→∞

1

t

Xt−1

τ=0





X

2

j=1





X

i∈W

E{U

j

i

(τ)}+

X

i∈M

E{U

j

i

(τ)+ U

j

i

′(τ)}+

X

i∈R

E{Ui(τ)}









≤

B+δ(h ̄ − h

∗

)

ǫ

(34)

lim inf

t→∞

h(μ ̄(t)) ≥ h

∗ −

B

δ

, (35)

where μ ̄(t) was defined in (3) and h ̄ is given by

h ̄ , lim sup

t→∞

1

t

Xt−1

τ=0

E{h(μ(τ))}.

The previous Theorem is a slight modification of Theorem 5.4 in [8] and for brevity the proof is omitted.

Remark 5.1: Note that since the flows μi, j(t) are upper bounded by μ

max

i

and the function h is continuous, there

exists hmax so that h ̄ − h

∗ ≤ hmax. In addition, let μmax , min{μmax

i

}.

The next Theorem describes the performance of the flow control algorithm.

Theorem 5.2: Let Assumptions 2.1 through 3.2 hold. For any positive parameter δ the flow control algorithm

stabilizes the (virtual) queues associated with the constraints of the optimization problem (14) and gives the following

upper bounds:

lim sup

t→∞

1

t

Xt−1

τ=0





X

2

j=1

X

i∈W

E{U

j

i

(τ)}+

X

i∈M

E{U

j

i

(τ)+ U

j

i

′(τ)}+

X

i∈R

E{Ui(τ)





≤

NB+δhmax

μmax

(36)

lim inf

t→∞

h(μ ̄(t)) ≥ h(μ

∗

)−

BN

δ

, (37)

where μ

∗

is the solution of (13) and where μ ̄(t) satisfies (3).

Proof: Let ǫ1 be a small quantity of product flow added to the inputs of queues U

1

i

(t) for all i ∈ M ∪W and

queues Ui(t), for i ∈ R. It follows that the dynamics of the aforementioned queues become

U

1

i

(t+1) = max







U

1

i

(t)−

X

b

μi,b(t),0







+μ1,i(t)+ǫ1,∀i ∈ M,

U

1

i

(t+1) = maxnU

1

i

(t)−μi,i

′(t),0

o

+

X

a

μa,i(t)+ǫ1,∀i ∈ W,

U

1

i

′ (t+1) = max







U

1

i

′(t)−

X

b

μi

′

,b(t),0







+μi,i

′ (t)+ǫ1,∀i ∈ W,

Ui(t+1) = max{Ui(t)−di(t),0}+

X

a

μa

′

,i(t)+ǫ1,∀i ∈ R,
16

and let Λǫ1denote the capacity region of the supply chain under the additional flow ǫ1, and  ̄μ

∗

(ǫ1) denote the

solution of (13), when Λ is replaced by Λǫ1

. Then, by Corollary 3.1 applied to the capacity region Λǫ1, we have

that there exists a stationary randomized flow control algorithm, that chooses the flows based on the current state

of the supply chain, and gives

E{μ

∗

1,i

(ǫ1)}+ǫ1 = E







X

b

μ

∗

i,b

(ǫ1)







,∀ i ∈ M,

X

a

E

n

μ

∗

a,i

(ǫ1)

o

+ǫ1 = E

n

μ

∗

i,i

′ (ǫ1)

o

,∀ i ∈ W,

E

n

μ

∗

i,i

′ (ǫ1)

o

+ǫ1 = E







X

b

μ

∗

i

′

,b

(ǫ1)







,∀ i ∈ W,

E







X

a

μ

∗

a

′

,i

(ǫ1)







+ǫ1 ≤ d ̄

i

,∀ i ∈ R,

where  ̄μ

∗

i, j

(ǫ1) = E

n

μ

∗

i, j

(ǫ1)

o

.

Similarly, assuming that a small flow ǫ2 is added to the inputs of queues Ui(t)

2

, their dynamics become

U

2

i

(t+1) = maxnU

2

i

(t)−μ1,i(t),0

o

+

X

b

μi,b(t)+ǫ2,∀i ∈ M,

U

2

i

(t+1) = max







U

2

i

(t)−

X

a

μa,i(t),0







+μi,i

′ (t)+ǫ2,∀i ∈ W,

U

2

i

′ (t+1) = maxn

U

2

i

′ (t)−μi,i

′(t),0

o

+

X

b

μi

′

,b(t)+ǫ2,∀i ∈ W.

Denoting by Λǫ2the capacity region under the additional flow ǫ2,  ̄μ

∗

(ǫ2) represents the solution of (13) when Λ is

replaced by Λǫ2.

As before, by Corollary 3.1 applied to the capacity region Λǫ2

, there exists a stationary randomized flow control

algorithm, that chooses the flows based on the current state of the supply chain, and gives

E







X

b

μ

∗

i,b

(ǫ2)







+ǫ2 = E{μ

∗

1,i

(ǫ2)},∀ i ∈ M,

E

n

μ

∗

i,i

′ (ǫ2)

o

+ǫ2 =

X

a

E

n

μ

∗

a,i

(ǫ2)

o

,∀ i ∈ W,

E







X

b

μ

∗

i

′

,b

(ǫ2)







+ǫ2 = E

n

μ

∗

i,i

′ (ǫ2)

o

,∀ i ∈ W,

where  ̄μ

∗

i, j

(ǫ2) = E

n

μ

∗

i, j

(ǫ2)

o

.

Note that by Assumption 3.2, such ǫ1 and ǫ2 do exist.

The flow control algorithm described in the previous section minimizes the right-hand side of inequality (32) for

all possible policies based on the current state of the supply chain. In particular, it does this against the previously

mentioned stationary policies, generated by adding the additional flows ǫ1 and ǫ2. Consequently, under the flow

control algorithm, it follows that

∆(U(t))−δE {h(μ(t))|U(t)} ≤ BN ̄ −

X

i∈M

U

1

i

(t)





X

b

μ ̄

∗

i,b

(ǫ1)−μ ̄

∗

1,i

(ǫ1)





−
17

−

X

i∈M

U

2

i

(t)





−

X

b

μ ̄

∗

i,b

(ǫ2)+μ ̄

∗

1,i

(ǫ2)





−

X

i∈W

U

1

i

(t)





μ ̄

∗

i,i

′ (ǫ1)−

X

a

μ ̄

∗

a,i

(ǫ1)



−

−

X

i∈W

U

2

i

(t)



−μ ̄

∗

i,i

′ (ǫ2)+

X

a

μ ̄

∗

a,i

(ǫ2)



−

X

i∈W

U

1

i

′ (t)





X

b

μ ̄

∗

i,b

(ǫ1)−μ ̄

∗

i,i

′ (ǫ1)





−

−

X

i∈W

U

2

i

′ (t)





−

X

b

μ ̄

∗

i,b

(ǫ2)+μ ̄

∗

i,i

′(ǫ2)





−

X

i∈R

Ui(t)





E{di(t)} −X

a

μ ̄

∗

a

′

,i

(ǫ1)



−

−δ

X

(i

′

, j)

fi

′

, j( ̄μ

∗

i

′

, j

(ǫ1))+δ

X

i∈M

g1,i( ̄μ

∗

1,i

(ǫ1))+δ

X

(i, j)

gi, j( ̄μ

∗

i, j

(ǫ1))+

+δ

X

(i,i

′

)

gi,i

′ ( ̄μ

∗

i,i

′ (ǫ1))+δ

X

(i

′

, j)

gi

′

, j( ̄μ

∗

i

′

, j

(ǫ1)),

Denoting ǫ = min{ǫ1, ǫ2}, the above inequality becomes

∆(U(t))−δE {h(μ(t))|U(t)} ≤

≤ BN ̄ −ǫ

X

2

j=1

X

i∈M

U

j

i

(t)−ǫ

X

2

j=1

X

i∈W

h

U

j

i

(t)+ U

j

i

′(t)

i

−ǫ

X

i∈R

Ui(t)−δh(  ̄μ

∗

(ǫ1)).

By Theorem 5.1 it follows that

lim sup

t→∞

1

t

Xt−1

τ=0





X

2

j=1





X

i∈W

E{U

j

i

(τ)}+

X

i∈M

E{U

j

i

(τ)+ U

j

i

′(τ)}+

X

i∈R

E{Ui(τ)}









≤

≤

BN +δ(h ̄ − h(μ

∗

(ǫ1)))

ǫ

≤

BN +δhmax

ǫ

(38)

and

lim inf

t→∞

h(μ ̄(t)) ≥ h(μ

∗

(ǫ1)))−

BN

δ

, (39)

The performance bounds in (38) and (39) hold for any values of ǫi such that 0 < ǫi ≤ μmax, for i = 1,2. However,

the particular values of ǫi only affect the values of the bounds and not the control algorithm. Therefore, the bounds

can be optimized separately over all possible values of ǫi, i = 1,2. Obviously, the bound (38) is minimized when ǫ

approaches μmax. It can be shown that the optimal solution of (13) when the capacity region is replaced by Λǫ1

, is

continuous in ǫ1. Consequently, as ǫ1 approaches zero, the capacity region Λǫ1

approaches Λ and μ

∗

(ǫ1) approaches

μ

∗

. Therefore, the bound (39) is minimized when ǫ1 goes to zero, and the result follows.

Remark 5.2: Note that inequality (36) shows that under the flow control algorithm, the queues remain stable,

i.e., the long-run flows are feasible. In addition, inequality (37) shows that the solution provided by the flow control

algorithm can get arbitrarily close to the optimal solution, by making δ arbitrarily large.
18

Fig. 3: Example of network topology: 2 branches, 2 retailers, up and downstream crossings

VI. Numerical example

The flow control algorithm described in the previous sections was implemented and tested on the supply chain

network shown in Figure 3. The cost function corresponding to each link (i, j) of the network has the form gi, j(μi, j) =

ai, jμ

2

i, j

+ bi, jμi, j, while the revenue function is given by f(μ) = cμ

1

p

i, j

+ d, where ai, j = 0.1, bi, j = 0.3, c = 3, d = 2,

p = 1.8. The maximum output rate at node i is assumed to be equal to μmax = 6×Li, where Liis the number of links

going out of i. This sets an “average” maximum rate of 6 for each link. A link has two states, ON and OFF. and the

links ON-OFF processes are assumed to be i.i.d. with an ’ON’ probability of 0.9. The demand processes are taken

to be independent and uniformly distributed between 0 and 3 at each time, with an average of 1.5. In addition, two

values for the parameter δ are considered, 0.1 and 0.9, respectively, to show its influence on the queues’ backlog.

The queues’ backlog over time as well as the running averages of the queues (Figures 4) for the two considered

values of δ are plotted. Both the plots for the forward U

(1)

i

and backward U

(2)

i

queues are shown. Recall that these

queues are virtual queues introduced as a consequence of modeling the inequality constraints as queues. However,

the forward queues can be interpreted as real queues at nodes of the network. Also, notice that there is no backward

queue defined for the queue at a retailer. The queues’ backlogs of each branch are shown in a 4-by-2 panel where

the left column corresponds to the froward queues (from top to bottom nodes) and the right column corresponds

to the backward queues.

Plots of the flow rates on the different links and their running averages (Figure 5) are also depicted. The link

flow rates are shown in a 2-column panel. The left column shows links (originating) in the left branch and the right
19

0 2000 4000 6000 8000 10000

0

20

40

F Queue 1

0 2000 4000 6000 8000 10000

0

50

B Queue 1

0 2000 4000 6000 8000 10000

0

50

F Queue 3

0 2000 4000 6000 8000 10000

0

20

40

B Queue 3

0 2000 4000 6000 8000 10000

0

20

40

F Queue 5

0 2000 4000 6000 8000 10000

0

20

40

B Queue 5

0 2000 4000 6000 8000 10000

0

50

F Queue 7

(a) (Queues branch 1 for δ = 0.1

0 2000 4000 6000 8000 10000

0

20

40

60

F Queue 2

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 2

0 2000 4000 6000 8000 10000

0

20

40

60

F Queue 4

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 4

0 2000 4000 6000 8000 10000

0

20

40

60

F Queue 6

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 6

0 2000 4000 6000 8000 10000

0

20

40

60

F Queue 8

(b) Queues branch 2 for δ = 0.1

Fig. 4: Queue levels
20

0 2000 4000 6000 8000 10000

0

20

40

F Queue 1

0 2000 4000 6000 8000 10000

0

50

B Queue 1

0 2000 4000 6000 8000 10000

0

50

100

F Queue 3

0 2000 4000 6000 8000 10000

0

20

40

B Queue 3

0 2000 4000 6000 8000 10000

0

20

40

F Queue 5

0 2000 4000 6000 8000 10000

0

20

40

B Queue 5

0 2000 4000 6000 8000 10000

0

50

F Queue 7

(c) Queues branch 1 for δ = 0.9

0 2000 4000 6000 8000 10000

0

20

40

F Queue 2

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 2

0 2000 4000 6000 8000 10000

0

20

40

60

F Queue 4

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 4

0 2000 4000 6000 8000 10000

0

20

40

F Queue 6

0 2000 4000 6000 8000 10000

0

20

40

60

B Queue 6

0 2000 4000 6000 8000 10000

0

20

40

F Queue 8

(d) Queues branch 2 for δ = 0.9

Fig. 4: Queue levels (cont.)
21

column displays the rates of links (originating) in the right branch. In Figure 5a, the left column shows (from top

to bottom) the rates in links L1,L3,L5,L7 of the topology on Figure 3; the right column shows the rate at links

L2,L4,L6,L8. The cross-branch links are shown in the bottom subplots (the rates at links L9,L11 are shown in the

bottom of the left column and links L10,L12’s rates in the bottom of the right column). Finally, to emphasize on

the convergence of the average rates, a zoom-in of the rate plots is presented to focus only on the [0,1.8]-range of

the y-axis (i.e., the rates). This is shown in Figure 6.

A. Discussion: Queues’ backlog

The queues’ backlogs are shown in Figures 4. It can be observed that the queues are oscillating but are not

growing unbounded. This is exactly the stability of queues predicted by the theory. In fact, the average rates also

satisfies the stability condition of Proposition 4.1, as well. The theory however, does not predict anything about the

convergence of the average queues’ backlog. Yet, it can be observed that the average backlog seems to converge

for all (forward and backward) queues. An interesting follow up of this study is to prove/disprove convergence of

average queue and to determine under which conditions convergence is guaranteed.

From the figures, it can be noticed that in general, the forward queues at the manufacturers and at the first (upper)

warehouses are in average more loaded than the queues at the second (lower) warehouses and at the retailers. This

is a consequence of the back-pressure algorithm, which forces upstream nodes to reduce their rate and consequently

build up their queues when downstream nodes are congested. Hence, in general, queues close to the destination

tend to have a smaller backlog. It can be also observed that the cross links serve to balance the load to reduce

the variations in each queue. Finally, it can be noticed that the queue fluctuations increase for higher values of

the parameter δ (Figures 4c-4d). Recall that setting δ large implies focusing on maximizing the profit (and getting

arbitrarily close to the optimal solution), but at the cost of increased product congestion in the queues.

B. Discussion: Link Rates

The rates at the different links are shown in Figures 5 and 6. A certain number of observations can be made

from the figures.

First, the rates are random due to the randomized control algorithm. However, for all runs of the simulation, the

average rate converges for each link. Furthermore, at each retailer, the value to which the average aggregate rate

converges is less than 1.5, the average demand at each market. This is a necessary condition for the stability of the

queues as was stated in Proposition 4.1. The average rate at the other links are such that the conservation of flow

principle is satisfied at each node (which is what was expected).

When there are two links departing from a node, traffic can either be split (when both links are up), or entirely

sent over one link (especially when the other link is down). To see which choice will be made at a given node,

one can analyze the cost function gi, j(μi, j) = aμ

2

i, j

+bμi, j. Assume that at a branching node, traffic is split such that

a rate of μ is sent over one link and 12 −μ over the other link 0 ≤ μ ≤ 12. The (local) total cost of such routing

is aμ

2 +bμ+a(12 −μ)2 +b(12 −μ) = a(2μ2 −24μ) +12b +144a. Analyzing this cost as a function of μ, it can be
22

observed that it is minimized when μ = 12, implying that at a branching node, when both links are up, the entire

traffic should be sent over one of the links. This is what is observed at nodes 2, 3, 6 and 7 for network in Figure 3,

where the traffic on the links departing from such nodes is (almost all the time) either 0 or equal to the maximum

rate of 12.

VII. Acknowledgements

The authors would like to express their gratitude to Allison Barnard Feeney, Mark Campanelli and Sharon J.

Kemmerer for their valuable comments and suggestions.

VIII. Conclusions

In this paper the management of flow product in a supply chain was addressed. Generally speaking, the main

contribution of the paper to the literature consists of the introduction of a distributed algorithm for the flow

management in a random and time-varying supply chain, that is not based on stochastic approximation. In more

detail, the contributions are as following. Motivated by recent events, a random and time-varying model for a supply

chain was proposed which induced a stochastic nature of the flows. A stochastic optimization problem aimed at

maximizing the profit function of a firm in terms of the time-averages of the flows and subject to constraints induced

by the supply chain was formulated. A distributed, dynamic algorithm for solving the aforementioned optimization

problem was proposed. Under this algorithm, at each time instant decisions are based only on the current state of

the supply chain. In addition, decisions do not need information on the probability distribution of the supply chain.

It was shown that the long-run, time-averages of the flows generated by the algorithm can get arbitrarily close to

the optimal solution of the stochastic optimization problem.

Another, indirect contribution of the paper is that it exposes the reader to new techniques for solving stochastic

optimization problems. This is beneficial to the operational research literature based in part on optimization theory.

References

[1] C.H. Aikens. Facility location models for distribution planning. European Journal of Operational Research, 22(3):263 – 279, 1985.

[2] B. Behdani. Risk management and paradoxical situations in modern supply chains. The CIP Report on Supply Chain, July 2011.

[3] R. Blasgen. Supply chain management: The driving force behind the global economy. The CIP Report on Global Supply Chain, July

2011.

[4] Canadian Competition Bureau. Competition bureau concludes examination into gasoline price spike following hurricane Katrina

(http://www.competitionbureau.gc.ca), 2006.

[5] S. Chopra and P. Meindl. Supply Chain Management (5th Edition). Prentice Hall, 2012.

[6] M.C. Chou, G.A. Chua, C.P. Teo, and H. Zheng. Process flexibility revisited: The graph expander and its applications. Operations Research,

59(5):1090 – 1110, 2003.

[7] J. Dong, D. Zhang, and A. Nagurney. Supply chain supernetworks with random demands. in Urban and Regional Transportation Modeling:

Essays in Honor of David E. Boyce, pages 289–313, 2004.

[8] L. Georgiadis, M.J. Neely, and L. Tassiulas. Resource allocation and cross layer control in wireless networks. Foundations and Trends in

Networking, 1(1):1–144, 2006.

[9] H.J. Kushner. Stochastic stability and control. Acad. Press, 1967.
23

0 2000 4000 6000 8000 10000

0

10

20

Link 1

0 2000 4000 6000 8000 10000

0

10

20

Link 2

0 2000 4000 6000 8000 10000

0

10

20

Link 3

0 2000 4000 6000 8000 10000

0

10

20

Link 4

0 2000 4000 6000 8000 10000

05

10

Link 5

0 2000 4000 6000 8000 10000

05

10

Link 6

0 2000 4000 6000 8000 10000

0

10

20

Link 7

0 2000 4000 6000 8000 10000

0

10

20

Link 8

0 2000 4000 6000 8000 10000

0

10

20

Link 9

0 2000 4000 6000 8000 10000

0

10

20

Link 10

0 2000 4000 6000 8000 10000

0

10

20

Link 11

0 2000 4000 6000 8000 10000

0

10

20

Link 12

(a) Rates forδ=0.1

0 2000 4000 6000 8000 10000

0

10

20

Link 1

0 2000 4000 6000 8000 10000

0

10

20

Link 2

0 2000 4000 6000 8000 10000

0

10

20

Link 3

0 2000 4000 6000 8000 10000

0

10

20

Link 4

0 2000 4000 6000 8000 10000

05

10

Link 5

0 2000 4000 6000 8000 10000

05

10

Link 6

0 2000 4000 6000 8000 10000

0

10

20

Link 7

0 2000 4000 6000 8000 10000

0

10

20

Link 8

0 2000 4000 6000 8000 10000

0

10

20

Link 9

0 2000 4000 6000 8000 10000

0

10

20

Link 10

0 2000 4000 6000 8000 10000

0

10

20

Link 11

0 2000 4000 6000 8000 10000

0

10

20

Link 12

(b) Rates forδ=0.9

Fig. 5: Link Rates
24

0 2000 4000 6000 8000 10000

0

1

Link 1

0 2000 4000 6000 8000 10000

0

1

Link 2

0 2000 4000 6000 8000 10000

0

1

Link 3

0 2000 4000 6000 8000 10000

0

1

Link 4

0 2000 4000 6000 8000 10000

0

1

Link 5

0 2000 4000 6000 8000 10000

0

1

Link 6

0 2000 4000 6000 8000 10000

0

1

Link 7

0 2000 4000 6000 8000 10000

0

1

Link 8

0 2000 4000 6000 8000 10000

0

1

Link 9

0 2000 4000 6000 8000 10000

0

1

Link 10

0 2000 4000 6000 8000 10000

0

1

Link 11

0 2000 4000 6000 8000 10000

0

1

Link 12

(a) Rates for δ = 0.1

0 2000 4000 6000 8000 10000

0

1

Link 1

0 2000 4000 6000 8000 10000

0

1

Link 2

0 2000 4000 6000 8000 10000

0

1

Link 3

0 2000 4000 6000 8000 10000

0

1

Link 4

0 2000 4000 6000 8000 10000

0

1

Link 5

0 2000 4000 6000 8000 10000

0

1

Link 6

0 2000 4000 6000 8000 10000

0

1

Link 7

0 2000 4000 6000 8000 10000

0

1

Link 8

0 2000 4000 6000 8000 10000

0

1

Link 9

0 2000 4000 6000 8000 10000

0

1

Link 10

0 2000 4000 6000 8000 10000

0

1

Link 11

0 2000 4000 6000 8000 10000

0

1

Link 12

(b) Rates for δ = 0.9

Fig. 6: Link Rates - zooming
25

[10] Y.K. Lin. A simple algorithm for reliability evaluation of a stochastic-flow network with node failure. Computers and Operations Research,

28(13):1277 – 1285, 2001.

[11] M. Nagarajan and G. Sosic. Game-theoretic analysis of cooperation among supply chain agents: Review and extensions. European Journal

of Operational Research, May 2006.

[12] A. Nagurney. Supply chain network design under profit maximization and oligopolistic competition. Transportation Research E, (46):281–

294, 2010.

[13] Q. Qiang and A. Nagurney. Modeling of supply chain risk under disruptions with performance measurement and robustness analysis. in

Managing Supply Chain Risk and Vulnerability: Tools and Methods for Supply Chain Decision Makers, Editors T. Wu and J. Blackhurst,

Springer, Berlin, Germany, 2009, pages 91–111, 2009.

[14] W.F. Samuelson and S.G. Marks. Managerial Economics, 4th ed. Wiley, 2003.

[15] T. Santoso, S. Ahmed, M. Goetschalckx, and A. Shapiro. A stochastic programming approach for supply chain network design under

uncertainty. European Journal of Operational Research, (167):96–115, 2005.

[16] R. Skulte and T. Wikerson. Critical infrastructure and the supply chain. The CIP Report on Global Supply Chain, July 2011.

[17] I. Varkony. Critical infrastructure protection: Volatility and risks in the global chain. The CIP Report on Global Supply Chain, July 2011.

[18] C.J. Vidala and M. Goetschalck. A global supply chain model with transfer pricing and transportation cost allocation. European Journal

of Operational Research, 129(1):134 – 158, 2001.

[19] D. Zhang, J. Dong, and A. Nagurney. A supply chain network economy: Modeling and qualitative analysis. in Innovations in Financial

and Economic Networks, 2003.