An efficient MILP framework for integrating nonlinear process dynamics and control in optimal production scheduling calculations

Computers and Chemical Engineering 110 (2018) 35–52

Contents lists available at ScienceDirect

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

An efficient MILP framework for integrating nonlinear process dynamics and control in optimal production scheduling calculations Morgan T. Kelley a , Richard C. Pattison a,1 , Ross Baldick b , Michael Baldea a,∗ a b

McKetta Department of Chemical Engineering, The University of Texas at Austin, Austin, TX 78712, United States Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, TX 78712, United States

a r t i c l e

i n f o

Article history: Received 18 August 2017 Received in revised form 22 November 2017 Accepted 22 November 2017 Available online 28 November 2017 Keywords: Integrated scheduling and control Nonlinear dynamics Surrogate models Hammerstein–Wiener models Finite step response models Lagrangian relaxation

a b s t r a c t The emphasis currently placed on enterprise-wide decision making and optimization has led to an increased need for methods of integrating nonlinear process dynamics and control information in scheduling calculations. The inevitable high dimensionality and nonlinearity of first-principles dynamic process models makes incorporating them in scheduling calculations challenging. In this work, we describe a general framework for deriving data-driven surrogate models of the closed-loop process dynamics. Focusing on Hammerstein–Wiener and finite step response (FSR) model forms, we show that these models can be (exactly) linearized and embedded in production scheduling calculations. The resulting scheduling problems are mixed-integer linear programs with a special structure, which we exploit in a novel and efficient solution strategy. A polymerization reactor case study is utilized to demonstrate the merits of this method. Our framework compares favorably to existing approaches that embed dynamics in scheduling calculations, showing considerable reductions in computational effort. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Competitive pressure from global market forces has placed a heightened emphasis on information exchange, coordination, and integration of all decision-making layers of the chemical supply chain. Significant developments in this area, supported by advances in computer hardware, data exchange, storage, and optimization algorithms, have already led to substantial economic benefits for chemical operations (Grossmann, 2005). This coordination often extends to inclusion of power grid and power supply networks such as distributed energy systems (Diangelakis and Pistikopoulos, 2017) into the operation and control of chemical systems. Two essential layers in the decision-making hierarchy of a chemical enterprise are production scheduling and process control. The interface between scheduling and control represents, in effect, an interaction between business-driven decisions (scheduling) with situation- and safety-driven decisions (control) (Baldea and Harjunkoski, 2014). Integrating these activities is therefore key in maximizing operational profits by meeting demand and ensuring that production targets (i.e., the setpoints transmitted to the

∗ Corresponding author. E-mail address: [email protected] (M. Baldea). 1 Current address: Apeel Sciences, Goleta, CA 93117, United States. https://doi.org/10.1016/j.compchemeng.2017.11.021 0098-1354/© 2017 Elsevier Ltd. All rights reserved.

control system) are met safely in the presence of disturbances and operational uncertainties. Integrating scheduling and control decisions has been shown to improve operations for several industrial entities, including, polymer and metal production, wastewater treatment, air separation, and energy storage systems (Engell and Harjunkoski, 2012; Touretzky and Baldea, 2014; Touretzky et al., 2016). Integrating scheduling and control becomes particularly important when a chemical plant operates in fast-changing markets (e.g., markets with real-time electricity pricing). In order to operate more profitably under such circumstances, the plant must be able to quickly change production rates or product grades, taking advantage of excess production capacity and product storage facilities when available. This is reflected in fast and frequent changes in scheduled production targets, often over time intervals shorter than the (closed-loop) time constant of the process. As a result, the process may permanently operate in a transient mode (as opposed to being predominantly at or close to a steady state). Under these circumstances, it is crucial that the scheduling calculations account explicitly for the dynamics of the process and the performance of its control system, such that the aforementioned scheduled transitions are feasible and economically optimal (Baldea and Harjunkoski, 2014). Any approach for integrating scheduling and control must combine the long scheduling time horizon with the frequent execution

36

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

Sets Sets i j j g k m

scheduling time slot dynamics time slot dynamics time slot product grade piecewise linear breakpoints Lagrangian relaxation iteration

Continuous variables t time cycle time tC ti length of scheduling time slot i tgP production time of each product P tgi production time of each product in slot i amount of product g produced pg steady state production time tiss titrans transition time in each slot ending time for slot i tie tis starting time for slot i error eij ui scale-bridging model input Hammerstein block output hi xij state variable y ij state-space output wij scale-bridging model output  ij vector of scheduling-relevant variables w  inv w vector of storage-relevant variables ij ijk sgij st tgi in/out

fst

gij

RIEij  i,m

special ordered sets of type 2 variable storage holdup time spent in storage by product g in slot i flow into/out of storage of product g reverse-integrated error complicating constraint violation

Integer variables zgi binary variable defining production schedule bij binary variable defining RIE special ordered sets of type 1 variable ijk Parameters and coefficients number of scheduling time slots NI TI maximum allotted time per scheduling slot NJ number of dynamics time slots length of each dynamics time slot TJ  dom dominant system time constant tset settling time timescale conversion factor ˛ Dg demand of product g for each i and j production rate of product g Pg Pricei associated operating costs for product g over time large number used in big-M linearization M pwijk value of the piecewise linear function at breakpoint k bpijk breakpoint k of piecewise linear function RIE tolerance   Lagrangian relaxation (LR) termination tolerance LDi,m Lagrange multiplier positive parameter for Lagrangian multiplier calcui,m lation Sj finite step response (FSR) model parameter

of the control system, thus accounting for both (longer term) economic performance and (short term, real-time) safety and stability requirements. The multiple time-scale nature of this problem results in stiff models that are computationally intensive (Baldea and Harjunkoski, 2014). Furthermore, first-principles dynamic process models are almost invariably nonlinear and high-dimensional (Wang and Shan, 2007). Together, these features represent a considerable challenge to the effective integration of scheduling and control, including solving the resulting problems in a practical amount of time. Motivated by the above, in this work, we propose a novel computationally-efficient scheduling formulation which integrates scheduling and process dynamics/control information. The key contributions of this paper are: • a general framework for developing (exact) linearizations of low-order, data-driven representations of the closed-loop process dynamics (introduced in our previous work by Pattison et al., 2016), based on commonly-employed model classes, such as Hammerstein–Wiener (HW) and finite step response/finite impulse response (FSR/FIR) models. • a production scheduling formulation for continuous processes, incorporating these models. Importantly, this scheduling problem is formulated as a mixed-integer linear program (MILP), and can be solved using powerful solvers currently available. • a Lagrangian relaxation (LR) scheme for increased computational efficiency in solving the aforementioned scheduling problem. A polymerization case study is included to demonstrate the proposed approach. The article is organized as follows: we begin with a review of the relevant literature concerning the integration of production scheduling and process control. The next section contains an overview of scheduling techniques and a presentation of the model classes (HW and FSR) considered. This review is followed by methods for linearization of the selected models. From here, the model dimensions are reduced using their unique structural properties, and the MILP scheduling problem is given with a computationallyefficient LR strategy. Following the theoretical content, a case study is presented on optimal scheduling in a polymerization reactor. 2. Literature review When presented with the need to account for process dynamics, conventional scheduling problem formulations typically circumvent computational efficiency issues by capturing dynamic information in terms of estimated transition times between products in a set product wheel, and/or constraints concerning the maximum rate of change for, e.g., production rate transitions (Zhang and Grossmann, 2016; Maravelias, 2012). However, these techniques inherently assume the process reaches steady-state prior to set-point changes, and that it operates (mostly) at steady state. These premises are no longer valid when set-point changes occur at frequencies comparable to or higher than the dominant time constant in the process, and may therefore result in dynamically infeasible transition sequences (Pattison et al., 2016; Chu and You, 2012). Further, transition time estimates and transition rate constraints are often chosen to be very conservative (to account for safety and equipment limitations), typically resulting in a suboptimal solution. An evolution of these basic concepts consists of separating a production time slot into two sub-intervals, transition from the previous slot and steady-state operation, estimating transition dynamics or calculating them offline by enumerating all possible pairwise switches (transitions) (Nyström et al., 2005; Tong et al., 2015). While this is relatively manageable for systems with

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

37

Fig. 1. Scale-bridging model (SBM) for integration of scheduling and control. Adapted from Du et al. (2015).

a small product wheel, obtaining such calculations for nonlinear systems that transition between a continuum of product grades or production rates is computationally intractable, and may still be quite inaccurate. Motivated by the above, numerous previous studies concerning the integration of process dynamics and control information in scheduling calculations approached the problem by embedding the full-order differential algebraic process model in the scheduling problem formulation. Zhuge and Ierapetritou (2016) introduced a novel decomposition scheme wherein the problem was decomposed into two subproblems, an MINLP for the optimal production sequence and an NLP for the production times based on product demand. While workable in principle, practical implementations of full-order models for scheduling calculations are limited to relatively small-scale systems whose models are low dimensional. This includes, e.g., polymerization reactors (Du et al., 2015; Baldea et al., 2015; Prata et al., 2008; Chu and You, 2012; Zhuge and Ierapetritou, 2016), batch reactor systems (Chu and You, 2014; Nie et al., 2012, 2014), and batch separations processes (Nie et al., 2012). Conversely, it has been shown that implementing full-order dynamic models into scheduling calculations is not practical (from the point of view of computation time) for problems of a meaningful size and complexity (Du et al., 2015; Pattison et al., 2016; Nyström et al., 2005; Flores-Tlacuahuac and Grossmann, 2010). It is also important to note that embedding only the dynamic model of the process in the scheduling problem does not explicitly account for the presence or performance of a control system. Rather, the onus of computing the control actions is placed on the scheduling problem, with the resulting solution being in effect an open-loop one from a control perspective. These shortcomings have provided the impetus for deriving low-order models of the closed-loop dynamics of the process, to be used in scheduling calculations. Such models can be derived

either from detailed process models, via model reduction, or from experimental or simulated operating data, using system identification techniques. Beal et al. (2017) proposed to integrate scheduling and control by representing nonlinear dynamics with linear relations. These linear relations were time-scaled using throughput rates and a residence time approximation. While computationally inexpensive, the linear model does not capture all aspects of the nonlinear process which may be optimized in a scheduling calculation, and amounts to representing nonlinear dynamics via approximated residence times (note that this is closely related to previous works which rely on estimated transition times to capture process dynamics). In our previous work, we defined a class of such models, which we refer to as scale bridging models (SBMs) (Fig. 1). Initially discussed by Park et al. (2014), SBMs have been successfully used to represent specific supervisory controllers (Du et al., 2015), implemented as soft constraints in model predictive control formulations (Baldea et al., 2015), and used to represent complex dynamic behavior of air separation units (Pattison et al., 2016, 2017) for optimal scheduling. While SBMs have proven to be computationally efficient and accurate representations of full-order dynamics, they are typically nonlinear, and therefore, their usage has led to mixed-integer nonlinear program (MINLP) formulations of the resulting integrated scheduling/control problems. Conversely, optimal scheduling problems are typically represented and solved using mixed-integer linear program (MILP) frameworks. This is motivated both theoretically (due to optimality guarantees that can be obtained for the solution) and practically, given that most available enterprise optimization software used in industry is geared towards the solution of integer/linear optimization problems. Based on the above, the objective of this paper is to provide a systematic means for deriving

Fig. 2. Hierarchy of decisions for process operations.

38

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

Fig. 3. Diagram showing time slot structure, with sampling points for the discretized process dynamics, j spaced uniformly throughout a scheduling slot, i, of maximum length TI . The grey region highlights the excess time length preallocated in the problem definition.

accurate linear approximations (and, in some cases, exact linearizations) of SBMs, using multiple potential model formats/classes.

Such that Profit is maximized for the given makespan/scheduling horizon.

1. 3. Problem definition We consider a continuous process system that is capable of producing multiple products and/or modulating its production rate to reflect market changes (e.g., changes in supply and demand quantity and/or prices). Facilities are assumed to exist for product storage, and can be used either to store product(s) until the delivery date, or as a means to modulate production and satisfy demand. In defining the operating strategy for this system, we begin from the well-established decision-making hierarchy shown in Fig. 2. A planning step is assumed to determine production targets based on contractual agreements and anticipated product demand, and inform the decisions of the scheduling layer. In turn, at the scheduling level, the optimal production sequence is determined, and passed on to the control layer in the form of targets/setpoints. The control system determines the values of the process manipulated variables in closed-loop, ensuring, (i) setpoint tracking and, (ii) closed-loop stability and disturbance rejection. Real-time measurements of process variables and fault and disturbance information can be fed back to the scheduling layer to inform rescheduling decisions (Touretzky et al., 2017). Based on this hierarchy, the following problem is to be addressed: Given 1. 1. Production information • product specifications • (time-sensitive) product demand (rates) • (time-sensitive) product prices • (time-sensitive) feedstock availability • (time-sensitive) feedstock prices 2. Process information • process dynamic model • control system configuration • equipment limitations and safety constraints 3. Storage system information • storage system dynamic model • capacity and flow constraints Determine 1. 1. Production timing and sequencing (for multiple products or single product) • assignment of product(s) to time slots • time slot durations and total makespan (if variable) • production rate (if variable) at each time instant

We discuss these elements in what follows, with a focus on the type and format of process and storage system information to be embedded in the scheduling problem. We distinguish between processes creating multiple products and those which make the same product(s) at time-varying production rates. 3.1. Production timing and sequencing 3.1.1. Time basis We utilize a slot-based formulation, with NI scheduling time slots, each of variable length, ti , such that the total cycle time, tC , is given by: tC = NI



ti

(1)

i

where ti is the difference between the end time of slot i, tie , and the start time of slot i, tis . We note that the ending and starting times, tis and tie , are based on the production time necessary for each product (as explained later) ti = tie − tis

(2a)

s ti+1

(2b)

=

tie

tneI

= tC

(2c)

t1s

=0

(2d)

Each time slot is allocated a maximum duration, TI ≥ ti . Within each scheduling slot, i, there exists a set of NJ time intervals of uniform length TJ , which are used to represent the discretized dynamics of the process (as will be discussed later). The maximum allowable scheduling slot length, TI , is an integer multiple of the dynamics time step length, TJ , i.e., NJ ∈ Z: TI =

1 NJ TJ ˛

(3)

Above, ˛ is a conversion factor used to convert between the time units of the scheduling level and the dynamics level. For the case where the scheduling time slots are of constant duration, we utilize the notation ti = TI to denote that the slot length is equal to its maximum value. We note here that our definition of the parameter TI is equivalent to “pre-allocating” a sufficient of time for production of each product. This approach leads to an increase in the number of variables in the problem formulation (since TI ≥ ti , the values of the state variables at times greater than tie are not meaningful). Nevertheless, this device is necessary to ensure that the problem is

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

39

indefinitely. The equations below give the product assignment constraints:



zgi = 1 ∀i

(6a)

zgi = 1 ∀g

(6b)

g

 i

The binary variable zgi is used to determine the sequence of production; zgi is 1 if product g is produced in scheduling slot i and 0 otherwise. Eq. (6a) ensures that only one product, g, is made in each time slot, i, Eq. (6b) ensures that a product is only made once during a cycle. To maintain the cyclical nature of the production sequence, a constraint is imposed such that the states at the beginning of the cycle are equal to their set-points at the end of the cycle:  x1,1 = h NI Fig. 4. Process evolution during production of product grade g in slot i. The transition periods are shaded in grey. The finer time grid used for discretizing the process dynamics is omitted for clarity.

formulated in the most general way and can capture the process dynamics in their entirety. Fig. 3 describes this construct. By contrast, previous works (see, e.g., Flores-Tlacuahuac and Grossmann, 2006) assume that each transition is complete in each scheduling time slot, and that the desired steady-state is reached every time. In turn, this means that the scheduling time slot can be explicitly separated into a transition and a production period, with the dynamics being accounted for only in the former. Remark 1. The length of TJ , or sampling time, is determined by the system dynamics (or inertia), using, e.g., heuristic relations of the type (see Braun et al., 1999; Seborg et al., 2010): 0.01 ≤

TJ ≤ 0.05 dom

tset tset ≤ TJ ≤ 15 6

(4) (5)

where  dom is the dominant time constant of the system, and tset is the system settling time (Seborg et al., 2010). 3.1.2. Scheduling constraints with multiple products We first consider the case where a product wheel with a finite number of products is defined. This can reflect a situation where key characteristics such as composition differ between products (this is the case, e.g., in the production of different grades of polymer (Du et al., 2015; Daoutidis et al., 1990). Under these circumstances, the scheduling calculation aims to determine the optimal sequence of production and time slot lengths (and the total makespan), such that product demand is met at minimum cost (or with maximum profit). Since in this case the product characteristics are specified, it is reasonable to assume that each production slot comprises two separate time intervals (Flores-Tlacuahuac and Grossmann, 2006). First, a transition interval, where the process dynamics are manifest in the transition characteristics between products; waste product is typically generated in this time interval. This is followed by a steady state interval, where the desired product characteristics have been reached and are maintained for as long as it is necessary to produce the required quantity of product. Fig. 4 shows this interval separation. Once the demand is met, a transition occurs to manufacture the next product in the production sequence. For the purpose of this discussion, we assume that each product is only manufactured once during a cycle and that the cycle repeats

(7)

The constraint is formulated as above (rather than setting initial and final states equal is adequate since in a cyclical process the states reach steady-state (set-point) values prior to the end of the scheduling time slot, i. The production time required for each product is given by variable tgP , which is calculated based on the production rate and demand. In order to determine tgP , constraints must be enforced such that the amount produced is greater than or equal to the demand but less than some maximum amount, pmax . g In this case, we assume that total demand, Dg , is to be met at the end of the cycle, where products are delivered all at once: Dg ≤ pg ≤ pmax g

(8a)

pg ≥ Pg tgP

(8b)

Here, the variable pg is the amount of product g produced, and Pg is the rate of production of product g. Eq. (9b) is only linear if Pg is a constant. If Pg is variable, special care should be taken to linearize the bilinear term Pg tgP . If instead demand is specified in terms of a demand rate, D˙ g , (8) becomes: D˙ g Tc ≤ pg ≤ pmax g pg ≥ Pg tgP

(9a) (9b) tiss ,

Within each slot, i, there exists a steady-state period, in which the system is at steady state and the product(s) meet(s) specP , of product g in slot i, required to ifications. The production time, tgi meet demand, lies within this steady-state period. We recall here the mathematical device described above (i.e., the “pre-allocation” of the number of sample times j to each slot i), and observe that P , with the samples situated to the tiss may in fact be larger than tgi e the right of ti not being accounted for in the optimization problem objective (see Fig. 4). We thus have: P tgi = tgP zgi

(10)

To linearize the bilinear term in (10), big-M linearization is used, where TI = M: P tgi ≤ zgi TI

(11a)

P tgi ≥ tgP − TI (1 − zgi )

(11b)

P tgi ≤ tgP + TI (1 − zgi )

(11c)

The total steady-state production time, tiss , is a subset of the P ≤ t ss maximum alotted time per scheduling slot, TI , such that tgi i ti ≤ TI . Lastly, titrans is the transition time in each slot. The equations below complete the time formulation for the problem: titrans = TI − tiss

(12a)

40

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

tie = tis + titrans +



P tgi

(12b)

g

One of the key issues in this framework is the determination of the durations of the transition and production intervals within a time slot (and the subsequent computation of the variable bij ), based on the dynamic response of the process, as represented by its model. The typical approach to this end (as proposed by, e.g., Flores-Tlacuahuac and Grossmann, 2006) is to use the duration of the transition interval as a decision variable, with an additional constraint that the system states reach values corresponding to the product assigned to the respective slot at the end of the transition. In practice, this is implemented in conjunction with a discretization scheme (e.g., collocation on finite elements) that allows for the time step to vary in proportion to the transition interval duration. This latter condition is not met by a discrete-time data-driven representation of the process dynamics, where the discretization time step/sample time is fixed. As a consequence, we propose the use of the Reverse Integrated Error (RIE) as introduced by Touretzky et al. (2016). The key idea behind RIE is to integrate (backwards in time) the difference between the variable trajectory and its setpoint. The main motivation for using RIE compared to other steady-state determination methods lies in the avoidance of erroneous transition times for the case of overshoot and/or oscillatory behavior (Touretzky et al., 2016). The following set of equations shows how RIE was determined. A set, j is used to index RIE over each dynamics discretization time-step: pos

neg

{eij , eij } ≥ 0 pos

eij

neg

− eij

pos

eij = eij

RIE ij =

(13a)

= ui − wij

(13b)

neg

+ eij



(13c)

eij

(13d)

j≥j

Relations (13a)–(13c) define the absolute value of the error between the setpoint ui and the scheduling-relevant variable trapos jectory, wij , in an MILP framework. Above, eij takes on the positive neg

values of ui − wij and eij

pos

the negative values. In order to ensure neg

that when summed, {eij , eij } give the absolute value of the error,

pos neg eij , {eij , eij }

are specified as part of a special ordered set of type 1 (SOS1), discussed in detail in Section 4.2.1. Eq. (13d) defines RIE ij as the backwards integral of eij . A binary variable bij is defined such that when the product grade is within specifications at steady-state, bij is 1, where  is a system-dependent tolerance defining the maximum allowable RIE ij (Fig. 5):

Fig. 5. Diagram showing Reverse Integrated Error (RIE) calculation, where on-spec production (within RIE tolerance) occurs in the shaded region and the grey-dashed line represents the RIE tolerance, .

3.1.3. Scheduling constraints assuming fixed product(s) and variable production rates In this case we assume that the properties of the product(s) made by the plant do not change. While the same products are made in a continuous stream at all times, the operating state of the plant can change as the production rates can be varied to accommodate (fast) changes in market conditions. To maximize economic benefits, the scheduling calculation can be set up such that the time slots allow for frequent changes (e.g., hourly) in the production targets (and the setpoints furnished to the control system). As such, the duration of the scheduling time slots may be shorter than the time constant of a typical chemical plant, with the consequence that plant operation may never reach steady state. Emphasis is therefore placed on ensuring that the scheduled plant operation is optimal and feasible at every time point dynamic constraints (i.e., abiding by safety limits and equipment limitations while meeting product quality specifications at all times). Remark 2. A key factor in formulating the scheduling problem in this case is the availability of forecasts for demand and feedstock/utility prices. Such information is included in the objective function and therefore should be given consideration when determining the length of the scheduling horizon, tC . Intuitively, tC (typically fixed in this situation) should be less than or equal to the shortest time span of these time series data. Furthermore, the length of the scheduling time slots should be correlated to the dynamic properties of the process and with the sampling time of the aforementioned economic data.

RIE ij ≤  + M(1 − bij )

(14a)

3.2. Inventory information

RIE ij ≥  − Mbij

(14b)

The time spent in storage and the amount of each product in storage can be used to determine the cost of storing material for use in the objective function. In addition, models of the storage system dynamics are necessary to ensure that storage capacity constraints are not violated, i.e., the inventory level, s, remains below a maximum value and is physically realizable:

The following equation determines the steady-state time, tiss , based on bij , such that tiss is the sum of the dynamics time steps where product g is within specifications: tiss = TJ



bij

(15)

j

The timing equations in this section can be solved simultaneously, where the production time is dependent on when the system reaches steady state and products are within specifications, as defined by RIE ij .

0 ≤ sgij < sgmax

(16)

3.2.1. Multiple products In the case where multiple products are manufactured over makespan tC , storage facilities are assumed to be used to store products until the end of the cycle, when the entire quantity of product is

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

delivered. We calculate the time spent in storage for each product, st , starting from the end of the production interval in the respectgi tive time slot and up to the end of the cycle, tC , following the model description given by Flores-Tlacuahuac and Grossmann (2006) and Pinto and Grossmann (1994): st tgi = zgi (tC − tie )

(17)

In order to linearize the bilinear term in (17) we define M = NI TI and write: st tgi ≥ tC − tie − M(1 − zgi )

(18a)

st tgi ≤ tC − tie + M(1 − zgi )

(18b)

st tgi

(18c)

≤ Mzgi

The stored amount of each product is assumed to be the amount of material produced, pg . The time spent in storage for each product is then added to the objective function and multiplied by the cost of storage (Cgst ) as an associated operating cost, cgst : cgst = Cgst



st tgi pgi

(19)

i

We note that the cost of storage is a recurring operating cost associated with storing product at the proper conditions until it is delivered to customers. 3.2.2. Fixed product(s) and variable production rates Analogous to the multi-product case, we assume that a storage facility/tank exists for each product. Part of the product stream can be diverted to the corresponding storage facility when the production rate exceeds demand and, conversely, the production rate can be supplemented with stored material when the economics of the system dictate it. As a consequence, the storage system may itself be in a transient state at all times, and an appropriate model must be included to capture the evolution of inventory levels: ds = fsin (t) − fsout (t) dt

(20a)

fsin (t) = fsout (t) + P(t) − D(t)

(20b)

fsout (t)

(20c)

= D(t) − P(t)

In discrete time using a forward Euler discretization scheme: sg,i,j = TJ (fsin − fsout ) + sg,i,j+1 g,i,j g,i,j

(21a)

fsin = fsout + Pgij − Dgij gij gij

(21b)

fsout = Dgij − Pgij gij

(21c)

where sgij is the inventory holdup, and fsout and fsin represent the gij gij flow leaving and entering the storage system, respectively. In the above expressions, demand is given as a rate. However, demand may also be a fixed value depending on the type of problem addressed. If demand is a value to be met at the end of the time horizon, (20) and (21) can instead be represented using: ds = P(t) dt

(22a)



P(t) ≥ D

(22b)

t

Or, in discrete time: sg,i,j = TJ Pg,i,j + sg,i,j+1

 i

j

Pg,i,j ≥ Dg

(23a) (23b)

41

3.3. Process information The arguments presented above make it critical that a representation of the process dynamics be included in the formulation of the scheduling problem. Either full order, typically first-principles, or low order models can be used, as discussed below. 3.3.1. Full-order models The dynamics of process systems can be captured via first-principles models comprising material and energy balance equations, as well as relevant constitutive relations. These models are typically described by systems of differential algebraic equations in the input-affine form (Baldea and Daoutidis, 2012): x˙ = f(x) + G(x)u

(24a)

0 = h(x)

(24b)

Initial efforts towards integrating production scheduling and control (Flores-Tlacuahuac and Grossmann, 2006; Zhuge and Ierapetritou, 2012) were based on embedding such models (discretized using, e.g., finite difference, finite element or Runge-Kutta schemes) in production scheduling calculations. While advantageous in principle (since it uses the most accurate process model without any further manipulation), this approach suffers from two significant shortcomings. First, the state vector x ∈ Dx ⊂ Rn can be very large (i.e., n = O(105 )), and f(x) : Dx → Rn , G(x) : Dx → Rn×m and h(x) : Dx → Rl are highly nonlinear functions. The resulting scheduling problems are in turn high-dimensional and nonlinear and, as a consequence, are difficult to solve. Second, the vector u of manipulated variables becomes a decision variable, and the solution of the scheduling problem is in effect an open-loop control solution, with the values u(t) being computed for the entire makespan/scheduling horizon. This solution is thus optimal in the nominal case, when the model (24) is a perfect representation of the plant dynamics (i.e., there is no “plant-model mismatch”) and no disturbances are present. Optimality is, however, lost when these conditions are not met, and system stability may be lost as well under certain circumstances. In order to address the second shortcoming mentioned above, it is possible to incorporate along with model (24), the representation of the control system of the process u = K(x, xsp )

(25)

where, for simplicity, we assume a state feedback control law. While the dimension m of the input u ∈ Du ⊂ Rm is typically much lower than n, and incorporating (25) in the scheduling model does not increase its size considerably, the dimensionality challenge mentioned previously remains. Additionally, it is possible that the control law in (25) cannot be obtained explicitly (this is the case with, e.g., model predictive control) and, as a consequence, the scheduling problem would become a bi-level optimization problem, with the associated solution challenges. 3.3.2. Low-order models The above arguments provide a strong motivation for the derivation of low-order models, with the state vector x¯ of the low-order ¯ dim(x). An additional argument in favor model having dim(x) of adopting a low-order modeling approach for capturing process dynamics in scheduling calculations can be made on an economic basis. Namely, not every process variable (and, equivalently, model state) is relevant to the schedule calculation, and the number of scheduling-relevant variables may in effect be small. As discussed in Pattison et al. (2016), scheduling-relevant variables either (i) affect operating cost, thus having an impact on the optimal schedule, and/or (ii) are at or near their bounds either during transitions

42

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

or during steady-state operation at any point in the process operating domain. The approach proposed by Pattison et al. (2016) for identifying scheduling-relevant variables is summarized below: 1. Identify the scheduling-relevant process and quality constraints (PCs and QCs): Scheduling-relevant constraints relate to process performance, safety, and efficiency and are at or near their bounds (active) during transition or steady-state operation. The PCs and QCs which are not near bounds during operation (inactive) do not hinder the process agility and therefore are not scheduling-relevant. 2. Determine the subset of process state variables related to scheduling-relevant PCs and QCs: These variables make up the list of scheduling-relevant variables which must be captured using low-order models. 3. Determine the subset of input/output process variables with scheduling objective-relevant dynamics: This list includes, e.g., raw feed flowrate and heating/cooling rates (on the input side), and variables which track the properties of the process output stream which factor into the objective function (on the output side). The scheduling-relevant variables can then be represented using low-order models. The general structure of such low-order models, which we refer to as scale-bridging models (SBMs) to acknowledge their role in bridging the time scales of scheduling and control decisions, is thus such that their inputs ui consist of targets set at the scheduling level, and their outputs wij are the dynamic responses of the scheduling-relevant variables described above. Such models are thus typically in a multiple-input, single-output (MISO) form since separate models can be identified for each relevant variable. The derivation of scale-bridging models (SBMs), can be carried out following two main paradigms: • model-based: in this case, the starting point is a full-order, detailed plant model, and techniques such as perturbation analysis (Baldea and Daoutidis, 2012) are used to obtain a lower dimensional representation of the process dynamics. The process typically entails an asymptotic analysis, followed by a nonlinear coordinate transformation that results in a (minimal order) state-space realization of the remaining (slower) dynamic modes. Another option is the use of nullspace projections, where an equivalent dynamical system with improved computational performance is derived (Nie et al., 2013; Yu et al., 2015). The advantage of model-based methods is that they result in physically meaningful low-dimensional models, whose states can be interpreted from the perspective of the physical and chemical phenomena captured by the original, first-principles representation. On the other hand, the theory for determining the perturbation parameters and the coordinate transformations required for the aforementioned approaches is in general not constructive, in the sense that, while conditions can be defined for these to exist, it is often not trivial to find appropriate expressions for a particular practical system. Empirical dimensionality reduction methods such as proper orthogonal decomposition (POD) (Lang et al., 2009) and the use of balanced empirical Gramians to capture nonlinear system behavior near an operating point (Hahn and Edgar, 2002) have also been proposed. These techniques produce models without physically meaningful states, but require less physical intuition about the system/model at hand, being hence more broadly applicable. • data-driven: in this case, the starting point is a set of process operating data, that are either collected from the operation of an existing plant, or generated via simulating a detailed dynamic process model. The data reflect the impact of deliberate excitations imposed on the process as well as the consequences of measurable disturbances, and should ideally cover the entire

operating range of the plant. System identification techniques (see the text by Ljung (1987) for a thorough treatment of the topic) can then be used to derive low-order models of the scheduling-relevant process dynamics. In this paper, we focus on two specific classes of data-driven SBMs, namely, models of the Hammerstein–Wiener (HW) form, and finite step response/finite impulse response (FSR/FIR) form, which will be described in Sections 4.1.1 and 4.1.2.

3.4. Scheduling problem formulation Based on the above, we define two general optimal scheduling problems that account for process dynamics. Eq. (26) gives the general formulation of the optimal scheduling problem using full-order models, where Pricegij represents any costs associated with operation such as utilities, feedstocks, and  inv is the vector of inventory variables such as flow storage costs, w gij  gij is the vector of into and out of the storage system, etc., and w scheduling-relevant variables: max J = ui

 i

j

inv  gij , w  gij (Pricegij , w , Dg )

(26)

g

s.t. Full-order Continuous time Process model (24) Inventory model (16), (18), (20), and (22) Initial Conditions Process Constraints Quality Constraints The optimal scheduling problem using low-order SBMs can be written analogously as: max J = ui

 i

j

inv  gij , w  gij (Pricegij , w , Dg )

(27)

g

s.t. SBMs Inventory model (16), (18), (21), and (23) Initial Conditions Continuity Conditions (30) Process Constraints Quality Constraints The full-order scheduling problem (26) is high-dimensional and nonlinear, expressed as a mixed integer dynamic optimization (MIDO), which is typically discretized as an MINLP. The aim of this work is to express (26) as an MILP with low dimensionality, leading to formulation (27). The rest of this work will provide SBM expressions and linear reformulations of these models, and a case study to show their value as low-order modeling techniques.

4. Data-driven scale-bridging models For simplicity, we focus on single input, single output (SISO) systems, noting that the arguments below can be readily extended to MISO systems.

4.1. Model classes considered 4.1.1. Hammerstein–Wiener (HW) models The HW structure comprises a linear dynamic model flanked by static nonlinear transformations at the input (Hammerstein) and output (Wiener) (Fig. 6). Below we show a single input, single out-

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

43

Fig. 6. General HW model form.

put HW structure whereby the linear dynamics are represented by a continuous-time state-space model: h = H(u) dx = Ax + Bh dt

(28)

y = C x w = W (y) where H and W are continuous input and output nonlinear functions, A, B, and C are the matrices of the linear state-space system, u are the inputs (setpoints), x are the states, and w is the system output. For the case of multiple variables represented by differ There ent SBMs, the variables in question are represented using w. exists a vast body of work on the identification of HW models (Wang and Ding, 2012; Narendra and Gallman, 1966; Bai, 2002; Wills et al., 2013), with practical applications supported, e.g., by the popular MATLAB System Identification Toolbox (MATLAB, 2016). An application of HW structures as SBMs representing (continuoustime) process dynamics in scheduling calculations was reported by Pattison et al. (2016). The continuous state space model can be converted to discrete time using, methods such as zero- or first-order hold, impulse invariant, and zero-pole matching (Ogunnaike and Harmon Ray, 1994). Owing to space limitations, below and in the remainder of the paper, we will rely on the zero-order hold method. However, the developments below are general and are not limited to the discretization approach. The discretization of H(u) and W(y) will be discussed later: hi = H(ui ) xij+1 = Axij + Bhi

(29)

yij = C xij wij = W (yij )

For the purpose of integrating (29) in a scheduling model (which we will discuss in more detail below), a continuity condition must be enforced for each state variable, x, such that the value of each state at the beginning of the scheduling time step/slot is equal to its value at the end of the previous scheduling time step/slot: xi,1 = xi−1,NJ

(30)

Note that the sample time corresponds to the value of TJ in the scheduling problem formulation. 4.1.2. Finite step response (FSR) models FSR models are a class of non-parametric models and are (along with their counterpart, finite impulse response (FIR) models) used extensively in practice, particularly when the order and time delay of the dynamics of the system under consideration are unknown (Wise and Ricker, 1993). SISO FSR models are of the general form (Seborg et al., 2010):

 J−1

wi,j+1 = wi,j=0 +

[Sn ui,j−n+1 ] + SJ ui,j−J+1

(31)

n=1

where Sj are the step response coefficients. A benefit of these models is that they are applicable to systems with unusual dynamic behavior. However, the identification of such models can be cumbersome as the number of coefficients increases (Wise and Ricker,

Fig. 7. Piecewise linear representation of the output function W(yij ), where the green circles represent the locations of the breakpoints (bpijk ) in the input space and their locations (pwijk ) in the output space. The index k indicates a breakpoint.

1993; Dayal and MacGregor, 1996). In particular, this formulation requires a sum over all previous time-steps, which can be computationally intensive for an MILP optimal scheduling framework due to coupling throughout the time-horizon. Note that FSR models are discrete by nature and their fidelity is highly dependent on sample time (TJ ) selection as well as proper assessment of system settling time. 4.2. Linear representations of Hammerstein–Wiener SBMs The choice of the functional form of the nonlinear functions H(u) and W(y) in Eq. (28) is dictated by the practical application under consideration. Typically, polynomials and piecewise linear functions give good results, and have also been the focus of significant theoretical development (Aryani et al., 2017; Bloemen et al., 2001; Patwardhan et al., 1998). In pursuing our objective of devising linear reformulations of data-driven SBMs, we note that piecewise linear approximations can be used to represent most nonlinear functions (including polynomials), and can be defined to represent these nonlinear transformations with the desired degree of closeness/accuracy. A discussion on how this is performed is beyond the scope of this paper, and we direct the reader to the text by Eriksson et al. (2004) for details. Motivated by the above, in this study we will focus on the case where H(u) and W(y) in (28) are piecewise linear functions, and emphasize the derivation of exact representations of the corresponding discrete-time HW scale-bridging models in scheduling calculations. In this section, we present three methods of exact linearization of piecewise linear functions and draw comparisons. Considering specifically this class of models, we represent the (piecewise linear) input and output functions H(u) and W(y) as follows: • The piecewise linear input function (Hammerstein block) is given as: hi = H(ui ) = +pwH ik

pwH − pwH ik+1 ik H

H

bpik+1 − bpik

H

(ui − bpik )

∀ bpHik < ui ≤ bpHik+1

(32) H

The index k represents “breakpoints,” bpik gives the values of the breakpoints/divisions of the input space, and pwH are the ik corresponding values in the output space (Fig. 7).

44

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

Fig. 8. Linearization strategy comparison, where active breakpoints are represented as red squares and the block input follows the dashed grey line.

• The output piecewise linear function (Wiener block) described in a similar manner, with the exception that the input variable yij (which is the output of the linear dynamic block) is indexed on both i and j to account for the fact that the discretization of the linear dynamics (index j) is finer than the discretization of the scheduling horizon, given by the index i:

wij = W (yij ) =



H =1 ik



W (yij − bpijk ) W W bpijk+1 − bpijk

(36a)

k

H ≤ bH ik ik

pwW − pwW ijk+1 ijk

W ∀ bpW ijk < yij ≤ bpijk+1

+pwW ijk

Constraints are placed on the SOS2 variables such that only two adjacent SOS2 variables are nonzero (active), imposing that they in effect become weighting variables for linear interpolation (Fig. 8):

bH =2 ik

(36b) (36c)

k

(33)

An example of a piecewise linear function is shown in Fig. 7, with the above variables identified on the diagram.

bH + bH ≤ 1 ∀k > k + 1 ik ik

(36d)

bH ∈ {0, 1} ik

(36e)

W =1 ijk

(37a)

 k

The above expressions lend themselves naturally to a mixedinteger linear reformulation, and we discuss and contrast potential approaches below.

W ≤ bW ijk ijk



bW =2 ijk

(37b) (37c)

k

4.2.1. Linearization using special ordered sets Special ordered sets (SOS) are variable sets defined such that only specific elements can be non-zero, thereby providing optimization algorithms with better means of branching to find the optimal values of integer variables. In special ordered sets of type 1 (SOS1), only one of the elements can be non-zero, while in special ordered sets of type 2 (SOS2), two consecutive elements can be non-zero. The introduction of SOS2 has been motivated by the need for efficiently representing nonlinear and piecewise linear functions (particularly objective functions) in linear optimization models (Beale and Tomlin, 1970) and are thus a natural choice for defining linear reformulations of (32) and (33). To this end, we define SOS2 variables H = {H } and W = {W } such that: ijk ik ui =



H

[H bpik ] ik

(34a)

[H pwH ] ik ik

(34b)

k

hi =

 k

yij =



W

[W bpijk ] ijk

(35a)

k

wij =

 k

[W pwW ] ijk ijk

(35b)

bW + bW ≤ 1 ∀k  > k + 1 ijk ijk

(37d)

bW ∈ {0, 1} ik

(37e)

We refer to this SOS2 reformulation as “Explicit SOS2.” Note that this reformulation of the Hammerstein and Wiener blocks in the SBMs requires a rather large number of binary variables, and may appear to be computationally taxing. Nevertheless, this formulation allows MILP solvers to intelligently explore the binary search tree, thereby considerably reducing computation time (Keha et al., 2004). Furthermore, several solvers (e.g., CPLEX (CPL, 2015), BDMLP (BDM, 2017), and XPRESS (XPR, 2017)) natively support the SOS2 variable type by enforcing the adjacency conditions (Eqs. (36b)–(36e) and (37b)–(37e)), in which case the user-input model equations reduce to (34), (35), (36a), and (37a). We refer to this reformulation as “implicit SOS2.” Switching now to discussing the use of SOS1 variables, we observe that they are able to increase computational efficiency (Beale and Tomlin, 1970) for the case when the state space model inputs hi are constant over the scheduling time slot. For the scheduling problems discussed earlier, this corresponds to the case of multiple, separate products. In this context, the SOS1 variables, H }, are defined to reformulate the piecewise constant expresK = {ik sion of the Hammerstein block (note that the Wiener block should still be described by a piecewise linear, rather than piecewise constant function):

 k

H ik =1

(38a)

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

ui =



H ik bpk

(38b)

H ik pwk

(38c)

k

hi =

 k

H ik = {0, 1}

(38d)

SOS1 structures consist of a set of continuous variables, where only one member, k, can be non-negative. When utilized to linearize piecewise constant functions, constraints (38a) and (38d) are imposed. Similar to SOS2, the SOS1 variable type is built-in to certain MILP solvers such as CPLEX (CPL, 2015), BDMLP (BDM, 2017), and XPRESS (XPR, 2017); a feature that can be exploited when implementing these models in specific software. 4.2.2. Linear reformulation using big-M The underlying if–then conditions of expressions (32) and (33) H and z W , which transcan be reformulated using binary variables zik ijk form the input, ui , in the case of the Hammerstein block, and, the output of the state-space model, yij , for the Wiener block into the state space input hi and model output wij , respectively. The corresponding expressions are given below: hi = H(ui )k=0 +



+



H [(H(ui )k − H(ui )k−1 )zik ] = H(ui )k=0

k

AH z H = H(ui )k=0 + ik ik

k

H Bik

(39)



 k

AW z W = W (yij )k=0 + ijk ijk



W Bijk

(40)

k

Big-M techniques can then be used to reformulate the bilinear terms in Eqs. (39) and (40):

H Bik

≥

AH ik

≤

AH ik

Method

Discrete variables

Equations

Big-M Explicit SOS2 Implicit SOS2 SOS1

3(NI NK )(NJ + 1) 3(NI NK )(NJ + 1) 1(NI NK )(NJ + 1) 1(NI NK )(NJ + 1)

7(NI NK )(NJ + 1) 5(NI NK )(NJ + 1) 2(NI NK )(NJ + 1) 2(NI NK )(NJ + 1)

4.2.3. Comparison of reformulation approaches This section will compare reformulation approaches in terms of their applications and impact on problem size. Fig. 8 shows how an input is transformed for each linearization method, where SOS1 is only applicable to the Hammerstein input nonlinearity H(u) and SOS2 and big-M can be applied to both input/output nonlinearities H(u) and W(y). Table 1 summarizes the model size changes that result from each of the reformulation approaches presented above. Note that the number of continuous variables in the HW model is independent of the reformulation method. The big-M and the explicit SOS2 formulations add the same number of variables, but the number of equations for the explicit SOS2 method is lower given the lack of a bilinear term. The implicit SOS2 reformulation requires the lowest number of user-defined variables and equations; nevertheless, the actual model dimensions are the same as in the explicit case, with the benefit that using SOS2 variables allows for a more efficient exploration of the binary search tree. 5. Reducing the dimensionality of linear reformulations of SBMs

W [(W (yij )k − W (yij )k−1 )zijk ] = W (yij )k=0

k

H Bik

Table 1 Model size (number of equations and discrete variables (binary and SOS)) for each linearization method for a HW model.

k

wij = W (yij )k=0 + +



45

−M

H

H (1 − zik )

(41a)

+M

H

H (1 − zik )

(41b)

H H Bik ≥ −M H zik

(41c)

H H Bik ≤ M H zik

(41d)

W W Bijk ≥ AW − M W (1 − zijk ) ijk

(42a)

W W Bijk ≤ AW + M W (1 − zijk ) ijk

(42b)

W W Bijk ≥ −M W zijk

(42c)

W W Bijk ≤ M W zijk

(42d)

where the values of MH and MW should be selected such that they are greater than both the difference between the highest and lowest values of the breakpoints, bp, and piecewise values, pw. In this H and z W take a value of 0 for all formulation, the binary variables zik ijk breakpoints greater than the input value and 1 for all breakpoints less than the input value (Fig. 8). This approach is straightforward to implement but becomes computationally intensive given that a separate binary variable must be defined at each breakpoint for each time step, with further H and BW . equations needed to linearize the bilinear terms Bik ijk

The linear model reformulations presented in the previous section can result in scheduling MILPs that are of considerable size. Nevertheless, the structure of the SBMs and their role in the formulation of the scheduling problem enables several strategies for reducing problem size; these are discussed below. 5.1. Hammerstein–Wiener models We begin by observing that the input and output variables of the SBMs are likely to be bounded between upper and lower limits defined by the scheduling problem as either process constraints or quality constraints. From the point of view of the HW representation of SBMs, these upper and lower bounds on the model output, i.e., the scheduling-relevant variable wij , must be contained within the codomain of function W(yij ) and, respectively, the domain of state-space output variable yij . Note, however, that these bounds do not necessarily span the entire codomain (Fig. 9). This naturally limits the number of breakpoints to consider in the representation of the piecewise function W(yij ), and reduces the problem size. Furthermore, if the output wij of W(yij ) is not included in the objective function (27), the exact value of wij is not needed (i.e., one is only interested to know whether this variable is within its bounds). In this case, the output piecewise linear function W(yij ) can be approximated by a linear function using the two breakpoints corresponding to the lower and upper bounds (Fig. 9). This observation can be further extended: up the upper (yij ) and lower (yijlo ) bounds on the state space output yij defined by the codomain of W(yij ) can be used in place of the output function. With this representation, variables which are bounded but are not in the objective function can be represented using Hammerstein models and bounds on the output of the linear dynamic up block during optimal scheduling calculations. That is, yijlo ≤ yij ≤ yij up

rather than wijlo ≤ wij ≤ wij , which substantially reduces the num-

46

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

in practical circumstances, the number of response coefficients required to accurately capture the system dynamics can be quite large (Eq. (31)). Thus, using FSR representations can increase model size and the computational cost of the scheduling problem. Here, we note that the FSR structure affords a useful simplification in the case where the system dynamics are quite fast and steady state is reached by the end of each scheduling time slot, i. In this situation (and recalling that the inputs to the SBMs are constant during each scheduling slot), the values of ui,j are zero for all 1 < j < NJ , and the expression (31) can be simplified to: wij = wi−1,NJ + Sj (ui − ui−1 )

(43)

We note that this simplification is suitable only when the system is “well-behaved” and the dynamic response does not exhibit large transients that must be captured and constrained. Fig. 9. Simplification of the output function W(yij ) using upper and lower bounds (process and/or quality constraints) on wij . Red squares represent the new breakpoints spanning the feasible region (non-shaded) and excluding the infeasible region (shaded). The red line (connecting the red squares) represents the linear approximation of the output function W(yij ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6. Solution strategy for integrated scheduling and control using linear SBM reformulations

ber of variables necessary given that the Wiener portion of the HW model is not necessary for representation of variables not in the objective function. Implementing these observations can have a significant impact on problem size, especially given that the values of W(yij ) must be calculated at every sample time j in each scheduling slot i, adding a substantial number of variables to the problem (Table 1) and increasing computation time.

• the scheduling variables change over every scheduling slot i and represent inputs to the SBMs that capture the schedulingrelevant dynamics of the system • the outputs of the SBMs are represented with a finer time discretization, with the output values predicted at every time point j in each slot i • the system states must meet a continuity condition (30) across scheduling slots

5.2. Finite step response models

The above features make this problem block-angular in structure, where the overarching “common problem” contains the objective function, spanning a series of subproblems containing relevant models and constraints (Elfeky et al., 2008) linked only by state continuity conditions (30) and the cyclic endpoint constraint (7) (Fig. 10).

Finite step response (FSR) models are useful for capturing complex linear process dynamics, and have been successfully used in practice for capturing the behavior of systems whose nonlinearity is not very strong in the operating region of interest. However,

The optimal scheduling problem defined in Eq. (27) has the following features:

Fig. 10. Block-angular problem structure, where subproblems generated for each scheduling slot i are divided by dynamic time slots j and are linked by state continuity conditions, making up the overarching common problem.

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

The above problem is given by (27) and will be referred to as PI. The structure of PI lends itself well a Lagrangian relaxation (LR)based solution strategy, which is discussed below. 6.1. Lagrangian relaxation of continuity constraints The LR framework consists of dualizing the complicating constraints: these constraints are multiplied with a Lagrange multiplier and included in the objective function. Then, the solution is computed by means of an iterative scheme, alternating between solving a series of simpler subproblems and updating the values of the Lagrange multipliers. The key to successfully applying LR lies in identifying the constraints to relax and defining an appropriate update strategy such that a lower bound is generated for a minimization problem (respectively, an upper bound is generated for a maximization problem) (Fisher, 2004). It is also helpful to have an upper bound value (in the case of a minimization) or a lower bound value (in the case of a maximization), both of which can reduce solution time considerably and are typically defined using heuristics. The continuity constraints (30) and (7) are ideal candidates for relaxation in a LR approach, since their relaxation leads to a decoupled series of sub-problems, each representing a time slot in the scheduling horizon. Note that the index g is used to denote the case where multiple products are produced. If only a single product is produced, g is of size 1 in the formulations below and can be neglected. The continuity constraints (30) and (7) are ideal candidates for relaxation in a LR approach, since their relaxation leads to a decoupled series of sub-problems, each representing a time slot in the scheduling horizon. Note that the index g is used to denote the case where multiple products are produced. If only a single product is produced, g is of size 1 in the formulations below and can be neglected. 1. Initialize: The problem is first solved to find initial values of the Lagrange multipliers, based on constraint violations. The complicating constraints (30) are removed in this initial solution: max Ko =



ui

i

j

inv  gij , w  gij ␺(Pricegij , w , Dg )

(44)

g

s.t. Process model (43 and/or 29) Inventory model (16), (18), (21), and (23) Initial Conditions Process Constraints Quality Constraints In some practical cases, it is beneficial to solve this problem as a Relaxed Mixed Integer Program (RMIP), where the problem is solved with integer variables relaxed to be continuous. 2. Calculate/update Lagrange multipliers: We begin by defining the parameter  i,m to describe the absolute value of the continuity constraint violation for a state xijm , where the m index represents the iteration number. The absolute value of the constraint violation is used to obtain certain guarantees about the optimality of the relaxed subproblem (PII) solution shown below, given in Theorem 1. Similar to the developments pos neg in (13a)–(13c), xi,m and xi,m represent the positive and negative (complicating) constraint violations, respectively: pos

neg

pos

neg

xi,m − xi,m = xi,m − xi−1,NJ ,m , xi,m − xi,m = xi,1,m − hNI ,m , pos

neg

xi,m , xi,m ≥ 0

i > 1 i=1

(45a) (45b) (45c)

pos

neg

i,m = xi,m + xi,m

47

(45d)

This value,  i,m , can then be used to calculate an initial stepsize (Stepsizei,0 ) and to update the values of the Lagrange multipliers at every iteration by a value denoted as Stepsizei,m : Stepsizei,0 =

Ko − LB

 + ||i,0 ||2

Stepsizei,m = m

Jm − Lm

 + ||i,m ||2

(46) (47)

where  is a small value which ensures the stepsize is finite for the case of no constraint violation, LB is the problem-dependent heuristically-set lower bound on the solution used to initialize a stepsize and shrink the feasible region, Lm is the value of the relaxed subproblem objective function, (49), and m is a tuning parameter (typically, a value m = 2 is chosen, and is reduced by half if there is no improvement between iterations). The above Stepsize calculation is based on the developments in Fisher (1985) and Bertsimas and Tsitsiklis (1997). The Lagrange multipliers are then updated within each iteration using: LDi,m = max(0, i,m−1 + Stepsizei,m ∗ i,m )

(48)

where i,m−1 is a parameter, typically the value of LDi,m−1 . A value i,0 > 0 is selected as an initial guess. 3. Solve the relaxed subproblem: The relaxed subproblem (PII) is formulated as: max Lm =



ui i

j

inv  gijm , w  gijm (Pricegijm , w , Dg ) − LDi,m i,m

(49)

g

s.t. Process model (43 and/or 29) Inventory model (16), (18), (21), and (23) Initial Conditions Process Constraints Quality Constraints 4. Check convergence: Several termination criteria are available for LR approaches. The below criterion takes into account both the optimality gap and the dualized constraint violations by measuring the change in the Lagrange multiplier (48) between iterations, m, where is the tolerance, typically on the order of 10−3 : |LDi,m−1 − LDi,m | ≤

(50)

If convergence is met, the iterations are terminated, if not, the algorithm returns to step 2 and the next iteration is started. Penalizing the absolute value of the continuity constraint violations through LR allows us to make claims about the optimality of the relaxed subproblem PII. While the below claims are supported by several previous works (Rush and Collins, 2014; Bertsekas, 1999), the theorem and proof below follow the results in Theorem 4.1 of Guignard (2003). Theorem 1. If PI is dynamically feasible, then the optimal solution to PII is optimal for PI. Proof. We assume a dynamically feasible solution exists to PI (i.e., state continuity constraints in the HW models can be met ( i,m = 0), and the production schedule is dynamically feasible). Then, when equality constraints  i,m = 0 are met, the minimum objectives of PI and PII are the same, L* = J*m , and the relaxed sub-problem PII is equivalent to the original problem PI. As such, the optimal solution to the relaxed sub-problem PII is optimal for the original problem PI.䊐

48

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52 Table 2 Summary of symbols and their definitions for the MMA case study full-order model (51) from Du et al. (2015). Symbol

Fig. 11. Depiction of the continuous MMA polymerization reactor.

7. Case study In this section, the previously derived methods are applied to a methyl methacrylate (MMA) process and the results are compared in terms of computational time required for scheduling with full-order and reduced-order models previously investigated by Du et al. (2015) and Daoutidis et al. (1990). The problems were implemented as an MILP in GAMS (GAMS, 2016) and solved to global optimality using CPLEX 12 to generate the optimal production sequence and timing. This optimal schedule was simulated on the full-order process model to generate state-trajectories and verify the accuracy and feasibility of the low-order scale-bridging models for optimal scheduling calculations. The computations were ® carried out on a PC running Windows 10 64 bit, having a 3.4 GHz ® Intel Core i7 processor with 16 GB of RAM and 8 threads. 7.1. Process description The continuous polymerization process presented by Daoutidis et al. (1990) is considered. Free-radical polymerization of methyl methacrylate (MMA) takes place in a continuously stirred tank reactor (CSTR) (Fig. 11), with azo-bis-isobutyronitrile as initiator and toluene solvent. A cooling jacket removes the heat generated by the exothermic reactions. The material and energy balance equations are given below (51), originally from Daoutidis et al. (1990) and rely on the mechanism presented by Congalidis et al. (1989): dC m =− dt





Zp exp

dC i = −ZI exp dt dD0 =− dt



dD1 = Mm dt

dt

P0 =

=



EI RT



 

−Ep RT



−Efm



Cm Po +

RT



V



+ ZTd exp



 + Zfm exp

−ETd RT

−Efm

 Cm P0 +

 Cm P0 −

RT

F D0 V

F D0 V

(51)

− Hp UA F (T − Tja ) + (Tin − T ) P0 − cp cp V V

UA F cw (Tw0 − Tja ) + (T − Tja ) V0 w cw V0





F in (C − Cm ) V m

F I CIin − FC I

−ETc RT

−Ep RT

Cm



+ Zfm exp

CI +

Zp exp







0.5ZTc exp

dT = Zp exp dt dT ja



−Ep RT

 

2aC I ZI exp −EI /RT

ZTd exp −ETd /RT



+ ZTC exp −ETC /RT

0.5 

The set of five coupled differential equations and P0 definition (51) comprise the full-order model for the MMA polymerization model. The variables in these equations are defined in Table 2. For a more in-depth discussion of the full-order model we direct the reader to the work by Daoutidis et al. (1990). Based on the operating conditions (namely, reactor temperature, T), the reactor can produce four product grades g = {A, B, C, D}, defined in terms of the molecular weight distribution . The respec-

Definition

Continuous variables Cm Monomer concentration in Inlet monomer concentration Cm Inlet initiator concentration Ciin CI Initiator concentration Molar concentration of dead polymer chains D0 D1 Mass concentration of dead polymer chains T Reactor temperature Jacket temperature Tja Inlet cooling water temperature Tw0 Inlet feedstream water temperature Tin Initiator flowrate FI Reactor feed flowrate F cw F Cooling water flowrate Molecular weight distribution  Parameters and coefficients Z Reaction rates Activation energies E Reactor volume V Cooling jacket volume V0 Heat transfer coefficient U A Heat transfer area between the reactor and the jacket Density of reaction mixture  w Density of the cooling water Concentration of live polymer chains P0 Heat capacity of the cooling water cw Heat capacity of the reaction mixture cp Symbols p Propagation I Initiation Chain transfer to monomer fm Termination via coupling TC Termination via disproportionation Td

tive setpoints for these variables are represented using T¯ and . ¯ The molecular weight distribution  is described below:  = D1 /D0

(52)

where D0 and D1 are the molar and mass concentrations of the dead polymer chains, respectively. The feedstock is the same for all products. Grade transitions are accomplished by manipulating the initiator (FI ) and cooling water (Fcw ) flow rates. The initiator and cooling water feed flow rates are constrained to be lower than 0.1 m3 /h and 4 m3 /h, respectively. The amount of polymer grade g produced (pg ) is based on demand Dg , which is assumed to be satisfied at the end of the time horizon and therefore is not time-based: Dg ≤ pg ≤ 1.1Dg

(53)

The aim of the scheduling problem is to determine the optimal (cyclic) production sequence, i.e., the order of grades g, the producP  ), and the transition time between grades (t trans ), tion time (tgP , tgi i such that profit is maximized and demand (Dg ) is met. Product demand and price information is provided in Table 3. Based on the above, the vector of relevant operating costs is: Pricegij = {Selling Price, Inventory Cost}. 7.1.1. Model identification In this section, we cast the scheduling problem for the MMA process in the form (27). Here, the scheduling relevant variables  = {, T, F I , F cw }. The variables  and T define the transition are: w and production times in each scheduling time slot, and therefore directly affect the objective function, while FI and Fcw are subject to constraints during normal operation. In order to generate the relevant SBMs, we proceeded as follows:

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

49

Table 3 Cost and production data based on Du et al. (2015). Grade, g

Demand (kg)

Price ($/kg)

Inventory cost ($/kg/h)

 ¯ (kg/mol)

T¯ (K)

A B C D

9550 12,220 13,220 14,150

0.9 1.2 1.3 1.5

0.010 0.011 0.012 0.05

44,032 54,798 76,835 95,820

329.8 328.8 327.5 326.9

Table 4 Sample times for HW model discretization and FSR model identification. Variable

TJ (min)

Fijcw FijI

5.96

Tij ij

2.68 2.68

Table 7 Summary of scheduling calculations introduced in this work, where P3–P6 are all formulated on an MILP framework. Label

Table 5 HW model statistics for MMA case study, where “output” specifies the variable modeled and “input” the model setpoint, ui . Input

Output

H(u) breakpoints

State space model order

W(y) breakpoints

% Fit

T¯i T¯i T¯i T¯i

Fcw FI T 

4 4 1 4

2 2 2 2

10 13 0 0

99.79 99.35 99.93 99.89

Fijcw

99.90

FijI

99.78

Tij ij

99.57 99.49

%Fit = 100

ref

wij −





ref wij

wij

F , F , T,  Fcw , FI , T,  T,  None

FSR models

Models with LR

None None Fcw , FI , Fcw , FI , T, 

Fcw , FI , T,  None T,  N/A

7.1.2. Problem formulations Problem (49) was solved using the data-driven SBMs described above; for comparison purposes, the results using the full-order model (i.e., problem (26)) and the results without LR (27) are included, where a combination of HW/FSR models was used to represent dynamics in the scheduling calculation. The six problems compared are referenced using the notation P1–P6.

max

zgi ,t s ,t e ,t P ,t trans ,ui i

(i) The full-order model of the process as described by Daoutidis et al. (1990) (51) was simulated under the input–output linearizing controller described by Du et al. (2015). A series of stepwise set-point changes were imposed to generate training and validation data sets. (ii) These data were used to identify the data driven SBMs using the sample times, TJ within the range specified in (4), given in Table 4. Hammerstein models were identified for T and  given the simple dynamics involved, and HW models were identified for Fcw and FI . Model size and fits are give in Table 5, where fit was calculated using (54):



I

• Problem P1 relies on a full-order process model embedded in scheduling calculations (Du et al., 2015) on an MINLP framework:

Table 6 FSR model statistics for MMA case study. % Fit

cw

P3 P4 P5 P6

5.96

Variable

HW models

(54)

The trapezoid rule was used to estimate integrals in (54) and ref wij are the model test data. (iii) The continuous-time state-space models were discretized with sample times, TJ , given in Table 4. The piecewise constant Hammerstein functions, H(u), were linearized and discretized using an SOS1 reformulation since the process features discrete products. The output nonlinearities (piecewise linear Wiener functions, W(y)) were linearized using SOS2 variables. (iv) FSR models were developed following the procedure given in Section 4.1.2 using the sampling times given in Table 4. Model fits to test data are given below, as calculated by (54) (Table 6). (v) The HW and FSR models were subsequently reformulated using the ideas introduced in Sections 5.1 and 5.2.

i

Profit

(55)

i

i

s.t. Scheduling constraints Full process model (Fijcw , FijI , Tij , ij ) Inventory model Initial Conditions Process Constraints Quality Constraints • Problem P2 relies on SBMs that reflect the linear behavior imposed by the input–output feedback linearizing controller described by Du et al. (2015), and is formulated as a MINLP: max

zgi ,t s ,t e ,t P ,t trans ,ui i

i

i

Profit

(56)

i

s.t. Scheduling constraints SBM (Fijcw , FijI , Tij , ij ) Inventory model Initial Conditions Process Constraints Quality Constraints • Problems P3–P6 are based on the data-driven SBMs introduced in this work (Table 7) and follow the format in (49). Remark 3. In implementing the proposed LR scheme, a lower bound (LB) on the profit was determined by using the profit associated with producing only the lowest grade of polymer, i.e., the least profitable, for the estimated (maximum) TC . 7.1.3. Results and discussion Problems P1–P6 were solved, and the optimal time-varying values of the decision variables were used to simulate the optimal schedule on the full-order model (51). Table 8 gives the model sizes and the results of P1–P6, showing the cost estimated by the optimization in comparison to the actual cost provided by the full-order model, in addition to the total cycle time.

50

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

Table 8 Computational statistics for the MMA case study, where P3-P6 solved with an optimality gap of 0.00%. Objective function values were re-calculated by simulating the optimal schedule on the full-order process model. Problem

Equations

Continuous variables

Discrete variables

CPU (s)

Obj ($)

Predicted obj ($)

TC (h)

P1 P2 P3 P4 P5 P6

5040 1100 25,083 25,050 22,842 17,514

4605 1025 22,376 22,279 20,071 14,695

16 16 1360 1360 1360 1360

89.77 6.970 12.28 34.39 1.53 1.58

55,407 56,188 56,832 56,832 56,832 56,015

57,627 55,677 57,300 57,300 57,300 57,482

47.45 46.72 48.64 48.64 48.64 47.93

Fig. 12. Variable trajectories for the MMA case study for P5.

The number of discrete variables in P3–P6 increased in comparison to P1 and P2 due to the use of RIE for determining transition time (additional binary variable, bij ). However, the solution strategy we propose lowered the computation times of P3–P6 below that of P1, P5–P6 were solved faster than P2. We note that there are discrepancies between the predicted profit (calculated by the optimization problem) versus actual profit (calculated by simulating the proposed scheduling using the full-order model). These can be attributed to several factors; in the case of P1, the equivalent of an optimization-based controller that outperforms the inversion-based controller implemented in the simulation is used in the schedule optimization problem. Further, for all problems (P1–P6), there is a difference in the discrete-time numerical representation of the model, where the schedule optimization problem utilizes collocation of finite elements (P1–P2) and uniform time steps (P3–P6), while and the simulation uses a variable-step time integration method. Finally, there are inherent discrepancies in the way transition times are estimated, whereby the percent offset

(with a threshold set at 1% of the setpoint) is used in all simulations, P3-P6 rely on the RIE method described above, and the approach described by Flores-Tlacuahuac and Grossmann (2006) is used for P1–P2. The schedules in P1–P6 are in agreement, suggesting the sequence g = {A, B, C, D}, with P3–P6 improving on the objectives in P1–P2. In comparing P3 to P4, LR reduces computation time by 64.3% and produces the same objective function value (Theorem 1). Furthermore, the solution to P5 is equivalent to P3–P4, which is reasonable given that in P3–P5, the same HW models were utilized to represent the two variables used to determine transition time (via RIE) in the objective function (ij , Tij ). The use of FSR models (P5, P6) clearly reduces computation time, with P5–P6 having the lowest model dimensions and shortest computation times. It should, however, be emphasized again that FSR models do not actually capture gain nonlinearity. While P5–P6 are solved faster than P3–P4, the FSR models in P5–P6 are only applicable for variables which reach steady-state before the end of the

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

scheduling time slot, and for which dynamics are linear. Therefore, in the case of transient operation, HW models are the more attractive option. P5 combines both HW and FSR models has the lowest solution time and highest objective function value, reducing computation time from the full-order method (P1) by 98.3%. This suggests that a judicious consideration of model forms results in both fast solution times and high accuracy. Plots of the scheduling for the optimal schedule in P5 are included relevant variables w below. Fig. 12 shows that the setpoints were tracked well by T and . In addition, no constraint violations occurred for Fcw or FI . The predicted variable trajectories are comparable to the actual trajectories obtained by simulating the full-order model, showing the ability of the SBMs to capture dynamics in a computationally parsimonious way.

8. Conclusions In this work, we present a novel method for optimal scheduling under dynamic constraints for continuous processes operating under highly variable conditions. Our approach is based on using low-order, data-driven models of the closed-loop process dynamics, which we refer to as “scale-bridging models (SBMs).” We relied on specific SBM formats selected from amongst models that are widely-used in the control community (namely, Hammerstein–Wiener and finite step response models). A key contribution of this work was to show that such models can be linearized (often exactly) for the purpose of incorporating them in optimal production scheduling problems. While counter-intuitive, this effort results in scheduling problems that are in a MILP format, and can be solved readily using currently available optimization solvers. These scheduling problems account explicitly for the process dynamics and the performance of the control system, and can therefore be counted on to preserve the closed-loop stability properties of the process. The theoretical developments were illustrated using a case study considering the cyclic scheduling of a polymerization process. Our approach was compared to previous results, evincing excellent computational performance improvements (reductions by over 98.3% of CPU time) while providing the same (or improved) economic results. Along with the fact that the process models can be identified from closed-loop operating data, this computational performance is extremely encouraging for carrying out industrial implementations of these concepts, which we will pursue in our future work.

Disclaimer This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

51

Acknowledgements The authors gratefully acknowledge partial financial support from the US Department of Energy, under award DE-OE0000841, and from the National Science Foundation (NSF) through the CAREER Award 1454433 and Award CBET-1512379.

References 2015. IBM ILOG CPLEX Optimization Studio CPLEX User’s Manual. IBM, Armonk, NY. BDMLP Solver, 2017. GAMS User Guide, Chapter Solvers, 24.7 ed, Washington, DC. 2017. XPRESS Solver Engine. FrontlineSolvers, Incline Village, NV. Aryani, D., Wang, L., Patikirikorala, T., 2017. On identification of Hammerstein and Wiener model with application to virtualised software system. Int. J. Syst. Sci. 48 (6), 1146–1161, http://dx.doi.org/10.1080/00207721.2016.1244303, ISSN 0020-7721. Bai, E.-W., 2002. A blind approach to the Hammerstein–Wiener model identification. Automatica 38 (6), 967–979, http://dx.doi.org/10.1016/S00051098(01)00292-8, ISSN 00051098. Baldea, M., Daoutidis, P., 2012. Dynamics and Nonlinear Control of Integrated Process Systems. Cambridge University Press, New York City, ISBN 978-0-521-19170-8. Baldea, M., Harjunkoski, I., 2014. Integrated production scheduling and process control: a systematic review. Comput. Chem. Eng. 71, 377–390, http://dx.doi. org/10.1016/j.compchemeng.2014.09.002, ISSN 00981354. Baldea, M., Du, J., Park, J., Harjunkoski, I., 2015. Integrated production scheduling and model predictive control of continuous processes. AIChE J. 61 (12), 4179–4190, http://dx.doi.org/10.1002/aic.14951. Beal, L.D., Park, J., Petersen, D., Warnick, S., Hedengren, J.D., 2017. Combined model predictive control and scheduling with dominant time constant compensation. Comput. Chem. Eng. 104, 271–282, http://dx.doi.org/10.1016/j.compchemeng. 2017.04.024, ISSN 00981354. Beale, E., Tomlin, J., 1970. Special facilities in a general mathematical programming system for nonconvex problems using ordered sets of variables. In: Proceedings of the 5th International Conference on Operational Research, London. Bertsekas, D., 1999. Nonlinear Programming. Athena Scientific, Belmont, MA. Bertsimas, D., Tsitsiklis, J.N., 1997. Introduction to Linear Optimization, vol. 30. John Wiley & Sons, Inc., Hoboken, NJ, USA, ISBN 1886529191. Bloemen, H.H.J., Van Den Boom, T.J.J., Verbruggen, H.B., 2001. Model-based predictive control for Hammerstein–Wiener systems. Int. J. Control 74 (5), 482–495, http://dx.doi.org/10.1080/00207170010014061, ISSN 0020-7179. Braun, M.W., Rivera, D.E., Stenman, A., Foslien, W., Hrenya, C., 1999. Multi-level pseudo-random signal design and “model-on-demand” estimation applied to nonlinear identification of a RTP wafer reactor. In: American Control Conference (ACC). IEEE, pp. 1573–1577. Chu, Y., You, F., 2012. Integration of scheduling and control with online closed-loop implementation: fast computational strategy and large-scale global optimization algorithm. Comput. Chem. Eng. 47, 248–268, http://dx.doi.org/10. 1016/j.compchemeng.2012.06.035, ISSN 00981354. Chu, Y., You, F., 2014. Integrated planning, scheduling, and dynamic optimization for batch processes: MINLP model formulation and efficient solution methods via surrogate modeling. Ind. Eng. Chem. Res. 53 (34), 13391–13411, http://dx. doi.org/10.1021/ie501986d, ISSN 0888-5885. Congalidis, J.P., Richards, J.R., Ray, W.H., 1989. Feedforward and feedback control of a solution copolymerization reactor. AIChE J. 35 (6), 891–907, http://dx.doi. org/10.1002/aic.690350603, ISSN 0001-1541. Daoutidis, P., Soroush, M., Kravaris, C., 1990. Feedforward/feedback control of multivariable nonlinear processes. AIChE J. 36 (10), 1471–1484, http://dx.doi. org/10.1002/aic.690361003, ISSN 0001-1541. Dayal, B.S., MacGregor, J.F., 1996. Identification of finite impulse response models: methods and robustness issues. Ind. Eng. Chem. Res. 35 (11), 4078–4090, http://dx.doi.org/10.1021/ie960180e, ISSN 0888-5885. Diangelakis, N.A., Pistikopoulos, E.N., 2017. A multi-scale energy systems engineering approach to residential combined heat and power systems. Comput. Chem. Eng. 102, 128–138, http://dx.doi.org/10.1016/j.compchemeng. 2016.10.015, ISSN 00981354. Du, J., Park, J., Harjunkoski, I., Baldea, M., 2015. A time scale-bridging approach for integrating production scheduling and process control. Comput. Chem. Eng. 79, 59–69, http://dx.doi.org/10.1016/j.compchemeng.2015.04.026, ISSN 00981354. Elfeky, E.Z., Sarker, R.A., Essam, D.L., 2008. Partial decomposition and parallel GA (PD-PGA) for constrained optimization. In: 2008 IEEE International Conference on Systems, Man and Cybernetics. IEEE, pp. 220–227, http://dx.doi.org/10. 1109/ICSMC.2008.4811278, ISBN 978-1-4244-2383-5. Engell, S., Harjunkoski, I., 2012. Optimal operation: scheduling, advanced control and their integration. Comput. Chem. Eng. 47, 121–133, http://dx.doi.org/10. 1016/j.compchemeng.2012.06.039, ISSN 00981354. Eriksson, K., Estep, D., Johnson, C., 2004. Applied Mathematics: Body and Soul. Springer Berlin Heidelberg, Berlin Heidelberg, ISBN 978-3-642-05658-1. Fisher, M.L., 1985. Applications oriented guide to Lagrangian relaxation. Interfaces 15 (2), 10–21, http://dx.doi.org/10.1287/inte.15.2.10, ISSN 00922102.

52

M.T. Kelley et al. / Computers and Chemical Engineering 110 (2018) 35–52

Fisher, M.L., 2004. The Lagrangian relaxation method for solving integer programming problems. Manag. Sci. 50 (12), 1861–1871, http://dx.doi.org/10. 1287/mnsc.1040.0263, ISSN 0025-1909. Flores-Tlacuahuac, A., Grossmann, I.E., 2006. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 45 (20), 6698–6712, http://dx.doi.org/10.1021/ie051293d, ISSN 0888-5885. Flores-Tlacuahuac, A., Grossmann, I.E., 2010. Simultaneous scheduling and control of multiproduct continuous parallel lines. Ind. Eng. Chem. Res. 49 (17), 7909–7921, http://dx.doi.org/10.1021/ie100024p, ISSN 0888-5885. GAMS, 2016. General Algebraic Modeling System (GAMS). Release 24.7.4. Grossmann, I., 2005. Enterprise-wide optimization: a new frontier in process systems engineering. AIChE J. 51 (7), 1846–1857, http://dx.doi.org/10.1002/aic. 10617, ISSN 0001-1541. Guignard, M., 2003. Lagrangian relaxation. Top 11 (2), 151–200, http://dx.doi.org/ 10.1007/BF02579036, ISSN 1134-5764. Hahn, J., Edgar, T.F., 2002. An improved method for nonlinear model reduction using balancing of empirical Gramians. Comput. Chem. Eng. 26 (10), 1379–1397, http://dx.doi.org/10.1016/S0098-1354(02)00120-5, ISSN 00981354. Keha, A.B., de Farias, I.R., Nemhauser, G.L., 2004. Models for representing piecewise linear cost functions. Oper. Res. Lett. 32 (1), 44–48, http://dx.doi.org/10.1016/ S0167-6377(03)00059-2, ISSN 01676377. Lang, Y.-d., Malacina, A., Biegler, L.T., Munteanu, S., Madsen, J.I., Zitney, S.E., 2009. Reduced order model based on principal component analysis for process simulation and optimization. Energy Fuels 23 (3), 1695–1706, http://dx.doi. org/10.1021/ef800984v, ISSN 0887-0624. Ljung, L., 1987. Ljung L System Identification Theory for User.pdf, ISSN 00051098. http://linkinghub.elsevier.com/retrieve/pii/0005109889900198. Maravelias, C.T., 2012. General framework and modeling approach classification for chemical production scheduling. AIChE J. 58 (6), 1812–1828, http://dx.doi. org/10.1002/aic.13801, ISSN 00011541. MATLAB, 2016. MATLAB 2016a. Narendra, K., Gallman, P., 1966. An iterative method for the identification of nonlinear systems using a Hammerstein model. IEEE Trans. Autom. Control 11 (3), 546–550, http://dx.doi.org/10.1109/TAC.1966.1098387, ISSN 0018-9286. Nie, Y., Biegler, L.T., Wassick, J.M., 2012. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 58 (11), 3416–3432, http://dx.doi.org/10.1002/aic.13738, ISSN 00011541. Nie, Y., Biegler, L.T., Villa, C.M., Wassick, J.M., 2013. Reactor modeling and recipe optimization of polyether polyol processes: polypropylene glycol. AIChE J. 59 (7), 2515–2529, http://dx.doi.org/10.1002/aic.14144, ISSN 00011541. Nie, Y., Biegler, L.T., Wassick, J.M., Villa, C.M., 2014. Extended discrete-time resource task network formulation for the reactive scheduling of a mixed batch/continuous process. Ind. Eng. Chem. Res. 53 (44), 17112–17123, http:// dx.doi.org/10.1021/ie500363p, ISSN 0888-5885. Nyström, R.H., Franke, R., Harjunkoski, I., Kroll, A., 2005. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 29 (10), 2163–2179, http://dx.doi.org/10.1016/j. compchemeng.2005.07.006, ISSN 00981354. Ogunnaike, B., Harmon Ray, W., 1994. Process Dynamics, Modeling, and Control. Oxford University Press, New York City. Park, J., Du, J., Harjunkoski, I., Baldea, M., 2014. Integration of scheduling and control using internal coupling models. In: 24th European Symposium on Computer Aided Process Engineering, Pts A and B, Budapest, pp. 529–534. Pattison, R.C., Touretzky, C.R., Johansson, T., Harjunkoski, I., Baldea, M., 2016. Optimal process operations in fast-changing electricity markets: framework for scheduling with low-order dynamic models and an air separation application. Ind. Eng. Chem. Res. 55 (16), 4562–4584, http://dx.doi.org/10. 1021/acs.iecr.5b03499, ISSN 0888-5885.

Pattison, R.C., Touretzky, C.R., Harjunkoski, I., Baldea, M., 2017. Moving horizon closed-loop production scheduling using dynamic process models. AIChE J. 63 (2), 639–651, http://dx.doi.org/10.1002/aic.15408, ISSN 00011541. Patwardhan, R.S., Lakshminarayanan, S., Shah, S.L., 1998. Constrained nonlinear MPC using Hammerstein and Wiener models: PLS framework. AIChE J. 44 (7), 1611–1622, http://dx.doi.org/10.1002/aic.690440713, ISSN 00011541. Pinto, J.M., Grossmann, I.E., 1994. Optimal cyclic scheduling of multistage continuous multiproduct plants. Comput. Chem. Eng. 18 (9), 797–816, http:// dx.doi.org/10.1016/0098-1354(93)E0021-Z, ISSN 00981354. Prata, A., Oldenburg, J., Kroll, A., Marquardt, W., 2008. Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor. Comput. Chem. Eng. 32 (3), 463–476, http://dx.doi.org/10.1016/j. compchemeng.2007.03.009, ISSN 00981354. Rush, A.M., Collins, M., 2014. A tutorial on dual decomposition and Lagrangian relaxation for inference in natural language processing. J. Artif. Intell. Res. 45, 305–362, http://dx.doi.org/10.1613/jair.3680, ISSN 10769757. Seborg, D.E., Mellichamp, D.A., Edgar, T.F., Doyle, F.J., 2010. Process Dynamics and Control. John Wiley & Sons, Inc., Hoboken, NJ, USA, http://dx.doi.org/10.1007/ s13398-014-0173-7.2, ISBN 9780470128671. Tong, C., Palazoglu, A., El-Farra, N.H., Yan, X., 2015. Energy demand management for process systems through production scheduling and control. AIChE J. 61 (11), 3756–3769, http://dx.doi.org/10.1002/aic.15033, ISSN 00011541. Touretzky, C.R., Baldea, M., 2014. Integrating scheduling and control for economic MPC of buildings with energy storage. J. Process Control 24 (8), 1292–1300, http://dx.doi.org/10.1016/j.jprocont.2014.04.015, ISSN 09591524. Touretzky, C.R., Harjunkoski, I., Baldea, M., 2016. A framework for integrating scheduling and control using discrete-time dynamic process models. ESCAPE. Touretzky, C.R., Harjunkoski, I., Baldea, M., 2017. Dynamic models and fault diagnosis-based triggers for closed-loop scheduling. AIChE J. 63 (6), 1959–1973, http://dx.doi.org/10.1002/aic.15564, ISSN 00011541. Wang, D.-Q., Ding, F., 2012. Hierarchical least squares estimation algorithm for Hammerstein–Wiener systems. IEEE Signal Process. Lett. 19 (12), 825–828, http://dx.doi.org/10.1109/LSP.2012.2221704, ISSN 1070-9908. Wang, G.G., Shan, S., 2007. Review of metamodeling techniques in support of engineering design optimization. J. Mech. Des. 129 (4), 370, http://dx.doi.org/ 10.1115/1.2429697, ISSN 10500472. Wills, A., Schön, T.B., Ljung, L., Ninness, B., 2013. Identification of Hammerstein–Wiener models. Automatica 49 (1), 70–81, http://dx.doi.org/10. 1016/j.automatica.2012.09.018, ISSN 00051098. Wise, B.M., Ricker, N.L., 1993. Identification of finite impulse response models with continuum regression. J. Chemometr. 7 (1), 1–14, http://dx.doi.org/10.1002/ cem.1180070102, ISSN 0886-9383. Yu, M., Miller, D.C., Biegler, L.T., 2015. Dynamic reduced order models for simulating bubbling fluidized bed adsorbers. Ind. Eng. Chem. Res. 54 (27), 6959–6974, http://dx.doi.org/10.1021/acs.iecr.5b01270, ISSN 0888-5885. Zhang, Q., Grossmann, I.E., 2016. Planning and Scheduling for Industrial Demand Side Management: Advances and Challenges. Springer International Publishing, Cham, http://dx.doi.org/10.1007/978-3-319-28752-2 14, ISBN 9783319287522. Zhuge, J., Ierapetritou, M.G., 2012. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 51 (25), 8550–8565, http:// dx.doi.org/10.1021/ie3002364, ISSN 0888-5885. Zhuge, J., Ierapetritou, M.G., 2016. A decomposition approach for the solution of scheduling including process dynamics of continuous processes. Ind. Eng. Chem. Res. 55 (5), 1266–1280, http://dx.doi.org/10.1021/acs.iecr.5b01916, ISSN 0888-5885.