A General State-Space Formulation for Online Scheduling

: We present a generalized state-space model formulation particularly motivated by an online scheduling perspective, which allows modeling (1) task-delays and unit breakdowns; (2) fractional delays and unit downtimes, when using discrete-time grid; (3) variable batch-sizes; (4) robust scheduling through the use of conservative yield estimates and processing times; (5) feedback on task-yield estimates before the task ﬁnishes; (6) task termination during its execution; (7) post-production storage of material in unit; and (8) unit capacity degradation and maintenance. Through these proposed generalizations, we enable a natural way to handle routinely encountered disturbances and a rich set of corresponding counter-decisions. Thereby, greatly simplifying and extending the possible application of mathematical programming based online scheduling solutions to diverse application settings. Finally, we demonstrate the effectiveness of this model on a case study from the ﬁeld of bio-manufacturing.


Introduction
Scheduling plays an important role in all industrial production facilities [1].Contingent on the scale of operation, optimization based scheduling methods can even achieve multi-million dollars increase in profits [2].Thus, considerable effort has been devoted towards developing optimization models that accurately represent the decision making flexibility in these facilities [3].Maravelias (2012) [4] provides a unified notation and a systematic framework for the description of chemical scheduling problems.Further, significant advances in solution methods, now enable us to solve small size scheduling problems.For example, a highly constrained scheduling instance over a network of 8 processing units, 19 tasks, and 26 materials, with a realistic scheduling horizon of 2 weeks, was shown to be solved to optimality in less than 1 min on an ordinary office computer [5].Thus, being able to generate and revise schedules in an online fashion, so as to account for new information and disturbances, is very much a reality now.The Dow Chemical Company has already adopted online scheduling in many of its production facilities [6][7][8].
Scheduling models, as have been developed till now, were not necessarily designed with an emphasis on being natively ready for implementation in an online scheduling setting.Thus, the online framework utilizing a model had to be tailored to that specific model, and required many ad-hoc (heuristic) adjustments to be able to represent and resolve a disturbance to the schedule [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].The introduction of the state-space idea to chemical production scheduling alleviated many of these issues that arose from having to make ad-hoc adjustments [24].
In this work, we present a generalized and extended state-space model which is suitable for implementation in an online scheduling setting.Further, we propose a new scheme for updating the state of the process, as well as an overall formulation to enforce constraints (through parameter/variable modifications), based on feedback information, on future decisions.Although here, we focus on expanding the modeling scope, it is important to point out that once a model has been adopted, there are still many other factors which influence the performance of an online scheduling method [25].These factors are: the online optimization horizon length, the re-computation trigger and its frequency if periodic, allowable changes from one online iteration to the next, any added constraints (e.g., terminal constraints), and the modeling of uncertainty (deterministic vs. stochastic optimization) [8].
This paper is structured as follows.In Section 2, we present a brief background on chemical production scheduling and discuss the state-space model of Subramanian et al. [24].In Section 3, we present a reformulated state-space model, based on a new convention, and showcase the generalizations on it one at a time.In Section 4, we present the final integrated model, with all generalizations present simultaneously, which requires more than simply concatenating all the individual generalizations together.Finally, in Section 5, we demonstrate the applicability of our proposed new model to a case study taken from the field of bio-manufacturing.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.

Background
In this section, we present the necessary background to be able to follow through the new general state-space model that we propose in this paper.Here, first, we layout the general problem statement, a standard problem representation framework, and briefly describe model classification, and solution methods.For a detailed discussion, the reader is referred to the following review papers [1,3,4,26].Second, we show a mixed integer linear programming (MILP) based widely adopted scheduling model.Third, we describe the typical state-space formulation adopted in model predictive control (MPC) technology.Finally, we provide a short overview of the state-space based scheduling model pioneered by Subramanian and co-workers [24,27,28].

General Problem Statement
The general scheduling problem can be stated as follows.Given: (i) Production facility data (e.g., unit capacities and connectivity), (ii) Production recipes (e.g., processing times and mixing rules), (iii) Production costs (e.g., material holding costs), (iv) Material availability (e.g., raw materials delivery amounts and dates), (v) Resource availability (e.g., maintenance schedule and utility levels), and (vi) Production targets or orders with due-times; scheduling seeks to find: Number and the associated processing-sizes of the needed tasks, (ii) Assignment of these tasks to processing units, and (iii) Timing (or just the sequence) of these tasks on the assigned units; so as to meet production targets at minimum cost, or to maximize profit if production beyond the given target is allowed.Apart from minimization of cost, or maximization of profit, the objective can also be minimization of makespan, or minimization of earliness, or any other suitable objective for the considered application.In general, several processing characteristics and constraints could also be present such as sequence dependent changeovers, setup times, storage constraints, time-varying utility costs, etc. [1,3].

Problem Representation
Before a scheduling problem can be solved, we need an abstract framework to represent the different elements of the problem, viz., the production facility, the associated production recipe, etc.The state task network (STN) enables this representation [29].Under this representation, tasks are carried out on units (equipment), and they transform materials (states) from one to another.Apart from the material to be processed and the equipment to process these materials on, these tasks can also require resources, such as, utilities, manpower etc.Another popular framework, is the resource task network (RTN) [30].In contrast to STN, in which materials, units, and utilities are treated as different from one another, in RTN, these are treated at par, all termed together as resources.We use the STN representation in this paper, but the general modeling ideas presented are also easily adaptable to the RTN representation.
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 + k /I − k ; task i consumes/produces material k equivalent to ρ ik / ρik mass fraction of 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 I j ; The processing time of task i, when executed on unit j, is denoted by τ ij .On any given unit, only one task can be performed at a time with its batch-size between lower (β min ij ) and upper capacities (β max ij ); the associated fixed and proportional production costs of carrying out task i on unit j are α F ij and α P ij respectively.Feed, intermediate, and product materials are denoted by k ∈ K F /K I /K P ; there are possible incoming deliveries (ζ kt ) and outgoing orders (ξ kt ) at certain times for selected materials; the selling price, inventory cost, and backlog cost of material k are γ k , γ I NV k , and γ BO k , respectively.In Figure 1, we see a process network's STN representation comprising of 4 material nodes (circular) labeled M0-M3 and 4 task nodes (rectangular) labeled T1-T4.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 capacities (β min /β max ) are also shown here.Material prices (γ) are shown adjacent to the material nodes.

Model Classification
Scheduling models can be classified on the basis of (i) optimization decisions; (ii) modeling elements; and (iii) modeling of time [4].Models that employ a time-grid are either continuous time or discrete time models.In discrete time models, the fixed time-grid spacing is denoted by δ.Events can take place only at these grid time-points.Thus, all time-related parameters are rounded in a conservative direction, such that the resulting schedule computed using these new parameter values is feasible even for the original parameter values.Hence, processing times and raw material delivery dates are rounded up, while due dates are rounded down so as to match with an integer multiple of δ.
Even though, having a discrete time-grid introduces the above approximation error, discrete time models have several advantages over continuous time-grid models.For example, accounting for utility consumption, inventory and backlog costs, time varying prices, or time-dependent resource availability introduces non-linearities in continuous time models, but not so in discrete time models [26].Furthermore, discrete time models are, in general, at least as effective as continuous time models, and in fact are better suited for large scale problems with several additional processing features [31].In this work, we employ a discrete time-grid for our state-space model.

Scheduling MILP Model
The discrete time STN MILP scheduling model modified from Shah et al. (1999) [59] comprises of Equations ( 1)- (6).Time is represented by index t ∈ T. Binary variable W ijt , when 1, implies task i is starting on unit j at time t.Variable B ijt ∈ [β min ij , β max ij ] denotes its batch-size.The assignment constraint (Equation ( 1)) ensures only one task can be executed on a unit at a time.
Equation ( 2) ensures that the batch-size of a task, if initiated, is within its upper and lower bounds.
S kt , which is the variable denoting inventory of material k during time-period (t − 1, t], is calculated in Equation (3) as a balance of production/consumption and outgoing (V kt )/incoming (ζ kt ) shipments.
Equation ( 4) couples the outgoing shipment variable V kt with demand, ξ kt , for material k at time t.Backlog variables, BO kt , denote pending demand during time-period (t − 1, t], and are penalized in the cost minimization objective function (Equation ( 5)).
Finally, the domain of all the variables is restricted via Equation ( 6):

Standard form of State-Space Models
State-space model formulations have been useful, alongside frequency domain models, in process control [60][61][62][63][64]. Now, as optimization based control and economic MPC are becoming the new standard, state-space models have become ubiquitous [65,66].In the most general form, a state-space based model can be written as dx dt = f (x, u, d); where x are the states, u are the manipulated inputs, and d are the disturbances.The function f (•) is not theoretically restricted to the class of linear functions, but is typically approximated as linear due to computational tractability considerations.The linear difference equation form for f (•) yields the model as: where, A, B, and B d are state-space matrices and t is the index for time.The states x need not be associated with a physically identifiable entity in the plant.Some can have a direct physical meaning, while others can be artificial (e.g., augmented) constructs so as to enable the modeling exercise.
The output (measurements y) is related to the states and inputs as y = h(x, u), where h(•) can be non-linear, but is typically linear (e.g., y(t) = Cx(t) + Du(t), where C and D are coefficient matrices).
The control optimization model has to follow the plant physical constraints and any other imposed constraints due to operational strategy (e.g., for environmental concerns) or those that enable better closed-loop properties (e.g., economics and stability).These constraints, when linear, can take the general form: where, E x , E u , and E d are the coefficient matrices of the states, inputs, and disturbances, respectively.If there are any equality constraints, these can also be represented as two opposite inequality constraints, so as to conform to the general form (Equation (8)).For example, the following constraints are equivalent: Thus, any equality constraints that we propose from here on, can be easily converted to the general inequality form through the use of the above trick.Finally, the objective function takes the form: ,...,u(N−1) V N (x(0), u(0), x(1), u(1), ..., x(N − 1), u(N − 1)) (10) where N is the number of discrete time-points in the online optimization horizon.
A wealth of literature focuses on the closed-loop properties of the aforementioned iterative control methods, with novel and most recent results, specifically, in presence of discrete inputs, discussed in Rawlings and Risbeck (2017) [67].

Scheduling State-Space Model
Motivated by process control approaches, Subramanian et al. (2012) [24] proposed a state-space model (Equations ( 5) and ( 11)-( 16)) for the chemical production scheduling problem.For brevity, we present the formulation for constant batch-sizes ( There are two distinct features of this model.First, the "complete status" of the plant can be interpreted solely from the variables (states) at that moment in time.This is made possible by lifting past actions/inputs (the task start binary variables, W ijt ) which have a lagged effect on the "current status" of the plant.Second, observed uncertainties are treated as disturbances, and represented as parameters in the model equations.These two features, together, allow for the model to be kept identical in each online scheduling iteration without any ad-hoc adjustments (due to observation of uncertainty).Thus, the model is in "online ready" form.In addition, due to the use of the state-space formulation, which is popular in process control models, this model also happens to be a very suitable candidate for integration of scheduling and control [68].
To enable lifting of inputs, new task-states (variables) Wn ijt are defined.Although this increases the number of variables in the model, it is matched by an equal increase in the number of equations (the lifting equations, Equations ( 11) and ( 12)).Thus, no new degrees of freedom are introduced.When the task starts, n is zero (n = 0, Equation ( 11)), and when the task finishes, n equals the processing time of the task (n = τ ij ).To express task delay and unit breakdown disturbances, new parameters Ŷn ijt and Ẑn ijt are defined, respectively.Ŷn ijt , when 1, denotes a delay of δ h in task i during time-period [t − δ, t), where δ, as defined in Section 2.1.3,is the granularity of the discrete time-grid.Ẑn ijt , when 1, denotes break-down of unit j while executing task i during time-period [t − δ, t).For ease of presentation, we assume from here on that δ = 1 h.
In the absence of delays or breakdowns, the lifting equations effectively represent the relation: The lifted variables are defined only till n = τ ij , because a "look-back" beyond that value of n is not needed.The effect of past inputs, for n > τ ij , is already, indirectly, contained in the inventory and backlog variables S kt and BO kt .The lifted states, Wn ijt , are augmented to the future states (see Figure 2).Here, τ ij for the tasks is assumed to be 3. Lifting of past inputs enables knowing the complete status of the plant by looking at the states (variables) only at that moment in time.In the absence of delays or breakdowns, the lifting equations effectively represent the relation: Arrows show which variables are equal due to the lifting equations (Equations (11) and (12), with no delays or breakdowns).Variables in green or red have a value of 1, rest have value 0. Information is carried over from one iteration to the next through the update step (Equations ( 17)-( 19)).
In the assignment constraint (Equation ( 13)), parameters Ŷτ ij ijt and Ẑτ ij ijt are included, to ensure that the unit appears to be busy, and no new tasks can be started, when there is a delay or breakdown observed at a time when a task is about to finish.Additionally, for multi-period breakdowns, parameter Ẑτ IT,j IT,jt is made 1, where IT is a fictitious "idle task", with τ IT,j = 1, that keeps the unit busy through the duration of the multi-period breakdown.
In inventory balance (Equation ( 14)), βC ijkt and βP ijkt are parameters that denote material handling loss during consumption and production of material k, respectively.When a delay or breakdown is observed at the end of a task, the terms Ŷτ ij ijt and Ẑτ ij ijt , which are subtracted from Wτ ij ijt , prevent erroneous multiple counting of the material amount produced by that task.
Finally, Equation (16) shows the bounds on the variables present in the model.
Next, we describe the online update step, i.e., how information is carried over from one online iteration to the next.Since the scheduling horizon is advanced by 1 h (the model is kept identical), the state at t = 0 (initial condition) for the next iteration is matched with the state at t = 1 of the previous iteration.This is shown in Figure 2, and achieved through the online "update equations" (Equations ( 17)-( 19)), in which σ denotes the iteration number.Variables σ S k(t=0) , σ BO k(t=0) , and σ Wn ij(t=0) for n ≥ 1 which represent lifted task-states, are assigned fixed values through the update step.But, σ W0 ij(t=0) , which represents degrees of freedom to start new tasks at t = 0, is not fixed.This is identical to how the online updates are performed for the no disturbance case [25].However, since, here we are dealing with the case where disturbances can be present, the disturbance parameters ( ξkt , Ŷn ijt , Ẑn ijt , βP ijkt , and βC ijkt ) are also assigned appropriate values to reflect the observed disturbances.However, these parameters do not participate in the update equations.These influence the prediction of states, for t ≥ 1, in the online iteration σ.
Figure 3A,B, show the evolution of task-states when a 2 h delay is observed right after a task starts and just before a task is about to finish, respectively.The 2 h duration of this multi-period delay is known immediately in iteration σ.However, the model formulation also does allow for representing the observation of consecutive, possibly independent, 1 h single-period delays, one at a time in succeeding iterations.These collectively, in hindsight, appear to be a single multi-period delay, but are actually not.
It is quite evident, that n, now in the presence of delays, loses its physical meaning of denoting how much progress has been made on the task.For example, for the task in Figure 3A (iteration σ), due to the 2 h delay, the task-states evolve as W1 All that can be said now is that when W0 ijt = 1, the task has just started at time t, and when Wτ ij ijt = 1 and Ŷτ ij ijt = 0 simultaneously, then the task has finished.As we will show in Section 3.1, we overcome this limitation by introducing a new convention to map observed disturbances to the disturbance parameters, and hence, are able to preserve the physical meaning of n, even when disturbances are present (see Figure 4).When a 2 h delay is observed, through the lifting equations, Wn ijt evolve over the green trajectory, leading to the task correctly finishing 2 h late in iteration σ.Here, τ for the task is 3. Arrows show which variables are enforced as equal by the lifting equations (Equations ( 11) and ( 12), with delays present).Variables and parameters in green have a value of 1, rest have value 0. (A) The task now finishes at t = 4, instead of at t = 2. (B) The task now finishes at t = 2, instead of at t = 0. Through Equation ( 13), the unit is kept busy at t = 0 and 1, by the inclusion of the terms Ŷτ=3 and Ŷτ=3 (t=1) , and hence a new task is prevented from starting at these times.In addition, these terms in Equation ( 14), prevent the task's produce from erroneously contributing to inventory (S k,(t=1) and S k,(t=2) ).
Figure 5A,B, show the evolution of task-states when a breakdown just before t = 0 is observed and is known to have a 2 h unit downtime (for repairs), right after a task starts and just before a task is about to finish, respectively.Given the observation of breakdown, we would expect intuitively, and unlike what is shown in Figure 5A,B, that none of the task-states Wn ijt are active for the green task at t = 0. We show how this is achieved through the new model discussed in Section 3.1 (see Figure 6).When a 2 h delay is observed, through the lifting equations, Wn ijt evolve over the green trajectory, leading to the task correctly finishing 2 h late in iteration σ.Here, τ for the task is 3. Arrows show which variables are enforced as equal by the lifting equations (Equations ( 22)-( 24), with delays present).Variables and parameters in green have a value of 1, rest have value 0. (A) The task now finishes at t = 4, instead of at t = 2. Due to the update step, parameters Ẏ, Ŷ0 0 and variable X 0 are 1, hence, in iteration σ, due to the optimization model, W0 0 , X 1 , and W0 1 are also 1. (B) The task now finishes at t = 2, instead of at t = 0.The update step ensures that the true progress, n, of the task is reflected in the task-states at t = 0, i.e., n = 2.

Modeling Generalizations
In Section 3.1, we present a new state-space model formulation that differs, from the state-space model of Subramanian et al. (2012) [24], in the convention that is followed for mapping observed disturbances to the disturbance parameters.Although both models are accurate, this new convention ensures that the task-states, in the presence of disturbances, follow a more intuitive notation.Specifically, the meaning of n as the progress of a task, is maintained.In addition, we define several new parameters to systematically account for disturbances.
In Section 3.2, we show how to handle fractional delays and unit downtimes (due to unit breakdowns).In Section 3.3, we expand the scope of the model to account for variable batch-sizes.Thereafter, in Sections 3.4-3.9we present generalizations that can be applied to the state-space model, one at a time.Afterwards, in Section 4, we present the final model equations with all generalizations present simultaneously.As we will see in that section, for all the generalizations to work in the presence of each other, a few more modifications are necessary.

New Basic Formulation
The new state-space model relies on a comprehensive update step of the task-states, in between the online iterations, to promptly reflect the delays and breakdowns in the task-states.The inventory and backorder update stay the same (Equations ( 18) and ( 19)) as in the model of Subramanian et al. (2012) [24].The task-states update is modified from Equation (17) to Equations ( 20) and (21).
The parameters Ẏn ij which, if 1, represent a 1 h delay in task with progress status n.Note the dot (•) instead of the hat (∧) on the symbols of these parameters.Since these parameters are exclusively for the update step, and do not directly participate in the optimization model, these need not be indexed by time-neither t (iteration σ − 1) nor t (iteration σ).Similarly, Żn ij denotes a breakdown of unit j on which task i with progress status n was running.X ijt is a new binary variable, defined for all time-points, which, when 1, captures the information about delays in a task with progress status n = 0, i.e., when the task gets delayed right after it starts.The use of this variable, in Equations ( 22) and ( 23), will become clear when we discuss the optimization model.We also define a new parameter Λjt , which, when 1, denotes the unit is unavailable for the time-period [t, t + 1).This parameter participates in Equation (25). Figure 5.When a breakdown is observed, further evolution of the task-states for the task (the green trajectory), running on the unit that broke, stops.Here, τ for the task is 3 and the unit downtime (blue) is 2 h.Arrows show which variables are enforced as equal by the lifting equations (Equations ( 11) and ( 12), with breakdown present).Variables and parameters in green, blue, or red have a value of 1, rest have value 0. The green task is suspended at t = 0.A new task (red) can only start at t = 2, once the unit downtime is over.(A) Through Equation ( 13), the unit is kept busy at t = 0 and 1, due to the terms W1 t=0 and Ẑ1 IT,1 , respectively.(B) Through Equation ( 13), the unit is kept busy at t = 0 and 1, due to the terms Zτ=3 t=0 and Ẑτ=1 IT,1 , respectively.Additionally, the term Ẑτ=3 (t=0) in Equation ( 14), prevents the green task-state ( Wτ ij =3 t=0 ) from erroneously contributing to the inventory.
Figure 6.When a breakdown is observed, further evolution of the task-states for the task (the green trajectory), running on the unit that broke, stops.Here, τ for the task is 3 and the unit downtime (blue) is 2 h.Arrows show which variables are enforced as equal by the lifting equations (Equations ( 22)-( 24)).
Variables and parameters in green, blue, or red have a value of 1, rest have value 0. The green task is suspended at t = 0.A new task (red) can only start at t = 2, once the unit downtime is over.Through Equation ( 25), the unit is kept busy at t = 0 and 1, by the inclusion of the terms Λt=0 and Λt=1 .(A) The parameter Ż0 , through Equation ( 20), prevents the task-state from evolving from The parameter Ż2 , through Equation ( 20), prevents the task-state from evolving from (σ−1) W2 0 to σ W3 0 .
For a multi-period delay of φ h, in task i on unit j, in addition to Ẏn ij , parameters Ŷn ijt are activated for t = 0, t = 1, ..., t = φ − 2. For unit breakdowns with downtime duration of φ h, in addition to Żn ij , parameters, Λjt are activated for t = 0, t = 1, ..., t = φ − 1.Thus, single-period delays do not result in activation of any Ŷn ijt parameters, but single-period breakdowns require activation of Λj(t=0) .Having described the update step, we now describe the optimization model.In this model, the lifting equations consist of Equations ( 22)- (24).
When there is a φ h multi-period delay in a task with progress n = 0, the update step assigns X ij(t=0) = 1 and Ŷ0 ijt = 1 ∀t ∈ {0, 1, ..., φ − 2}.This ensures that W0 ijt stays activated for next (φ−1) h, but with W ijt = 0, i.e., the task is not erroneously interpreted as a new task start.If there are no delays, then, through Equations ( 21) and ( 22), X ijt = 0 and any new task that starts with W ijt = 1, through Equation (23), results in W0 ijt = 1.Equation ( 23) is a constraint that we impose on the inputs (W ijt ) given the states (X ijt and W0 ijt ), and if needed can be converted to inequality form through use of Equation (9).Variables X ijt are either fixed (t = 0) by the update step or are equated to the delay parameters in the optimization model, hence, can be declared as free variables with no explicit bounds.The assignment constraint (Equation ( 25)) includes the parameter Λjt to account for unit downtime.Additionally, it contains the variable W0 ijt on the left-hand side, and not variable W ijt , to correctly account for the unit being busy, specifically, when a delay in a task with progress n = 0 is observed.
The inventory balance, Equation (26), in contrast to Equation ( 14), does not require any corrective delay or breakdown terms.This is because, for any task, the states W ijt (task-start) and Wτ ij ijt (task-end), even if delays or breakdowns are observed, are active only at most once.

Fractional Delays and Unit Downtimes
In Figures 4 and 6, we showed the cases where delays and unit downtime are integer multiples of time-grid spacing δ.Additionally, the unit breakdown was assumed to take place at almost the time-point t, i.e., very close to an integer multiple of δ.However, if δ is not very small, then these assumptions may not be good.Given any fractional delays (π delay ), downtimes (π down ), or unit breakdown time (π break ), we need an appropriate scheme for the (online iterations) update step, to ensure realistic rounding of these to integer values, so as to keep the task-finish and unit-availability times, in sync with the discrete time-grid.A single task can have multiple separate delays, hence, we index the delay time with index r (recurrence), i.e., π r delay .A breakdown, however, can occur only once, at π break , following which, the unit downtime, π down , starts.
For the first delay, a rounded up value is applied in the update steps, i.e., the delay is assumed to be π 1 delay /δ .For every additional ψ th delay, the difference, φ = (∑ r=1 π r delay )/δ , dictates how much additional, integer φ, delay is applied in the update steps.Figure 7A shows a numerical example for fractional delays.
When a unit breaks down, the parameter Żn is always activated, so as to suspend the running task.The key challenge is to identify, for how many next time-points the unit would be unavailable.This dictates, if, and how many, Λt parameters are activated.This is done as follows.On breakdown, the unit becomes unavailable from π break to π break + π down (in iteration σ, π break < 0).Hence, all Λt that span integer multiple of δ, t ∈ (π break , π break + π down ] are activated.This also means, if (π break − π break /δ δ + π down ) < δ, none of the Λt are activated, i.e., the unit breaks down and comes back online before the immediate next time-point.This is illustrated in Figure 7B.

Variable Batch-Sizes
To account for variable batch-sizes, we define variables, B ijt which denotes the batch-size of the task that just starts, Bn ijt for lifted task batch-size states, and B X ijt to represent batch-size of task that is delayed with progress status n = 0. We define parameters B Ẏij and B Żn ij that participate in the update steps and these denote the batch-size of the task delayed and suspended due to unit break-down, respectively.Further, we define parameter B Ŷn ijt for the optimization model.Here, for ease of discussion, the time-grid is global, i.e., it is not reset for each iteration.(A) A delay (shown as oblique green pattern) of 0.66 h in time-period (0, 1) is observed.Since, 0.66 = 1, a delay of 1 h is applied at t = 1.This would ensure that the task finish is aligned with t = 4, even if the task actually ends at t = 3.66.The horizontal green pattern represents the fictitious extra task runtime to align with the discrete time-grid.Next, another delay of 0.2 h is observed in time-period (1, 2).Since, 0.66 + 0.2 − 0.66 = 0, no additional delay is applied in the update steps.This makes sense, because, the task finishes at t = 3.86 in reality.Since, the previous delay was applied as 1 h, the task is now still thought to finish at t = 4, in alignment with the time-grid.Finally, when another delay of 0.66 h is observed in time-period (2,3).Since, 0.66 + 0.2 + 0.66 − 0.66 + 0.2 = 1, a 1 h delay is applied in the update step.This correctly ensures that the task is now thought to end at t = 5, which is the round up of the true end time of t = 4.52 h.(B) A break-down is observed at t = 0.2.Thus, Ż0 = 1 for the update step between the iterations starting at t = 0 and t = 1.If the downtime (blue) is 0.66 h, then the unit actually becomes available at t = 0.86.Thus, Λ(t=1) is not activated.This is indeed the case from our mathematical procedure as well, since, t ∈ (0.The additional update steps, Equations ( 27) and ( 28), due to variable batch-sizes, are as follows: The optimization model now requires Equations ( 29)- (31) for lifting the batch-size: It might appear in Equation ( 30) that when B X ijt > 0, nothing prevents B ijt from also erroneously taking on a positive value.This was not an issue in Equation ( 23) because the W ijt and W0 ijt variables there were binary.However, Equation (2) ensures that B ijt can only take a non-zero value when W ijt = 1.Since, through Equation ( 23), W ijt = 0, whenever X ijt = 1, B ijt also takes value 0. The update steps ensure that X ijt and B X ijt can only be non-zero simultaneously.
The inventory balance (Equation ( 32)) now incorporates the new batch-size variables, B ijt and Bn ijt , rather than the task-state binary variables (W ijt and Wn ijt ) which was the case in Equation (26).
Remark: In principle, we can completely avoid defining the new parameters B Ẏij , B Żn ij , B Ŷijt , and variable B X ijt by reformulating Equations ( 27)- (31), so as to only use parameters Ẏij , Żn ij , Ŷijt , and variable X ijt .For example, Equation ( 31) can be reformulated to Equation (34).
Since, Equation ( 34) entails the multiplication of variables with parameters, by itself, it is an acceptable alternate linear formulation.However, when task termination is allowed as a scheduling decision (Section 3.7), this reformulation results in bi-linear terms which are undesirable.Thus, we indeed define the new parameters B Ẏij , B Żn ij , B Ŷijt , and variable B X ijt , and use Equations ( 27)-( 31) in their native form without the simplifying reformulation discussed in this remark.

Robust Scheduling: Batch-Sizes
In many applications, it can be prudent to schedule batches bigger than what are needed to just satisfy the nominal demand.This can be, for example, due to the possibility of seeing a demand spike, or to pro-actively compensate for typical material handling losses when a batch finishes.To do so, the parameter ρik in material inventory balance (Equation (35b)) can be substituted by a scaled down value ( ρr ik ), where ρr ik < ρik .This results in bigger batches starting, since the model now under-predicts the yield of materials from any given batch-size.In order to, however, correctly account for the actual inventory resulting from the finishing of a task, the nominal value of ρik is used at t = 0, along with any yield-loss or material handling loss disturbance (Equation (35a)).As it can be seen in Figure 8, which is a simple illustration of this modeling generalization, as the iterations progress, a task-finish-state eventually hits t = 0, yielding the large yield proportionate to the true (nominal) value of ρik .Now, if there are any material handling losses ( βP ijk(t=0) ), they can be subtracted from the true yield (in Equation (35a)).It is worth noting that, although βP ijkt and βC ijkt are defined for all time-points, they are possibly active only at t = 0, if the corresponding uncertainty is observed.Hence, these parameters can be, in principle, dropped from Equation (35b).
We can write the above two equations, compactly together, as follows: . In iteration σ − 1, the green task (τ = 2) with a "large batch" is finishing at t = 1, but, due to the use of ρr k , is anticipated to produce only half of what the demand is.Thus, another identical green task is scheduled to start at t = 4, to satisfy the demand.With the use of ρr k , if the demand could be still be satisfied with a single "large" batch, a second batch wouldn't be scheduled.In iteration σ, the earlier green task yields a large amount of material, in line with its large batch-size due to the true value of ρik used at t = 0. (A) Since, here, there was no material handling loss, thus the second green task need not run now, as the demand was satisfied by the first batch itself.(B) Although the green task results in a large yield at t = 0, a small material handling loss (β P k(t=0) < 0) at t = 0, requires the second green task to still be scheduled in order to meet the demand, but now with a smaller batch-size.If there are no further material handling losses, there would be a small excess inventory of material produced by the second green task, since its batch-size was decided assuming the yield to be lower ( ρr k ) but would in reality by higher ( ρk ).

Robust Scheduling: Processing Times
Uncertainty in the processing times is very common in scheduling [69].A popular approach to proactively manage this uncertainty is to robustify the schedule by adding a delay buffer to each task's processing time [70].Once this robust schedule has been computed, it is advantageous to adjust it online, by taking into account the feedback on actual finish times of the tasks [71].In discrete-time models, this has been typically done using ad-hoc adjustments in between the online iterations.To the best knowledge of the authors, there is not yet a systematic way to be able to naturally handle this adjustment within an optimization model.
We show here how we can extend the state-space model to produce robust schedules, from a processing time point of view, and yet seamlessly allow for tasks to finish, after they have been running for their nominal processing times plus the delays.We define a new parameter τ r ij , which denotes the conservative processing time of the tasks (τ r ij > τ ij ).Thereafter, we modify the lifting (Equation ( 24) modified to Equations (38a), (38b) and ( 31) modified to Equations (39a) and (39b)), assignment (Equation ( 25) modified to Equations (40a) and (40b)), and inventory balance equations (Equation (32) modified to Equations (41a) and (41b)), such that the nominal value of processing times (τ ij ) is employed at t = 0, and the conservative value (τ r ij ) is employed for t > 0. No other model or update equations are modified.An illustration is given in Figure 9.For the green task (τ = 2), a delay buffer of 1 h is chosen, i.e., τ r ij = τ ij + 1. (A) In iteration σ − 1, the model predicts that the green task will finish at t = 2.In the next iteration σ, when there is no delay in the green task, having run for 2 h, it correctly finishes at t = 0.No ad-hoc model changes are required in between the iterations to make this possible.Consequently, to take advantage of this nominal finish time of the green task, the red task, now starts 1 h earlier than what it was scheduled for in the previous iteration.(B) The green task, with conservative processing time, is scheduled to finish at t = 2 in iteration σ − 2. In iteration σ − 1, a delay of 1 h is observed.The processing time delay buffer is maintained, and now the task is anticipated to finish at t = 2.In the next iteration (σ), the task correctly finishes at t = 0, accounting for the 1 h delay.The red task can now start on time at t = 0, or equivalently t = 2. Thus, having a delay buffer was useful, since it predicted a realistic start time for the red task in iteration σ − 2 itself.
The lifting equations, Equations (38b) and (39b), contain variables W n ij(t=1) , B n ij(t=1) ∀j, i ∈ I j , n ∈ {τ ij + 1, τ ij + 2, ..., τ r ij }, which are not coupled back to any variables at t = 0.These variables must be fixed to zero, otherwise, the optimization can assign these variables a spurious value so as to erroneously generate inventory (in Equation (41b) through variables B τ r ij ijt ∀j, i ∈ I j , t ∈ {1, 2, ..., τ r ij − τ ij }).We can write the above equations, compactly together, as follows: where Finally, please note that in this approach, as can be seen in Figure 9B (iteration σ − 1), an a priori fixed buffer time (τ r ij − τ ij ) is added to the task duration, irrespective of whether delays have already been observed (during the execution of the current task).It can be argued that the buffer was meant to absorb the actual delays, and hence, should be cut back for tasks that actually get delayed.This is a fair critique.However, owing to feedback, the true task-finish is accounted for (see Figure 9A iteration σ), and there is no wasted equipment time due to the unused buffer, which otherwise results in idle time in static robust scheduling approaches.

Feedback on Yield Estimates
The material handling loss ( βC ijkt / βP ijkt ) disturbance parameters can be used to represent yield losses as well, but these parameters are assigned a value only at t = 0, i.e., when a task actually ends and a material handling loss is observed.In many applications, the actual yield of a task can, in fact, be estimated during the task's execution, and does not always come as a surprise when the batch finishes.To incorporate the information about anticipated yield loss in future, these parameters can be assigned values for t > 0. However, this information has to be then carried over from one iteration to the next, with corresponding decrement in the t index of these parameters to reflect the shifting time-grid, till the task finishes.In addition, if the task is delayed, these parameter values also have to be delayed.This requires cumbersome mechanisms to handle this information in between the online iterations.
Adapting the state-space model to lift the yield-loss information forward provides a much more natural way to handle this feedback.Thus, we define a new free variable Ln ijt , which is analogous to the batch-size variable Bn ijt , but, as we can see in Equation ( 47), when the task finishes, instead of adding to, it subtracts from the inventory.
Thus, in a way, this new variable can be thought to be the negative counterpart of the batch-size variable.
The yield loss can also be delayed by use of the parameters L Ẏn ij and L Ŷn ijt so as to stay in sync with the task-finish time, or nullified in case of a unit breakdown through the parameter L Żn ij .This is achieved using the update steps (Equations ( 48) and ( 49)) and the model lifting equations (Equations ( 50)-( 52)).
λn ij is a parameter the denotes the "additional" yield loss observed (anticipated) in that iteration (σ).When delays and yield losses are observed simultaneously, then the parameters L Ẏn ij and L Ŷn ij take the value (σ−1) Ln−1 ij(t =0) + λn−1 ij , i.e., the total yield loss up to and including in iteration σ.
Please note that there are no L ijt variables in Equation (51).If these variables were to be present, they would have to be fixed to zero.Otherwise the optimization itself can spuriously assign non-zero values to these variables, such that the intended true batch-size is B ijt − L ijt .

Task Termination
In the update step (Equations ( 20), ( 21), ( 27) and ( 28)) we made an implicit assumption that past decisions (task-states with n > 0 at t = 0) are fixed.In general, more decisions from the previous iteration can be considered fixed or the deviation from them penalized in the current iteration.This has been suggested in the literature to reduce schedule nervousness [23,72,73].
However, to the best knowledge of the authors, no model to date considers canceling/terminating tasks already underway, as an optimization decision.This is surprising, given that whenever is needed, this would be a natural decision for a human scheduler.For example, to prioritize processing for a new rush order, an unfinished process unrelated to this rush order may need to be terminated.Other routine possibilities for task termination are excessive delays or yield losses, as is commonly the case in bio-manufacturing [74,75].This decision to terminate should be an outcome of the online optimization so as to best react to the observed disturbances or new information.We define task termination as a "willful" decision to discontinue a task.This is in contrast with task suspension, which is a "forced" discontinuation of a task as a result of a unit breakdown or loss of utility support.Since, preemption is not customary in chemical processes, we assume a total loss of output of a task that is terminated.This is in agreement with how the output of task suspensions is treated.
This termination of tasks can be achieved by "softening" the initial conditions (task-states at t = 0) in an online iteration.We introduce a new binary variable, T n ij , which, when 1, denotes termination of task i, on unit j, which has run-index n.Since, all variables values for iteration σ − 1 are now a parameter for iteration σ, we can write the following linear equations to achieve this softening: Hence, the update equations (Equations ( 20), ( 21), ( 27) and ( 28) modified to ( 53)-( 56)) now become part of the model.The update equations for inventory (Equation ( 18)) and backlog (Equation ( 19)) stay unmodified, and are not softened.
If we had no disturbances, just softening the initial task-states would have sufficed.However, when we have delays, we have to also ensure that we appropriately nullify the effect of delay parameters, which have been already assigned a value at the start of an iteration (optimization).Thus, wherever the delay parameters Ŷn ijt appear, we multiply these with (1 − T n ij ).Since the coefficients of the (1 − T n ij ) terms are the delay parameters, when there is no delay, these terms are also, consequently, absent.Parameters Ẏτ ij ij , B Ẏτ ij ij , Ŷτ ij ijt , and B Ŷτ ij ijt are always zero, hence, variables T τ ij ij do not participate in the model.A unit breakdown implicitly disrupts a task, hence, in such a situation, the question of termination does not arise.Overall, the lifting equations are modified to Equations ( 57)- (60) with Equations ( 23) and (30) remaining unchanged.
We do not index variable T n ij with time (t), because it serves no purpose to terminate a task in future (t > 0).If a task is already under execution and has to be terminated in the future for some reason, it is always better to terminate it right away (t = 0).
In addition to including a cost associated with task termination, ∑ j ∑ i∈I j ∑ in the objective, we can also enforce a pre-specified unit downtime (τ T j ) following every task termination, by including a summation term in the assignment constraints for that many time-points: We can write the above two equations, compactly together, as follows: where An added advantage of the compact form, is that the unit downtime can now be a function of the task that was terminated, i.e., τ T ij is indexed by i as well, and not just j.To systematically account for unit downtime resulting from task-termination in a previous iteration, i.e., if (σ−1) T n ij = 1 ∀j, i ∈ I j , n ∈ {1, 2, 3, ..., τ ij − 1}, the parameter Λjt has to be activated for t ∈ {0, 1, 2, ..., τ T j − 2} in iteration σ.Thereafter, in each subsequent iterations, the downtime is decremented by 1, and the parameter Λjt appropriately activated for the corresponding time-points.Alternatively, we can define a new binary variable and lift it, to keep the unit deactivated for the remaining downtime in subsequent iterations.The new variable, which would be subtracted on the right hand side of Equations (61a) and (61b), in lieu of Λjt , can be thought of as the unit unavailability variable.Kondili et al. (1993) [29] proposed a formulation for "hold" tasks, which are dummy tasks that can be used to model storage of materials in a processing unit, while waiting to be unloaded.This is especially important for production facilities that follow a no intermediate storage (NIS) policy for certain materials.In the network shown in Figure 1, task T4 is a hold task.The purpose of this task is to keep material M3 residing in unit U3.So as to ensure that the hold task can only store material in a unit that it was originally produced in, we write the state-space version of the original constraint proposed by Kondili et al. (1993) [29] in Equation (64).

Post-Production Storage in Unit
where I HOLD is the set of hold tasks.When multiple materials are produced in a unit, here we assume that the corresponding hold task emulates the simultaneous storage of all these materials in the unit.The ρHOLD,k mass-coefficient then dictates in what proportion are the materials released from the unit.If only a certain individual material is held, and the others unloaded, the term Bτ ij ijt , and the constraint is written only for that material k which is held.Further, for a perishable material, ρHOLD,k < ρ HOLD,k ; that is, a fraction of the material is lost (perishes or deactivates) when held.

Unit Capacity Degradation and Maintenance
In many processes, such as polymerization reactions, or purification processes, it is not uncommon for the unit capacity to "degrade" after a task has been processed on that unit [76][77][78][79][80][81].This can be, for example, due to residue formation (e.g., scaling) or impurity accumulation (e.g., membrane pore blockages).Some degradation is gradual and predictable, while some degradation may occur suddenly and unexpectedly.
To model unit degradation, we define a new non-negative variable, C ijt , which denotes the capacity of the unit j, to perform task i that starts at time t.For an un-degraded unit, C ijt is initialized to the value β max ij .This variable value is passed over from one iteration to another, through the update equation: A new disturbance parameter, μij , which when negative, represents extent of sudden (unexpected) partial loss in unit capacity.Conversely, a positive value represents renewal of unit capacity.This positive value could be, for example, a result of installing a new repaired unit in place of an old unit that broke down.
Through Equation ( 66), we define a balance on the unit capacity.ρ C ii j is a parameter that denotes the gradual degradation in capacity of unit j to perform task i, due to execution of task i on that unit.

ρ C
ii j is either negative or zero.
The degraded unit can be typically restored to its full capacity through a maintenance task (e.g., cleaning), which we denote with the abbreviation MT.We define I MT as a set of maintenance tasks, and J MT as a set of units which can degrade, and consequently, need a corresponding maintenance task.Further, we add the maintenance task (MT) to the set of tasks I j that can be performed on unit j.The value of Mτ MT,j ijt , the variable which we define below, dictates the restored capacity due to completion of a maintenance task.
Like any conventional task, the maintenance task has a start binary (W MT,jt ) associated with it, which is appropriately lifted.The assignment equation (Equation ( 25)) ensures that the maintenance task can only run when the unit is not running any other task.Since the maintenance task does not consume or produce materials, the conventional batch-size of this maintenance task, B MT,jt is fixed to zero.Instead, we define a new type of batch-size, Mn ijt , specific to the purpose of this task, which we term as the maintenance-size.Since we assume that only one kind of a maintenance task exists for every unit, the index i in this maintenance-size variable is not MT.This variable denotes, how much capacity to perform task i is restored, when maintenance is performed on unit j.This variable is also lifted, similar to the batch-size variable, and can be delayed or suspended due to breakdowns.Similarly, the parameters M Ẏn ij , M Ŷn ijt , and M Żn ij denote the maintenance-sizes of the maintenance task corresponding to capacity restored to perform task i.The update equations associated with this variable are: The model equations, for lifting, are: The batch-size of new tasks is upper bounded by the unit capacity variable (Equation ( 72)).This ensures that only smaller batches can now be processed if C ijt < β max ij .If a task just finishes, the restored or degraded capacity due to that is also accounted for while upper bounding the task-batchsize B ijt .
We define the maintenance task with fixed processing time (τ MT,j ), and assume that whenever performed, the unit is restored to its full capacity (i.e., C ijt = β max ij ).This requires the maintenance-size, M ijt , to be the difference between the deteriorated unit capacity (C ijt ), before the maintenance starts, and the upper capacity limit (β max ij ), to which the unit has to be restored to.This is achieved using Equation ( 73): Finally, we specify the lower and upper bounds for variable C ijt : and define a new parameter α M ij , which denotes the proportional cost of maintenance of a unit.The cost term, in the objective, for a maintenance task is

Integrated Model
In this section, we present the complete model with all generalizations present simultaneously.For brevity, we write index σ in only those equations, where σ − 1 is also present.Everywhere else, all variables are those of the current iteration σ.
Variables X ijt , B X ijt , L X ijt , M X ijt are free variables, however, through the update and model equations, they are always equated to a parameter value, hence, are not degrees of freedom.The decision variables, or in other words-the inputs from a systems perspective [24], are W ijt , B ijt , T n ij , M ijt , and V kt .The objective is:

Case Study
Bio-manufacturing is a type of manufacturing in which molecules of interest, such as, metabolites, drugs, enzymes, etc., are produced through the use of biological systems such as living micro-organisms [82].Several commercial sectors, rely on these, for example, pharmaceuticals, food and beverage processing, agriculture, waste treatment, etc.Furthermore, there is an increasing thrust towards finding biological routes for production of bulk chemicals [83].The use of live systems in bio-manufacturing, however, introduces several operational challenges.These include batch-to-batch variability, parallel growth of both, the desired product as well as undesired toxic byproducts in the same batch, and possible random shocks that can lead to complete failure of a batch [75].This makes it an interesting area for application of scheduling methods.In this section, we present an example, motivated from bio-manufacturing, to demonstrate all modeling generalizations discussed in Section 3, using the integrated model equations outlined in Section 4.
In general, bio-manufacturing processes can be divided into an upstream bio-reaction (e.g., fermentation) stage and a downstream purification stage [74].The upstream stage typically consists of two steps: cell culture preparation in the lab and the bio-reaction.These two steps are task T1 and T2, respectively, in the network in Figure 1.The downstream purification stage typically consists of three steps: centrifugation, chromatography, and filtration.Among these three steps, chromatography takes the longest and the chromatograph columns are prone to unpredictable failures.Hence, we assume that chromatography, being the dominant step, is representative of the complete purification stage.This is task T3 in the network in Figure 1.Overall, we choose this simplified system (network) for an easier illustration of the capabilities of our general state-space model.
In our example, as features, we allow for possible delays in the cell culture preparation (T1), and small yield losses in the bio-reaction (T2), including possible substantial yield losses due to sudden cell death.Thus, we carry out robust scheduling, using a conservative processing time (τ r ) for task T1 and a conservative mass-conversion coefficient ( ρr ), against small yield losses, for task T2.Further, in the downstream stage, we assume that the chromatograph column (U3) loses capacity with usage, part of which is predictable.Executing a maintenance task can restore capacity on the chromatographs.To enforce the no-intermediate storage (NIS) policy for material M2, the inventory variable S M2,t is fixed to zero for all time-points.Raw material, M0, is assumed to be available in an unrestricted supply, as needed.Selected instance parameter values, other than the ones already shown in Figure 1, are outlined in Table 1.An order for 15 kg of M3 is due at the 14th hour of the day.To meet this order, online scheduling is carried out, with a horizon of 16 h, and re-optimization every 1 h, starting at the 0th hour.All optimizations are solved to optimality using default solver options in CPLEX 12.6.1 (IBM Corporation, North Castle, NY, USA) via GAMS 24.4.3 (GAMS Development Corporation, Fairfax, VA, USA), installed on an Intel Xeon (E5520, 2.27 GHz, 8 core processor) machine (Intel Corporation, Santa Clara, CA, USA), with 16 GB of RAM and Linux CentOS 7 operating system (Red Hat Inc., Raleigh, NC, USA).The schedules obtained in selected online scheduling iteration are shown in Figure 10.For the remaining online iterations, the predicted schedule is identical to the respective previous iteration (but with time-grid shift).
The nominal makespan, without any disturbances or robustification, to meet this order is 13 h.However, in iteration 0 (see Figure 10), T1 is started at t = 0, instead of t = 1, since a conservative processing time, τ r = 3 h is in use.Further, the batch-sizes for T1 and T2 are 16.67 kg, since a conservative yield parameter ( ρr = 0.9) is in use for T2.This predicts production of 15 kg of M2, which through the two T3 batches, of 6.25 kg and 8.75 kg, makes 15 kg of M3.
In the next iteration, a fractional delay of 0.5 h is observed.Due to the use of a discrete time-grid with granularity δ = 1 h, and being the first delay for this task, this is rounded up to 1 h.Since, now the order is predicted to be late, the batch-sizes of T3 are revised so as to meet as much of the order as possible on time (β max = 10).Going forward, no more delays are observed in T1, hence, it finishes at t = 0 in iteration 3.This is because the nominal τ is in use at t = 0 in Equations ( 85) and (86).Consequently, the downstream tasks are all scheduled earlier now, matching up with the initial predicted schedule in iteration 0. Thus, the conservative processing time was useful, in making T1 start earlier at 0th hour.
In iteration 4, due to sudden cell death, 90% yield loss in T2 is observed (anticipated at task finish).Hence, T2 is terminated.T1 is restarted (with conservative delays are after nominal processing time of 2 h.Thus the start times of the downstream tasks are pulled forward by 1 h. In iteration 9, 10% yield loss is observed (anticipated at task finish).Since, this is not substantial, T2 is not terminated.Instead, a new train of tasks is scheduled to start at t = 2 to compensate for the lost yield.In iteration 11, when T2 finishes, it results in nominal yield ( ρ) minus the 10% anticipated yield loss.Hence, 15 kg of M2 is produced.Thus, the new train of tasks previously scheduled, but not yet started, in iteration 9 are canceled.
In iteration 11, in addition to the gradual decline in chromatograph capacity due to usage (ρ C T3,T3,U3 < 0), which does not affect starting the next 5 kg T3 task, a sudden loss in capacity is observed ( μT3,U3 = −8).Hence, a maintenance task (MT) is scheduled with maintenance-size M T3,T3,t=1 = 9.34.Once this maintenance is over, the pending task T3 can start.No further disturbances are observed.Consequently, the order is fully met at the 19th hour.
If task termination, conservative processing time, and conservative yield were not used, the order would have been fully met only at the 27th hour (Gantt chart not shown).Hence, using the general state-space model enabled a richer set of decision making, resulting in an overall better schedule.

Conclusions
We developed a general state-space model, particularly motivated by an online scheduling perspective, that allows modeling (1) task-delays and unit breakdowns with a new, more intuitive convention over that of Subramanian et al. (2012) [24]; (2) fractional delays and unit downtimes, when using discrete-time grid; (3) variable batch-sizes; (4) robust scheduling through the use of conservative yield estimates and processing times; (5) feedback on task-yield estimates before the task finishes; (6) task termination during its execution; (7) post-production storage of material in unit; and (8) unit capacity degradation and maintenance.Further, we propose a new scheme for updating the state of the process, as well as an overall formulation to enforce constraints (through parameter/variable modifications), based on feedback information, on future decisions.We demonstrate the effectiveness of this model on a case study from the field of bio-manufacturing.Through this new state-space model, we have enabled a natural way to handle routinely encountered processing features and disturbance information in online scheduling.The general features that we address are found in several industrial sectors, namely, pharmaceuticals, fine chemicals, pulp and paper, agriculture, steel production, oil and gas, food processing, bio-manufacturing, etc.The proposed model, therefore, greatly extends and enables the possible application of mathematical programming based online scheduling solutions to diverse application settings.Finally, it is important to note, that although here we presented the model using STN based representation, these generalizations can also be adapted to RTN based representation.

Nomenclature
Indices/sets i ∈ I tasks j ∈ J units (equipment) k ∈ K materials t ∈ T time-points/periods I ⊇ I j tasks that can be carried out in unit j I ⊇ I +

S kt
inventory level of material k during period (t − 1, t] T n ij binary variable, when 1, denotes termination of task i, with run-status n, on unit j V kt outgoing shipment to meet demand for material k at time t W ijt binary variable, when 1, denotes task i starts on unit j at time-point t Wn ijt lifted task-start variables X ijt when 1, captures the information about delays in a task with progress status n = 0 B X ijt the batch-size of delayed task with progress status n = 0 L X ijt yield-loss of delayed task with progress status n = 0 M X ijt maintenance-size of delayed maintenance task with progress status n = 0

Figure 1 .
Figure 1.STN representation of a process network.This network, for a bio-manufacturing process, is described in detail in Section 5.It has three steps (tasks) for production of a pharmaceutical ingredient (material M3).T1 is the task of preparing the cell cultures in lab-size beakers.T2 denotes the task of having the cell culture grow and produce the pharmaceutical active ingredient, on a feed of sugars, in a bio-reactor.T3 is a purification task, which is carried out in chromatograph columns.Finally, T4 is a dummy task to model storage of material M2 inside unit U2. ρ 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.

Figure 2 .
Figure2.Task-states are shown for two online iterations -numbered σ − 1 and σ.Each iteration uses its own local time-grid which is reset to start from 0. Here, τ ij for the tasks is assumed to be 3. Lifting of past inputs enables knowing the complete status of the plant by looking at the states (variables) only at that moment in time.In the absence of delays or breakdowns, the lifting equations effectively represent the relation:Wn ijt = W ij(t−n) ∀j, i ∈ I j , n.Arrows show which variables are equal due to the lifting equations (Equations (11) and (12), with no delays or breakdowns).Variables in green or red have a value of 1, rest have value 0. Information is carried over from one iteration to the next through the update step (Equations (17)-(19)).

Figure 3 .
Figure 3.When a 2 h delay is observed, through the lifting equations, Wnijt evolve over the green trajectory, leading to the task correctly finishing 2 h late in iteration σ.Here, τ for the task is 3. Arrows show which variables are enforced as equal by the lifting equations (Equations (11) and (12), with delays present).Variables and parameters in green have a value of 1, rest have value 0. (A) The task now finishes at t = 4, instead of at t = 2. (B) The task now finishes at t = 2, instead of at t = 0. Through Equation (13), the unit is kept busy at t = 0 and 1, by the inclusion of the terms Ŷτ=3

Figure 4 .
Figure 4.When a 2 h delay is observed, through the lifting equations, Wnijt evolve over the green trajectory, leading to the task correctly finishing 2 h late in iteration σ.Here, τ for the task is 3. Arrows show which variables are enforced as equal by the lifting equations (Equations (22)-(24), with delays present).Variables and parameters in green have a value of 1, rest have value 0. (A) The task now finishes at t = 4, instead of at t = 2. Due to the update step, parameters Ẏ, Ŷ0 0 and variable X 0 are 1, hence, in iteration σ, due to the optimization model, W0 0 , X 1 , and W0 1 are also 1. (B) The task now finishes at t = 2, instead of at t = 0.The update step ensures that the true progress, n, of the task is reflected in the task-states at t = 0, i.e., n = 2.

Figure 7 .
Figure7.The task has a nominal processing time τ = 3. Time-grid spacing is δ = 1 h.Here, for ease of discussion, the time-grid is global, i.e., it is not reset for each iteration.(A) A delay (shown as oblique green pattern) of 0.66 h in time-period (0, 1) is observed.Since, 0.66 = 1, a delay of 1 h is applied at t = 1.This would ensure that the task finish is aligned with t = 4, even if the task actually ends at t = 3.66.The horizontal green pattern represents the fictitious extra task runtime to align with the discrete time-grid.Next, another delay of 0.2 h is observed in time-period (1, 2).Since, 0.66 + 0.2 − 0.66 = 0, no additional delay is applied in the update steps.This makes sense, because, the task finishes at t = 3.86 in reality.Since, the previous delay was applied as 1 h, the task is now still thought to finish at t = 4, in alignment with the time-grid.Finally, when another delay of 0.66 h is observed in time-period(2,3).Since, 0.66 + 0.2 + 0.66 − 0.66 + 0.2 = 1, a 1 h delay is applied in the update step.This correctly ensures that the task is now thought to end at t = 5, which is the round up of the true end time of t = 4.52 h.(B) A break-down is observed at t = 0.2.Thus, Ż0 = 1 for the update step between the iterations starting at t = 0 and t = 1.If the downtime (blue) is 0.66 h, then the unit actually becomes available at t = 0.86.Thus, Λ(t=1) is not activated.This is indeed the case from our mathematical procedure as well, since, t ∈ (0.2, 0.86] does not include any integer time-point.If π down = 1.5 h, then t ∈ (0.2, 1.7] does span t = 1, and consequently Λ(t=1) = 1.Finally, if π down = 2.25 h, then t ∈ (0.2, 2.45] spans t = 1 and t = 2, which results in Λ(t=1) = 1 and Λ(t=2) = 1.

Figure 9 .
Figure 9.For the green task (τ = 2), a delay buffer of 1 h is chosen, i.e., τ r ij = τ ij + 1. (A) In iteration σ − 1, the model predicts that the green task will finish at t = 2.In the next iteration σ, when there is no delay in the green task, having run for 2 h, it correctly finishes at t = 0.No ad-hoc model changes are required in between the iterations to make this possible.Consequently, to take advantage of this nominal finish time of the green task, the red task, now starts 1 h earlier than what it was scheduled for in the previous iteration.(B) The green task, with conservative processing time, is scheduled to finish at t = 2 in iteration σ − 2. In iteration σ − 1, a delay of 1 h is observed.The processing time delay buffer is maintained, and now the task is anticipated to finish at t = 2.In the next iteration (σ), the task correctly finishes at t = 0, accounting for the 1 h delay.The red task can now start on time at t = 0, or equivalently t = 2. Thus, having a delay buffer was useful, since it predicted a realistic start time for the red task in iteration σ − 2 itself.

k
tasks producing material k I ⊇ I − k tasks consuming material k I ⊇ I HOLD hold (storage) tasks I ⊇ I MT maintenance tasks J ⊇ J i units suitable for carrying out task i J ⊇ J MT units which can degrade, and consequently, need a corresponding maintenance taskK ⊇ K F feed (raw) materials K ⊇ K I intermediates K ⊇ K P final products Parameters α F ij fixed costof running task i on unit j α P ij proportional cost of running task i on unit j α T ij cost of terminating task i on unit j α M ij proportional cost of maintenance task i on unit j β ij fixed batch-size of task i executed on unit j β min ij /β max ij min/max capacity on batch-size of task i executed on unit j βP ijkt / βC ijkt material unloading/loading loss during production/consumption of material k γ k selling price of material k γ I NV k inventory cost of material k γ BO k backlog cost of material k δ discretization of time-grid; length of time-periods ζ kt incoming shipment of material k at time t Żn ij disturbance parameter denoting unit breakdown B Żn ij batch-size of task suspended due to unit breakdown L Żn ij yield-loss size of task suspended due to unit breakdown M Żn ij maintenance-size of maintenance task suspended due to unit breakdown θ τ ijt dummy parameter as defined in Equation (46) θ ρ ikt dummy parameter as defined in Equation (37) θ T ijt dummy parameter as defined in Equation (63) λn ij yield loss in task i running on unit j, with run status n Λjt binary parameter, when 1, denotes unit j unavailable during time [t, t + 1) μij when negative, represents extent of sudden partial loss in unit capacity ξ kt demand for material k at time t ξkt demand disturbance for material k at time t π break actual (fractional) time at which a unit breaks down π down fractional downtime in a unit π delay fractional delay in a task π r delay r th delay in a task ρ ik mass-conversion coefficient (material consumption) ρik mass-conversion coefficient (material production) ρr ik conservative mass-conversion coefficient for production ( ρr ik < ρik ) ρ C ii j deterioration in unit capacity to perform task i, due to performing task i on that unit j. σ online iteration number τ ij processing time of task i on unit j τ r ij conservative processing time of task i (τ r ij < τ ij ) on unit j τ T j task independent unit j downtime after terminating a task τ T ij task dependent unit j downtime after terminating task i Ẏn ij , Ŷn ijt single-/multi-period disturbance parameters denoting delay B Ẏn ij , B Ŷn ijt single-/multi-period disturbance parameters denoting batch-size of a delayed task L Ẏn ij , L Ŷn ijt single-/multi-period disturbance parameters denoting yield-loss size of a delayed task M Ẏn ij , M Ŷn ijt single-/multi-period disturbance parameters denoting maintenance-size of delayed maintenance task φ duration of delay or breakdown, in multiples of δ ψ recurrence count of delay for a task Variables B ijt batch-size of task i on unit j Bn ijt lifted batch-size BO kt backlog level of material k during period (t − 1, t] C ijt capacity of unit j to perform task i during period (t − 1, t] Ln ijt lifted yield-loss variables M ijt maintenance-size of the maintenance task Mn ijt lifted maintenance-size

Table 1 .
Parameter values for the case study.