Robust supply chain performance via Model Predictive Control

Computers and Chemical Engineering 33 (2009) 2134–2143

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Robust supply chain performance via Model Predictive Control Xiang Li, Thomas E. Marlin ∗ Department of Chemical Engineering, McMaster University, 1280 Main Street West, Hamilton, ON, L8S 4L7, Canada

a r t i c l e

i n f o

Article history: Received 31 October 2008 Received in revised form 9 June 2009 Accepted 22 June 2009 Available online 4 July 2009 Keywords: Supply chain Robust MPC Stochastic optimization

a b s t r a c t This paper presents a novel robust Model Predictive Control (MPC) method for real-time supply chain optimization under uncertainties. This method optimizes the closed-loop economic performance of supply chain systems and addresses different sources of uncertainties located external to and within the feedback loop. The future system behavior is predicted by a closed-loop model, which includes both the open-loop system model and a controller model described by an optimization problem. The robust MPC formulation involves the solution of a constrained, bi-level stochastic optimization problem, which is transformed into a tractable problem involving a limited number of deterministic conic optimization problems solved reliably using an interior point method. The robust controller is applied to a real industrial multi-echelon supply chain optimization problem, and its performance is shown to reduce stock-outs without excessive inventories. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Supply chain optimization, or supply chain management, is “a set of approaches utilized to efficiently integrate suppliers, manufacturers, warehouses and stores, so that merchandise is produced and distributed at the right quantities, to the right locations, and at the right time, in order to minimize system-wide costs while satisfying service level requirements” (Simchi-Levi, Kaminski, & Simchi-Levi, 1999). With large incentives to simultaneously improve customer service and reduce inventories, supply chain optimization has become essential for industry, and it has attracted intense interest in the process systems engineering community (Shah, 2005). Traditionally, supply chain management has employed heuristics or mathematical programming techniques for simplified representations of real systems, e.g., methods ignoring capacity constraints. The Materials Resource Planning (MRP) and Enterprise Resource Planning (ERP) methods are widely available in commercial software products. Recently, ERP models have been extended and solved using more accurate models and rigorous mathematical programming methods (e.g., Pinedo, 2000). Even these recent approaches have not considered model uncertainty or the integration of real-time measurements. As indicated by the thorough review by Sarimveis, Patrinos, Tarantilis, and Kiranoudis (2008), control technology can be applied to regulate the dynamics of the supply chain system with periodical feedback information. Model Predictive Control (MPC) is a

∗ Corresponding author. Tel.: +1 905 525 9140x27125; fax: +1 905 521 1350. E-mail address: [email protected] (T.E. Marlin). 0098-1354/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2009.06.029

multivariable control method that has been successfully applied to industry for decades (Maciejowski, 2002; Qin & Badgewell, 2003). Due to its ability to control multivariable systems with input and output constraints, MPC is a natural choice to control or optimize supply chain systems (e.g., Perea-Lopez, Ydstie, & Grossmann, 2003; Tzafestas, Kapsiotis, & Kyriannakis, 1997). However, uncertainties can degrade MPC performance, especially for the supply chain system where uncertainties are usually significant. To prevent significant performance degradation caused by uncertainties, Wang, Rivera, and Kempf (2007) introduce an upper level stochastic optimizer that executes infrequently and provides constraint back-off parameters to the lower level MPC controller. However, this approach cannot respond quickly to changes in the control structure, such as when a manipulated variable is temporarily placed on manual or is placed in operation after having been in manual. A better way to address uncertainties is to consider the uncertainties explicitly in every MPC calculation. There are basically three types of robust MPC algorithms: (a) Adopt the idea from the traditional optimal control theory and basically aim at robust stability (e.g. Kothare, Balakrishnan, & Morari, 1996). This type of robust MPC cannot explicitly address disturbance uncertainty (which is typically the demand uncertainties in supply chain optimization) and sacrifices some dynamic performance for robust stability. (b) Parameterize the future control actions with linear affine output or state feedback laws (Goulart, Kerrigan, & Maciejowski, 2006; Van Hessem, 2004). This method successfully addresses disturbance uncertainty, but it becomes intractable when applied to problems with uncertainty in the feedback model.

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

Nomenclature Variables in supply chain modeling Bo matrix in the state lower bounds (15) CI cost of the inventories penalty on back orders in objective (11) CIO CIP–SKU constant convert IP unit to SKU unit CS cost of the SKU in transportation penalty on back orders in objective (10) Co Cu cost of the system decisions cost of the system states Cx D customer demand I1 IP inventory at plant IP storage I2 SKU inventory at plant distribution center I3 SKU inventory at regional distribution center I3 -O I3∗ P IP manufacturing rate S SKU in transportation F1 IP flow from IPM to IPS IP flow from IPS to SKUM F2 F3 SKU flow from SKUM to DC F4 SKU shipping flow from DC to RDC F4,max transportation capacity F5 arrival SKU shipment at RDC SKU flow from RDC to costumer F6 Rs SKU manufacturing rate T discrete model time interval Tc measurement and controller execution period TP time between two P decisions SKU manufacturing machine running time Ts Ts,max upper bound on SKU production time  ration of F4 entering different virtual routes  SKU transportation time Variables in robust MPC formulation d disturbance d˜ disturbance over future n time steps dˆ disturbance over future 2n − 1 time steps e feedback error in MPC formulations n horizon for MPC calculation t degrees of freedom of robust MPC (30) u decision variables umin lower bounds on decision variables umax upper bounds on decision variables ur reference for the decision variables ur in future n time steps u˜ r x state variables xmin lower bounds on state variables reference for the state variables xr x˜ r xr in future n time steps A, B, Bd coefficients of open-loop model Lup , Mup , Nu coefficients of closed-loop model for decision variables Lxp , Mxp , Nx coefficients of closed-loop model for state variables Kx , Ke , Kd , Kxr , Kur coefficients of control law (27) ˛ confidence of each chance constraint Subscripts i time step number in QP problem (26) j indicates different virtual routes in the disjunctive model (22)–(24) k time step number l indicates different upper bounds on decision variables (33)

m1 , m2 p

2135

special time step numbers in Eqs. (2) and (3), respectively indicates the uncertain value of the variable with presence of uncertainty

(c) Approximate the future controller action with unconstrained nominal MPC (Warren & Marlin, 2003), which cannot handle the hard bounds on the decision variables explicitly. While uncertain feedback systems could be analyzed using dynamic programming, the computations are intractable for largescale systems. Therefore, approximate dynamic programming has been developed to address uncertain dynamic systems such as supply chains (Lee & Lee, 2004; Powell & Van Roy, 2004). Again, this method requires knowledge of not only uncertainties, but also all control structures that might occur. In this paper, robust MPC is employed for supply chain optimization with significant uncertainties. A novel robust MPC method is developed which (1) accounts for uncertainty in both disturbances and feedback model mismatch, (2) accounts for correlations in the uncertainties, (3) allows hard bounds on the decision variables, and (4) involves the real-time solution of a tractable, convex, second order conic optimization problem. The method is introduced and applied to a real industrial multi-echelon supply chain optimization problem in this paper. Section 2 introduces the industrial supply chain problem. Section 3 describes the nominal MPC formulation. Section 4 presents the robust MPC formulation and its reformulation for fast calculation in real-time applications. Section 5 shows the advantage of the robust MPC over nominal MPC through simulation results. Section 6 gives a summary and suggests future research. 2. The industrial supply chain optimization problem Fig. 1 shows a simplified diagram of the multi-echelon supply chain system. In this system, unlimited raw materials are processed in plant IPM to produce the intermediate product (denoted by IP) at rate P (unit IP/84 h). The IP is stored in the storage IPS with inventory I1 (unit IP) and then processed at a rate F2 (unit IP/h) in plant SKUM into the final product, called Stock Keeping Unit (denoted by SKU). SKU is then sent to the plant distribution center DC with inventory I2 (unit SKU). SKU is shipped at a rate F4 (unit SKU/12 h) to the regional distribution center RDC with inventory I3 (unit SKU). The unit cost of SKU shipment is constant, because if the SKUs do not fill up a truck, other products can fill in. The average customer demands D (unit SKU/h) to RDC can be estimated from historical data, but a large variability exists in demands. If a demand cannot be satisfied immediately, a stock-out occurs and the unfilled part of the order is recorded as a back order O (unit SKU), which must be satisfied by later shipments before new demands will be satisfied. In all cases considered, plants IPM and SKUM have sufficient capacity to satisfy the total customer demands over the time horizon. Also, the maximum storage capacities are unlimited. The following periods between control actions have been taken from the industrial application and could be further investigated. The decision on IP manufacturing is made once every 84 h (3.5 days). The decision on SKU manufacturing is made once every day. The decision on SKU shipment is made once every 12 h. The goal of the supply chain optimization is to minimize the system cost of the supply chain (including inventory cost, manufacturing cost and transportation cost) while satisfying customer

2136

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

Fig. 1. The simplified diagram of the industrial supply chain system.

demands (if possible) by making decisions on the IP manufacturing rate P, SKU manufacturing machine running time Ts (h) as well as the SKU transportation quantity F4 . The uncertainties in the system include the SKU manufacturing rate Rs (unit SKU/h), the product transportation time  (h) and the customer demands D. In this paper, we only consider one type of material and product and one regional distribution center in the supply chain system, but the method can address multiple materials, products, and regional distribution centers. We make the following assumptions based on the real industrial problem and needs for real-time computing: (1) In each manufacturing decision interval, the IP or SKU manufacturing is continuous and the manufacturing rate is constant. The production scheduling problems (if needed) are assumed to be solved locally, which is out of the scope of this supply chain optimization. This decomposition is typical (Pinedo, 2000). (2) The SKU shipments to RDC only occur at predetermined times each day. (3) The daily quantities of customer demands are assumed to be continuous, and the demand rate is assumed to be constant within each day. (4) Fractional numbers are allowed in the solution of the supply chain optimizer and are rounded to an integer for implementation in practice. (Note that the plant simulation used in this study to evaluate closed-loop supply chain performance enforces all integer values where required.)

include robust performance. Here, nominal MPC means the classical MPC method that uses a deterministic model to predict the future system behavior, the deterministic prediction model is called nominal model. The nominal MPC is composed of three parts: the nominal state-space dynamic model of the supply chain system, the economic objective function and the constraints on the variables. 3.1. The nominal state-space dynamic model We develop the nominal model in the state-space form in the discrete time domain. The time interval of the discrete time model is constant and chosen to be the smallest time between optimizer decisions. As shown in Fig. 2, decisions are implemented at three different periods, all integer multiples of 12 h. Therefore, we set the time interval, T, to be 12 h. Also, all the inventories in this supply chain system are measured once per 24 h, and we let the supply chain optimizer execute only when new measured information on inventories are available. So supply chain optimizer execution period, Tc , is 24 h. The nominal model can be built based on the mass balance of each inventory in the supply chain system. 3.1.1. IPS inventory I1,k+1 = I1,k + F1,k T − F2,k T ,

where k in the subscript denotes the sequence number of sampled time steps. The IP flow out of plant IPM, F1,k , is determined by the IP manufacturing decision P with,

3. The nominal MPC formulation F1,k = This section gives the detailed information on the nominal MPC formulation of this supply chain optimization problem; the nominal MPC will be used for comparison with the robust MPC and the nominal formulation will be extended in the next section to

(1)

Pm1 , TP

(2)

where Pm1 is the decision variable determined at time step m1 that produces F1,k ; TP (=84 h) is the decision time interval of P. The IP flow sent to plant SKUM, F2,k , is determined by the SKU manufac-

Fig. 2. Relationship between discrete model time interval and the decision time intervals.

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

turing decision Ts,m2 with,

3.3. The constraints on the variables

CIP –SKU Rs,k Ts,m2 , Tc

F2,k =

(3)

where Ts,m2 denotes the decision variable that determines the SKU manufacturing time at time step m2 that producesF2,k ; Rs,k is the (uncertain) production rate at time step k; CIP-SKU converts the units of intermediate product to SKU. 3.1.2. DC inventory I2,k+1 = I2,k + F3,k T − F4,k T,

(4)

where the SKU flow F3,k is equal to the IP flow F2,k because the manufacturing equipment has no inventory, so F2,k Rs,k Ts,m2 = CIP –SKU Tc

F3,k =

2137

First, all the inventories and the SKUs in transportation should be nonnegative, so I1,k+1 , I2,k+1 , I3,k+1 , S ≥ 0

(14)

According to the definition of xk , (12) can be written into xk+1 ≥ xmin − Bo Ok+1

(15)

where the diagonal matrix

⎡

⎤

0

⎢ ⎣

..

Bo = ⎢

⎥ ⎥. ⎦

. 0

(5)

1 Second, all the back orders should be nonnegative, so,

3.1.3. RDC inventory I3,k+1 = I3,k + F5,k T − F6,k T

(6)

Ok+1 ≥ 0

Ok+1 = Ok + Dk T − F6,k T

(7)

Third, the decision variables are all nonnegative, and they are also subject to upper limits.

∗ =I Eq. (7) denotes the back order balance. Let I3,k 3,k − Ok , then Eqs. (6)–(7) can be written into ∗ ∗ I3,k+1 = I3,k + F5,k − Dk

(8)

∗ where I3,k+1 can be negative (when back orders exist), zero, or positive. Also, the arrival shipment to RDC, F5,k , is from the plant distribution center DC, so

F5,k = F4,k−(/T )

(9)

The delayed time between the departure shipment F4 and arrival shipment F5 in Eq. (9) is caused by transportation time . We can express Eq. (9) in another form by introducing additional variables S = (S1 ,. . .,S/T )T that denote the SKU in the transportation during a time step:

⎛

F5,k s1,k .. .

⎞

⎛

⎞

s1,k−1 ⎜ ⎟ ⎜ s2,k−1 ⎟ ⎜ ⎟=⎜ . ⎟ ⎝ ⎠ ⎝ .. ⎠ s/T,k F4,k

(10)

Eqs. (1)–(5), (8) and (10) can be combined into the following state-space model xk+1 = Axk + Buk + Bd dk

(11)

T (S T , I1 , I2 , I3∗ )

where x = contains the state variables, u = (P, Ts , F4 )T contains all decision variables and dk = Dk denotes the disturbances. 3.2. The economic objective function According to the goal of the supply chain optimization, the economic objective of the nominal MPC can be written as min

uk ,Ok

 k

T + CS,k+1 Sk+1 + CO,k+1 Ok+1 ]

(12)

where C1,k+1 contains the inventory costs, Cu,k contains the production and transportation costs, CS,k contains the inventory cost of SKUs in transportation, CO,k denotes the penalty on the backorders. According to the definition of the state vector xk , the objective function can be also written in the following form min

uk ,Ok

k

(17)

Ts,k ≤ Ts,max

(18)

F4,k ≤ F4,max

(19)

where the bound (18) denotes the total manufacturing time of a day cannot exceed Ts,max (24 h), and the bound (19) denotes the quantity of SKU in a shipment cannot exceed the capacity of the trucks available. Since IP manufacturing capacity is unlimited, there is no upper bound on P, F1 . However, since the IP manufacturing decision P is made once every 84 h and not reconsidered until 84 h have elapsed, F1 in the subsequent several time steps must be set to the values determined every 84 h. Therefore, we need to force these F1 to be the value that has been determined in a previous optimization, which can be realized by forcing the upper and lower bounds on these F1 to be the determined F1 value. Therefore, the F1 variables are optimized less frequently than other variables, but they can be nonzero between the 84 h based on the values determined when they are allowed to be degrees of freedom in the optimization. So these bounds are always active and the way to handle them can be found in Section 4. All of these the bounds on the decision variables can be written in the following form: umin ≤ uk ≤ umax

T Cx,k+1 xk+1

+

 k

T Cu,k uk

+

 k

CIO,k+1 Ok+1

(13)

(20)

with the minimum and maximum values determined prior to each optimization as required to achieve the above discussed strategy. Therefore, the objective function (13), the nominal model (11) and the constraints (15), (16) and (20) construct the nominal MPC formulation, min

s.t.

T + Cu,k uk



uk ≥ 0

uk ,Ok

T [CI,k+1 (I1,k+1 , I2,k+1 , I3,k+1 )T

(16)



T Cx,k+1 xk+1 +

k



T Cu,k uk +

k

xk+1 = Axk + Buk + Bd dk + e0



T CIo,k+1 Ok+1

(21a)

k

(21b)

umin ≤ uk ≤ umax

(21c)

xk+1 ≥ xmin − Bo Ok+1

(21d)

Ok+1 ≥ 0

(21e)

k = 0, . . . , n − 1 where n denotes the same horizon for both the decision and state variables (i.e., the same control and prediction horizons), the

2138

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

Fig. 3. Approximate disjunctive model for transportation time uncertainty.

additional variable e0 in Eq. (21b) denotes the feedback information which contains the difference between the predicted and the measured inventories. Note that for this supply chain system, we assume that all the inventories can be measured at each time step, but the SKUs in the transportation (i.e., S in state vector x) cannot. Therefore, we set the feedback elements in e0 , which correspond to S, to be zero at each time step. 4. The robust MPC formulation 4.1. Description of the uncertainties with uncertain parameters The new robust MPC formulation addresses the uncertainties by using a prediction model that has the same structure as the nominal model (21b) but uncertain parameters. The uncertainty in customer demand D can be directly described by the uncertainty in parameter dk . The uncertainty in SKU manufacturing rate Rs can be described by the uncertain elements in matrix B that depend on Rs . The uncertainty in SKU transportation time, however, makes the structure of the model uncertain. Our solution method accounts for parametric uncertainty; therefore, we propose to use an alternative disjunctive model with uncertain parameters to approximate this structural uncertainty. The concept is illustrated in Fig. 3. The shipment from DC is assumed to reach the RDC through several (in this work three) different virtual routes with different but deterministic transportation times,  n (nominal transportation time),  max (maximum transportation time) and  min (minimum transportation time). Then, we have F4,j = j F4 ,

0 ≤ j ≤ 1,



j = 1

(22)

j

F5,j,k = F4,j,k−(j /T )

(23)

F5 =

(24)



troller on the closed-loop behavior (e.g. Kothare et al., 1996; Goulart et al., 2006). Generally, we can address the effect of the controller on the closed-loop behavior by including a controller model that accounts for the future MPC actions, which results in the closed-loop robust MPC formulation, as shown in Fig. 4(b). In the closed-loop formulation, both the (uncertain) process dynamic model and the controller model are solved at every time step in the future horizon, so the future controlled and manipulated variables (uk ) are uncertain. The new robust MPC method approximates the future MPC algorithm (which is actually the robust MPC algorithm itself) with a nominal MPC algorithm. This nominal MPC optimizes its objective function by adjusting the trajectory of manipulated variables. As a result, the robust MPC is given in the following formulation. min



x˜ r,k+1 ,u˜ r,k

s.t.

T Cx,k+1 xk+1 +

k



T Cu,k uk +

k



T CIo,k+1 Ok+1

(25a)

k

up,k = NMPC(xp,k , ep,k , d˜ k, x˜ r,k+1 , u˜ r,k )

(25b)

xp,k+1 = Axp,k + Bp,k up,k + Bd dp,k + e0

(25c)

ep,k+1 = xp,k+1 − (Axp,k + Bup,k + Bd,k dk )

(25d)

umin ≤ up,k ≤ umax

(25e)

xmin − Bo Ok+1 ≤ xp,k+1

(25f)

Ok+1 ≥ 0

(25g)

For all Bd,k , dp,k in uncertainty region, k = 0, . . . , n − 1 In this formulation, the objective of the optimization is to minimize the nominal economic objective function. The values of the manipulated variables are determined by the simulated, nominal controller; therefore, the degrees of freedom for optimization T

F5,j

T T , . . . , xr,k+n ) and the become the reference states x˜ r,k+1 = (xr,k+1

j

where Eq. (23) describes the time delay between the departure shipment F4,j and arrival shipment F5,j due to the transportation time  j , and it can be transformed into the form of Eq. (10) with additional variables. The ratios  j are uncertain. So, the uncertainty in SKU transportation time can be approximately described by the uncertain elements in B that depend on  j . 4.2. Optimization of the uncertain closed-loop dynamics With the presence of the uncertainties, the closed-loop behavior of the system in the future is uncertain. A simple (but inaccurate) way to describe the future uncertain closed-loop dynamics is to enhance the dynamic model (21b) by including uncertain parameters for the future time steps and solve for the future decisions uk . This “open-loop” approach to robust MPC, shown in Fig. 4(a), is not satisfactory because it omits the future correcting actions by the MPC itself, and these future controller actions will partially compensate for the model uncertainty and unmeasured disturbances. Therefore, open-loop robust MPC may give very conservative control and poor dynamic performance. As is widely recognized, robust MPC methods must address the effect of the con-

Fig. 4. (a) Open-loop and (b) closed-loop robust MPC.

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143 T

reference decisions u˜ r,k = (uTr,k , . . . , uTr,k+n−1 ) in the future horizon (time step k = 0,. . .,n − 1). Eq. (25b) means the decision variables at the future kth time step are determined by the nominal MPC control T law NMPC(xp,k , ep,k , d˜ k, x˜ r,k+1 , u˜ r,k ), where d˜ k = (dT , . . . , dT ) k

k+n−1

contains the forecast future disturbances from the kth to the (k + n − 1)th time step. We use up,k to represent the uncertain decision variables in the future kth time step so as to differentiate them from the their nominal value uk . Similarly, xp,k , ep,k denote the uncertain state and feedback variables in the future kth time step. Eqs. (25c)–(25d) denote the uncertain process model and Bp,k , dp,k denote the uncertain parameters of the model at time step k. The simulated nominal MPC solves the following Quadratic Program (QP) at the future kth time step, min ui



T

(xi+1 − xr,i+1 ) Q (xi+1 − xr,i+1 ) +

i



(26a)

s.t. xi+1 = Axi + Bui + Bd di + ep,k

(26b)

umin ≤ ui ≤ umax

(26c)

i = k, . . . , k + n − 1 where Q, R are the weighting matrices for state and decision variables respectively. Since Q, R are usually selected to be diagonal matrices with nonzero diagonal elements, the QP problems (26a)–(26c) are strictly convex. Note that for the optimal solution to QP (26), only the decisions at the current time step (uk ) will be considered to be implemented in the closed-loop model. In the nominal controller, we must consider the bounds (26c) on the decision variables explicitly, because these bounds are “hard bounds” which can never be violated during implementation due to strict limits in the real problem. However, we do not include the bounds on the state variables explicitly because they are addressed with the constraints (25f) when solving the robust MPC formulation (25). As shown in Fig. 4, the robust MPC formulation (25) is a bilevel stochastic optimization problem, which is very difficult to solve especially in real time. The next several subsections discuss the reformulation of this difficult problem so that solutions can be obtained for real-time applications. 4.3. Reformulation with known active bounds To reformulate the bi-level problem (25), we make the following assumption: A bound for a manipulated variable at each time step is either active or inactive for all the realizations of the system. With this assumption, the bi-level problem (25) can be transformed into a single level problem using a method introduced by Clark and Westerberg (1983). We can remove all the inactive bounds, which will not affect the solution, and transform all the active bounds into equality constraints. The solution of a convex QP problem with equality constraints can be obtained by solving its first order optimality condition (Nocedal & Wright, 1999), from which up,k can be expressed as a linear function of xp,k , ep,k , d˜ k, x˜ r,k+1 , u˜ r,k , say, up,k = Kx xp,k + Ke ep,k + Kd d˜ k + Kxr x˜ r,k+1 + Kur u˜ r,k

(27)

where Kx , Ke , Kd , Kxr , Kur are constant matrices. Details of these matrices are not shown here due to the space limitations; see Li (2009) for more details. Let us define new variables tk that combine all decision variables into one vector. tk = Kxr x˜ r,k+1 + Kur u˜ r,k

then Eq. (27) becomes up,k = Kx xp,k + Ke ep,k + Kd d˜ k + tk

(29)

According to the control law (29) and the process model (25c)–(25d), we can write the closed-loop model of the system into

⎛

up,k = Lup,k ⎝

⎞

t0 .. . tn−1

⎛

⎞

t0 .. .

xp,k+1 = Lxp,k ⎝



⎠ + Mup,k

x0 e0 dˆ



⎠ + Mxp,k

tn−1

 + Nu,k d˜ p

x0 e0 dˆ

(30)

 + Nx,k d˜ p

(31)

T

(ui − ur,i ) R(ui − ur,i )

i

2139

(28)

k = 0, . . . , n − 1. T T ) contains forecast future disturbances, where dˆ = (d0T , . . . , d2n−1 Lup , Mxup , Lxp , Mxp , and d˜ p contains uncertain (and certain) param-

eters. Again, the details of the matrices Lup , Mxup , Nup , Lxp , Mxp , Nxp are not shown here due to space limitations; see Li (2009) for more details. With closed-loop model (30)–(31), we can transform robust MPC formulation (25) into min

tk ,Ok



T Cx,k+1 xk+1 +

k

k

⎛

t0 .. .

s.t. up,k = Lup,k ⎝

⎛ xp,k+1 = Lxp,k ⎝



T Cu,k uk +



⎞



⎠ + Mup,k

⎞

⎠ + Mxp,k

tn−1

(32a)

k

tn−1 t0 .. .

T CIo,k+1 Ok+1



x0 e0 dˆ

x0 e0 dˆ

 + Nu,k d˜ p

(32b)

 + Nx,k d˜ p

(32c)

umin ≤ up,k ≤ umax

(32d)

xmin − Bo Ok+1 ≤ xp,k+1

(32e)

Ok+1 ≥ 0

(32f)

For all Lup , Mxup , Lxp , Mxp , d˜ p in uncertainty region, k = 0, . . . , n − 1 where the degrees of freedom (˜xr,k+1 , u˜ r,k ) appear in tk . This problem is a single level stochastic Linear Program (LP), and its solution is explained in Section 4.5. The next section explains how the active set of the original problem (25) is determined. 4.4. The “DMC” heuristic to determine active bounds We use a heuristic to obtain the active bounds on the decision variables in an iterative way, which includes the following three steps. (1) Assume no bounds are active in the future and solve problem (32). (2) Problem (32) is solved. If one (or more) decision bound, which is not already set active, is active, go to step (3); otherwise, end the iterative procedure. (3) Select the decision variables that reach their bounds at the earliest time step to their bound value by adding appropriate equality constraints to the “inner” optimization problem (26). Then, we get a new closed-loop model (30)–(31) and so that a new problem (32). Go to step (2).

2140

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

Fig. 5. The heuristic used to obtain the active bounds on decision variables.

Note that during the above iterative procedure, any decision variable at any time step that has been set to its bound value will be so constrained for the remainder of the iteration. Fig. 5 illustrates the heuristic. At the first iteration, the problem is solved with the assumption that no bounds will be active in the future. However, the solution gives the uncertain trajectory of the manipulated variable whose boundary at the 3rd time step encounters its upper bound. According to the active set heuristic, the manipulated variable at the 3rd time step is set to its upper bound, and the problem is solved again. The procedure is repeated until the manipulated variables at all time steps are forced to their bounds or are unbounded at the solution. A similar heuristic has been successfully applied for constrained Dynamic Matrix Control (DMC) (Prett & Gillette, 1979) in industry for more than 20 years. This “DMC heuristic” addresses the hard bounds on the manipulated variables iteratively so that at every iteration only a least square problem is to be solved in the MPC calculation. Due to its success in the deterministic MPC formulation along with good results to be shown, we believe the idea of the heuristic is appropriate for the stochastic MPC formulation introduced here. Note that the number of iterations in the heuristic is only proportional to the length of the input horizon, which is a good property for real-time application. However, the heuristic does not guarantee the “global optimum” of the solution, i.e. there may be another active set that is better than the one found by the heuristic. If the heuristic recovers the correct active set, the solution will be the global optimum. 4.5. Solution with the chance-constrained program The stochastic LP problem (32) contains linear inequalities with uncertain parameters, which can be transformed into deterministic constraints. Let us consider an upper bound of an uncertain decision variable, which can be written (according to (32b) and (32d)) as,

⎛

Lup,l ⎝

t0 .. . tn−1

⎞

⎠ + Mup,l



x0 e0 dˆ

 + Nu,l d˜ p ≤ umax,l

(33)

where the subscript l denotes the lth upper bounds on the decision variables. The idea of a chance-constrained program is to guarantee the feasibility of the constraint at the given confidence level ˛, i.e.

transform bound (33) into the following form,

⎛

⎛

t0 .. .

Pr ⎝Lup,l ⎝

⎞



⎠ + Mup,l

tn−1

x0 e0 dˆ

⎞



+ Nu,l d˜ p ≤ umax,l ⎠ ≥ ˛

(34)

If we assume Lup,l , Mup,l and d˜ p obey normal distribution, then the constraint (34) is equivalent to the following deterministic constraint (Lobo, Vandenberghe, Boyd, & Lebret, 1998):

⎛ E(Lup,l ) ⎝

t0 .. . tn−1

⎞ ⎠ + E(Mup,l )



x0 e0 dˆ

 + Nu,l E(d˜ p ) − umax,l + ˚−1 (˛)

  T  1/2 T ˆ Nu,l , 1)  , x0T , e0T , d, × Vu,l (t0T , ..., tn−1  ≤0

(35)

2

where E(·) denotes the expected value of the parameters in the brackets, ˚−1 (˛) denotes the inverse cumulative probability function of normal distribution, Vu,l denotes the covariance matrix of (Lup,l , Mup,l , d˜ pT ), which reflects the effects of the correlated uncertain parameters on the closed-loop system behavior. Vu,l is obtained through Monte Carlo sampling as follows: (1) Randomly select a sample of the uncertain parameters of the open-loop system. Rs and  are assumed to be normally distributed with the known variance from industrial data (see Table 1 in Section 5). Customer demand sample is selected from historical data directly (which gives d˜ p ). (2) Calculate Lup,l , Mup,l using the sample of Rs , , d˜ p and other certain parameters. (3) Repeat procedures (1–2) for 100 samples of the uncertain parameters and obtain different groups of Lup,l , Mup,l , d˜ p , which are then be used to calculate Vu,l according to the standard technique (Box, Jenkins, & Reinsel, 2008). Substantial computing for the Monte Carlo sampling is performed off-line as part of the controller design and tuning, and therefore, it does not affect the tractability of the real-time solution. The deterministic constraint (35) (for ˛ > 0.5), called second order conic constraint, is a convex inequality. If we transform all the uncertain linear inequalities in the same way, the problem (32) will become a deterministic Second Order Cone Program (SOCP)

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

2141

Table 1 The value of the parameters used in the simulation study. Parameter

Value

Parameter

Value

Control and prediction horizon n (day) Discrete model time interval T (day) Measurement and controller execution period Tc (day) Nominal Rs (unit SKU/h) Rs range with 90% confidence (unit SKU/h) Nominal  (h)  range with 90% confidence (h) D range with 90% confidencea (unit SKU/h)

14 0.5 1 16.7 13.3–22.2 144 132–156 0–1.6

SKU Inventory cost, in CI ($/unit SKU) IP Inventory cost, in CI ($/unit IP) SKU Transportation cost, in Cu ($/unit SKU) IP manufacturing cost, in Cu ($/unit IP) SKU manufacturing cost, in Cu ($/unit SKU) Back order cost CO ($/unit SKU) SKU transportation capacity F4,max (unit SKI/12 h) Confidence of each chance constraint ˛

1.2 7.2 0.01 1 0.1 1200 40 90%

a

Historical data of day-to-day customer demand D is known and samples from historical data are used to characterize the uncertainty in D.

(Ben-Tal & Nemirovski, 1999) that can be solved efficiently and reliably with the state-of-the-art interior point optimizer, e.g. CPLEX or IPOPT (Wächter & Biegler, 2006). By combining the active set heuristic with the conic formulation of the stochastic program, we can solve the original robust MPC formulation (25) (approximately, because of the active set heuristic and the confidence level ˛ in the probabilistic reformulation) by solving a limited series of SOCPs, which is tractable for real-time implementation. 5. Simulation results The robust MPC has been applied to the real multi-echelon industrial supply chain problem introduced in Section 2. Table 1

shows the values of the parameters used in the simulation study. Two scenario cases were considered: (Case a) the forecast prediction of customer demand, SKU manufacturing rate and SKU transportation time are exactly correct and (Case b) the forecast prediction of customer demand per day is the average demand during the 28 days which is different from the actual demand. Also, in Case b the nominal SKU manufacturing rates and SKU transportation times differ from their nominal values at each time period in a random manner; the sampled values could be outside of the 90% limits reported in Table 1. We emphasize that Case a is a very idealized scenario because it assumes that the demand forecast is perfect for the controller’s entire prediction horizon at every execution.

Fig. 6. The supply chain simulation results with the nominal MPC. (a) Nominal Rs , t, D are correct. (b) Nominal Rs , t, D are incorrect.

2142

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143

Fig. 7. The supply chain simulation results with the robust MPC. (a) Nominal Rs , t, D are correct. (b) Nominal Rs , t, D are incorrect.

To establish the performance for a standard technology, the nominal MPC method (described by formulation (21)) was used to optimize the supply chain system. In addition, the same scenario (with the same realizations of random variables) were optimized in closed-loop using the robust MPC method described by formulation (32). Fig. 6 shows the simulation results with the nominal MPC in the two cases. We observe in Fig. 6(a) that when the nominal prediction of the uncertain parameters is correct, the nominal MPC performs very well—it controls the inventory near zero while satisfying customer demands essentially all the time. (Note that the small back orders occurring around the 28th day were caused by the modeling approximations we made when building the nominal model.) However, Fig. 6(b) shows that when model and forecast mismatch occurs, the nominal MPC performed unsatisfactory because of the large numbers of back orders occurring after the 20th day. The nominal controller reduces inventories and does not have a “safety stock” to enable it to respond well to mismatch. Fig. 7 presents the simulation results with the robust MPC in the two cases, in which we can find there were no back orders in either case. In Case a without mismatch the robust controller maintains more inventory than absolutely required because it is maintaining safety stock for realizations that include higher demands and longer manufacturing and transportation times than their expected values. In Case b, the safety stock prevented backorders.

We emphasize that the robust MPC obtains the minimum safety stock through optimization according to the known uncertainty characterization. Therefore, the safety stock is the minimum for the uncertainty and scenarios experienced, which is better than setting constant safety stocks based on past experience. The significance level of the probability of constraint violation can be tuned by the engineer to reflect the importance of backorders versus the cost of inventory. The figures represent the behavior of nominal and robust MPC for a specific mismatch realization. The performance for a wider range of realizations is summarized in Table 2. To generate the results for Table 2, we ran the simulation with nominal and robust MPC for eight different cases, in which the parameters Rs , , D are certain or uncertain. The results in the table report the average behavior for 30 realizations of the uncertain parameters. The uncertainty realizations were selected randomly from their probability distributions and therefore, were allowed to exceed the ranges used when designing the controller. The results in Table 2 clearly demonstrate the reduced occurrence of backorders when the plant is controlled by the robust MPC. Some small number of backorders occurred with robust MPC because of (few) realizations with demand variation larger than the design distribution. We could increase the significance level to reduce the potential back orders; but then we have to maintain more “safety stock”. In contrast, the average back order was nonzero

X. Li, T.E. Marlin / Computers and Chemical Engineering 33 (2009) 2134–2143 Table 2 Summary of the simulation results with nominal and robust MPC in the eight cases. Case number

Rs



D

1 2 3 4 5 6 7 8

C C C C U U U U

C C U U C C U U

C U C U C U C U

Average  O (unit SKU) of 30 samples Nominal MPC

Robust MPC

2.2 527.0 32.8 582.8 2.3 530.2 33.0 599.6

0 54.9 0 65.0 0 55.5 0 65.9

Note: C = certain; U = uncertain. (Note that the small back order value in case 1 with nominal MPC was due to the modeling approximations.)

in all the seven mismatch cases when nominal MPC was used, and the increase in backorders was substantial in most cases. Table 2 also shows that the uncertainty in customer demand D has the dominant effect on the system performance, and this uncertainty must be addressed in the MPC calculation to reduce the back orders. In addition, uncertainties in multiple parameters increase the number of backorders. While the effects of the uncertainties in the SKU manufacturing rate Rs and transportation time  are less significant, addressing them in the robust MPC calculation is still essential to reduce the back orders. Since Rs ,  appear in the feedback model and D is the disturbance to the system, both the feedback model mismatch and the disturbance uncertainty must be addressed to achieve the performance for the robust MPC in Table 2. The simulation was performed on a PC with Intel Core 2 Duo, 4GB memory and Windows Vista. The solution for the plant simulation is programmed in MATLAB and the controller SOCP optimization problems are solved in GAMS with the interior point (barrier) solver of CPLEX 10.2.0. Each SOCP problem has 70 decision variables, 204 linear constraints and 116 second order cones. The CPU time for solving the SOCP problems for one controller execution is less than 2 s. The method can be extended to multiple materials (types of products) and multiple regional distribution centers (Li, 2009). The number of variables and computer memory scales linearly with the problem size, and the SOCP solution time is polynomial in the size; for example, extension to five material types in the supply chain increases the computing time to about 1 min. 6. Summary and future work This paper presents a novel robust MPC method that addresses models with uncertainty in both the feedback process and the disturbance forecast, with correlated uncertainty descriptions. The robust MPC handles hard bounds on decision variables through a heuristic, and it only requires solving a limited number of second order cone programs online. The simulation study shows that the robust MPC can determine the optimal safety stock with the known information on uncertainties, which is a key advantage of the robust MPC method over nominal MPC for supply chain optimization. The simulation study also shows the ability of the robust MPC to address both the model mismatch and disturbance uncertainty for this supply chain optimization problem. The robust MPC formulation provides robust performance by preventing constraint violations. Approximations implemented to yield tractable real-time computation prevent a guarantee of optimality. Neither nominal nor robust stability is mathematically

2143

guaranteed, but nominal stability could be included by adding stability constraints (e.g., Mayne, Rawlings, Rao, & Scokaert, 2000). Although not shown here, the method can solve larger problems with more distributions centers and multiple products. In addition, it can be applied to process control applications where significant model mismatch exists and maintaining feasibility during transients is critical. Acknowledgements This research was supported by the McMaster Advanced Control Consortium (MACC). Collaboration of an industrial sponsor was of great assistance. References Ben-Tal, A., & Nemirovski, A. (1999). Robust solutions of uncertain linear programs. OR Letters, 25, 1–13. Box, G. E. P., Jenkins, G. M., & Reinsel, G. C. (2008). Time series analysis: Forecasting and control. Hoboken, NJ: John Wiley. Clark, P. A., & Westerberg, A. W. (1983). Optimization for design problems having more than one objective. Computers and Chemical Engineering, 7, 259–278. Goulart, P. J., Kerrigan, E. C., & Maciejowski, J. M. (2006). Optimization over state feedback policies for robust control with constraints. Automatica, 42, 523–533. Kothare, M. V., Balakrishnan, V., & Morari, M. (1996). Robust constrained model predictive control using linear matrix inequalities. Automatica, 32, 1361–1379. Lee, J. M., & Lee, J. H. (2004). Approximate dynamic programming strategies and their applicability for process control: A review and future directions. International Journal of Process Control, Automation and Systems, 2, 263–278. Li, X. (2009). Robust model predictive control for process control and supply chain optimization. Ph.D. Thesis. McMaster University, in press. Lobo, M. S., Vandenberghe, L., Boyd, S., & Lebret, H. (1998). Applications of secondorder cone programming. Linear Algebra & Its Applications, 284, 193. Maciejowski, J. M. (2002). Predictive control: With constraints. England: Prentice Hall. Mayne, D. Q., Rawlings, J. B., Rao, C. V., & Scokaert, P. O. M. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36, 789–814. Nocedal, J., & Wright, S. J. (1999). Numerical optimization. New York: Springer. Perea-Lopez, E., Ydstie, B. E., & Grossmann, I. (2003). A model predictive control strategy for supply chain management. Computers and Chemical Engineering, 27, 1201–1218. Pinedo, M. L. (2000). Planning and scheduling in manufacturing and services. New York: Springer. Powell, W. B., & Van Roy, B. (2004). Approximate dynamic programming for high dimensional resource allocation problems. In J. Si, A. Barto, W. B. Powell, & D. Wunsch (Eds.), Learning and approximate dynamic programming: scaling up to the real world. New York: Wiley. Prett, D. M., & Gillette, R. D. (1979). Optimization and constrained multivariable control of a catalytic cracking unit. In AIChE National Meeting Houston, TX. Qin, S. J., & Badgwell, T. A. (2003). A survey of industrial model predictive control technology. Control Engineering Practice, 11, 733–764. Sarimveis, H., Patrinos, P., Tarantilis, C. D., & Kiranoudis, C. T. (2008). Dynamic modeling and control of supply chain systems: A review. Computers and Operations Research, 35, 3530–3561. Simchi-Levi, D., Kaminski, P., & Simchi-Levi, E. (1999). Designing and managing the supply chain: Concepts strategies and case studies. New York: McGraw-Hill/Irwin. Shah, N. (2005). Process industry supply chains: Advances and challenges. Computers and Chemical Engineering, 29, 1225–1235. Tzafestas, S., Kapsiotis, G., & Kyriannakis, E. (1997). Model-based predictive control for generalized production planning problems. Computers in Industry, 34, 201–210. Van Hessem, D. H. (2004). Closed-loop model predictive control with application to stochastic model-based chemical process operations. Ph.D. Thesis. Delft University of Technology. Warren, A. L., & Marlin, T. E. (2003). Improved output constraint-handling for MPC with disturbance uncertainty. In Proceeding of the American Control Conference (pp. 4573–4578). Wächter, A., & Biegler, L. T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106, 25. Wang, W., Rivera, D. E., & Kempf, K. G. (2007). Model predictive control strategies for supply chain management in semiconductor manufacturing. International Journal of Production Economics, 107, 56–77.