Framework for studying online production scheduling under endogenous uncertainty

Framework for Studying Online Production Scheduling under Endogenous Uncertainty

Journal Pre-proof

Framework for Studying Online Production Scheduling under Endogenous Uncertainty Dhruv Gupta, Christos T. Maravelias PII: DOI: Reference:

S0098-1354(19)30607-6 https://doi.org/10.1016/j.compchemeng.2019.106670 CACE 106670

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

9 June 2019 19 October 2019 29 November 2019

Please cite this article as: Dhruv Gupta, Christos T. Maravelias, Framework for Studying Online Production Scheduling under Endogenous Uncertainty, Computers and Chemical Engineering (2019), doi: https://doi.org/10.1016/j.compchemeng.2019.106670

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Framework for Studying Online Production Scheduling under Endogenous Uncertainty Dhruv Gupta, Christos T. Maravelias∗ Department of Chemical and Biological Engineering, University of Wisconsin-Madison 1415 Engineering Drive, Madison, WI 53706, USA

Abstract We propose a framework for studying online production scheduling in the presence of endogenous uncertainties. We address uncertainties in (i) processing times; (ii) batch yields; and (iii) unit operating status. First, we illustrate how uncertainty can result in infeasibilities in the incumbent schedule and propose a model for systematic schedule adjustment to restore feasibility in the absence of new scheduling inputs. In this model, we define variables to track and penalize changes between the new and old schedule. Second, we discuss the different probability distributions for the three uncertainties that we consider in this work and how the parameters for these distributions change with sampling frequency. Third, we present a formal procedure for carrying out closed-loop simulations and evaluating closedloop performance in the presence of these uncertainties. Finally, using this framework we draw useful insights for the design of online scheduling algorithm in the presence of the above three uncertainties. Keywords: Chemical production scheduling; process uncertainty; disturbances; re-optimization; model predictive control. 1. Introduction Scheduling plays an important role in all industrial production facilities (Harjunkoski et al., 2014). Contingent on the scale of operation, optimization based scheduling methods can even achieve multi-million dollars increase in profits (Kelly and Mann, 2003). Once ∗

Corresponding author: [email protected]

Preprint submitted to Elsevier

December 9, 2019

a schedule has been obtained, disruptions or arrival of new information can make the incumbent schedule suboptimal or even infeasible, thus motivating the need for online (re)scheduling (Gupta et al., 2016). On the basis of decision-dependence, uncertainty can be classified into two types: (i) exogenous; and (ii) endogenous. In the exogenous case, the probability distribution and the time of observation of the uncertainty is independent of any decisions. For example, the demand for a product. In contrast, in the endogenous case, either the time of observation or the underlying probability distribution of the uncertainty parameters is influenced by the decisions (Grossmann et al., 2016). For example, the success of a drug trial can be known only when the trial is actually conducted (Colvin and Maravelias, 2011). Similarly, the choice of the drilling technology, and the decision to drill, together, affect the accessible size of an oil reserve, and the time to observe that size (Goel and Grossmann, 2006). This makes modeling and testing of approaches in the presence of endogenous uncertainty challenging. An important consideration when calculating schedules is the incorporation of uncertainty in the online model. For this several approaches exist: (i) chance constrained; (ii) stochastic; (iii) robust; and (iv) deterministic optimization. Chance constrained optimization is aimed towards maintaining feasibility (satisfying constraints) with a target probability. Stochastic optimization approaches produce a recourse strategy by maximizing or minimizing an expectation function comprising of the decisions and uncertainty represented through random variables. Although, stochastic optimization approaches, if solved to near optimality, may find good (open-loop) solutions fully utilizing the underlying probability distribution of the uncertainty, their online computation, even for relatively small problems, is typically highly expensive. Robust optimization approaches focus on maintaining feasibility and have the advantage of retaining similar computational tractability as the deterministic model. Further, they are applicable even in the absence of probability distribution of the uncertain parameters, as long as the bounds on the parameters are available. However, a disadvantage of robust optimization is that decisions are over-conservative and incur worst-case costs

2

irrespective of the actual uncertainty observation. When uncertainty is present, online scheduling methods need not always explicitly model it; these methods can employ a deterministic model, and still mitigate the effect of uncertainties by treating them as disturbances, and relying solely on recourse, through feedback (Subramanian et al., 2012). Gupta and Maravelias (2019) proposed a simulation based design framework for online scheduling in which they gave recommendations on how to choose an appropriate horizon, re-optimization time-step, optimality gap, objective function modifications, and added constraints to achieve good closed-loop performance in the presence of exogenous uncertainty. The presence of endogenous uncertainty makes the simulation of online schedules challenging. This is because of two reasons: (i) the uncertainty is dependent on the decisions, which complicates the sampling of uncertainty parameters; and (ii) the majority of these uncertainties are observed in the batches already under execution which can result in infeasibilities. In this work, we address these challenges enabling researches to evaluate different online scheduling algorithms as well as carry out systematic simulation based design (as recommended in Gupta and Maravelias (2019)). We propose a framework for studying online scheduling in the presence of uncertainty in (i) processing times; (ii) batch yields; and (iii) unit operating status. The paper is structured as follows: In Section 2, we provide necessary background. In Section 3, we present a motivating example to understand how uncertainty can result in infeasibilities and emphasize the importance of a good resolution mechanism for these infeasibilities. Further, we discuss the implication of the time at which an uncertainty is observed, what makes it a disturbance, and how to generate perfect information solutions for benchmarking. In Section 4, we propose a framework in which we introduce a model for adjusting scheduling decisions for resolving infeasibilities, we discuss appropriate probability distributions for the modeling of the three uncertainties, their scaling with sampling frequency, and a simulation procedure to evaluate closed-loop performance. Finally, in Section 5, we

3

analyze closed-loop simulation results and draw useful insights for the design of the online scheduling algorithm. Throughout the text, we use lower case Latin characters for indices, uppercase Latin bold letters for sets, uppercase Latin characters for variables, and Greek letters for parameters. 2. Background In this section, first, we outline the general chemical production scheduling problem statement, problem representation framework, model classification, and solution methods. Second, we describe what online scheduling is. Third, we present the state-space scheduling model used in this work. Fourth, we discuss the types of endogenous uncertainty. Finally, we review related research. 2.1. Chemical production scheduling 2.1.1. Problem statement The chemical production scheduling problem can be defined as follows. Given: (i) (ii) (iii) (iv) (v) (vi)

Production facility data (e.g., unit capacities and connectivity); Production recipes (e.g., processing times and mixing rules); Production costs (e.g., material holding costs and unit startup costs); Material availability (e.g., raw materials delivery amounts and dates); Resource availability (e.g., maintenance schedule and utility levels); and Production targets or orders with due-times;

scheduling seeks to find: (i) Selection and sizing of batches; (ii) Assignment of these batches to units; and (iii) Timing of these batches; so as to meet the production targets at minimum cost, maximum profit, minimum makespan, or while employing any other suitable objective for the given application. In general, several processing constraints or problem features such as time-varying utility costs, sequence 4

dependent changeovers, storage constraints, etc., can be present (M´endez et al., 2006; Harjunkoski et al., 2014). In this work, we consider problems based on network structures, but our proposed methods are general. This is because network based representation can also be utilized to address problems in non-network based environments (Sundaramoorthy and Maravelias, 2011a). 2.1.2. Problem representation A popular representation for chemical production scheduling problem is the State Task Network (STN) (Kondili et al., 1993). The STN representation primarily comprises of tasks i ∈ I, units j ∈ J, and materials k ∈ K. The set of tasks producing/consuming material k are − denoted by I+ ρik mass fraction of k /Ik ; task i consumes/produces material k equivalent to ρik /¯

its batch-size (ρik < 0 for consumption and ρ¯ik > 0 for production). The subset of tasks that can be carried out on unit j are denoted by Ij ; the processing time of task i, when executed on unit j, is denoted by τij . On any given unit, the batch-size of the task can be between the lower (βijM IN ) and upper (βijM AX ) bounds. The fixed and proportional production costs P F respectively. Feed, intermediate, and product and αij of carrying out task i on unit j are αij

materials are denoted by k ∈ KF /KI /KP ; there are possible incoming deliveries (ζkt ) and outgoing demand (ξkt ) at certain times for selected materials; the selling price, inventory cost, and backlog cost of material k are γk , γkIN V , and γkBO , respectively. In Fig. 1, we see the STN representation of a process network with 7 tasks (T1-T7), 8 materials (M0-M7), and 5 units (U1-U5). 2.1.3. Model Classification Models that employ a time-grid are either continuous time or discrete time models (Maravelias, 2012). In discrete time models, the fixed time-grid spacing is denoted by δ. For the purpose of all numerical examples in this work, we use δ = 1 h. Tasks are allowed to take place only at these grid time-points. Thus, all time-related parameters are rounded such that the resulting schedule computed using these new parameter values is feasible even for the

5

Figure 1: STN representation of a process network. The task nodes are denoted by rectangles and the material nodes by circles. Arcs connect task nodes with corresponding input/output material nodes. Tasks can be carried out in compatible units and could require utilities. Task-unit mapping and task batch-size M IN M AX capacities (βij /βij ) are also shown here. Material prices (γk ) are shown adjacent to the material nodes. ρik /¯ ρik values (not shown in the figure) are either −1 or +1 depending on whether the task is consuming the material or producing it.

original parameter values. Even though having a discrete time-grid introduces an approximation error, discrete time models have several advantages over continuous time models. For example, accounting for inventory and backlog costs, utility consumption, time varying prices, or time-dependent resource availability introduces non-linearities in continuous time models, but not so in discrete time models (Velez and Maravelias, 2014). Furthermore, discrete time models are, in general, better suited for large scale problems with several additional processing features (Sundaramoorthy and Maravelias, 2011b; Lagzi et al., 2017). 2.1.4. Solution methods To tackle the computational challenge of mixed-integer linear programming (MILP) scheduling models, several solution methods have been proposed. Specifically, for problems in network environment, researchers have developed: (i) tightening methods based on preprocessing algorithms and valid inequalities (Blomer and Gunther, 2000; Velez et al., 2013; Burkard and Hatzl, 2005); (ii) reformulations (Sahinidis and Grossmann, 1991; Yee and Shah, 1998; Velez and Maravelias, 2013b); (iii) decomposition methods (Papageorgiou and Pantelides, 1996; Bassett et al., 1996; Kelly and Zyngier, 2008; Wu and Ierapetritou, 6

2003; Castro et al., 2011; Calfa et al., 2013); (iv) heuristics (M´endez and Cerd´a, 2003; Relvas et al., 2009); (v) hybrid methods (Maravelias and Grossmann, 2004; Roe et al., 2005); and (vi) heirarchical methods based on different grids (Lee and Maravelias, 2018a,b). Finally, parallel computing has been utilized to obtain faster solutions (Subrahmanyam et al., 1996; Velez and Maravelias, 2013a). 2.2. Online scheduling Online scheduling is a generalization of rescheduling as an ongoing process. In each online scheduling iteration an open-loop (optimization) problem is solved (see Fig. 2). The implemented schedule which can also be called the closed-loop schedule (or closed-loop solution) results from solving several open-loop problems that account for feedback. We denote the scheduling horizon by H and the time between two periodic online scheduling iterations, the re-optimization time-step, by ∆. The re-optimization time-step should not be confused with δ, which denotes the time-grid spacing employed in discrete time scheduling models. In the online model, a finite rolling horizon is used because it is otherwise computationally expensive to solve large long-term scheduling models and, with reasonable (un)certainty, information only for finite future is known. As time passes by, feedback and new information gradually becomes available; these should be accounted for as soon as possible to determine new decisions. Engell (2009) considered production scheduling as a control problem under uncertainty and provided a good overview on the role of feedback. Analogous to the control input move penalties in MPC (Rawlings and Mayne, 2009), there can be penalties (nervousness costs) imposed in the objective function to minimize the differences between schedules computed in adjacent iterations (Lee et al., 2019). Online scheduling can also be carried out as rule based rather than optimization based. Sabuncuoglu and Karabuk (1999) showed for a flexible machine shop environment that although optimum-seeking approaches use tight resources effectively, dispatch rules can be superior when information rapidly becomes obsolete. Dispatch rules based scheduling approaches, however, are not popular in chemical production environments. 7

Figure 2: In online scheduling, schedules are computed iteratively. Here, to make a schedule for 5 days, each 7 day iteration is computed on the start of each day. The schedule for that day is implemented, and the next day, a new schedule is computed for the next 7 day. At the end of each day, feedback, such as, any delayed batches, is incorporated in while making the schedule for the next 7 day. Further, the demand for the new day (new information) at the end of the horizon due to the advancement of the horizon is accounted for in the current iteration. In this example, H = 7 days and ∆ = 1 day.

2.3. Model In this work, we employ a discrete time-grid state-space model. This is because statespace is a natural (re)formulation for models used in online scheduling. Although, we use the STN discrete time state-space model of Gupta and Maravelias (2017), the ideas and conclusions presented in this work are general and are not model specific. They are also applicable to resource task network (RTN) (Pantelides, 1994) discrete time models as well as continuous time models in STN and RTN representations (please see supplementary material for details). Binary variable Wijt , when 1, implies a batch of task i is starting on unit j at time t; n n ¯ ijt ¯ijt variable Bijt denotes its batch-size. W and B are the lifted task-states. The index n

is the progress status in a batch of a task. For a batch at the start, n = 0, and at finish, n = τij . Skt and BOkt are the inventory and backlog level, respectively, of material k during time-interval (t − 1, t]. We carryover the scheduling decisions and the material inventory and backlog levels from the previous (open-loop) iteration (σ − 1) to the current iteration (σ) using Eqs. 1 - 6. In these equations, we write the prescripts σ and σ − 1 for the variables of the current and the previous iteration, respectively. In other equations, when these prescripts are absent, they are σ. For iteration σ − 1, we denote the time index by t00 and for iteration σ by t; and in each iteration, we reset the time-grid to begin from 0. Binary parameters Y˙ ijn is 1 when there is a delay in a batch of task i on unit j with progress status n. 8

B

Y˙ ijn is

the batch-size of the delayed batch. Xijt and B Xijt are variables which are equated to the 0 B ˆ0 appropriate delay parameters (Y˙ ij0 , B Y˙ ij0 , Yˆijt , Yijt ; Eqs. 4, 6, 7, 10); as explained in Gupta

and Maravelias (2017), when delays are observed, these variables are necessary to preserve the physical meaning of n, the progress status of the batch. Binary parameter Z˙ ijn is 1, when there is a breakdown on unit j while executing a batch of task i with progress status n. The corresponding batch-size parameter is B Z˙ ijn . σ Sk(t=0) σ BOk(t=0)

= (σ−1) Sk(t00 =1) ∀k

(1)

= (σ−1) BOk(t00 =1) ∀k

(2)

¯n σ Wij(t=0)

¯ n−100 − Y˙ n−1 + Y˙ n − Z˙ n−1 = (σ−1) W ij ij ij ij(t =0)

σ Xij(t=0)

= Y˙ ij0

¯n σ Bij(t=0)

¯ n−100 − B Y˙ n−1 + B Y˙ n − B Z˙ n−1 = (σ−1) B ij ij ij ij(t =0)

B σ Xij(t=0)

= B Y˙ ij0

∀j, i ∈ Ij , n ∈ {1, 2, ..., τij }

∀j, i ∈ Ij

(3) (4)

∀j, i ∈ Ij , n ∈ {1, 2, ..., τij }

∀j, i ∈ Ij

(5) (6)

¯ n ) with ¯ n and B Eqs. 7-12 are the lifting equations through which the task-states (W ijt ijt the correct progress status (n) are activated. For a multi-period delay of φdelay h, in a batch n are activated for t = 0, t = 1, ..., t = of task i on unit j, in addition to Y˙ ijn , parameters Yˆijt

φdelay − 2. 0 Xij(t+1) = Yˆijt

∀j, i ∈ Ij , t

¯ 0 = Wijt + Xijt W ijt

(7)

∀j, i ∈ Ij , t

n n ¯ ij(t+1) ¯ n−1 − Yˆ n−1 + Yˆijt W =W ijt ijt B

0 Xij(t+1) = B Yˆijt

(8)

∀j, i ∈ Ij , t, n ∈ {1, 2, ..., τij }

∀j, i ∈ Ij , t

0 ¯ijt B = Bijt + B Xijt

(9) (10)

∀j, i ∈ Ij , t

n n ¯ij(t+1) ¯ n−1 − B Yˆ n−1 + B Yˆijt B =B ijt ijt

∀j, i ∈ Ij , t, n ∈ {1, 2, ..., τij }

(11) (12)

Eq. 13 is the assignment constraint, which ensures that only one batch of any task can 9

execute at a time on a production unit. For unit breakdowns with downtime duration of φbreak h, in addition to Z˙ ijn , parameters Zˆjt are activated for t = 0, t = 1, ..., t = φbreak − 1. ij −1 X τX

i∈Ij n=0

n ¯ ijt ≤ 1 − Zˆjt W

∀ j, t

(13)

Eq. 14 enforces lower and upper bounds on the batch-size of the tasks. βijM IN Wijt ≤ Bijt ≤ βijM AX Wijt ∀ i, j, t

(14)

Eqs. 15 and 16 are the inventory balances for all the materials, where, SkT erminal is the terminal inventory level. Vkt is the outgoing shipment to meet demand for material k at time t. βˆijkt is the observed yield-loss.

Sk(t+1) = Skt +

X X j

SkT erminal = Skt +

i∈Ij ∩I+ k

¯ τij − βˆijkt ) + (¯ ρik B ijt

X X j

i∈Ij ∩I+ k

X X j

¯ τij − βˆijkt ) + (¯ ρik B ijt

i∈Ij ∩I− k

ρik Bijt − Vkt + ζkt ∀ k, t < |T| (15)

X X j

i∈Ij ∩I− k

ρik Bijt − Vkt + ζkt ∀ k, t = |T| (16)

Eqs. 17 and 18 are the backlog balance for all product materials, where, BOkT erminal is the terminal backlog level. BOk(t+1) = BOkt + ξkt − Vkt ∀ k ∈ KP , t < |T|

(17)

BOkT erminal = BOkt + ξkt − Vkt ∀ k ∈ KP , t = |T|

(18)

Eq. 19 enforces the correct domain ranges for the variables. n n ¯ ijt ¯ijt Wijt , W ∈ {0, 1}; Bijt , B , Skt , BOkt , SkT erminal , BOkT erminal , Vkt , Xijt , B Xijt ≥ 0

10

(19)

We choose the scheduling objective (Eq. 20) to be the minimization of cost which is composed of inventory costs (γkIN V = 0.15γk ), backlog costs (γkBO = 1.5γk ), and unit startup costs (αF = $100).

min

X k

γkIN V (

X

Skt + SkT erminal ) + γkBO (

t

X t

!

BOkt + BOkT erminal )

+

XXX j

i∈Ij

F αij Wijt

t

(20) 2.4. Endogenous uncertainty Endogenous uncertainty is of three types: (i) Type I - the underlying probability distribution is affected by the decisions, but not the time of observation, e.g., maintenance can restore equipment health and thus affects the probability of breakdown but a breakdown can still be observed at any time the unit is in operation; (ii) Type II - time of observation is affected by the decisions, but not the underlying probability distribution (i.e., exogenous uncertainty with endogenous observation), e.g., in clinical testing of pharmaceutical drugs, only the observation time of drug pass/fail is affected by testing decisions, but not the probability of pass/fail; and (iii) Type III - the time of observation as well as the underlying probability distribution of uncertainty is affected by the decisions, e.g., choice of drilling technology, and the decision to drill, together, affects the accessible size of an oil reserve, and the time to observe that size. Delays, yield-losses, and breakdowns are endogenous uncertainties of type II. 2.5. Related research Rescheduling strategies have been suggested in different contexts (Vieira et al., 2003; Li and Ierapetritou, 2008; Gupta et al., 2016). However, the strategies proposed till date are quite restrictive, either in terms of the scope of the problems that can be addressed, or in terms of the allowable corrections (e.g, only re-timing is allowed but not re-assignment). Kanakamedala et al. (1994) investigated a bi-level least impact heuristic for rescheduling. Huercio et al. (1995) presented a heuristic for rescheduling through batch start time 11

shifting and unit reallocation. Kim and Lee (1997) proposed rule-based reactive scheduling, triggered by a monitoring algorithm, for multi-product batch processes. Huang and Chung (2003) applied heuristics to reschedule pipeless plants on occurrence of unexpected events, such as the failure of a processing unit. Elkamel and Mohindra (1999) addressed rescheduling as a multi-objective optimization problem to account for delays, relative customer importance, and production costs. Vin and Ierapetritou (2000) demonstrated optimization based reactive scheduling using a continuous time MILP model to handle two kinds of disturbances: unit breakdowns and rush order arrivals. M´endez and Cerd´a (2003) presented a user controlled flexible MILP systematic tool to repair and generate updated schedules to react to the occurrence of unforeseen events. Janak et al. (2006) proposed partial rescheduling by identifying and fixing batches not directly affected by an observed disturbance. Qiao et al. (2018) proposed a schedule repair method when machine breakdowns are observed in dynamic seminconductor manufacturing systems. Letsios and Misener (2018) proposed a lexicographic scheduling and approximation algorithm. Honkomp et al. (1999) developed a simulator to evaluate and compare reactive scheduling policies for handling processing time uncertainty. Cui and Engell (2010) applied a two-stage stochastic model in a moving horizon formulation to an expandable polystyrene manufacturing plant. Chu and You (2014) solved a reduced problem online, achieving computational tractability and schedule stability. Kopanos and Pistikopoulos (2014) used a state-space model in a multiparametric optimization framework. Nie et al. (2015) formulated a production scheduling state-space model for RTN representation. Silvente et al. (2015) used a rolling horizon framework for the simultaneous optimization of energy supply and demand planning in microgrids. Touretzky et al. (2016) proposed a methodology to detect faults in the process that act as trigger for rescheduling. Pattison et al. (2016) updated schedules to respond to process and market disturbances. Lindholm et al. (2013) formulated a scheduling model for continuous production sites and, given updates of actual production each day and incoming orders, solved it every day in a rolling horizon fashion.

12

Multiple works have proposed stochastic programming formulations to address endogeneous uncertainty starting with Jonsbr˚ aten et al. (1998). Goel and Grossmann (2006) addressed endogenous uncertainty through stochastic optimization when optimization decisions influence the time of information discovery for a subset of the uncertain parameters. Gupta and Grossmann (2011) proposed solution methods to address the computational challenge associated with these stochastic optimization models. Colvin and Maravelias (2010) studied scheduling of R&D projects which involve endogenous uncertainty. Christian and Cremaschi (2018) proposed a branch and bound algorithm to solve large-scale multistage stochastic programs with endogenous uncertainty. Recently, Lappas and Gounaris (2016) used adjustable robust optimization to generate schedules with a recourse strategy, in which they address uncertainty in processing times (an endogenous uncertainty). Grossmann et al. (2017) reviewed recent advances in methods for optimization under uncertainty and their application to process systems engineering. 3. Preliminaries We first present a motivating example to illustrate how uncertainty can render a schedule infeasible necessitating an immediate adjustment (correction) in order to continue its execution. Second, we discuss the implications of the time when uncertainty is observed, what makes it a disturbance, and how to generate perfect information solutions for benchmarking. Third, we summarize selected insights on the design of online scheduling algorithm from Gupta and Maravelias (2019). 3.1. Motivating example For the network in Fig. 1 an order for 10 tons of M6, and 5 tons each of M5 and M7 is due at t = 11 h. Given zero starting inventory, we find a schedule and see that the order can be met on time (Fig. 3A). A delay is observed in the second batch of T1 at t = 4, which makes the schedule for t ≥ 4 infeasible (Fig. 3B). This is because the batches of T4 and T2 cannot start unless they receive M1 from the (delayed) batch of T1. There are no slack variables 13

Figure 3: (A) Initial schedule calculated at t = 0 to meet an order due at t = 11. (B) When a delay of 1 h (fine horizontal white lines) in T1 is observed at t = 4 (thick vertical red line), the remaining schedule becomes infeasible. (C) The infeasibility can be resolved through a right-shift in the schedule starting at t = 4, but the order at t = 11 is then met late by 1 h. (D) By swapping the batches of T2 and T4, and other related adjustments, feasibility is restored and the order due at t = 11 is met on time.

(e.g., backlog) which can be utilized to remain feasible. Thus, an adjustment is necessary to move forward with the execution of the schedule. Further, simply right-shifting the batches of T4 and T2 back by 1 h does not resolve the infeasibility because then materials M4 and M2 are not available at the start time for the batches of T7 and T5. An 1 h right-shift in the complete schedule from t = 4 resolves the infeasibility (Fig. 3C), but there exists a much better schedule adjustment (Fig. 3D) to mitigate the infeasibility and still meet the order due at t = 11 on time. Hence, we see how uncertainty can result in infeasibilities. These are not trivial to resolve, e.g., by modifying only the start times of the directly affected batches (here T4 and T2) because of the dependence of other indirectly related batches (here T7 and T5). Further, the standard scheduling model (see Section 2.3) does not directly contain variables amenable for shifting, swapping, etc. of batches, which are the natural actions taken on the production floor to restore feasibility. Finally, it is important to understand that an uncertainty observation need not always result in an infeasibility; e.g., in the case when a delay is observed but there is sufficient time already built into the schedule to absorb it (say, through robust scheduling approaches). 14

However, to save on costs and to keep intermediates inventory low, schedules are often “packed” with little to no room for absorbing any observed uncertainties. 3.2. Observation of uncertainty 3.2.1. Observation time Given any uncertain parameters in a model, the uncertainty observation time (η) is the time ahead at which the true value of the parameters gets known. Naturally, knowing the true value as far ahead as possible (e.g., η = H, i.e., no uncertainty) helps in making better decisions. For instance, orders are not allowed to change (deterministically known) too close to the current time; e.g., if orders within the next two weeks are finalized, then η = 2 weeks. However, delays, yield-losses, and unit breakdowns are observed after they have already occured, i.e., η = 0. In online scheduling, given that decisions made betweeen t = 0 and t = ∆ are fixed, any uncertainty observed at η ∈ [0, ∆) can potentially cause an infeasibility. When an uncertainty is observed at η ≥ ∆, it does not result in infeasibility, since, decisions for t ≥ ∆ are tentative and re-computed anyways in the next online iteration. 3.2.2. Disturbances When the model solved online is deterministic, the uncertain parameters are assigned an estimated (nominal) value which could be, e.g., the most likely value or the mean value. When the true value of an uncertain parameter is observed as a surprise (η = 0) and is different than the estimated value, the observed uncertainty is termed as a disturbance. This disturbance can (but not always) render the current solution suboptimal or infeasible. For 0 example, the nominal value for the parameters Y˙ ij0 and Yˆijt is 0; however, when a delay

disturbance is observed the parameters are assigned a value of 1. 3.2.3. Perfect information benchmark Obtaining a perfect information benchmark helps in evaluating the effectiveness of alternate methods for addressing uncertainty. In the case of exogenous uncertainty (e.g., demand),

15

obtaining the (expected) perfect information cost is straightforward. It is the mean of offline solutions, each obtained for an individual exogenous uncertainty sample drawn from a distribution prior to computing the decisions. Obtaining this benchmark, however, is challenging for endogenous uncertainty. Specifically, for endogenous uncertainty of type II, the time at which the uncertain parameters are observed is dependent on the decisions which are themselves dependent on knowing the sample value of the uncertain parameters. For delays, breakdowns, and yield-losses, we can obtain the benchmark as follows: for every task i, assuming a stream of practically infinite batches of this task, we can sample, from the distribution of delays, yield-losses, and unit breakdowns, for the batches in (and the time within the batch at) which these uncertainties are observed. Having obtained these samples, we can now calculate an offline schedule, for each sample, with added constraints to enforce the sample value of uncertainties in the corresponding batches. On obtaining sufficient samples and corresponding schedules, we can then calculate a distribution for the perfect information objective. 3.3. Insights on design of online scheduling algorithm Gupta and Maravelias (2019) discussed a systematic way to design online scheduling algorithms to achieve better closed-loop performance. They first outlined and quantitatively defined the key production system characteristics that are pertinent to the design of online scheduling algorithms: (i) production capacity and load (Λ); (ii) order size max-mean relative difference (); and (iii) time-constants including the demand uncertainty observation time (η; explained in Section 3.2.1). The production capacity of a network is the maximum rate at which it can produce a fixed composition of products. The load is a fraction between 0 and 1, which represents how much of the capacity is utilized for meeting a given demand. It is calculated as the ratio of averaged (over time) demand and the production capacity of the network. The order size max-mean relative difference is the ratio of the difference of the max and mean to the mean size of the order (assuming that the order sizes are bounded and the max exists). 16

Second, they discussed how these system characteristics affect the choice of the reoptimization time-step and the horizon length. They concluded that the re-optimization time-step should be either a factor or a multiple of the time between orders, and re-optimizing more frequently than the greatest common factor of time-related data does not yield any additional benefit towards closed-loop performance. Further, horizon is a strong increasing function of load, time between orders, and order size max-mean relative difference. In the deterministic case, a longer horizon can be used to overcome infrequent re-optimization. Third, they evaluated the role of demand uncertainty and explored how it can be mitigated through better tuning of the online scheduling algorithm. They showed that a horizon sufficient for deterministic case is also sufficient for the case under demand uncertainty. The re-optimization time-step needs to be smaller than a threshold (dependent on the system characteristics, especially, η), else the least possible closed-loop cost cannot be achieved. To mitigate the effect of uncertainty, prompt re-optimization is beneficial and is the correct strategy as opposed to increasing the horizon length which results in no benefit. Finally, when the necessary re-optimization time-step and horizon length identified after tuning are impractical to achieve, they demonstrated how non-zero optimality gap, objective function modifications, and added constraints allow the use of a shorter horizon or a larger re-optimization time-step. 4. Simulation framework First, we discuss an adjustment model to restore feasibility when the observation of uncertainty makes the current schedule infeasible. Second, we discuss the appropriate probability distributions for the three uncertainties we focus on in this work. Finally, we propose a systematic closed-loop evaluation framework which can be used to better understand the effect of endogenous uncertainties and test different online scheduling strategies to mitigate them.

17

4.1. Adjustment model The adjustment model for resolving schedule infeasibilties is based on re-assignment on units and re-timing (start times) of the batches along with modifications to the their batch-sizes. This model is solved in a rolling horizon fashion at time-points in between the re-optimizations, and on observing infeasibility, its purpose is to simulate the limited adjustment moves that would be potentially carried out by a human on the production floor in lieu of a full re-optimization. There are three input decision variables in the scheduling model: (i) Wijt ; (ii) Bijt ; and (iii) Vkt . The others are state variables dependent on these input variables. Hence, it is sufficient to keep track of adjustments in these input variables and not the state variables V B (e.g., inventory (Skt )). We first define parameters, ΩW ijt , Ωijt , and Ωkt , to store the values of

variables Wijt , Bijt , and Vkt in the incumbent schedule, which is infeasible. Next, we define a binary variable Aijj 0 tt0 which, when 1, denotes that a batch of task i which was assigned to unit j and was supposed to start at time t, is now re-assigned to unit j 0 and its new start W which, when 1, denotes that a batch of task time is t0 . Further, we define binary variable Xijt

i, on unit j, with start time t will not be executed anymore. Every batch in the incumbent 0 0 W schedule (ΩW ijt = 1), which is executed (Xijt = 0) now on unit j and at time t , activates

the corresponding Aijj 0 tt0 variable through Eqs. 21 and 22 (see Fig. 4), while Wij 0 t0 is the activated task start binary for this adjusted batch. X

t0 ,j 0 ∈Ji

X

t,j∈Ji

W Aijj 0 tt0 = ΩW ijt − Xijt ∀i, j ∈ Ji , t

Aijj 0 tt0 = Wij 0 t0 ∀i, j 0 ∈ Ji , t0

(21) (22)

Similarly, we keep track of any change to the batch-sizes using Eqn. 23, where PijB0 t0 ∈ R is equal to the change in the batchsize (positive when increased and negative when decreased).

18

Figure 4: Illustration of re-assignment and re-timing of a batch through the adjustment model.

We find the absolute value of the change in batch-size (QB ij 0 t0 ), through Bij 0 t0 −

X

t,j∈Ji

B 0 0 ΩB ijt Aijj 0 tt0 = Pij 0 t0 ∀i, j ∈ Ji , t

(23)

B 0 0 QB ij 0 t0 ≥ Pij 0 t0 ∀i, j ∈ Ji , t

(24)

B 0 0 QB ij 0 t0 ≥ −Pij 0 t0 ∀i, j ∈ Ji , t

(25)

Finally, we measure the changes in shipments V Vkt − ΩVkt = Pkt ∀k, t

(26)

V QVkt ≥ Pkt ∀k, t

(27)

V QVkt ≥ −Pkt ∀k, t

(28)

The bounds on these new variables are W V B V Aijj 0 tt0 , Xijt ∈ {0, 1}, QB kt , Qkt ≥ 0, Pkt , Pkt ∈ {−∞, ∞}

(29)

In order to choose feasible new inputs (i.e., Wij 0 t0 , Bij 0 t0 , and Vkt ), we also make the original model equations (Eqs. 1 - 19) part of the adjustment model. Hence, one way to think of this model is to consider it equivalent to the original scheduling model but now with the input decisions constrained to execute only the already scheduled batches with (penalized) changes to their start-time, unit assignment, and batch-sizes. Note that only the

19

batches which are yet to start can be adjusted. The new objective function with added penalties for changes is

min

X k

! X X XXX F γkIN V ( Skt + SkT erminal ) + γkBO ( BOkt + BOkT erminal ) + αij Wijt t

+

X

t

J γjj 0 Aijj 0 tt0 +

ijj 0 tt0 :(j6=j 0 ,j,j 0 ∈Ji )

+

X

γijB0 QB ij 0 t0 +

i,j 0 ∈Ji ,t0

X

X

ijj 0 tt0 :(t6=t0 ,j,j 0 ∈J

j

γttT 0 Aijj 0 tt0

i∈Ij

t

(30)

i)

γkV QVkt

kt

J T 0 V where, in this work, for illustrative purposes, we use γjj 0 = $200, γtt0 = $10|t − t |, γk = $1,

and γijB0 = $10. For further details and guidance on how these penalty values are chosen, the reader is referred to Lee et al. (2019). 4.2. Uncertainty sampling We discuss the probability distributions appropriate for the three uncertainties we address in this work along with scaling of the distribution parameters with the sampling frequency for the purpose of simulations. In the discrete-time scheduling model, uncertainty can be sampled at every time-point. From the motivating examples, we learnt that infeasibilities can arise at the resolution of δ. Hence, we provide scaling relations for uncertainty sampled every δ h. 4.2.1. Task delays The occurrence of task delays can be modeled using the poisson distribution (Ross, 2009). This is because of the following practical and reasonable assumptions: (i) probability of single-period (δ) delay is finite and small, and probability of multi-period delays is much smaller than the probability of single-period delay; (ii) the occurrence and duration of non consecutive delays is independent of each other; and (iii) the probability distribution of a delay occurring at any time-point during the task is independent of the run status (n) of the task i.e. same probability of delay at the start, middle or end of the task. 20

The poisson distribution is: φdelay ∼ P(φdelay = l, λT D ) =

TD

e−λ

(λT D )l l!

where, l is an

integer, and λT D is the expected value of l. From this distribution, given λT D , we can sample for the delay (random variable φdelay ) at the current time-point. If φdelay = 0, there is no delay, when φdelay = 1, there is a delay of 1 time-interval, and when φdelay ≥ 2, the observed delay is a multi-period delay. To ascertain a value for λT D , we can use historical plant data to obtain the mean delays (µD i ). The sum of poisson random variables is also a poisson random variable with its mean as the sum of the means of the summed poisson random variables. Thus, sampling at every time-point for the duration of a batch, λT D =

δ D µ . τij i

The obtained delay from the poisson

distribution is then φdelay δ, or φdelay time-intervals. 4.2.2. Unit breakdowns A unit breakdown is typically observed when a batch is being executed on it, and results in suspension of the batch with the partially processed material inside it declared a loss. In general, for the purpose of short-term scheduling, it is a good assumption that breakdowns have a small and constant failure rate. Thus, these can be modeled using the exponential UBt

probability distribution: P(t, λU B ) = λU B exp−λ

(Lipson and Sheth, 1973). From past

operational data, we would typically know the mean time between failures (MTBF) for any particular unit, which is equal to 1/λU B (mean of exponential distribution). Given that we sample the unit breakdown random variable every δ h (when a batch is executed on it), the probability of breakdown within a time-interval is given by the cumulative exponential distribution function from 0 to δ: 1−exp−λ

UBδ

. In other words, the occurrence of a breakdown

follows the Bernoulli distribution with its parameter p = 1 − exp−λ

UBδ

.

When we observe a unit breakdown, there is typically a downtime (e.g., for repairs) associated with it. For simplicity, we assume a constant breakdown time of φbreak = 2δ in this work. In general, however, the value of φbreak can be sampled from an appropriate probability distribution.

21

4.2.3. Yield losses We assume that yield loss is rare, and its occurrence can be modeled using an exponential distribution (similar to unit breakdowns). The parameter λY L for the exponential distribution can be calculated as the reciprocal of the mean number of batches between yield losses YL

(MNBYL). The probability of observing a yield loss, at the end of a batch is 1 − exp−λ , irrespective of δ. In general, yield disturbance can be known even before the batch finishes. The corresponding model is discussed in detail in Gupta and Maravelias (2017). For simplicity, here we assume that yield disturbance is observed only at the end of a batch. When a yield loss is observed, we also need to find out how much loss occurred. In general, historical data can be approximated by an appropriate distribution for this. For simplicity, we model the yield loss magnitude as uniformly distributed between 0 and 0.1 fraction of the batch-size. 4.3. Closed-loop implementation We discuss how we conduct closed-loop simulations so as to draw conclusions on how the uncertain parameters and the online scheduling algorithm affect the closed-loop quality. We consider three sources of variation contributing to a change in the closed-loop quality: (i) change in the online scheduling algorithm; (ii) replications to account for uncertainty in breakdown, delay, yield-loss; and (iii) randomness in demand which is a key characteristic of short-term scheduling. Although, there is randomness, there is no uncertainty in demand. When varying the online scheduling algorithm or the uncertainty distribution, keeping demand samples unchanged reduces the inherent variation in closed-loop quality due to the randomness of demand. We discern statistical significance of the change in mean closed-loop quality using analysis of variance (ANOVA) (Wonnacott and Wonnacott, 1972), and decide the number of demand samples and the number of replications a priori. When the needed statistical significance is not achieved, we run additional simulations (replications). The individual steps are shown in Fig. 5 and described next. Once we choose the online 22

Figure 5: Steps for carrying out closed-loop simulations and the associated replications.

algorithm and a demand sample, every replication is a closed-loop run composed of multiple open-loop iterations (up to σmax ). In iteration σ = 0, we solve the scheduling model. Then we sample endogenous uncertainty which is observed as disturbances for the batches under execution at t = 0. For all units which are executing a batch, we first sample for unit breakdown. If a breakdown is not observed and there is no ongoing delay, we sample for delay. If no delay is observed and the batch is scheduled to finish, we finally sample for yield-loss. Thereafter, we advance the horizon and update the initial state (Eqs. 1 - 6) and, if observed, the disturbance parameters in the (new) current iteration from the previous iteration. If the start of the current iteration σ coincides with a time-point which is a 23

multiple of the re-optimization time-step (∆), we solve the scheduling model. Otherwise, we just allow the schedule to execute without revising it. If the schedule turns to be infeasible (due to disturbances), we solve the adjustment model to restore feasibility. When the closedloop run is finished, we move to the next replication (κ), and when we are finished with the replications (κmax ) for a given demand sample we start the replications for the next demand sample. We stop when we are finished with the replications for all demand samples. This gives us the closed-loop cost for the chosen uncertainty distribution and the online scheduling algorithm. Next, we change the uncertainty distribution or the online algorithm and repeat the steps to obtain new closed-loop cost. For studying online scheduling under only exogenous uncertainty (e.g., demand, with η < H), we can choose the uncertain parameter samples apriori, and need not carry out replications (i.e., κmax = 1). The need for replications is to account for decision dependence of the uncertain parameters. Interestingly, for type II endogenous uncertainty, in which the observation time is dependent on the decisions but not the distribution, we can choose the replication samples a priori as well. This can be achieved by using the sampling algorithm for the perfect information (η = H) benchmark discussed in Section 3.2.3 but with imperfect information (η < H). 5. Results Using the proposed simulation framework, we study the effect of the different endogenous uncertainties on closed-loop solution quality and how this effect can be potentially mitigated through better design of the online scheduling algorithm. We carry out these studies on the network shown in Fig. 1. We corroborate the results through a study on more networks and found qualitatively the same trends (see supplementary material). We model the randomness in demand through a symmetric triangular probability distribution with max-mean relative difference  = 10%. In the network of Fig. 1, the bottleneck unit is U3, with maximum possible production for each product material (assuming iden24

tical demand) equal to 10/(2 + 2 + 3) = 1.43 ton/h. This is the production capacity of the network. We consider the following three different demands for all product materials: 0.33, 0.66, 1 ton/h or equivalently (when divided by 1.43 ton/h) Λ = 0.23, 0.46, and 0.70. We model the uncertainty observation time (η), for task-delays, unit breakdowns, and yield-losses, as zero. We find, for the deterministic case, that using horizon H = 24, 30, and 36 h result in identical closed-loop cost. Therefore, for subsequent closed-loop simulations (for the uncertainty case) we use H = 24 h (for reasoning refer to Sec. 3.3). We verify that increasing the horizon, in the presence of task delays, yield-losses, and unit breakdowns, does not result in improvement in closed-loop cost. We decide on using 10 demand samples of each demand pattern and 10 replications for each demand sample. These overall 100 samples for each closed-loop setting give us statistically significant results (p-value = 0.05 for ANOVA). A total of 16,200 (1-week long closed-loop, δ = 1 h) simulations are carried out in parallel in a cluster of 24 Intel Xeon machines running CentOS 7 operating system. The needed 2,721,600 (open-loop) optimizations for the closed-loop simulations are solved to within 0.1% optimality using default solver options in CPLEX 12.6.1 via GAMS 24.4.3. In every panel of Figs. 6-8, the costs are an average over all the demand samples and the replications, and are scaled by a single chosen mean value (described in the caption) within each panel. Further, note that in Figs. 6-8, the schedule costs are different across different loads (keeping other parameters the same). Hence, a cost of 1, for example, in Fig. 6A is not the same as a cost of 1 in Fig. 6B or Fig. 6C. We also mark 95% confidence intervals for each (cost) data-point shown. First, we study the effect of task delays and the re-optimization time-step on closed-loop cost (see Fig. 6). For each task, we choose the Poisson distribution parameter λT D , i.e., the mean length of delay, as a fraction of its processing time τ . We see that as we increase the mean length of delay, the closed-loop cost increases. The increase in cost is substantial for large loads (Λ ≥ 0.7). This is because of the larger backlog costs associated with not meeting the due-date. As expected, when re-optimizing more often, it results in smaller increase in

25

A

3

28

Cost

Cost

2

5

1

1

1 0

0.1 0.2 0.3 0.4 0.5

=1 =3 =6

19 10

3

1.5

C

37

7

2.5

Cost

B

9

0

0.1 0.2 0.3 0.4 0.5

TD

0

TD

0.1 0.2 0.3 0.4 0.5 TD

Figure 6: Effect of task delays on closed-loop cost for different loads. Each data point is an average of 100 samples. The data is scaled, within each panel, by the corresponding data-point to ∆ = 1 h and λT D = 0. (A) Λ = 0.23; (B) Λ = 0.46; and (C) Λ = 0.70.

the cost due to uncertainty. Importantly, we observe that the difference in the costs for ∆ = 1 (full-optimization at each time-point) and ∆ = 3 is small. This implies that the adjustment algorithm is working well, and that suboptimal (due to the adjustment constraints) scheduling inputs are still acceptable as long as we re-optimize frequently. However, when ∆ = 6, the suboptimal inputs from the adjustment model are not revised for a longer time and this results in substantial increase in the closed-loop costs. A

2.5

11

Cost

Cost

1.5

3 2 1

1 0

0.1 0.2 0.3 0.4 0.5

C

16

4

2

Cost

B

5

=1 =3 =6

6 1

0

0.1 0.2 0.3 0.4 0.5

YL

YL

0

0.1 0.2 0.3 0.4 0.5 YL

Figure 7: Effect of yield losses on closed-loop cost for different loads. Each data point is an average of 100 samples. The data is scaled, within each panel, by the corresponding data-point to ∆ = 1 h and λY L = 0. (A) Λ = 0.23; (B) Λ = 0.46; and (C) Λ = 0.70.

Second, we study the effect of yield-loss on closed-loop cost (see Fig. 7). On the x-axis for these plots, we vary, λY L , the inverse of the mean number of batches between yield-loss (MNBYL) from infinity (no yield-loss) to yield-loss at every 2 batches (λY L = 1/2). Similar to the case of task delays, we find that larger loads magnify the effect of yield-loss. Because yield-loss is calculated as a fraction of the batch-size, for larger batches (when Λ ≥ 0.7), the absolute amount of yield-loss is larger, resulting in bigger backlog costs. Closed-loop

26

performance for ∆ = 1 h and ∆ = 3 h is close, but not so for ∆ = 6 h. A

5

61

13

2

7

1

1 0

2

4

6

100 UB

8

10

Cost

Cost

3

C

76

19

4

Cost

B

25

=1 =3 =6

46 31 16 1

0

2

4

6

8

10

100 UB

0

2

4

6

8

10

100 UB

Figure 8: Effect of unit breakdowns on closed-loop cost for different loads. Each data point is an average of 100 samples. The data is scaled, within each panel, by the corresponding data-point to ∆ = 1 h and λU B = 0. (A) Λ = 0.23; (B) Λ = 0.46; and (C) Λ = 0.70.

Third, we study how unit breakdowns affect closed-loop cost (see Fig. 8). We vary, λU B , the inverse of the mean time between failures (MTBF) from infinity (no breakdowns), down to every 10 h (λU B = 1/10 h-1 ). We can see that having more frequent breakdowns increases closed-loop cost and re-optimizing fast (∆ = 1 h) helps partially alleviate the increase in cost in comparison to slower re-optimization (∆ = 6 h). An important common observation in the above results is that the cost due to uncertainty rises sharply with load. This is because of two reasons: (i) the observation of these endogenous uncertainties is proportional to the production utilization; and (ii) higher demand, when not met on time due to uncertainty, results in larger backlog costs. Finally, we study what role horizon length has on mitigating costs associated with the above disturbances. We find, expectedly, that a longer horizon does not result in any benefit in reducing costs associated with disturbances that are observed at the start of the horizon (as a surprise). This is similar to how it is for demand uncertainty as showed in Gupta and Maravelias (2019). In general, choosing an appropriate re-optimization time-step (∆) and horizon length (H) for an application is a very important but non-trivial decision. As summarized in Section 3.3, Gupta and Maravelias (2019) discuss a systematic method to choose the necessary re-optimization time-step and horizon length.

27

6. Conclusions In this work, we proposed a framework for studying online scheduling in the presence of uncertainties whose observation time depends on optimization decisions; a common, but less studied type of uncertainty in scheduling. We specifically considered uncertainties in (i) processing times; (ii) batch yields; and (iii) unit operating status. First, we showed why uncertainty can result in infeasibilities and emphasized the importance of a good resolution mechanism for these infeasibilities. We illustrated how a systematic schedule adjustment, even if sub-optimal as compared to a full re-optimization, is more beneficial than simple heuristics, to stay feasible while waiting for new (optimal) scheduling inputs. Second, we proposed a new model for adjustment of the schedule against the observation of uncertainty and the consequent infeasibility in between the re-optimizations. In this model, we defined new variables to track and penalize input changes, so as to stay as close as possible to the incumbent infeasible schedule while obtaining an adjusted feasible schedule. Third, we discussed the different probability distributions for the three uncertainties that we consider in this work and the scaling of the parameters of these distributions with the sampling frequency. Fourth, we formalized the procedure for carrying out closed-loop simulations and evaluating closed-loop performance in the presence of these uncertainties. Finally, we drew insights from the closed-loop simulation results: (i) the observation of these endogenous uncertainties, and consequently the associated costs, increase when production is higher; (ii) higher demand (load) when not met on time due to uncertainty, results in larger backlog costs; (iii) re-optimization is beneficial in tackling these uncertainties and becomes more and more important as the uncertainty increases; and (iv) horizon longer than what is necessary to achieve good performance in the deterministic case, does not result in any benefit when the above uncertainties are present. To our knowledge, the work presented herein is the first systematic framework for studying online scheduling under endogenous uncertainties, and thus takes us one-step forward, towards synthesizing and testing general strategies to obtain high-quality closed-loop sched28

ules.

29

Nomenclature Indices/sets i∈I

tasks

j, j 0 ∈ J

units (equipment)

k∈K

materials

t, t0 , t00 ∈ T

time-points

n∈T

task status

I ⊇ Ij

tasks that can be carried out in unit j

I ⊇ I+ k

tasks producing material k

I ⊇ I− k

tasks consuming material k

J ⊇ Ji

units suitable for carrying out task i

K ⊇ KF

feed (raw) materials

K ⊇ KI

intermediates

K ⊇ KP

final products

Parameters Scheduling model F αij

fixed cost of running task i on unit j

P αij

proportional cost of running task i on unit j

M IN /β M AX βij ij

min/max capacity on batch-size of task i executed on unit j

γk

selling price of material k

γkIN V

inventory cost of material k

γkBO

backlog cost of material k

δ

discretization of time-grid; length of time-intervals

ζkt

incoming shipment of material k at time t

ξkt

demand for material k at time t

ρik /¯ ρik

mass-conversion coefficient (material consumption/production)

τij

processing time of task i on unit j

30

System characteristics η

uncertainty observation time



order size max-mean relative difference

Λ

load

Online algorithm ∆

re-optimization time-step

H

scheduling horizon

σ, σ max

online iteration number, maximum number of iterations

Disturbances n Z˙ ij

disturbance parameter denoting unit breakdown

B Z˙ n ij

batch-size of task suspended due to unit breakdown

Zˆjt

binary parameter, when 1, denotes unit j unavailable during time [t, t + 1)

n Y˙ ijn ,Yˆijt

single-/multi-period disturbance parameters denoting delay

B Y˙ n ,B Y ˆn ij ijt

single-/multi-period disturbance parameters denoting batch-size of a delayed task

βˆijkt

yield-loss of material k from task i that finished on unit j, at time t

Simulation framework κ, κmax

replication, maximum number of replications

µD i

mean delay in task i

φdelay , φbreak

delay, breakdown duration

λT D , λY L , λU B

probability distribution parameter for delays, yield loss, and breakdowns, respectively

Adjustment model ΩW ijt

values of Wijt variables in the incumbent schedule

ΩB ijt

values of Bijt variables in the incumbent schedule

ΩVkt

values of Vkt variables in the incumbent schedule

J γjj 0

penalty for re-assignment of a batch

31

T γtt 0

penalty for change in batch start time

B γij 0

penalty for change in batch-size

γkV

penalty on change in outgoing shipment

Continuous variables Bijt

batch-size of a batch of task i on unit j that started at time-point t

¯n B ijt

lifted batch-size

BOkt

backlog level of material k during time-interval (t − 1, t]

BOkT erminal

the terminal backlog level of material k, i.e., during time-interval (|T|, |T| + 1]

Skt

inventory level of material k during time-interval (t − 1, t]

SkT erminal

the terminal inventory level of material k, i.e., during time-interval (|T|, |T| + 1]

Vkt

outgoing shipment to meet demand for material k at time t

BX ijt

the batch-size of delayed batch with progress status n = 0

PijB0 t0

change in batch-size of adjusted task

QB ij 0 t0

absolute value of change in batch-size of adjusted task

V Pkt

change in outgoing shipment

QVkt

absolute value of change in outgoing shipment

Binary variables Wijt

when 1, denotes a batch of task i starts on unit j at time-point t

¯n W ijt

lifted task-start variables

Xijt

is 1, when there is a delay in a batch with progress status n = 0

W Xijt

equals 1, when an incumbent batch of task i on unit j at time t will not be executed

Aijj 0 tt0

when 1, denotes that a batch of task i which was assigned to unit j and was supposed to start at time t, is now re-assigned to unit j 0 and its new start time is t0

32

References Bassett, M. H., Pekny, J. F., Reklaitis, G. V., dec 1996. Decomposition techniques for the solution of large-scale scheduling problems. AIChE Journal 42 (12), 3373–3387. Blomer, F., Gunther, H.-O., mar 2000. LP-based heuristics for scheduling chemical batch processes. International Journal of Production Research 38 (5), 1029–1051. Burkard, R., Hatzl, J., jul 2005. Review, extensions and computational comparison of MILP formulations for scheduling of batch processes. Computers & Chemical Engineering 29 (8), 1752–1769. Calfa, B. A., Agarwal, A., Grossmann, I. E., Wassick, J. M., feb 2013. Hybrid BilevelLagrangean Decomposition Scheme for the Integration of Planning and Scheduling of a Network of Batch Plants. Industrial & Engineering Chemistry Research 52 (5), 2152–2167. Castro, P. M., Harjunkoski, I., Grossmann, I. E., feb 2011. Greedy algorithm for scheduling batch plants with sequence-dependent changeovers. AIChE Journal 57 (2), 373–387. Christian, B., Cremaschi, S., apr 2018. A branch and bound algorithm to solve large-scale multistage stochastic programs with endogenous uncertainty. AIChE Journal 64 (4), 1262– 1271. Chu, Y., You, F., 2014. Moving horizon approach of integrating scheduling and control for sequential batch processes. AIChE Journal 60 (5), 1654–1671. Colvin, M., Maravelias, C. T., 2010. Modeling methods and a branch and cut algorithm for pharmaceutical clinical trial planning using stochastic programming. European Journal of Operational Research 203 (1), 205–215. Colvin, M., Maravelias, C. T., dec 2011. R&D pipeline management: Task interdependencies and risk management. European Journal of Operational Research 215 (3), 616–628.

33

Cui, J., Engell, S., 2010. Medium-term planning of a multiproduct batch plant under evolving multi-period multi-uncertainty by means of a moving horizon strategy. Computers and Chemical Engineering 34 (5), 598–619. Elkamel, A., Mohindra, A., aug 1999. A rolling horizon heuristic for reactive scheduling of batch process operations. Engineering Optimization 31 (6), 763–792. Engell, S., 2009. Uncertainty, decomposition and feedback in batch production scheduling Sebastian. In: Computer-Aided Chemical Engineering. Vol. 26. pp. 43–62. Goel, V., Grossmann, I. E., sep 2006. A Class of stochastic programs with decision dependent uncertainty. Mathematical Programming 108 (2-3), 355–394. Grossmann, I. E., Apap, R. M., Calfa, B. A., Garc´ıa-Herreros, P., Zhang, Q., aug 2016. Recent advances in mathematical programming techniques for the optimization of process systems under uncertainty. Computers & Chemical Engineering 91, 3–14. Grossmann, I. E., Apap, R. M., Calfa, B. A., Garcia-Herreros, P., Zhang, Q., nov 2017. Mathematical Programming Techniques for Optimization under Uncertainty and Their Application in Process Systems Engineering. Theoretical Foundations of Chemical Engineering 51 (6), 893–909. Gupta, D., Maravelias, C. T., 2017. A General State-Space Formulation for Online Scheduling. Processes 5 (4), 69. Gupta, D., Maravelias, C. T., 2019. On the design of online production scheduling algorithms. Computers & Chemical Engineering submitted. Gupta, D., Maravelias, C. T., Wassick, J. M., 2016. From rescheduling to online scheduling. Chemical Engineering Research and Design 116, 83–97.

34

Gupta, V., Grossmann, I. E., nov 2011. Solution strategies for multistage stochastic programming with endogenous uncertainties. Computers & Chemical Engineering 35 (11), 2235–2247. Harjunkoski, I., Maravelias, C. T., Bongers, P., Castro, P. M., Engell, S., Grossmann, I. E., Hooker, J., M´endez, C. A., Sand, G., Wassick, J. M., 2014. Scope for industrial applications of production scheduling models and solution methods. Computers and Chemical Engineering 62, 161–193. Honkomp, S., Mockus, L., Reklaitis, G. V., 1999. A framework for schedule evaluation with processing uncertainty. Computers and Chemical Engineering 23 (4-5), 595–609. Huang, W., Chung, P. W. H., 2003. A constraint approach for rescheduling batch processing plants including pipeless plants. Computer Aided Chemical Engineering 14 (C), 161–166. Huercio, A., Espu˜ na, A., Puigjaner, L., 1995. Incorporating on-line scheduling strategies in integrated batch productioncontrol. Computers & Chemical Engineering 19 (95), 609–614. Janak, S. L., Floudas, C. A., Kallrath, J., Vormbrock, N., 2006. Production scheduling of a large-scale industrial batch plant. II. Reactive scheduling. Industrial and Engineering Chemistry Research 45 (25), 8253–8269. Jonsbr˚ aten, T. W., Wets, R. J.-B., Woodruff, D. L., 1998. A class of stochastic programs withdecision dependent random elements. Annals of Operations Research 82 (0), 83–106. Kanakamedala, K. B., Reklaitis, G. V., Venkatasubramanian, V., 1994. Reactive schedule modification in multipurpose batch chemical plants. Industrial & Engineering Chemistry Research 33, 77–90. Kelly, J. D., Mann, J., 2003. Crude oil blend scheduling optimization : an application with multimillion dollar benefits. Hydrocarbon Processing 82 (July), 47–54.

35

Kelly, J. D., Zyngier, D., nov 2008. Hierarchical decomposition heuristic for scheduling: Coordinated reasoning for decentralized and distributed decision-making problems. Computers & Chemical Engineering 32 (11), 2684–2705. Kim, M., Lee, I.-B., 1997. Rule-based reactive rescheduling system for multi-purpose batch processes. Computers & Chemical Engineering 21, S1197–S1202. Kondili, E., Pantelides, C. C., Sargent, R. W. H., 1993. A general algorithm for short-term scheduling of batch operations-I. MILP formulation. Computers and Chemical Engineering 17 (2), 211–227. Kopanos, G. M., Pistikopoulos, E. N., 2014. Reactive scheduling by a multiparametric programming rolling horizon framework: A case of a network of combined heat and power units. Industrial and Engineering Chemistry Research 53 (11), 4366–4386. Lagzi, S., Lee, D. Y., Fukasawa, R., Ricardez-Sandoval, L., aug 2017. A Computational Study of Continuous and Discrete Time Formulations for a Class of Short-Term Scheduling Problems for Multipurpose Plants. Industrial & Engineering Chemistry Research 56 (31), 8940–8953. Lappas, N. H., Gounaris, C. E., 2016. Multi-stage adjustable robust optimization for process scheduling under uncertainty. AIChE Journal 62 (5), 1646–1667. Lee, H., Gupta, D., Maravelias, C. T., 2019. Systematic Methods for Generating Alternative Production Schedules. Computers & Chemical Engineering (submitted). Lee, H., Maravelias, C. T., aug 2018a. Combining the advantages of discrete- and continuoustime scheduling models: Part 1. Framework and mathematical formulations. Computers & Chemical Engineering 116, 176–190. Lee, H., Maravelias, C. T., nov 2018b. Combining the advantages of discrete- and continuous-

36

time scheduling models: Part 2. systematic methods for determining model parameters. Computers & Chemical Engineering. Letsios, D., Misener, R., 2018. Exact Lexicographic Scheduling and Approximate Rescheduling. Arxiv preprint. Li, Z., Ierapetritou, M. G., 2008. Process scheduling under uncertainty: Review and challenges. Computers and Chemical Engineering 32 (4-5), 715–727. Lindholm, A., Giselsson, P., Quttineh, N.-H., Lidestam, H., Johnsson, C., Forsman, K., 2013. Production scheduling in the process industry. In: 22nd International Conference on Production Research. Lipson, C., Sheth, N., 1973. Statistical design and analysis of engineering experiments. McGraw-Hill, New York. Maravelias, C. T., 2012. General framework and modeling approach classification for chemical production scheduling. AIChE Journal 58 (6), 1812–1828. Maravelias, C. T., Grossmann, I. E., sep 2004. A hybrid MILP/CP decomposition approach for the continuous time scheduling of multipurpose batch plants. Computers & Chemical Engineering 28 (10), 1921–1949. M´endez, C. A., Cerd´a, J., 2003. Dynamic scheduling in multiproduct batch plants. Computers and Chemical Engineering 27 (8-9), 1247–1259. M´endez, C. A., Cerd´a, J., Grossmann, I. E., Harjunkoski, I., Fahl, M., 2006. State-of-the-art review of optimization methods for short-term scheduling of batch processes. Computers and Chemical Engineering 30 (6-7), 913–946. Nie, Y., Biegler, L. T., Villa, C. M., Wassick, J. M., apr 2015. Discrete Time Formulation for the Integration of Scheduling and Dynamic Optimization. Industrial & Engineering Chemistry Research 54 (16), 4303–4315. 37

Pantelides, C. C., 1994. Unified frameworks for optimal process planning and scheduling. In: Proceedings on the second conference on foundations of computer aided operations. Publications, Cache, New York, pp. 253—-274. Papageorgiou, L. G., Pantelides, C. C., 1996. Optimal campaign planning/scheduling of multipurpose batch/semicontinuous plants. 2. A mathematical decomposition approach. Industrial and Engineering Chemistry Research 35 (2), 510–529. Pattison, R. C., Touretzky, C. R., Harjunkoski, I., Baldea, M., 2016. Moving horizon closedloop production scheduling using dynamic process models. AIChE Journal. Qiao, F., Ma, Y., Zhou, M., Wu, Q., 2018. A Novel Rescheduling Method for Dynamic Semiconductor Manufacturing Systems. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 1–11. Rawlings, J. B., Mayne, D., 2009. Model Predictive Control : Theory and Design. Nob Hill Pub, Madison, WI, USA. Relvas, S., Barbosa-P´ovoa, A. P. F., Matos, H. A., mar 2009. Heuristic batch sequencing on a multiproduct oil distribution system. Computers & Chemical Engineering 33 (3), 712–730. Roe, B., Papageorgiou, L. G., Shah, N., may 2005. A hybrid MILP/CLP algorithm for multipurpose batch process scheduling. Computers & Chemical Engineering 29 (6), 1277– 1291. Ross, S. M., 2009. Introduction to Probability Models: Tenth Edition. Sabuncuoglu, I., Karabuk, S., 1999. Rescheduling frequency in an fms with uncertain processing times and unreliable machines. Journal of Manufacturing Systems 18 (4), 268–283.

38

Sahinidis, N., Grossmann, I., apr 1991. Reformulation of multiperiod MILP models for planning and scheduling of chemical processes. Computers & Chemical Engineering 15 (4), 255–272. Silvente, J., Kopanos, G. M., Pistikopoulos, E. N., Espu˜ na, A., oct 2015. A rolling horizon optimization framework for the simultaneous energy supply and demand planning in microgrids. Applied Energy 155, 485–501. Subrahmanyam, S., Kudva, G. K., Bassett, M. H., Pekny, J. F., jun 1996. Application of distributed computing to batch plant design and scheduling. AIChE Journal 42 (6), 1648–1661. Subramanian, K., Maravelias, C. T., Rawlings, J. B., 2012. A state-space model for chemical production scheduling. Computers & Chemical Engineering 47, 97–110. Sundaramoorthy, A., Maravelias, C. T., mar 2011a. A general framework for process scheduling. AIChE Journal 57 (3), 695–710. Sundaramoorthy, A., Maravelias, C. T., 2011b. Computational Study of Network-Based Mixed-Integer Programming Approaches for Chemical Production Scheduling. Industrial and Engineering Chemistry Research 50, 5023–5040. Touretzky, C. R., Harjunkoski, I., Baldea, M., dec 2016. Dynamic models and fault diagnosisbased triggers for closed-loop scheduling. AIChE Journal. Velez, S., Maravelias, C. T., 2013a. A branch-and-bound algorithm for the solution of chemical production scheduling MIP models using parallel computing. Computers & Chemical Engineering 55, 28–39. Velez, S., Maravelias, C. T., 2013b. Reformulations and branching methods for mixed-integer programming chemical production scheduling models. Industrial and Engineering Chemistry Research 52 (10), 3832–3841. 39

Velez, S., Maravelias, C. T., 2014. Advances in Mixed-Integer Programming Methods for Chemical Production Scheduling. Annual Review of Chemical and Biomolecular Engineering 5, 97–121. Velez, S., Sundaramoorthy, A., Maravelias, C. T., mar 2013. Valid Inequalities Based on Demand Propagation for Chemical Production Scheduling MIP Models. AIChE Journal 59 (3), 872–887. Vieira, G. E., Herrmann, J. W., Lin, E., 2003. Rescheduling manufacturing systems: a framework of strategies, policies, and methods. Journal of Scheduling 6, 39–62. Vin, J., Ierapetritou, M. G., 2000. A new approach for efficient rescheduling of multiproduct batch plants. Industrial and Engineering Chemistry Research 39 (11), 4228–4238. Wonnacott, T. H., Wonnacott, R. J., 1972. Introductory statistics. John Wiley and Sons. Wu, D., Ierapetritou, M. G., sep 2003. Decomposition approaches for the efficient solution of short-term scheduling problems. Computers & Chemical Engineering 27 (8-9), 1261–1276. Yee, K., Shah, N., mar 1998. Improving the efficiency of discrete time scheduling formulation. Computers & Chemical Engineering 22, S403–S410.

40

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: