Optimisation approaches for supply chain planning and scheduling under demand uncertainty

This work presents efficient MILP-based approaches for the planning and scheduling of multiproduct multistage continuous plants with sequence-dependent changeovers in a supply chain network under demand uncertainty and price elasticity of demand. This problem considers multiproduct plants, where several products must be produced and delivered to supply the distribution centres (DCs), while DCs are in charge of storing and delivering these products to the final markets to be sold. A hybrid discrete/continuous model is proposed for this problem by using the ideas of the Travelling Salesman Problem (TSP) and global precedence representation. In order to deal with the uncertainty, we proposed a Hierarchical Model Predictive Control (HMPC) approach for this particular problem. Despite of its efficiency, the final solution reported still could be far from the global optimum. Due to this, Local Search (LS) algorithms are developed to improve the solution of HMPC by rescheduling successive products in the current schedule. The effectiveness of the proposed solution techniques is demonstrated by solving a large-scale instance and comparing the solution with the original MPC and a classic Cutting Plane approach adapted for this work. © 2018 The Authors. Published by Elsevier B.V. on behalf of Institution of Chemical Engineers. This is an open access article under the CC BY license (http://creativecommons. org/licenses/by/4.0/).


Introduction
In this paper, we address the integration of planning and scheduling on multiproduct multistage continuous plants in a supply chain (SC) network with price elasticity of demands under uncertainty. For this, decisions at tactical and operational levels should be done together to improve the efficiency of the supply chain network.
Tactical decisions in SC are focused on different echelon nodes including plants, distribution centres (DCs), markets. Here, planning decisions have to also consider the interrelation between different echelon nodes in a complex SC network. These include assignment and transportation decisions on products; for example from which plant DCs must be supplied and from which DCs individual markets are served at each time period. Other important issue is the sequence-dependent changeovers along consecutive time periods.
This limitation forces the system to take into account the previous limitation usually becomes more difficult when the resolution of the problem increases, due to the increasing number of decisions that have to be made. Another issue relies in the combinatorial complexity when many products must be sequenced simultaneously in the system. Finally, the assignment of products to processing units could increase the problem complexity in terms of decision variables.
In this particular problem, multiproduct multistage continuous plants, with a single production unit per stage where sequencedependent changeovers occur between different products in each processing unit, are considered.
Common applications in multiproduct multistage batch plants with a single-unit stage can be found in automated manufacturing systems (Geiger et al., 1997;Karimi, 2003, 2004;Karimi et al., 2004;Aguirre et al., 2011Aguirre et al., , 2012Castro et al., 2011bCastro et al., , 2012Zeballos et al., 2011;Novas and Henning, 2012), surface-treatment process of metal components in aircraft manufacturing plants (Paul et al., 2007), hoist scheduling problem in electronics industry (Phillips and Unger, 1976;Shapiro and Nuttle, 1988;Crama 1997) and many other areas like car-part manufacturing (Lu et al., 2017), health-care applications (Rahman et al., 2015), clothing industry (Chang et al., 2015). For continuous plants, Grossmann (1994, 1998) addressed the cycle scheduling problem of multiple products in continuous plants considering a sequence of stages with a single production line. Then, Pinto (2002, 2004) and  proposed different formulations for a similar problem considering a single-unit stage in a flowshop plant with intermediate storage policies. Some similar applications for hybrid flowshop plants with dedicated units can be found in Birewar and Grossmann (1989); Moon et al. (1996); Munawar and Gudi (2005); Bhattacharya and Bose (2007) and Bhattacharya (2008, 2009), etc.
In this work, we are going to focus our attention on multiproduct multistage continuous plants in a supply chain network with demand under uncertainty. Efficient solution approaches that properly combine the integration of planning and scheduling decisions for the process industry in supply chain networks and also allow properly managing the uncertainty are still required. A brief literature review of main approaches in this area are presented below.

Literature review
There are some works that try to solve a similar problem by using different kind of solution methods (see Table 1). Bose and Pekny (2000) proposed a general Model Predictive Control (MPC) framework for planning and scheduling problems under uncertainty. Afterward, Perea et al., (2000) developed a dynamic approach for supply chain management mixing the ideas of process dynamics and control. Then, some of the authors, in Perea-Lopez et al. (2001), studied the decentralised control of supply chain systems by using classical control laws for the total cost, customer satisfaction and demand variations. Later on, they extend that work in Perea-Lopez et al. (2003) for the optimisation of a centralised supply chain by using a MPC and mixed-integer programming models. Braun et al. (2003) presented a robust and flexible decision framework for dynamically managing inventories and customer requirements in demand supply chain networks. Mestan et al. (2006) addressed the optimal operation of multiproduct supply chain systems by using MPC strategy and continuous/discrete dynamics and logic rules, for switching between different operating conditions. More recently, Kaplan et al. (2011) studied a centralised multi-period multiechelon multi product supply chain network with variations in demand and prices with the elasticity in price demand by using MINLP models.
Terrazas- Moreno et al. (2011) presented a Lagrangian method to study the behavior of the temporal and spatial decomposition in a multi-site multi-period production planning problem with sequence dependent changeovers. Liu et al. (2012) presented a Hierarchical Model Predictive Control (HMPC) approach proposed for complex planning and scheduling problems with demand under uncertainty. Then, Calfa et al. (2013) introduced a hybrid bi-level Lagrangian decomposition approach for tackling large-scale industrial case study. Subramanian et al. (2013) proposed a dynamic model and MPC approach to deal with multiproduct multi-echelon supply chain problem. Nasiri et al. (2014) developed a Lagrangian relaxation approach for the integrated supply chain plan-ning and distribution with stochastic demands. Muñoz et al. (2015) proposed a hierarchical-based strategy based on Lagrangian decomposition for the multi-period multisite supply chain planning and scheduling problem of multipurpose batch plant. More recently, Felfel et al. (2016) introduced a two stages stochastic linear programming approach for a multi-period, multi-product, multi-site, multi-stage supply chain planning problem under demand uncertainty, and Aqlan and Lam (2016) proposed a simulation-optimisation based tool for supply chain optimisation under risk and uncertainty.
Some other works for the integration of scheduling and control can be found in Zhuge and Ierapetritou (2014); Baldea and Harjunkoski (2014); Chu and You (2014); Fu et al. (2014Fu et al. ( , 2015, etc. Main approaches developed for the integration of production planning and transport can be seen in Mula et al. (2010) and Díaz-Madroñero et al. (2015). Strategic and operational planning approaches, including handling of uncertainties and multi-objective formulations, can be seen in Guillén et al. (2005); You and Grossmann (2008); Guillén-Gosálbez and Grossmann (2009);Papageorgiou (2009), etc. Finally, suitable reviews about integration of planning and scheduling in the process industry can be found in Maravelias and Sung (2009);Safaei et al. (2010); Shah and Ierapetritou (2012); Engell and Harjunkoski (2012).

Motivation
To the best of our knowledge, and despite of all previous works mentioned, there are limited robust solutions that properly solve this kind of problems by integrating constructive and improvement approaches in an efficient manner, especially for the supply chain planning problem with multiple production stages at plants under uncertainty.
In this work, we proposed a mixer-integer linear programming  Liu et al. (2012) to deal with the complex features of this problem including the uncertainty on product demand. The product demand is uncertain and follows a uniform distribution between specific upper and lower bounds as it was proposed in Liu et al. (2012). This well-known HMPC approach has been broadly used in the past for solving industrial size case studies in an iterative way. But, as many of the decision variables are successively fixed iteration by iteration, the final solution reported by this method could still be far from the global optimum. To overcome this issue, alternative solution techniques based on LS strategies, such as the one presented in Castro et al. (2011a), are proposed. LS strategies for the temporal, production and spatial decomposition are developed and implemented as an inner loop inside the HMPC in order to improve the solution obtained in an iterative way. Thus, for each iteration of the HMPC, a selected LS approach is applied by rescheduling successive products from the current schedule (Méndez et al., 2006). Thus, the proposed approach could be able to compensate the main drawback of the classic MPC approach and provide a better quality solution of the system for a realisation of uncertainty in shorter CPU time.
The remaining of this paper is structured as follows. In Section 2, we present the multiproduct multistage continuous plant planning and scheduling problem in a supply chain network study. Section 3 presents hybrid discrete/continuous time MILP formulations for this problem.
In order to solve the uncertainty, an efficient solution methodology is developed in Section 4, based on the main ideas of the HMPC and LS approaches. A Cutting Plane approach is also implemented for comparison. Computational analysis and results of a large-scale case study proposed are summarised in Section 5. Final conclusions and future work directions are provided at the end in Section 6.

Problem statement
Given a set of products i and echelon nodes n (plants s, distribution centres c and markets m) of the supply chain network and their possible connections and transportation times defined by T inn and inn , and given also the information about initial demands (D 0 int ), production rate (r iln ), production yields (y iln ), changeover times ( c ijln ) at each production stage l, minimum and maximum inventory levels (I min int , I max int ), initial selling prices (P 0 int ) and production cost (CP in ), changeover cost (CC ijn ) and unmet demand cost (CU in ) at each echelon node, we aim to find the total amount of products to produce (PR ilnt ) at each plant, to stores (I int ) at each echelon node and to sell (SV int ) or to unmet the demand (U int ) to each market at each time period t of the planning horizon H, as well as the price to charge (P int ), in order to maximise the total profit TP of the whole supply chain, as illustrated in Fig. 1 For this problem, important operational decisions related with the production assignment (L int , E int , F int ), sequencing (Z ijnt & ZF ijnt ) and timing (TS ilnt & T ilnt ) must be made in multiple stages l on a multiproduct continuous plant s to produce, store and deliver the products by the end of every time period t (see Fig. 2).
As can be seen in Fig. 2, E int represents the assignment of product i to plant n at time period t, while F int /L int , takes value 1 if this product is the first/last in the production sequence. The production sequence inside each time period is defined by Z ijnt which indicates whether product i is processed immediately before product j in the production sequence. This variable allows the execution of changeovers between different products inside each time period defined by ( c ijln ). According to this, initial processing time (TS ilnt ) and production time (PT ilnt ) in each stage l are defined to accomplish within bound limitations (LB, UB).
As the plants operate without any break, sequencing production decisions between consecutive time periods t-1 and t must be taken into account via binary variable ZF ijnt . This decision allows managing the production process of each stage l between consecutive time periods, adopting value 1 if product i is the last product of time period t-1 and product j is the first product at time period t.
All those planning decisions at each echelon node and scheduling decisions at each production plant result in a complex and challenging supply chain problem to solve. In the following section, different solution methodologies for this problem are discussed and analysed in detail.

Hybrid time formulations for Supply Chain Planning and Scheduling
In this work, we extend the previous MILP formulation developed by Liu et al. (2012) to consider multiple production stages in multiproduct continuous plants in a supply chain network. The proposed formulation is based on the main ideas of the travelling salesman problem (TSP) and the precedence-based relationships summarised in Méndez et al. (2006). In Fig. 3, we introduce the analogy from the TSP and the precedence-based ideas applied to our particular problem with the information of the main decision variables presented above.
If we do a parallelism between these two problems, we can interpret both as a unique problem. The TSP can be seen as a planning problem where each time period t represents a round-trip of the salesman where each node represents a product to produce and the arcs are the time needed to changeover from one product to the other. Following this idea, each round-trip can be seen as a scheduling problem, as it is presented in the right-hand side, where the time spent by the salesman at each node can be seen as the production time of this product in the production stage l. Thus, the TSP formulations can be used to model the planning and scheduling problem presented above.
Using the strength of the TSP and the precedence-based concepts, novel hybrid formulations are developed to tackle the supply chain planning and scheduling problem in multiproduct multistage continuous plants. This hybridisation refers to the concept of merging discrete and continuous time

Fig. 3 -Analogy between TSP and precedence-based ideas for a multiproduct multistage continuous plant.
representations for planning and scheduling decisions respectively in a single robust formulation.
First, we introduce the main notation used in these models and then in [M1] we extend the MILP model proposed by Liu et al. (2012) for a single-stage plant, to couple with multiple stages by using the main ideas of immediate-precedence concepts (see also Liu et al., 2008Liu et al., , 2009Liu et al., , 2010a. Then, in [M2], we present a TSP/general-precedence representation for the multiproduct multistage continuous plant in the supply chain network.

Notation
Time periods Set of markets PI n Set of products stored at echelon node n PM n Set of products sold at echelon node n PS n Set of products produce at echelon node n S Set of plants

Parameters CC ijn
Unit changeover cost from product i to j at echelon node n ($) CP in Unit production cost of product i at echelon node n ($) Changeover time from product i to j in stage l at echelon node n (hours) Binary variables E int 1 if product i is processed at echelon node n in time period t F int 1 if product i is the first one processed at echelon node n in time period t L int 1 if product i is the last one processed at echelon node n in time period t X ijnt 1 if product i is processed before product j at echelon node n in time period t Y inkt 1 if product i is sold at price level k at the echelon node n in time period t Z ijnt 1 if product i is processed just before product j at echelon node n in time period t ZF ijnt 1 if product i in period t-1 is processed just before product j in period t at echelon node n Positive variables CET1 lnt Changeover time in stage l at the beginning of time period t at echelon node n (hours)

TSP-continuous time formulation using immediate-precedence sequencing variables [M1]
Here, we introduce an optimisation model based on Liu et al.

Timing constraints
Timing constraints are presented in Eqs.
(1)-(4). Notice that Eqs. (1)-(2) define the lower (LB) and upper bounds (UB) of the processing time of each product i at production stage l at echelon node n in time period t (PT ilnt ), while Eqs. (3)-(4) determine the changeover time constraints and the upper bound constraint for each time period. Here, the total changeover time between two consecutive time periods c ijln is split in two consecutive time periods t and t-1 by variables CET1 lnt and CET2 ln,t-1 .

Sequencing constraints in the same unit
Eqs. (5)-(7) define sequencing decisions in the each processing unit. Eq. (5) is proposed based on the immediate-precedence concepts, while Eqs. (6)-(7) are used to bound the initial time of each task TS ilnt based on the information of the processing time PT ilnt and changeover times CET1 lnt and CET2 lnt .

Sequencing constraints in consecutive units
Then, Eqs. (8)- (9) are introduced in this model to take into account sequencing decisions between consecutive stages l and l + 1.

Production constraints
Eqs. (10)-(11) are derived to calculate the production amount PR ilnt determined by the unit production rate r iln and production time PT ilnt in each time period t and affected by the production yield at that stage, y iln , as well as the production amount from the previous stage.

TSP-continuous time formulation using general-precedence sequencing variables [M2]
This model is developed to consider sequence-depending changeover issues with general-precedence variables instead of the immediate-precedence variables proposed in [M1]. In the proposed model, Eq. (5) is replaced by Eqs. (12)-(15) after the introduction of the new variables. The main advantage of using general-precedence concepts relies on the possibility to obtain tight formulation in terms of binary variables. Thus, significant improvement can be achieved when we try to decompose the entire model in reduced sub-models to be solved sequentially.

Sequencing constraints in the same unit
Here, we introduce a new general-precedence variable X ijnt to consider sequencing decisions between different products in the same unit in Eqs. (12)-(13). This new variable only runs for product combinations when i > j, but not for all i / = j as the immediate-precedence formulation in [M1], reducing the number of binaries to a half. Here, variable Z ijnt can be set as a continuous variable. Then, we have to introduce additional constraints, Eqs. (14)-(15), to force this variable to takes 0-1 values. Despite of the reduced number of binary variables in this formulation, we need to create much more continuous variables and equations that also affect the convergence of the model.

Efficient MILP-based solution approaches
The solution approach developed in this work emphasises the strengths of MPC approach for dealing with the uncertainty and a tailored LS algorithm for this problem. The aim of this solution methodology is to build an initial solution quick using an HMPC approach developed in Liu et al. (2012) in a constructive phase, and then enhance this solution, in an improvement phase, by using a LS algorithm. Similar ideas of LS strategies can be found in Liu et al. (2010b); Méndez and Cerdá (2003); Kopanos et al. (2010); Castro et al. (2011b); Aguirre et al. (2014Aguirre et al. ( , 2017. More information about efficient solution approaches for MILP-based solution approaches for the process industry can be found in Kallrath (2000); Méndez et al. (2006) and Harjunkoski et al. (2014). The main idea of the improvement phase is to enhance the local solution obtained by the HMPC evaluating the impact of releasing a number of products (NP) and/or time periods (NT) from the schedule in the current horizon for each selected number of plants (NS). Thus, the HMPC tends to deal with the uncertain behaviour of the system, when actual information is revealed, providing an initial feasible and good quality result,

Fig. 4 -Basic ideas of the HMPC + LS solution approach.
while the LS aims to improve this solution at each iteration in a quick way.
One of the main differences from existing LS solutions strategies in the literature is here the possibility to decompose the problem in several ways, releasing decision variables for time periods, products or echelon nodes per iteration. We called that as temporal (LS T ), production (LS P ) and spatial (LS S ) decomposition respectively. We will demonstrate that a total improvement of 10% can be reached by using LS in an inner loop inside HMPC.
A basic notion of the HMPC + LS algorithm is presented in Fig. 4 Here it can be seen how these two algorithms interact at each iteration. Thus, given a control horizon length of the problem L CH and a planning horizon H, in each iteration iter, the algorithm starts running the HMPC for periods iter ≤ t < L CH + iter to obtain an initial local solution and then runs the LS algorithm for the same period in order to obtain an improved result. A detailed explanation of the different steps of HMPC approach and the different strategies for the LS are provided below in Sections 4.1 and 4.2, respectively. Once an initial local solution is obtained, all decision variables for period t = iter, and the iteration is updated by iter = iter + 1 to consider the following control horizon. The algorithm ends when the iteration reaches the horizon length of the problem iter = H.
In order to compare the efficiency of our solution strategy, a Cutting Plane approach is implemented inside the MPC. The basic idea of this algorithm is presented below in Fig. 5. As it can be seen, the algorithm starts similarly than the one described above defining a control horizon length L CH and a planning horizon H. At each iteration iter, a reduced model is solved to obtain an initial relaxed solution of the problem. Due to sequencing decisions are removed from the original model [M1], relaxed solutions may violate some restrictions of the original problem generating undesirables subtours. These subtours can be eliminated by including subtour breaking constraints. Thus, the following step of this algorithm consists to add integer cuts one-by-one in order to remove subtours iteratively. The addition ends when no more subtours are found, and a feasible integer solution is reported. This procedure is repeated till the number of iterations iter reach the planning horizon H. A detailed description of these steps is presented in Section 4.3.
The main features of the original MPC and the proposed HMPC, LS and Cutting Plane approaches are summarised in Table 2. Further details about equations needed and variables released and fixed per step are presented in the following Sections 4.1-4.3.

Hierarchical Model Predictive Control (HMPC)
The HMPC used in this work is an improved version of the one proposed in Liu et al. (2012). Our HMPC has three consecutive steps: the first one to determine which products have to be processed at each time period, the second one for sequencing decisions of selected products at each time period and the third one to check solution reliability. If the solution is not reliable, the addition of integer cuts at the end is possible, removing the current solution and giving the possibility of finding a new feasible solution. Infeasible solutions may happen for different reasons due to timing limitations, production bounds, and inventory and backlog levels. The idea of decomposing the problem in only two steps, solving first the assignment decisions and subsequently adding sequencing decisions, may turn the problem to an infeasible solution. In order to avoid that, additional timing constraints can be added to the first step to overestimate the production time required for the production of the selected products in a particular time period, as in Liu et al. (2012). The equations have the limitation to not consider in detail sequence-dependent changeovers, whose times, in some cases, are larger than the production time for a particular product, when triangular inequality constraints are violated. For this particular cases, the idea of maintaining the information about sequencing decisions, without explicitly considering timing constraints, will give us a better overestimation of timing decisions, and in many of the cases better bounds and consequently better results. Despite of that, possible sub-tours may occur, and additional sub-tour elimination constraints must be added in the second step, as timing constraints, in order to avoid infeasible solutions in the system. In this second step, sequencing decisions are optimised by fixing the information of the selected products performed at each time period. The result after this step may still return an infeasible solution, which is partially solved adding a new integer cut of the current solution and returning to the first step. This sequential procedure is performed till a feasible solution for assignment and sequencing decisions is found. Finally, assignment and sequencing decisions are fixed and timing, production, inventory, demand and the pricing decisions are optimised by solving a reduced model in the third step. A suitable description of each step is presented below while a detailed algorithm is shown in Fig. 6.

Initialisation: dealing with uncertainty
Before starting, the control horizon length L CH , forecast error int and iteration parameter iter must be initialised. Here we consider uncertainty in the demand forecast where D 0 int represents the initial demand at time period t, corresponding to a given initial price, P 0 in . This uncertain initial forecast demand D 0 int follows a uniform distribution, with an expected value of DF 0 int and a forecast error of ˛i nt ∈ (0, 1). In general, the initial demand is revealed at the beginning of the current time period t, allowing the possibility to react with the current schedule by changing some operational decisions to adapt to the new state of the system (reactive scheduling). Thus, at each iteration, iter, the uncertainty is revealed and the demand D 0 int is updated by a uniform distribution according to the following expression. For the remaining periods, t > iter, we assume that D 0 int = DF 0 int . After that, the algorithm starts by solving the problem in a hierarchical way.
∀n ∈ M, i ∈ PM n , t

First Step: solve the production assignment problem
In the first step, a reduced model [M1 1] is solved without sequencing constraints. To ensure that in every time period some products must be produced we introduce Eq. (16). Eq. (17) forces Z ijnt and ZF ijnt to adopt positive values to avoid possible infeasibilities at timing constraints in Eqs.  Step 1

Additional cuts to remove infeasible solutions
In case that an infeasible solution is reported due to timing constraints in Eqs.
(1)-(4), additional integer cuts can be applied by Eq. (20) to remove the previous solution (Z old ijnt , ZF old ijnt , E old int ). This equation should be added sequentially into model [M1 1] after each iteration (iter) returning to the first step. Thus, infeasible sequences are removed sequentially until a feasible solution is reached. When a feasible solution is reached, the algorithm proceeds to follow to the next step, updating the initial solution and going through LS algorithm.

Local Search algorithm (LS)
The LS approach is described in detail in Fig. 6. This LS algorithm starts by selecting a proper number of products (NP), time periods (NT) and plant (NS) to release at each iteration, denoted by iter1 for iterations of time periods, and iter2 for iterations of products. These parameters are set at the beginning, defining the search strategy of the algorithm in the selection procedure. The selection procedure defines a set of products i and periods t to release at each production plant where a nt and b int are binary parameters to set whether time period t, product i at selected plant n are released for re-optimisation. Thus, if the time period t is between iter1 and iter1 + NT, parameter a nt is set to one to release this time period t for each production plant n. In the same way, if the product i is between iter2 and iter2 + NP then this product is selected for reschedule by setting b int = 1 for each selected plant defined by iter3 and iter3 + NS. So, when a nt and b int are defined, some decision variables are fixed {X ijnt ,  [0,1] values due to the binary value of X ijnt . Based on that, we can define two types of decomposition, temporal decompositions and production decomposition. The first one is defined when releasing time periods defined by NT. The second one is done when a set of products are selected to reschedule by NP. These two techniques can be combined in order to create hybrid strategies. Also, those steps can be run for each plant n separately or for all the plants simultaneously, allowing the possibility of spatial decomposition aside from temporal and product decomposition as explained above. If spatial decomposition is applied, these steps are run in a closed loop for each plant n by fixing the partial solution after each iteration (dotted line in Fig. 6). Finally, when the total number of products iter2 + NP = |I| and time periods iter1 + NT = TH are reached, the LS algorithm finishes. All decision variables for the current time period t = iter are fixed, and Z ijnt variable returns to a dis- crete mode by setting Z.prior = 1, and then the upper loop goes back to the next MPC iteration. A simple example considering products i1, i2 and i3 for a plant n and time period t is presented in Fig. 7 showing which decisions variables are released and which remain fixed when a single product is selected to be optimised (a nt & b int ) by the algorithm. In this case, by selecting i2 (a nt = 1 & b i2,nt = 1), assignment decisions E i2,nt , for the plant n and time period t are released and optmised by [M2], as well as the sequencing decisions X i2,jnt , while the assignment decisions of the other two products, i1 and i3, (E i1,nt and E i3,nt ) and the generalprecedence sequencing decision between these two (X i3,i1,nt ) remain fixed, with b i1,nt = 0 and b i3,nt = 0. The rest of decisions of all products (L int , F int Z ijnt , ZF ijnt ) are reoptimised by [M2].

Decomposition techniques
In this section we summarise the main decomposition techniques (temporal, production and spatial decompositions) used in this work:

Cutting Plane approach
This Cutting Plane approach, specially adapted for this problem, has been embedded into the original MPC described above in order to solve the problem considering forecast demand under uncertainty.

First step: solve the relaxed problem without sequencing constraints
The first step provides a feasible solution of the reduced model

Second step: solve the relaxed problem by adding integer cuts
In the second step, sub-tours are identified and additional cuts are added to the reduced model [rM1] by Eq. (20). This repre- Table 3 -Set of products produced, stored and sold at echelon node n.

Example
In this section, we are going to study a large-scale version of the example provided by Liu et al. (2012). The original example proposed for a single-stage problem with 3 production plants (s1-s3), 8 distribution centres (c1-c8) and 16 markets (m1-m16) is extended to consider 3 single-unit stages (l1-l3) and 20 products (i1-i20) by using a similar set of data (see Supplementary  Material). The structure of the supply chain was presented in Fig. 1 -m12} and {c7-c8} with {m13-m16}. The set of products produced, stored and sold at each echelon node n are defined by PS n , PI n , PM n in Table 3. Important information about problem parameters are also summarised in Table 4.
In the analysis, we compare the results provided by the original MPC, in which model [M1] is solved by iterations, the HMPC given in Section 4.1, and HMPC + LS approaches proposed in this work by using temporal, production and spatial decompositions, considering price elasticity of demand. We also evaluate the behaviour of these approaches in comparison with a commonly known cutting plane approach at different control horizon lengths (L CH = 6,8,10). This cutting plane approach was embedded into the original MPC.
The time limit imposed for the LS algoritms is 60 s per iteration and the optimality gap is set to 1%. For HMPC, a CPU limit of 60 s per model at each iteration and the relative termination tolerance is set to 5%. Also, in the objective function, we use a weight for inventory w i = 2.5 and for the price change w p = 10. Finally, for a fair comparison with HMPC, MPC, and Cutting Plane is solved within a time limit of 180 s per iteration and an optimality gap of 5%.

Computational results
This example considers an error of 20% in the demand forecast for all periods, products and markets. In order to keep this uncertainty under control, the MPC is solved for the whole scheduling horizon with the real demand revealed for the current period, and the initial demand forecast for the following unrevealed periods. The results about Profit ($), CPU time (s) performance and Improvement percentage IP (%) for different methods considering price elasticity of demand are reported in Table 5. Results in Table 5 indicate that the best solution has been found by LS T ($1,149,871) after 2232 s for a horizon length of L CH = 10, while the worst solution has been reached by Cutting Plane approach ($835,707) for L CH = 6 after 1111 s. It is worth to notice that in general MPC could find good quality results (IP% = 3.4% in average in comparison with HMPC) but at expenses of longer CPU times. Also here it can be seen that the improvement percentage IP(%) of LS T after HMPC is in average 10.27%, while the IP(%) is 4.23% for LS P and 5.26% for LS S . All of HMPC + LS approaches perform better than the MPC and the Cutting Plane (average IP% = −1.35%) in about a quarter of the CPU time, which demostrates the effectiveness of integrating LS inside the HMPC for the large scale example considered in this work. Note that for this example, when L CH = 10 the MPC and Cutting Plane achieve lower profits than when L CH = 8, due to the more computational difficulties in each iteration with a longer control horizon, and the resulting further suboptimisality of the obtained solutions within the CPU limits.
A ranking profile for (a) profit and (b) computational time and (c) overall performance for the six different approaches under study is shown in Fig. 8. Here the overall ranking is calculated based on the average profit but using the average CPU time as a tiebreaker.
Ranking the profits from the worse (1) to the best (6), it can be seen in Fig. 8a, that the best solution for each horizon length L CH is obtained by LS T , followed by LS S and LS P , while MPC only could have a very good solution for L CH = 6 and worse than LS algorithms for the other L CH . HMPC and Cutting Plane are considerably worse than the others. On the other hand, the ranking for the computational time in Fig. 8b shows that MPC and Cutting Plane approaches report worse CPU time followed by LS S , LS T and LS P , respectively, while

Fig. 8 -Ranking profile for (a) Profit (b) CPU time and (c) Overall (Profit & CPU time) for different approaches.
HMPC provides the best CPUs for the L CH cases analysed. Thus, it is easy to demonstrate that HMPC + LS methods can take the advantages of the decomposition improving the MPC and reducing the total CPU time of the Cutting Plane approach. For this example, it can be seen that the temporal decomposition in LS T is more effective than other methods followed by spatial decomposition in LS S and production decomposition in LS P by considering an overall ranking in Fig. 8c. In general, Cutting Plane approach provides similar solutions to HMPC, but spend a similar CPU time than MPC, which demonstrates that this is not a promising approach for this particular problem. One of the main reasons is that Cutting Plane spends a lot of time adding cuts to the relaxed problem and running many iterations to find a feasible solution, while HMPC takes the advantages of the sequential procedure solving the model in a hierarchical way on a single iteration.
Finally, Fig. 9 shows the behaviour of the approaches presented above for this problem, analysing the estimated profit ($) and the improvement percentage (IP%) obtained at each iteration along the planning horizon considering a horizon length L CH = 10. Here it can be seen how HMPC takes the advantage of solving the problem in a hierarquical way providing better results than the original MPC and the Cutting Plane approach at each iteration providing good results with less variability. The corresponding data of Fig. 9 is summarised in Table 6. Here it can be seen that the average improvement percentage IP% of MPC against HMPC is around −7.5% with a standard deviation of 3.9% while the Cutting Plane is in average −1.4% with a standard deviation of 1.8%.
The results of the MPC and Cutting Plane approaches vary, alternating between good and bad solutions at each iteration in comparison with HMPC. Bad solutions maybe occur when solver is interrupted, leaving a huge relative gap. This may happen if full MILP model [M1] and/or reduced model [rM1] can't provide a solution with a small relative gap (<5%) after a time limit imposed for each iteration. However, in some iterations, the Cutting Plane approach has been able to find good quality results, even better than the one provided by the HMPC. On the other hand, LS methods (LS S , LS T and LS P ) behave similarly providing improved results (Avg.: 2.4%, 1.9%, 1.8%) than HMPC with less deviations (Stdev.: 0.8%, 0.9%, 1.0%).

Discussion
Three scenarios considering time horizon L CH = 6, 8, and 10 (weeks) of an extended problem from the literature were presented in this paper and solved by using six different approaches; a well-known MPC and a Cutting Plane approach  developed for this problem and proposed hierarchical MPC called HMPC and LS approaches for temporal, production and spatial decomposition. The integrated framework HMPC + LS has obtained better results than other existing approaches for the extended case study proposed here. In addition, further discussions on the proposed approaches and computational results are summarised as follow: • Novelty of HMPC: In order to deal with the uncertainty in demand forecast, a hierarchical MPC has been proposed, allowing solving the problem in a sequential way obtaining good initial results in short CPU time. The main contribution of our HMPC is the possibility of dealing with no reliable solutions, with the addition of additional integer cuts, returning always a feasible solution of the system in a few iterations. As we have mentioned before, infeasibilities may happen when then minimum production amount cannot be produced in a specific time period or due to the inventory or backlog restrictions. • HMPC vs MPC: One of the main drawbacks of the original MPC in comparison with HMPC relies on the CPU time needed to converge to a reasonable gap (5%). For that reason, we imposed a CPU time limit (3 min) at each iteration of the MPC, as well as the HMPC. According to this, we could demonstrate that the HMPC has obtained similar results than the original MPC, 3% worse in average, reducing the CPU time in one order of magnitude. • CPU limit selection: The solution of MPC and HMPC do not improve much by increasing the CPU time limit. Based on the study, a reasonable computational time of 3 min per iteration is sufficient enough to obtain good quality solutions per iteration. For the practical reasons, a CPU limit should be defined according to the model behaviour and the total computational time imposed for the comparison (1 h). • Novelty of LS: Taking the advantages of the solution reported by HMPC, LS decomposition has been integrated as an inner loop inside the HMPC finding improved results. Three LS approaches have been developed in order to decompose the problem in temporal, production and spatial ways allowing the possibility of improving the initial solution in between 4-10% average within an hour of CPU time.
• Cutting Plane approach: In order to compare with other existing decomposition methods, a Cutting Plane approach has been implemented. In order to find good results, this Cutting Plane approach has been embedded inside the MPC and solved iteratively. The average result reported by this method was similar to the one reported by the HMPC but with a higher computational effort. • HMPC + LS vs. Cutting Plane: The efficiency of the Cutting Plane approach relies on the number of cuts added to the MPC to find a feasible solution. The addition of these cuts increases the model size taking much more CPU time to find relatively good results. On the other hand, the integrated HMPC + LS is able to find a feasible solution very fast while spending the rest of the time trying to improve it.

Concluding remarks
An MILP-based model for the supply chain planning and scheduling of process plants under uncertainty has been presented. Planning decisions, such production amount, inventory, sales, backlog and pricing, and scheduling decisions for each time period, such as production assignment, sequencing and timing, have been considered. In this study, an extended version of a supply chain planning and scheduling problem from the literature has been considered. Here a multiple products multistage continuous plant with sequence-dependent changeovers and demand under uncertainty is addressed. To deal with this complex problem, different hybrid decomposition approaches have been introduced and tested. The main contribution relies on the integration of well known efficient solution techniques, like MPC and LSapproaches, into a unique framework in order to deal with production planning and scheduling under uncertainty in supply chain networks.
Results have demonstrated that the proposed framework has been able to obtain good quality results of a supply chain case study proposed in reasonable computational time. The solutions obtained were compared favorably with other decomposition techniques implemented for this problem. Finally, new features of novel techniques embedded inside the MPC must be studied further in the future in order to create a robust solution strategy for different types of uncertainties. In addition, the investigation of the solution nervousness of the proposed framework on a wider range of examples is an interesting future work direction. et al.

Assignment constraints
Assignment constraints based on the ideas of the TSP problem are defined in Eqs. (A.1)-(A.5). Here, binary variables F int and L int are considered to ensure that only one product is assigned as first or last at echelon node n in time period t by Eqs. (A.1)-(A.2). Then, variable Z ijnt is derived to take into account the immediate-precedence relationship between products i and j at echelon node n in time period t while ZF ijnt defines if product i in time period t−1 is processed just before product j in time period t at echelon node n. Thus, Eqs. (A.3)-(A.5) are derived to ensure only one predecessor and successor of every product i assigned to echelon node n in time period t by binary variable E int . Note that, Eqs. (A.3)-(A.5) are slightly modified from the original version of Liu et al. (2012) in order to have a reduced number of equations of the asymmetric TSP model (see Aguirre et al., 2017).

Inventory constraints
Inventory level I int of each product i at each echelon node n in time period t is stated in Eqs. (A.6)-(A.8), by relating production amounts PR ilnt , sales SV int and product flows Q inn t . Here the initial inventory is defined by I 0 in .

Demand constraints
Eq. (A.9) is provided to estimate the unmet demand U int of products at echelon node n at the end of each time period.
U int = D int − SV int ∀n ∈ M, i ∈ PM n , t (A.9)

Price elasticity of demand
When the initial demand D 0 int , and price P 0 in are known, the manager can determine the price P int to affect the actual demand in the market D int . These two continuous variables are related to each other by the price elasticity of demand coefficient PE int , as stated in Eq. (A.10).
in P 0 in ∀n ∈ M, i ∈ PM n , t (A.10)

Pricing constraints
In order to incorporate pricing decisions Eqs. (A.11)-(A.12) are introduced. Eq. (A.11) selects a certain price at level k from a set of available prices P L ink by incorporating a new binary variable Y ink . Moreover, Eq. (A.12) enforces that only one price could be selected for each product i in each market n at time period t. P int = k P L ink · Y inkt ∀n ∈ M, i ∈ PM n , t (A.11) k Y inkt = 1 ∀n ∈ M, i ∈ PM n , t (A.12)

Linearisation constraints
Eqs. (A.13)-(A.14) are considered to linearise the product of price Y inkt and sales SV int in the profit calculation. An auxiliary variable is introduced and can be activated only once by Y ink variable in Eq. (A.12).

Inventory deviation
Eqs. (A.15)-(A.16) are introduced in order to maintain the inventory deviations I D int in a low level along the entire planning horizon. This new variable is optimised with the total profit stated in Eq. (A.20) to reduce the differences between the inventory target I T int and the real inventory level I int .
I D int ≥ I T int − I int ∀n, i ∈ PI n , t (A.15) I D int ≥ I int − I T int ∀n, i ∈ PI n , t (A.16)

Price change
Similarly, price changes PC int are considered by Eqs. (A.17)-(A.18).This is calculated as the difference between prices at consecutive time periods t−1 and t. For the first period an initial price is considered by P 0 in .

Total profit
The total profit is provided in Eq. (A.19) by considering sales revenues, production costs, changeover costs, transportation costs and unmet demand costs. TP = t,k,n ∈ M,i ∈ PMn P L ink · SY inkt − t,n ∈ M,i ∈ PMn,l=|L| CP in · PR ilnt − t,n ∈ S,i,j ∈ PSn,i / = j CC ijn · Z ijnt − t>1,n ∈ S,i,j ∈ PSn CC ijn · ZF ijnt − t,n ∈ S,n ∈ T inn' ,i ∈ PIn CT inn · Q inn t − t,n ∈ M,i ∈ PMn CU in · U int (A.19)

Objective function
Finally, the objective function OF is defined by Eq. (A.20) introducing a weight for the total price change w p and a weight for the total inventory deviation w i .