CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 52, 2016 

A publication of 

 

The Italian Association 
of Chemical Engineering 
Online at www.aidic.it/cet 

Guest Editors: Petar Sabev Varbanov, Peng-Yen Liew, Jun-Yow Yong, Jiří Jaromír Klemeš, Hon Loong Lam 
Copyright © 2016, AIDIC Servizi S.r.l., 

ISBN 978-88-95608-42-6; ISSN 2283-9216 
 

An Economic Model Predictive Approach for Inventory 

Routing and Control with Time Windows Constraints: 

Application in the Distribution of Industrial Gases 

Athanassios Nikolakopoulos 

National Technical University of Athens, Heroon Polytechneiou 9, P.C. 15780, Athens, Greece  

nikolako@mail.ntua.gr 

In the process of industrial gas distribution, the uncertainty of future demand presents an important challenge 

for the production of efficient delivery plans, where the overall cost is minimized and product shortages must 

be avoided. The objective of this paper is to present a new method for producing future delivery plans that 

eliminate stock-out risk and minimize transportation and inventory holding costs, in the face of time window 

constraints. The proposed method follows a Economic Model Predictive Control (EMPC) approach, which 

makes use of a new mathematical programming model of the distribution process to predict future inventory 

levels. The method is illustrated using an example with a simple forecasting approach. It is shown that the 

stock-outs risk can be eliminated, but still balanced decisions between cost saving and risk aversion remain at 

the discretion of the planner. It is also shown that coordinated inventory routing and control through the EMPC 

approach brings economic benefits and hedges effectively against the market volatility, while opting for 

minimum overall costs may incur excessive shortage risk. 

1. Introduction 

In the present study the Vendor Managed Inventory (VMI) strategy is considered for the distribution of 

liquefied industrial gases, with a combined backhauling of the empty tanks. Following a VMI strategy a tool 

that optimizes dispatches produces short-term plans by solving Inventory Routing Problems (IRP) to decide 

when, how much to deliver to each customer, and what will be the delivery routes. Real problems involve 

inherently time constraints that appear in the form of hard or soft time windows, within which vehicles must 

reach the customers’ sites and start the delivery. While deterministic planning models are useful tools to help 

reduce costs by taking into account customer synergies in the inventory-distribution planning, uncertain 

demand fluctuations may significantly affect the decision-making across the industrial gas supply chain. Thus, 

it is necessary to extend the deterministic planning models to address these uncertainties, and develop 

effective optimization and control algorithms for these problems (You et al, 2011b). Supply chain managers 

also wish to balance profit maximization versus minimization of the stock-out risk. One way to minimize risk is 

to hold some buffer stock so that the supply chain can respond to rush orders or demand spikes. The present 

work deals with demand uncertainties by employing an EMPC approach to control customer inventories along 

forecasted trajectories of safety stock levels. The literature on the application of MPC on distribution systems 

involving reverse logistics, time constraints, routing and inventory control at the customers’ sites is absent. An 

MPC application to similar supply chain systems is that of Seferlis and Gianellos (2004), however, there were 

no routing decisions involved. Here, a multi-period deterministic model is first presented for the Inventory 

Routing with Simultaneous Pickup and Delivery and Time Windows problem (IRPSPDTW), which captures all 

attributes of reverse logistics, time constraints and inventory routing, as an extension of the single period 

model proposed by Nikolakopoulos (2015). Next, following a centralized EMPC strategy, the IRPSPDTW is 

solved with a rolling horizon. A solution is obtained over the entire time horizon but only the part of the solution 

related to the first time period is actually implemented before moving to the next period. Then the procedure is 

repeated with updated sequences. The methodology is illustrated using an example where an inventory 

                                

 
 

 

 
   

                                                  
DOI: 10.3303/CET1652155 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Nikolakopoulos A., 2016, An economic model predictive approach for inventory routing and control with time 
windows constraints: application in the distribution of industrial gases, Chemical Engineering Transactions, 52, 925-930  
DOI:10.3303/CET1652155   

925



routing systems for the distribution of industrial gases is stabilized in the face of hard time constraints, and 

balanced decisions are produced between the minimization of the overall cost and the minimization of stock-

out risk. 

2. Mathematical programming model of the IRPSPDTW 

In the multi-period IRPSPDTW an homogeneous fleet V with Kavl = V  available vehicles of equal capacity Q 

must serve a set N of n customers, N = {1, 2…, n}. The planning horizon considers T periods, t = 1, …, T. 

There is a demand di,t for every customer i N. Each customer i maintains its own inventory up to capacity Ci 
pays for inventory carrying cost of hi units per period per stock keeping unit. The complete directed graph 

induced by the customers is G = (N+, A), where N+ = N  {0, n+1} is the set of vertices. The depot is 
represented by vertices 0 and n + 1. A is the set of arcs connecting any pair of vertices, while no arc 

terminates in vertex 0, and no arc originates from vertex n + 1. A travelling cost ci,j  is assigned to each arc (i, j)

 A and  max ,max i jc c . Also, a traveling time tri,j corresponds to each arc (i, j) A, while a time window [ei, 

li] is assigned to each vertex i N+. This is in fact a time restriction, meaning that the time Tri,t at which service 
to customer i in period t begins, must lie within the respective time window. If in period t a vehicle reaches 

customer i at time Tri,t < ei then it waits for a time period equal to ei – Tri,t until the service begins. The time 

window, in which a vehicle can leave from or arrive at the depot is [e0, l0] = [en+1, ln+1] = [E, L]. Zero delivery 

demands and pickup demands are assigned to nodes {0, n+1}, that is, d0 = dn+1 = p0 = pn+1 = 0. By allowing 

some of the vehicles to remain at the depot there is no cost contributed. These vehicles do not serve any 

customer and their routes are represented by arc (0, n+1) to which zero travelling cost and time are assigned: 

c0,n+1 = tr0,n+1 = 0. The objective is to minimize the overall transportation and inventory carrying cost incurred 

over a specific planning horizon. Each customer is serviced at most once in each period t, ( ,i ip d Q ,  i

N), capacity constraints must be satisfied, and no split delivery or backlogging is permitted. 

The next five sets of variables are considered for the problem: 

zk,i,,j,t = 
1 if vehicle drives from to in period

0

     

 otherwise

k i j tìïï
í
ï
ïî

,     (i,j) A, k  V and t = 1, ..., T 

Yk,i,j,t and Pk,i,j,t are the vehicle loads of full and empty tanks on that trip respectively, and Ii,t and EIi,t are the 

inventories of full and empty tanks at customer i, at the end of period t respectively. The IRPSPDTW can then 

be formally described as the following multi-period commodity network flow model with capacity and time 

window constraints: 

, , , , ,

1 ( , )

T

i j k i j t i i t

t k V i j A i N

Min c z h I
   

  
 

  
     (1) 

s.t.      

, , ,

0

and1 , , 1,...,
N

k i j t

j
j i

z i N k V t T



      (2) 
 

{0} { 1}

, , , , , ,

0 0

0 , , 1,...,
N N n

k i m t k l i t

m l
m i l i

z z k V i N t T
  

 
 

         (3) 

, , , , , ,
0 , , , , 1,...,

k i j t k i j t
Y z Q k V i j N i j t T         

(4) 

 
{0} { 1}

, , , , , ,

0 0

0 , , 1,...,
N N n

k l i t k i m t

l m
l i m i

Y Y k V i N t T
  

 
 

         (5) 
 

, , , , , ,
0 , , , , 1,...,

k i j t k i j t
P z Q k V i j N i j t T         

(6) 

 
{ 1} {0}

, , , , , ,

0 0

0, , , 1,...,
N n N

k i m t k l i t

m l
m i l i

P P k V i N t T
  

 
 

         (7) 
 

, 1 , , , , , , , ,

1 0 0

, 1,...,
V N N

i t i t k l i t k i m t i t

k l m
l i m i

I I Y Y d i N t T


  
 

 
 

      
 
 
 

    
(8) 

 

926



, , 1 , , , , , , ,

1 0 0

, 1,...,
V N N

i t i t k i m t k l i t i t

k m l
m i l i

EI EI P P d i N t T


  
 

 
 

      
 
 
 

    
(9) 

 

, , , , , ,
, , , 1,...,

k i j t k i j t
Y P Q i N k V t T


        

(10) 

 

, ,
, , 1,...,

i t i t i
I EI C i N t T      

(11) 

 

) - -
*

, , , , , , , , , , ,
( (1 ) , ( , ) , , 1,...,

k i t k i j t i i j k j t k i j t
Tr z st t Tr z T i j A k V t T         

(12) 

 

, , , , , , , ,
, , , 1,...,

i k i j t k i t i k i j t

j N j N

e z Tr l z i N k V t T
 



 

   
         

   
   
   

(13) 

 

, ,
, {0, 1}, , 1,...,

k i t
E Tr L i n k V t T         

(14) 

 

, ,{ 1},
0, , , 1,...,

k i n t
Y i N k V t T


       

(15) 

 

,{0}, ,
0, , , 1,...,

k j t
P j N k V t T       (16) 

, , , , , , , ,
0, 0, 0, 0, , , , 1,...,

k i j t k i j t i t i t
Y P I EI i j N k V t T          

(17) 

 

{0, 1}
, , , , , , , , , , ,

, , , , , ( , ) , , 1,...,
k i j t k i j t k i j t i t i t

z Y P I EI Z i j A k V t T        
(18) 

 

The objective function (1) minimizes the inventory and the transportation costs. Constraints set (2) states that 

a vehicle will visit a location no more than once in a time period, and constraint set (3) ensures route 

continuity. Constraints (4) and (6) ensure that the number of tanks transported between two locations will 

always be zero whenever there is no vehicle moving between these locations, and that the amount 

transported is less than or equal to the vehicle’s capacity. Constraint sets (5) and (7) along with the other 

elements of the model ensure that efficient solutions will not contain sub-tours. Constraint sets (8) and (9) are 

the balances of full and empty tanks at customer i at the end of period t. According to the capacity constraint 

set (10), when a vehicle leaves a customer, the total load must not exceed the capacity Q. Constraint set (11) 

ensures that the amount of tanks (empty and full) at a customer cannot exceed its’ capacity to hold inventory 

Ci. Constraint set (12) forces continuity for the travel times and the times that vehicles reach the customers. 

The constant value 
*
T must be at least equal to (L - E), otherwise feasible solutions could be potentially cut 

off. Constraint set (13) restrict the times at which vehicles reach the customers within their time windows. 

Constraint set (14) forbids vehicles to leave from or return at the depot in times outside the depot’s operating 

time window. According to constraint sets (15) and (16), the vehicles must return to the depot empty of full 

tanks, and must leave the depot without empty tanks. Constraints (17) and (18) are the variable positivity and 

integrity constraints. 

3. The Economic Model Predictive Inventory Control Strategy 

In this work, the MPC strategy described in (Camacho and Bordons 1999) is adjusted for solving the inventory 

routing and control problem in the face of unknown future demand and the presence of time window 

constraints. According to the approach, the future outputs for the prediction horizon Np, are predicted at each 

instant t using the process model (2 - 18). The set of future control signals is calculated by optimizing an 

objective function, which balances the economic cost criterion of Eq(1) and a set-point tracking criterion, which 

is a quadratic function. If the control horizon is decided to be Nc = 1, the control signal u(t|t) (deliveries) is 

applied by the distribution system for only one period ahead, whilst the future calculated deliveries are 

rejected, because at the next sampling instant the output of the system (i.e. the demand and the inventories 

one period ahead) y(t+1) are already known and the procedure repeats itself with this new value and all the 

sequences are brought up to date. Thus u(t+1|t+1) is calculated using the receding horizon concept. 

The basic structure of the adopted strategy is presented in Figure 1. The desired trajectories of safety 
inventory levels are compared to the predicted outputs of the method and the errors are fed to the optimization 
model. The optimization model takes into account:  

- The future customers’ demands predicted by the forecasting model, 
- The capacity and routing constraints and  
- The cost function, where the future safety inventory tracking and the economic cost error are 

balanced. 

927



The future deliveries plan is produced by the optimization model (Eqs(2) - (18) and (20) - (23)) and used as 
future input signals for model (2 - 18) which predicts the future inventory levels. The procedure iterates to 
cover the whole planning period, by using the control rule of applying to the real system only the first input 
result for the deliveries: u(t|t) at each iteration and using the receding horizon concept to update forecasts and 

decision variables. 
 

 

Figure 1:  Structure of the Model Predictive strategy for Inventory Routing and Control 

We use a state space representation of the system: 

d
x Ax Bu B d

    (19) 

in which 
N

x Z are the system states i.e. the customer inventories 
1 2

...
N

I I I    ;  , ,i t i tx I . The manipulated 

input is 
N

u Z , i.e. the quantities to be delivered 
1 2

...
N

u u u    , where , , , ,
1 0 0

V N N

i k l i k i m

k l m
l i m i

u Y Y
  

 

 
 

 
 
 
 

   . The 

disturbance to the system is the customers’ demand 
N

d R , or 
1 2

...
N

d d d    . The outputs of the system are 

the states y  x. 

3.1 The EMPC model 
For the nominal disturbance (demand) ds, assuming stability of the system (Amrit and Rawlings, 2011) (A is 

the identity matrix), we adopt the model of the EMPC problem (20).  

   

 

1

0

min ( ) (1 ) ( )

. . ( 1) ( ) ( ) ( ), 0,..., 1

Np
c tr

t

d

ag x t a g x t

s t x t Ax t Bu t B d t t Np

u






    


     



 (20) 

The transportation and inventory cost are formalized by: 

  , , , , ,
( , )

( )
c

i j k i j t i i t

k V i j A i N

g x t c z h x
  

     (21) 

The safety inventory tracking cost is: 

     
'

, ,

1
( )

2

tr

t t s t j s
g x t x x P x x    (22) 

where matrix P introduces the cost of tracking the safety stock xs. A, B and Bd are identity matrices, and , , ,k i j tz , 

x(t) and u(t) derive from Eq(8), with respect to Eqs(2-7) and (9-18) for period t = 1,..., Np. Coefficients a and (1 
- a), where a ∈ [0, 1], are relative weightings assigned to the economic and the tracking costs respectively.  
Linearization of the quadratic term of the objective function 

Since P is diagonal, the non-linear term    
'

, ,

1

2
t t s t j s

x x P x x  , can be replaced by 

 2 2, , , , , ,
1

1
( ) 2

2

N

i i t i t i j s i j s

i

diag P x x x x


  
  . 

Constraints 

Eqs(2-18) 

Future  

errors 

Constraints 

Eqs(2-18) 

Objective function 

Eqs (20 - 23) 

Historical 

demand data 

Control 

Rule 

Optimization 

Model   

Forecasting model 

Eq(26) 

Predicted outputs 

(Future inventories) 

Set point 

(Trajectories of safety 

inventories) 

Process model  

Eqs(2-18) 

+ 
- 

Future input signals 

(Future deliveries) 

928



Then by substituting 
2

1,t t
w x  2,t tw x  Eq(20) becomes:  

   

 

2

, 1, , 2, , , ,

1

2

1, 1, 2, ,

2, 1, 2, ,

, ,

1

1
( ) ( ) 2

2

also let    4 ... and

2 ... ( )

1, 0,1

N
tr

j t j t j t j s j s

j

t t t Q t

t t t Q t t

Q

k t k t

k

g x t diag P w w x x

w y y Q y

w y y Qy x Q

y y






     


    


     

 






 (23) 

Model (20) with the linearization constraints of (23) is a mixed integer linear model, which can be solved to 
optimality within reasonable computational times for moderate instances. 

4. Illustrative Example 

The methodology is illustrated using an example that involves 1 supplier, 3 customers, and 3 vehicles. The 

vehicles’ capacity is Q = 20. The initial inventories of full and empty tanks are I0 = [2 15 20] and P0 = [11 5 2] 

respectively and the unit inventory holding cost is h = 1. The planning horizon spans 12 time periods. The unit 

transportation costs of node pairs are: c0, 1 = 15, c0, 2 = 18, c0, 3 = 22, c1, 2 = 32, c1, 3 = 14, c2, 3 = 34. The 

travelling times between node pairs are: t0, 1 = 1.8, t0, 2 = 2.1, t0, 3 = 2.6, t1, 2 = 3.8, t1, 3 = 1.6, t2, 3 = 4. Both cost 

and travelling time tables are symmetric. The average service time at each customer is set to st = 0.5. The 

earliest and latest times at which the service can begin at each customer are: [ei, li] = [6, 8] for i = 1, 2 and [6, 

9] for i = 3. The time window for the depot is [e0, l0] = [e4, l4] = [5, 14]; no vehicle can start before 5 and return 

latter than 14 time units (hours in a day). The targeting cost of the safety stock is a diagonal matrix P, where 

Pi,j = 2 for i = j.  

  

  

Figure 2: Inventory responses for a = 0, 0.2, 0.3 and 0.5 

The following demand forecast is adopted for predicting the future demands as disturbances of the dynamical 

system: perfect demand information is assumed for one period ahead, thus Nc = 1, (Np = 12). For the 

0 2 4 6 8 10 12
0

5

10

15

20

25

Periods

L
e
v
e
l

Inventory responses

 Customer 1

 Customer 2

 Customer 3

a = 0

0 2 4 6 8 10 12
0

5

10

15

20

25

Periods

L
e
v
e
l

Inventory responses

 Customer 1

 Customer 2

 Customer 3

a = 0.2

0 2 4 6 8 10 12
0

5

10

15

20

25

Periods

L
e
v
e
l

Inventory responses

 Customer 1

 Customer 2

 Customer 3

a = 0.3

0 2 4 6 8 10 12
0

5

10

15

20

25

Periods

L
e
v
e
l

Inventory responses

 Customer 1

 Customer 2

 Customer 3

a = 0.5

929



remainder of the horizon, the demand forecast is set equal to the nominal demand ds = [10 7 13]. It is 

assumed that the customers actual demand, reported in Table 1, is normally distributed around the nominal 

demand with variance σ2 = 0.5. The targeted safety stock is the solution of ( 1) ( ) ( ) ( )
d

x t Ax t Bu t B d t     for one 

period delivery delay: xs = [10 7 13]. Figure 2, presents the response of the inventories with α varying from 0 

(i.e. strategy for maximum shortage risk aversion) to α = 0.5 (i.e. opt for the economic criterion). In the first 

case, the results show that the safety stock levels are reached promptly and maintained throughout the whole 

planning horizon. For a = 0.2, the response of the closed-loop system, immediately settles to the safety steady 

state. For increasing values of a (0.3 and 0.5) the deviation of the inventory levels from the safety trajectory 

grows larger. This translates into reduced transportation and inventory holding costs. 

Table 3 reports the distribution of costs when a varies within [0, 1], where tradeoffs can be observed between 
the cost of tracking the safety inventories and the total, transportation and inventory costs. It is shown that 
balanced decisions can be made based on the change of a single parameter a; additionally the desired 
inventory tracking can be matched with the best possible routing schedules. The tracking cost grows 
exponentially with a, while the cost reduction does not follow the same pattern. For values of a near 1 the 
supply chain is exposed to excessive stock-out risk without equivalent economic rewards. 

Table 1:  Cost distribution with varying a 

a 0 0.2 0.3 0.4 0.5 0.8 0.95 1 

Total cost 1,680 1,562 1,542 1,449 1,333 1,199 1,080 1,045 

Transp. cost 1,320 1,205 1,182 1,100 979 906 878 855 

Inv. cost 360 357 360 349 354 293 202 190 

Track. cost 0 25 28 93 192 481 1,654 1,902 

5. Conclusions 
The paper presents a linear deterministic model for the multi-period IRPSPDTW and a dynamical predictive 
inventory control strategy for the industrial gas distribution business. It is shown that by using the proposed 
economic MPC strategy, the supply chain can be stabilized. The method offers a useful tool for producing 
balanced decisions between economic costs and tracking of the safety inventory while the produced 
schedules are always feasible and with optimal routing. It is also shown that balanced decisions can be made 
based on the results of the EMPC model by manipulating a single parameter. However, the supply chain is 
exposed to excessive stock-out risk without equivalent economic rewards for extreme values of this 
parameter. Whichever is the case, the model produces operationally feasible plans at minimum possible 
transportation cost. Future research may consider the development of custom algorithms for the solution of 
problem instances with higher cardinality that will be required since the IRP is a problem known to be NP-
hard.  

Reference 

Amrit R., Rawlings J.B., Angeli D., 2011, Economic optimization using model predictive control with a terminal 

cost, Annu. Rev. Control 35(2) 178-186. 

Camacho E.F., Bordons C., 1999, Model Predictive Control. Springer, Berlin, Germany. 

Nikolakopoulos A., 2015, A Metaheuristic Reconstruction Algorithm for Solving Bi-level Vehicle Routing 

Problems with Backhauls for Army Rapid Fielding, Springer Operations Research/ Computer Science 

Interfaces Series 56, 141-157. 

Seferlis P., Giannelos N.F., 2004, A two-layered optimization-based control strategy for multi-echelon supply 

chain networks, Comp. Chem. Eng. 28, 799-809. 

You F., Pinto J.F., Capon E., Grossmann I.E., Arora N., Megan L., 2011, Optimal Distribution-Inventory 

Planning of Industrial Gases: I. Fast Computational Strategies for Large-Scale Problems, Ind. Eng. Chem. 

Res. 50, 2910-2927. 

 

 

 

930