Production optimization for continuously operated processes with optimal operation and scheduling of multiple units

Computers and Chemical Engineering 30 (2006) 392–406

Production optimization for continuously operated processes with optimal operation and scheduling of multiple units Rasmus H. Nystr¨om∗ , Iiro Harjunkoski, Andreas Kroll Subdepartment of Plant Performance Improvement, ABB Corporate Research, D-68526 Ladenburg, Germany Received 1 June 2004; received in revised form 23 September 2005; accepted 27 September 2005 Available online 29 November 2005

Abstract A production optimization problem for continuously operated processes is presented and a solution strategy is proposed. The strategy consists of decoupling the complete problem into a mixed-integer linear programming (MILP) scheduling problem, including sequencing and allocation, and a multi-stage dynamic optimization (DO) problem, including the determination of optimal trajectories and setpoints. The scheduling problem is treated as a master problem and the DO problem as a primal problem, and the complete problem is solved through iteration between the two. The approach is similar in nature to standard methods for solving mixed-integer nonlinear (MINLP) problems, such as Outer Approximation and Benders Decomposition, but more adapted to the specific problem, by permitting more freedom in choosing the binary representation. The decomposition strategy implies flexibility in choosing the optimization tools required, and enables the treatment of larger problems. The approach is a generalization of a previously reported one, where only the single-unit case was discussed. Splitting the primal problem into smaller DO subproblems, parameter estimation from the DO subproblems, termination criteria and other topics are discussed. The target process for demonstration of the method is an industrial polymerization process. © 2005 Elsevier Ltd. All rights reserved. Keywords: Dynamic model-based scheduling; Multi-stage optimization; Product sequencing; Optimal manufacturing; Multiple-unit production optimization

1. Introduction Significant cost savings can be achieved and/or production efficiency can be increased through appropriate production optimization and scheduling for a production plant. A scheduling application for a production plant should provide optimal or near-optimal solutions to sequencing and allocation problems, taking the typically nonlinear characteristics of the plant into account. For a continuously operated production plant, the key characteristics for a scheduling application consist of parameters such as transition times and costs between products, as well as costs and throughputs related to the production of a given product on a given piece of equipment. A production optimization Abbreviations: BD, Benders Decomposition; DO, dynamic optimization; LB, lower bound; MIDO, mixed-integer dynamic optimization; MILP, mixedinteger linear programming/problem; MINLP, mixed-integer nonlinear programming/problem; MV1, MV2, MV3, names of manipulated process inputs in example; OA, Outer Approximation; QV1, QV2, names of quality variables in example; UB, upper bound ∗ Corresponding author. Tel.: +49 6203 71 6354; fax: +49 6203 71 6253. E-mail addresses: [email protected] (R.H. Nystr¨om), [email protected] (I. Harjunkoski), [email protected] (A. Kroll). 0098-1354/$ – see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2005.09.009

application should also consider higher-level information such as amounts ordered, supply of raw material, due dates, release dates, etc. The higher-level information is typically specified by customer orders, whereas the lower-level information is connected with the process dynamics and product specifications. The trend in today’s production industry is towards increased flexibility and order-driven production, see for instance studies by Sinclair (1987) and Backx, Bosgra, and Marquardt (1998). In response to the higher demand for flexibility, detailed product specifications may be part of the customer orders. Since the key characteristics needed for scheduling depend on the (possibly time varying) process, as well as the varying and possibly very numerous customer orders, storing all values required in a database might be infeasible. An alternative to this is a modelcentric approach, where a dynamic model is used as a central source of the data needed, and dynamic simulation is used to extract this data. Performing dynamic simulation is only possible when the inputs to the process are known. In this case, it means that the recipes for changing between products, as well as producing the products are known beforehand. In a scenario, where the overall production is to be as economical as possible, these recipes

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

Nomenclature Indices and constants indicating number m, n unit (machine) indices M constant, number of units (machines) p, b product indices P constant, number of products q quality variable index Q number of quality variables s, z stage indices Sm constant, number of stages for unit m Parameters and variables in overall/primal problem Dp constant, due date of product p fm nonlinear function describing dynamics of unit m gm nonlinear function describing outputs of unit m hm nonlinear function describing cost functions on unit m Js,m cost contribution of stage s on unit m Lm constant vector, general output lower limits for unit m Lmin,q constant, general lower bound of quality variable q Lp,q constant, lower bound of quality variable q for product p rp,m continuous auxiliary variable calculated from master problem solution and used in primal problem, minimum amount to be produced of product p on unit m Rp constant, required amount to be produced of product p t time variable Tmin constant, minimum value of make span Ts,m decision variable in primal problem, length (duration) of stage s on unit m Um constant vector, general output upper limits for unit m Umax,q constant, general upper bound of quality variable q us,m process input of unit m in stage s Up,q constant, upper bound of quality variable q for product p xs,m process state of unit m in stage s ys,m process output of unit m in stage s ys,m,gen general process output of unit m in stage s ys,m,q quality variable q of unit m in stage s
constant, termination tolerance of iterative optimization ϕs,m amount of product produced in stage s on unit m ϕ˙ s,m production rate in stage s on unit m ξs,m,p binary time-slot variable used to introduce overall problem, product p is/is not produced in stage s on unit m τs normalized integration time variable in stage s

393

Parameters and variables in master problem cp,m continuous variable in master problem, completion time of product p on unit m dp,m continuous variable in master problem, duration of product p on unit m, also used as auxiliary variable in primal problem dmin constant, minimum duration of a job I set of invalid or undesirable transitions M1 , M2,p constants, upper bounds chosen as tight as possible (“big-M” constants) Mp set of units m that can be used to produce product p tp,m continuous variable in master problem, starting time of product p on unit m T continuous variable in master problem, make span αp,m binary variable in master problem, indicates whether product p is produced on unit m or not βp,b,m binary variable (sequencing variable) in master problem, indicates whether product p directly precedes product b on unit m, β0,p,m = 1 means that product p is first on unit m, βp,P+1,m = 1 means that product p is last on unit m γp,b cost of a transition from product p to product b (assumed same for all units), constant master problem parameter estimated from the primal problem solution γ0,p,m cost of a transition from initial state on unit m to product p, constant master problem parameter estimated from the primal problem solution γ˙ p cost per time for production of product p, constant master problem parameter estimated from the primal problem solution δp,b duration of a transition from product p to product b (assumed same for all units), constant master problem parameter estimated from the primal problem solution δ0,p,m duration of a transition from initial state on unit m to product p, constant master problem parameter estimated from the primal problem solution ϕ˙ p produced amount per time of product p, constant master problem parameter estimated from primal problem solution

are part of the total optimization problem. This means that, for extracting data from the dynamic model, mere dynamic simulation will not suffice. Instead, dynamic optimization (DO) problems can be solved to obtain the information required, see e.g. (Bhatia & Biegler, 1996). These optimization problems can be expressed as so-called multi-stage dynamic optimization problems, studied e.g. by Vassiliadis, Sargent, and Pantelides (1994). From the above considerations, it is clear that the scheduling problem and the dynamic optimization problems are coupled. The results from the scheduling level have an impact on the dynamic optimization level, and vice versa. Thus, the rigorous

394

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

approach to production optimization is to treat the two as a complete problem. If the dynamic behaviour of the production plant is not decisive to the outcome of the production planning, nonlinear static models can be used together with discrete-time modeling to express the production campaign as a mixed-integer nonlinear problem (MINLP), see e.g. (Smania & Pinto, 2003). However, in many cases the transient behaviour is of high importance, e.g. in the polymer industry. Then, the complete production optimization problem, including sequencing, allocation and optimal operation decisions is a mixed-integer dynamic optimization (MIDO) problem. A subclass of the MIDO problems can be solved e.g. using the mixed-logic dynamic (MLD) optimization framework, see work e.g. by Bemporad and Morari (1999) or Gallestey, Stothert, Castagnoli, Ferrari-Trecate, and Morari (2003). The MLD framework assumes that the problem is linear, linear time varying, or piecewise linear/affine. If the problem is highly nonlinear and cannot readily be approximated by piecewise affine functions, an alternative is to use conventional methods for solving MINLP problems, most typically Benders Decomposition (BD) and Outer Approximation (OA), described in overviews e.g. by Grossmann (2002) and Floudas (1995). In these methods, the total problem is solved through iteration between a master problem providing a lower bound on the optimal objective, and a primal problem providing an upper bound on the optimal objective. The tool DICOPT (Viswanathan & Grossmann, 1990) can be used for solving MINLP methods using these methods. Studies using BD for the MIDO problem class have been done e.g. by Mohideen, Perkins, and Pistikopoulos (1996). OA for MIDO has been investigated e.g. by Oldenburg, Marquardt, Heinz, and Leineweber (2003) and Prata, Oldenburg, Marquardt, Nystr¨om, and Kroll (2005) using the tool DyOS (Brendel, Oldenburg, Schlegel, & Stockmann, 2003). An alternative approach to solving MIDO problems using BD has been presented by Bansal, Sakizlis, Ross, Perkins, and Pistikopoulos (2003). A study on using the MINLP method Extended Cutting Plane (ECP) for a dynamic problem including integer decisions has been done by Emet and Westerlund (2004). In the same paper, a comparison is made between the performances of some common methods for the problem in question. The example process of this study is a continuously operated polymerization process. Grade transition optimization in polymerization processes has been a topic of frequent studies, see for instance those by Van Brempt et al. (2003), McAuley and MacGregor (1992), and Cervantes, Tonelli, Brandolin, Bandoni, and Biegler (2002). Feather, Harrell, Lieberman, and Doyle (2004) have approached grade transition control using model predictive control (MPC) based on mixed-integer linear programming (MILP). The plant nonlinearity is taken into account through regime-based switching between linear models, and other discrete events are addressed as well. Mahadevan, Doyle, and Allcock (2002) studied how control-relevant objectives and accounting for model uncertainty affect the selection of an optimal production sequence (production wheel), when the plant is assumed to be controlled by a linear controller. Another study on solving economically optimal grade transition problems together with sequencing problems for a polymer reactor has been done by Tousain (2002). Here the parameters required

for scheduling are optimized once and for all and stored in a database for future use. Chatzidoukas, Perkins, Pistikopoulos, and Kiparissides (2003a) solved a MIDO problem for a polymerization process, containing combined transition optimization and optimal control structure selection. In (Chatzidoukas et al., 2003b)Chatzidoukas, Kiparissides, Perkins, & Pistikopoulos, 2003b), the same authors consider the production sequencing as part of the MIDO problem. This paper presents a production optimization MIDO problem including sequencing and allocation for parallel reactor lines, and an approach to solving it. The overall MIDO problem is separated into two levels, in likeness with OA and BD. On one hand, there is a scheduling level, making the sequencing and allocation decisions, and on the other hand, there is a multi-stage dynamic optimization level, calculating optimal transitions and setpoints. The scheduling level comprises the master problem and the multi-stage dynamic optimization level is the primal problem. The approach differs from that by (Chatzidoukas et al., 2003b), and from OA and BD, mainly in that a decoupled approach is used. The information passed between the two levels is not a linearization of the primal problem (as in OA), or Lagrange multiplicators from the primal problem (as in BD). Instead, the information consists of higher-level key parameters that can be used in a scheduling framework in a more straightforward manner. The method is inspired by the general guidelines for solving MIDO problems given by Allgor and Barton (1999), where the formulations of the master and primal problems are not specifically given.This study is a sequel to a single-unit study by Nystr¨om, Franke, Harjunkoski, and Kroll (2005), the results of which are generalized to include multiple units running in parallel, thus including an allocation problem in addition to the sequencing problem. The presented method also supports further decomposition of the problem, particularly, the dynamic optimization can be split up to span shorter time ranges, or into separate optimizations for different units. In the present study this is illustrated by performing the dynamic optimizations separately for the separate units operated in parallel. 2. Statement of the production optimization problem Assuming we have M parallel reactors, units, or machines with indices m = 1, . . . , M, we introduce so-called processing stages s = 1, . . . , Sm , where Sm is the amount of stages defined for unit m. The processing stages correspond to two different production modes of continuous operation, transition and production, as will be explained in more detail below. We can write the process dynamics of the units as the dynamic equations d xs,m (τs ) = fm (xs,m (τs ), us,m (τs )) · Ts,m , dτs ys,m (τs ) = gm (xs,m (τs ), us,m (τs ))

(1)

Above, xs,m (τs ) is the process state within a stage s on unit m, us,m (τs ) is the input to unit m, and ys,m (τs ) is the process output from unit m.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

The variable τs is the “normalized integration time” for stage s, for which holds

P


ξs,m,p = 1,

395

if stage s on unit m is a production stage (8)

p=1

τs = 0

at the beginning of processing stage s,

τs = 1

at the end of processing stage s

Since each product p is produced at most once on each unit m, we have (2)

Further, the variable Ts,m is the duration of stage s on unit m. In this dynamic representation, the time derivatives within a stage s are thus scaled according to dt = Ts,m dτs

Sm


ξs,m,p ≤ 1,

Sm M



where t is the actual time. The stages are coupled to each other through stage coupling constraints

m=1 s=1

xs,m (0) = xini,m ,

s = 2, . . . , Sm

(4)

s=1

1. The production plan is cost optimal. 2. Within a processing stage s producing p, the quality variable with index q is always within the bounds Lp,q and Up,q . 3. At the end of the production campaign, the produced amount of p is at least Rp . 4. The general operational constraints are met. To treat the constraints from the production requests, i.e., the constraints on quality variables and amounts, we partition the output ys,m according to ys,m (τs ) = [ ys,m,gen (τs

)T

ys,m,q (τs

)T

T

ϕ˙ s,m (τs ) ]

(6)

where ys,m,gen are general process variables, ys,m,q are the quality variables, and ϕ˙ s,m is the rate of production of material on unit m. Further, we define the binary time-slot variable ξs,m,p , see e.g. (Pinto & Grossmann, 1998). If ξs,m,p = 1, it means that product p is produced within stage s on unit m, i.e., (s, m) is a production stage where on-spec material is produced. In a transition stage, on the other hand, unwanted off-spec material is produced, i.e. P
p=1

ξs,m,p = 0,

ξs,m,p ≥ 1,

if stage s on unit m is a transition stage

(7)

∀p

(10)

We can now make a linear mapping between the productdependent quality variable constraints, and stage-dependent path constraints on the quality variables

(5)

We now introduce a number of product orders p = 1, . . . , P. The products are characterized by upper and lower bounds on so-called “quality variables”; for polymers typically melt index, etc. Let the number of quality variables be Q. The lower bound of a quality variable with index q = 1, . . . , Q for product p is denoted Lp,q and the upper bound is denoted Up,q . In addition to the quality variable specifications, the product orders specify the required amount to be produced of each product, Rp , as well as the due date Dp of the product, i.e., the time by which the production of the product should be complete. Given the set of orders and the initial state of the process, the task is to calculate a production plan such that

(9)

and since each product p is produced at least in one stage on one unit

(3)

xs,m (0) = xs−1,m (1),

∀m, p

s=1

ys,m,q (τs ) ≤

P


 ξs,m,p Up,q + 1 −

p=1

P


 ξs,m,p  Umax,q ,

p=1

∀τs , s, m, q

ys,m,q (τs ) ≥

(11) P


 ξs,m,p Lp,q + 1 −

p=1

P


 ξs,m,p  Lmin,q ,

p=1

∀τs , s, m, q

(12)

Thus, only one of the upper or lower bounds on a quality variable is active at a time. If no product (off-spec) is produced, the default values Umax,q and Lmin,q become active. For the due date of a product p, the following (bilinear) constraint can be formulated   s
Tz,m ξs,m,p ≤ Dp , ∀s, m (13) z=1

A similar mapping can be made between the product-dependent production amount constraints and the produced amounts of a product. First, we have the amount produced within a stage s on unit m as the dynamic function d ϕs,m (τs ) = ϕ˙ s,m (τs ) · Ts,m , dτs

ϕs,m (τs = 0) = 0

(14)

Assuming we produce each product order p only in one stage s on one unit m, we can express the amount to be produced of m as the stage end-point constraint ϕs,m (τs = 1) ≥

P


ξs,m,p Rp ,

∀s, m

(15)

p=1

Alternatively, one wants to enable the optimizer to split up the production between units, yielding the (bilinear) constraints Sm M

m=1 s=1

ξs,m,p ϕs,m (τs = 1) ≥ Rp ,

∀p

(16)

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

396

In some cases (see Section 3.3) one wants to force the total make span of each unit m, i.e., the ending time of the last stage on m to have at least a given value, i.e. Sm


Ts,m ≥ Tmin ,

∀m

(17)

s=1

We usually also have general constraints related to operational safety, model feasibility, etc. These are usually universal path constraints and can be expressed as ys,m (τs ) ≤ Um ,

all τs

(18)

ys,m (τs ) ≥ Lm ,

all τs

(19)

For notational simplicity, any inputs and states which are to be constrained are here taken to be included in the vector of outputs ys,m . Introduce the cost contribution of stage s on unit m d Js,m (τs ) = hm (us,m (τs ), ys,m (τs ), ξs,m,p ), dτ Js,m (τs = 0) = 0

(20)

This is supposed to be a growing function of τs , typically an integral consisting of raw material consumption, off-spec material production penalty, etc. The complete production optimization problem can now be stated as min

us,m (τs ),Ts,m ,ξs,m,p

Sm M



Js,m (τs = 1)

(21)

m=1 s=1

Outer Approximation (OA): slack variables were needed to a high degree even without having due dates, and the optimal result was not always obtained. Nonconvexities are introduced e.g. by the bilinear constraints (13) and (16). • Using standard methods such as Benders Decomposition (BD) or OA have the limitation that in-depth information about the dynamic optimization problem (Lagrange multipliers and/or linearizations) is required. This makes the problem hard to implement and debug, and limits the choice of dynamic optimizers. • The size of the problem may grow formidable with an increasing number of machines and stages. As a remedy to the problems stated above, an alternative iterative procedure to solve the single-unit problem has been proposed by Nystr¨om et al. (2005). In the following, a similar method for the multiple-machine case will be presented. This includes a mixed-integer linear problem (MILP), which is the master problem, and a dynamic optimization (DO) problem, which is the primal problem. In each iteration, key parameters are extracted from the primal problem and passed to the master problem. Correspondingly, results from the master problem are passed to the primal problem. In Section 3, the MILP scheduling problem for parallel production units is presented. In Section 4, the dynamic optimization problem is described. The correspondence between the two is discussed in Section 5. How parameters for the scheduling problem are estimated from the dynamic optimization is discussed in Section 5.1, and how the parameters for the dynamic optimization are calculated from the scheduling problem in Section 5.2. The complete iteration scheme is discussed in Section 6, and in Section 7, a case study is presented.

subject to the constraints given by (1), (4), (5), (7), (8), (9), (10), (11), (12), (13), (14), (16), (17), (18), (19), (20). The time-scaling variables Ts,m are auxiliary decision variables to the optimization problem, used in addition to the process input us,m (τs ), and causing the optimizer to determine the length of each processing stage as a part of the total optimization problem. The binary sequencing (time-slot) variables ξs,m,p cause the sequencing to be a subproblem of the complete problem. The stated problem belongs to the MIDO class. A study by Prata et al. (2005) on solving this type of problem using OA is currently ongoing. It is difficult to solve the complete problem due to the following reasons:

The master problem formulation presented here is an extension of a previous study by Nystr¨om et al. (2005) to the multipleunit case. It is related to the standard formulations for solving travelling salesman problems, see e.g. (Winston, 1994), (Pekny & Miller, 1991), or (Pinto & Grossmann, 1998). In the formulation, the following assumptions are made:

• The problem has varying structure. A given product can be produced either on a multitude of units, or in one go on a single-unit, requiring a variable number of production stages. The production stages may either require a transition stage between them (if they do not have overlapping specifications) or not. The first stage in a campaign may either need to be a transition stage or not, depending on the initial state and the product to be produced first. This leads to situations where one or several time-scaling variables Ts,m may have the value zero. Each stage where this occurs makes the optimization ill conditioned, since the decision variables us,m (τs ) then have no impact and the problem gets too many degrees of freedom. • The problem is highly nonconvex. This has been shown by solving the single-unit problem in (Nystr¨om et al., 2005) by

• The production of a product can be distributed on many units, but all production of said product must be (partly) overlapping in time. As a consequence, the production of a given product cannot be split up so that it is produced twice on the same unit, it can only be split up across units. (This is an assumption in the present formulation only, not a fundamental one). • Due to the above assumption, the production of a product can be split up into at most M production instances, where M is the number of units. Such instances are referred to as “jobs”. If one should like to split up the production of a product into more jobs than M, the current formulation requires that dummy products are introduced. • At least one product must be produced on each unit during the planned campaign. (This restriction is not necessary; it

3. A MILP for scheduling multiple continuously operated units (master problem)

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

can be relaxed through slight reformulations of the problem or by the introduction of dummy products). • The units are identical, i.e., the transition time between two given products is the same on each unit, as is the production cost per time and the produced amount per time. (This restriction can also be removed according to need). • The units may have different initial conditions at the start of the planned campaign. Thus, the initial transition to a given product is likely to have different duration and cost depending on the unit m. • The make span is to be included in the cost function in a manner which involves actual production costs. This is done by forcing the productions on the individual units to all have identical make spans, thus the units in some cases produce more material than necessary. For the sake of brevity, all mathematical symbols are not explained as they appear. Instead, it is referred to the nomenclature in the beginning. Notice that the binary time-slot variables ξs,m,p , used in Section 2 for introducing the overall problem, are replaced by a binary representation more suitable for handling sequence-dependent changeover. 3.1. Optimized objective

P
M


β0,p,m γ0,p,m +

p=1 m=1

+

P
M


P
P
M


βp,b,m γp,b

P+1


βp,b,m = αp,m ,

∀p, m

(25)

b=1

If product b is produced on unit m, b comes after one and only one product p, or b is first P


βp,b,m = αb,m ,

∀b, m

(26)

p=0

Total sum constraint, enforcing that the number of ones in βp,b,m is correct P P+1



βp,b,m =

p=0 b=1

P


αp,m + 1,

∀m

(27)

p=1

Note: this constraint requires/enforces that there is at least one product being produced on each unit. Only units m which can/are allowed to produce p should do so αp,m = 0,

∀m ∈ Mp ,

∀p

(28)

(29)

Relation between starting time of product b and completion time of product p on unit m

p=1 b=1 m=1

dp,m γ˙ p

If product p is produced on unit m, p comes before one and only one product b, or p is last

Auxiliary constraints to make the problem tighter  αp,m ∀p, b, m|p = b βp,b,m ≤ αb,m ∀p, b, m|p = b

The following optimization problem is to be solved:

min

397

(22)

p=1 m=1

where the first term is the cost of the transitions from initial state, the second term is the cost of the transitions between products, and the third term is the cost of the production stages. The transition (changeover) costs referred to above involve costs incurred by producing unwanted transition (off-spec) material during the transition from one product to another, as well as from the initial state. The cost (22) does not contain directly time-related costs such as penalty for lateness and earliness. An option is to introduce such costs instead of hard constraints such as the due date constraints given below. The cost (22) also does not directly contain a term penalizing the total make span, sum of completion times or such. Instead, the total time is indirectly optimized through optimization of the cost of production.

tb,m ≥ cp,m + δp,b − M1 (1 − βp,b,m ),

∀p, b, m|p = b (30)

where M1 is an estimated upper bound on the time horizon of the problem, obtained from the parameters δp,b , δ0,p,m , ϕ˙ p , and Rp . Relation between starting time of product p and transition time of product p from initial state on unit m tp,m ≥ δ0,p β0,p,m ,

∀p, m

(31)

Relation between starting time and completion time of product p on unit m cp,m = tp,m + dp,m ,

∀p, m

(32)

Relation between duration (continuous-valued decision variable) of product p on unit m and the binary decision to produce p on unit m

3.2. MILP constraints

dp,m ≤ M2,p αp,m ,

The optimization problem (22) is to be solved subject to the following constraints: Assignment/elimination of meaningless variables

where M2,p is an upper bound on the job duration for product p. Minimum duration of product p on unit m

β0,P+1,m = 0, βp,p,m = 0,

∀m

(23)

∀p, m

(24)

dp,m ≥ dmin αp,m ,

∀p, m

∀p, m

(33)

(34)

where dmin is a chosen minimum value for the time a job is allowed to take.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

398

The total amount of product p produced should at least equal the required amount of product p M


dp,m ϕ˙ p ≥ Rp ,

∀p

(35)

m=1

Due date constraints cp,m ≤ Dp ,

∀p, m|∃Dp

(36)

Requirement that all production of product p is overlapping in time across different units (optional) tp,m ≤ cp,n + M1 (2 − αp,m − αp,n ),

∀p, m, n|n = m

(37)

The make span is greater than or equal to all completion times T ≥ cp,m ,

∀p, m

(38)

Constraint forcing the completion time on each unit m to be equal to the make span T and thus equal to the completion time on all other units (see Section 3.4)

T =

P


β0,p,m δ0,p,m +

p=1

P
P


βp,b,m δp,b +

p=1 b=1

∀m

P


dp,m ,

p=1

(39)

Additional constraints can be added to the problem as required. For instance, if it is known a priori that a transition from product p to product b is undesirable or impossible, it can be prohibited by assigning the corresponding binary variable to equal zero βp,b,m = 0,

∀m,

∀(p, b) ∈ I

(40)

Such assignments generally lead to reduced degrees of freedom in the binary decisions and thus make the MILP easier to solve. Also, they will lead to a simplification of the complete production optimization problem in ensuring that unnecessary optimization tasks are not carried out in the primal problem. 3.3. Fixing the make span to enforce idling During periods of a low amount of orders, the optimization task of interest may not be to maximize the throughput, but to minimize the production cost with a given make span instead. This may be the case if the throughput cannot be reduced to the level required due to operational constraints, and a shut-down of the plant is infeasible. Such an optimization can be achieved by requiring that the make span have at least a given value T ≥ Tmin

(41)

The solution of the MILP problem then strives to calculate a production such that idling/overproduction occurs where it is most economical. The problem could also be expanded to determine whether it is profitable to shut down production lines by introducing dummy products and shutdown/startup costs.

3.4. Remarks Constraint (39), together with, (35), forces the make spans of all production units to be equal. Note that with a given sequencing, i.e., with βp,b,m given, the durations dp,m are forced to match the make span, since dp,m are the only free variables in the Eq. (39). This in turn forces the start times and end times to exactly match the transition times tb,m = cp,m + δp,b ,

if

βp,b,m = 1

(42)

In other words, the solution cannot be one with empty spaces between jobs that do not match up with the given transition times; the only option is to force production units to idle (to produce more product than required). This causes actual production cost to be reflected in the make span, if the make span is not constrained according to (41). Due to the structure of the problem, where the binary variables βp,b,m are dependent on the binary variables αp,m , it is beneficial to define the branching priority such that the variables αp,m are branched upon first. Note that some of the binary variables can be set to continuous in order to reduce the total number of binary variables. These are the ones which, due to constraints, have no alternative but to assume values of zero or one. In some cases, this can make the problem solve somewhat faster. One possibility is to relax the variables αp,m (to continuous between 0 and 1), which will always take the value zero or one due to (25) and (26), if βp,b,m are binary variables. Then, however, the advantage with the branching priority mentioned above is lost. Another alternative is to relax e.g. βp,P+1,m and β0,p,m for all p and all m. 4. Multi-stage DO for optimal production (primal problem) The MILP master problem solves the production optimization problem, including sequencing, allocation, and durations, based on the customer orders. This takes place on a high level, where the key product- and process-related parameters γp,b , γ0,p,m , γ˙ p , δp,b , δ0,p,m , and ϕ˙ p are assumed given. On a lower level, one needs to address how to operate the process, i.e. what process inputs to have in order to fulfil this scheme in an optimal fashion. This part of the production optimization problem is the primal problem. The primal problem is obtained by fixing the binary variables ξs,m,p of the overall MIDO problem in Section 2. This yields the multi-stage DO problem min

us,m (τs ),Ts,m

Sm M



Js,m (τs = 1)

(43)

m=1 s=1

subject to the dynamic multi-stage process representation (1), stage coupling constraints (4) and (5), quality variable path constraints (11) and (12), additional dynamic equations (14) for the product amount produced within a stage, general operational constraints as in (18) and (19), and dynamic cost equations (20). The primal problem consists of the determination of optimal transitions and optimal production set points as subproblems.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

The decision variables are the process inputs us,m (τs ) that have an impact on process operation, product quality, throughput, etc. The durations of the production stages, Ts,m , are also decision variables. Typically, the optimized objective in (43) is a function consisting of cost components for raw material, energy, and negative profit stemming from off-spec material production. The primal problem is highly nonlinear. This is caused by the scaling of the dynamic equations by the stage durations Ts,m . Additionally, nonlinearity is typically contained within the process model. Ideally, the primal problem would be treated as a complete dynamic optimization over all units and the whole time span. However, with an increased amount of reactors, increased dynamic model complexity, and increased planning time horizon, the size of the primal problem may easily grow to be overwhelming. In the present study, the primal problem is split up into a separate DO subproblem for each parallel unit. In order to enable treating the DO subproblems as separate problems, the auxiliary continuous-valued parameters rp,m are used to set up the DO problems from the master problem solution. These parameters give the minimum amounts to be produced of each product p on each unit m in the primal problem ϕs,m (τs = 1) ≥

P


ξs,m,p rp,m ,

∀m

(44)

p=1

The reason for introducing rp,m is to be able to split up the primal problem into unit-wise DO subproblems. How rp,m are obtained from the MILP is explained in section 5.2. The constraints (44) replace the constraints (15) or (16) posed on the amount of produced product in the total problem. Another set of additional constraints used in the primal problem in this study require production stages to have a minimum duration dp,m Ts,m ≥

P


ξs,m,p dp,m ,

∀m

(45)

p=1

These constraints replace the time-related constraint (17) in the total problem. This issue is also discussed in more detail in Section 5.2. Remark: Including time constraints such as the due date constraints (13)/(36), the overlapping constraints (37), or the make span constraints (17)/(38) in the DO problems will lead to DO problems, which may be unsolvable due to infeasibility. The approach in this article is to treat such constraints only in the master problem, and to check the primal problem for feasibility a posteriori. A solution, where the due date is violated, is discarded, i.e., it does not contribute to the upper bound. 4.1. DO preprocessor The stage setup in the DO problem depends on the result of the preceding master problem, i.e., on αp,m (allocation), and βp,b,m (sequencing). The work of setting up the DO problems from the solution of the master problem is done by a DO preprocessor. From the binary variables αp,m and βp,b,m , the preprocessor de-

399

termines the number of stages for each unit m, whether a stage s is a transition stage or a production stage, and which product is to be produced within each production stage. From this information, it in turn generates constraints for the DO problems accordingly. The DO preprocessor should exclude unnecessary processing stages, i.e., those whose length will be zero, since they make the DO problem highly ill-conditioned (see Section 2). The possibility of thus modifying the structure of the DO problem according to the result of the master problem, i.e., the possibility of decoupling the DO problem from the master problem, is a highly useful property differentiating the present method from the standard methods BD and OA. 5. Correspondence between primal and master problem In each iteration, the results of the DO subproblems for the separate units m, m = 1, . . . , M, are combined to form a total cost of the production. This represents the objective of the primal problem in that iteration. If the primal problem is feasible and the objective is smaller than for all the feasible primal problems so far, the objective will be used as new value for the upper bound (UB). The primal and master problems minimize the same objectives, although subject to different constraints and using different levels of process knowledge. In the primal problem, operational and safety constraints, etc. are handled, see above. In the master problem, high-level constraints are handled, such as which unit can be used to produce which product. Additionally, the constraint (39) forces the time horizon of all units to be the same. Note that the time constraints such as the due date constraints (36), the overlapping constraints (37), or the make span constraints (38) need not be forced in the DO subproblems of the primal problem. Instead, they will be approximately satisfied through the constraints (45). Thus, primal problems which are insolvable due to infeasibility can be avoided. An important distinction between the master and primal problem is that the primal problem is a discrete-time problem (with scaled time representation) and the master problem is a continuous-time problem. Thus, constraints posed on inventory levels, etc., cannot readily be addressed in the master problem. Instead, such constraints can be accounted for in the primal problem. The different properties of the continuous-time and the discrete-time formulation, and hybrid versions of the both has been a topic of previous studies, see (Schulz, Engell, & Rudolf, 1998) and references therein. 5.1. MILP parameter estimation using results from DO Solving the DO problems yields new information on the MILP parameters γp,b , γ0,p,m , γ˙ p , δp,b , δ0,p,m , and ϕ˙ p . The approach for updating the parameter estimates is as follows: 1. If there is no other estimate of an element, a default value is used; an estimated upper bound for ϕ˙ p , and a small one, typically zero, for the other parameters.

400

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

2. If a cost element in a parameter matrix is obtained for the first time, the value obtained is used as an estimate: • A first estimate for the cost contribution parameter γp,b is obtained as Js,m (τs = 1) in a transition stage s from p to b, a first estimate for γ0,p,m is obtained as Js,m (τs = 1) in an initial transition to p on unit m, and a first estimate for γ˙ p is obtained as Js,m (τs = 1)/Ts,m in a production stage s producing p. • A first estimate for the time parameter δp,b is obtained as Ts,m in a transition stage s from p to b, and a first estimate for δ0,p,m is obtained as Ts,m in an initial transition to p on unit m. • A first estimate for the production rate parameter ϕ˙ p is obtained as ϕs,m (τs = 1)/Ts,m in a production stage s producing p. 3. If a first estimate for a parameter has already been obtained previously, the estimate is updated using a suitable rule. The aim is that the master problem provides a reliable lower bound, which is not too loose. The aim is also that the master problem produces relevant solutions, which lead to a decrease of the upper bound. The “suitable rule” mentioned above thus should provide an estimate, which is as accurate as possible, otherwise the upper bound will suffer. However, it is also of interest that it produces an underestimate, in order for the master problem to provide a lower bound. The lower bound is clearly valid when most of the estimates are conservative default values, i.e., when a lot of DO information has not yet been obtained. The remaining question is what to do, when several values are available for forming an estimate. The possibilities include 1. Use an underestimate consisting of a minimum of all values obtained. This approach was used in the single-unit study by Nystr¨om et al. (2005), and provides an underestimate, which is quite rigorous, but may be too loose. 2. Use some average value of the values obtained. This is a trade-off between having realistic values and underestimates. 3. Use a minimum of the last few values obtained. This approach maintains an underestimate, which is not conservative and allows the quality of the estimate to improve as the iterations proceed. 4. Use only the most recent value obtained and discard old values. This is an extreme form of the previous method and was used in the example of this paper, since it turned out to give a good result. Selecting the estimation rule above provides a possibility of tuning the problem. For instance, one can start with 1. If the lower bound does not reach the upper bound after a comfortable amount of iterations, 2. is tried, etc. In the end, the iterations should look something like the ones in Figs. 3 and 6. Having to use 3. or 4. implies that the parameters are highly contextdependent, i.e., their values depend on the overall sequence and allocation. This may be due to varying production rate, or to the possibility of highly transient (short) production stages.

The issue of which transitions are tried out in dynamic optimizations and which are not, is determined by the MILP. If a transition is not used in the solution of the MILP, even if its cost estimate consists of a conservative default value (zero), it means that that transition is not feasible or interesting for the problem in question. Subsequently, this transition does not appear in the primal problems and unnecessary optimization effort by the dynamic optimizer can be saved. 5.2. Passing information from the master problem to the primal problem A solution from the master problem is characterized by the binary variables αp,m (allocation), and βp,b,m (sequencing). If the whole primal problem is solvable as a single problem including all units, these variables are enough to characterize the DO problem. Constraints corresponding to (39) would then be introduced into the primary problem as well, forcing all units to have the same make span, and the total production of all units to correspond to the required amount. However, if the primal problem is to be solved as a separate DO subproblem for each unit m, this approach cannot be used, since it involves constraints, which should be valid for all the units as a whole. A workaround to this problem is to introduce the auxiliary continuous-valued variables rp,m . These variables give the minimum amount of material p to be produced on each unit m, and are used in the stage end-point constraints (44). In order to give the DO subproblems enough degrees of freedom, the results from the master problem are interpreted as ratios of the total amount and rp,m are calculated according to ϕ˙ p dp,m rp,m = M Rp ˙ p dp,m m=1 ϕ

(46)

This guarantees that at least the required total amount according to the orders is produced within the primal problem, but still gives the DO subproblems flexibility in manipulating the throughput. In addition to the variables rp,m , the duration variables dp,m from the master problem are used in the decomposed primal problem. These variables are used in the constraints (45) to force the production of product p on unit m to take a certaint amount of time. This causes the time constraints (due date constraints

Fig. 1. Diagram of information exchange between master problem and primal problem.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

401

and production overlapping constraints) to be (approximately) fulfilled in the primal solution, although it is composed of DO subproblems. The information exchange between the master problem and primal problem is summarized in Fig. 1. 6. Iterative approach to solving the production optimization problem The complete production optimization problem can be solved in an iterative approach as follows. See also Fig. 2. 1. An initial sequencing and allocation is first obtained, for instance by solving the master problem of Section 3 with initial values of the parameters γp,b , γ0,p,m , γ˙ p , δp,b , δ0,p,m , and ϕ˙ p . 2. The values of the auxiliary variables rp,m and dp,m are obtained from the solution of the master problem according to Section 5.2, and passed to the primal problem together with the binary sequencing and allocation variables αp,m and βp,b,m . 3. The DO subproblems of the primal problem described in Section 4 are set up by the DO preprocessor and solved, and the solutions are combined to form the resulting objective of the primal problem. The upper bound (UB) is set equal to the smallest value of all feasible primal objectives so far. 4. The parameters γp,b , γ0,p,m , γ˙ p , δp,b , δ0,p,m , and ϕ˙ p are estimated from the solution of the primal problem as described in Section 5.1. 5. The master problem described in Section 3 is set up and solved. Integer cuts can be used here to eliminate previously tested binaries (see remark below). The lower bound (LB) is set equal to the solution of the last master problem. 6. The termination criterion is checked. The iterations can be terminated when the lower bound is sufficiently close to the upper bound LB ≥ UB −


(47)

where
gives the tolerance. If (47) is fulfilled, a termination flag is set, which terminates the procedure at a suitable point in the next iteration. 7. Back to step 2. Remark: Since the master problem solution is characterized not only by the binary sequencing/allocation variables βp,b,m and αp,m , but also by the continuous-valued variables rp,m and dp,m , the concept of using integer cuts in the master to eliminate sequences/allocations already tested is questionable. The solution is not completely described by binary variables, and it is possible that a sequence/allocation which has been previously tested will yield a better result with a new set of rp,m and dp,m . If the same sequencing/allocation is obtained in two successive iterations, however, it is likely that the primal problem has not changed much, and that the same sequencing/allocation would just be tried out indefinitely. For this reason, the approach used in this study is to check if the same sequencing/allocation is obtained from the MILP as in the last iteration. If this is the case,

Fig. 2. Schematic diagram of iterative production optimization.

this sequencing/allocation is added to a list of binary cuts which eliminates it from the master problem. A new master problem is then solved to yield a new result. 7. Case study: three parallel polymerization lines The method presented so far will now be demonstrated on a case where three parallel production units are used to produce polymer. The lines are considered to be decoupled, i.e., the raw material supply, actual production and finishing occur independently for all three lines. The scenario is that nine different orders of nine products, respectively, are to be produced. The product quality is characterized by two different product quality variables, QV1 and QV2. The order/product specifications can be found in Table 1. The amount of product to be produced varies according to the product p. One product (number 5) has a due date, before which its production should be completed. The product quality variables are controlled mainly by two manipulated variables, MV1 and MV2, which represent flow rates of additives to the process. The throughput of the process is mainly controlled by a third manipulated variable, MV3, which represents the raw material feed rate to the process. All manipulated variables are constrained with respect to the minimum/maximum values they can attain. The minimum job length was set to dmin = 4 h. The value of the termination tolerance
in (47) was set to zero, i.e., the lower bound of the objective function is required to cross the upper bound for the iterations to terminate. The problem is the same as in the single-unit study by Nystr¨om et al. (2005), except that the make span is not included in the cost function. Instead, the cost function (20) contains only a raw-material cost component and an offspec material penalty cost component. Further, the throughput is manipulated in this study through an additional input variable. The DO subproblems of the primal problem were solved using the publicly available solver HQP, see e.g. (Franke & Arnold, 1997). HQP is a multiple-stage (multiple-shooting) solver, whose S-function interface (MathWorks, 2002) supports features such as passive and active inputs, decimation of decision variables, and multiple samples within optimization intervals. The MILP master problem was solved using the commercial solver ILOG CPLEX (ILOG Inc., 2002).

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

402 Table 1 Order specifications for products 1, . . . , 9 Product p

Lower bound QV1 (Lp,1 )

Upper bound QV1 (Up,1 )

Lower bound QV2 (Lp,2 )

Upper bound QV2 (Up,2 )

Required amount (Rp ) [kg]

Due date (Dp ) [h]

Production possible on units (Mp )

1 2 3 4 5 6 7 8 9

1.0 6.0 11.0 6.0 11.0 6.0 1.0 11.0 1.0

3.0 8.0 13.0 8.0 13.0 8.0 3.0 13.0 3.0

0.5 1.5 2.5 0.5 1.5 2.5 1.5 0.5 2.5

0.7 1.7 2.7 0.7 1.7 2.7 1.7 0.7 2.7

60000 60000 60000 60000 100000 60000 60000 60000 200000

– – – – 10.0 – – – –

{1} {2} {3} {1, 2} {2, 3} {1, 3} {1, 2, 3} {1, 2, 3} {1, 2, 3}

The production optimization problem was solved for the two cases of operation at maximum capacity and operation at underfull capacity described in Section 5.2. In the present study, the parameter estimates for the primal problem were updated by the method of picking the last value obtained (see Section 5.1), and this approach gave a very good result. 7.1. Operation in full order scenario In the full order scenario, it is of interest to maximize throughput and minimize make span. Thus, in this case the constraint (41) is excluded from the master problem, and the production occurs as time-efficiently as possible. Fig. 3 depicts how the iterative procedure evolves for the example problem. The lower curve is the lower bound of the objective function obtained from the master problem solutions. The initializing MILP solution (with objective zero) is not shown in the figure. The upper curve depicts the upper bound of the objective obtained from the primal problems. A primal solution, which is infeasible (due to violated due dates) is shown with a marker differing from that of a feasible solution. An infeasible solu-

Fig. 3. Illustration of how iterative solution proceeds for the example production optimization problem.

tion does not contribute to the upper bound, nor does a solution, which has a larger value than the upper bound attained so far. The evolution of the iterative procedure can roughly be divided into three parts as is visible in the figure. In the first part, there is a large amount of zero estimates for the transition costs and times, thus the lower bound increases slowly. In the same part, the master problem has no real information about the parameter values and produces primal problems with arbitrarily bad values of the objective function. In the middle part, the zero estimates are being replaced by realistic estimates, and the lower bound rapidly increases. The values of the objective functions of the primal problems also get better and better on average. In the last part, almost all estimates used for the solution are realistic ones, and the lower bound increases mainly due to the integer cuts, which eliminate already tested sequences. The primal problem solutions now closely correspond to the master problem solutions. The best objective function of the primal problem (86.4) is obtained in the 9th iteration. The lower bound crosses the upper bound in the 11th iteration, and the iterative procedure is subsequently terminated. In Fig. 4 are shown Gantt charts depicting each master solution (left column) and Gantt charts depicting the corresponding primal solutions (right column). The charts are ordered from top to bottom, so that the uppermost chart to the left corresponds to the initializing master solution. In each chart is also shown the obtained value of the objective function, cf. Fig. 3. The time range of the Gantt charts is 0–40 h. The following issues are demonstrated in Fig. 4: • In the beginning of the iterative procedure (top of figure), there is no information on transition times and costs, thus the master solutions (to the left) have zero-length transitions. As more information on the transitions is obtained (top to bottom), the master solutions are forced to include actual transition times and costs. • The master solutions (to the left) always have equal make spans for all three units. The production is split up across units or overproduction is enforced to accomplish this. It turns out that the overproduction is negligible; splitting up the production over units is enough to cause the units to have equal makespan.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

403

Fig. 4. Gantt charts of the production plans obtained in the iterative procedure, with start of iteration at top. The left column contains charts obtained from the master problems (MILPs) and the right column charts obtained from the primal problems (DOs). In the top middle of each chart is shown the obtained objective function.

• The primal solutions (to the right) do not always have equal make spans for all three units. This is because the decoupling into DO subproblems does not allow constraints over units. However, as more information on the process is gathered (top to bottom), the three make spans from the subproblems better and better match each other and that of the master solution. • The best solution is obtained in the ninth primal solution (with objective function 86.4). The same sequencing/allocation is tested again in iteration 11, but with a slightly worse result. • The production of product 5 must be split up between two units because of the tight due date. The reason that some primals are infeasible (cf. Fig. 3) is that the due date of product 5 is violated on one of these units. • The initial state of the system is such that product 1 can be produced on the first unit without an initial transition. Fig. 5 shows an example result from a dynamic optimization. This is the first DO subproblem of the first primal problem of the iterative procedure, where the production sequence 6-8-4-1 is optimized for unit 1 (see Fig. 4, first row, right column). The figure contains the two quality variables (QV1 and QV2, to the left), and the two manipulated variables having the largest direct effect on these (MV1 and MV2, to the right). The vertical dashed lines in the QV plots indicate boundaries of stages. The narrow stages are transition stages and the broad ones are production stages. Within the production stages, the quality variables are constrained to be within a specific band, cf. Table 1. Within the transition stages, the quality variables may obtain values outside of the bands. The bands are indicated by

horizontal dashed lines in the QV plots. The horizontal dashed lines in the MV plots show the minimum/maximum values of the manipulated variables. 7.2. Operation in underfull order scenario In the underfull order scenario, there are not enough orders to allow production at full throughput. Thus, the throughput

Fig. 5. Quality variables and manipulated variables directly related to the product quality in an example sequence (6-8-4-1) for unit 1.

404

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

An illustration of how the iterations proceed for an example case, where the minimum makespan has been set to Tmin = 40 can be found in Fig. 6. The Gantt charts depicting the different solutions can be found in Fig. 7. The following issues are demonstrated in Fig. 7. • The master problem now has a fixed make span. Within this make span, it strives to plan the production such that the overproduction is as economical as possible. From the solutions it is evident that product 9 is the cheapest one to overproduce. • The lowest value of the objective function, 94.0, is found in the 10th iteration. However, this solution is infeasible (see Fig. 6). The best feasible solution, 95.6, is found in iteration 11, and the lower bound crosses the upper bound directly after that. Fig. 6. Illustration of how iterative solution proceeds for the example production optimization problem in underfull order scenario.

has to be reduced and/or overproduction has to be permitted. This is achieved by including the constraint (41) in the master problem, for the other part the production task is the same as in subsection 7.1. The constraint (41) forces the make span to be longer than it needs to be when operating at full capacity. The optimization then strives to decrease the throughput as much as the constraints allow, and introduces overproduction where this is most economical.

7.3. Discussion The case study presented shows that solving a complete production optimization involving both scheduling (sequencing and allocation) and dynamic optimization (operation) is feasible. For the problem presented here, a DO subproblem typically takes a couple of minutes to solve on a standard PC, and the MILPs seconds to minutes. With an increasing amount of products, the DO problem algebra size will grow linearly, whereas the amount of binaries in the MILP will grow quadratically. Thus, with increasing problem size, the master problem will at some point

Fig. 7. Gantt charts of the production plans obtained in the iterative procedure for optimization in underfull order scenario, with start of iteration at top. The left column contains charts obtained from the master problems (MILPs) and the right column charts obtained from the primal problems (DOs). In the top middle of each chart is shown the obtained objective function.

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

405

be the limiting factor, and alternative formulations or heuristics might be required to obtain solutions for it. The state of the art in scheduling is to use prespecified values or values stored in a database for the key characteristics of the process, whereas dynamic optimization is used to extract this information in the approach presented. Assuming that the information could be available in a database, the presented method will naturally be much slower than a conventional scheduling solution, since it is then recalculating information already available. In a real-world application, the method would probably have to be a hybrid one, where a database is cross-checked for the required information, and dynamic optimization used where the required information is absent. The scheduling problem solver might also need to address rescheduling issues, such as locking jobs which have already started or are about to start, generating new schedules which are similar to previous ones, etc. Tests with Outer Approximation for solving a single-unit version of the production optimization problem indicate that the problem is heavily nonconvex. This may be largely due to the fact that a discrete-time formulation was used in the Outer Approximation solution, whereas a continuous-time MILP formulation such as the one used in the present study is better suited for expressing sequence-dependent changeovers. It is not clear how this binary representation could be incorporated in an Outer Approximation approach. This demonstrates the benefits of the decoupled approach presented here—one is able to separate the complete problem into a discrete-time problem and a continuous-time problem (see Section 5), to formulate the scheduling problem in the way which is best suited for the problem at hand, as well as to include heuristics to speed up the scheduling problem.

conditions, such as changes in customer orders, priority changes, operation policies, or equipment failures. The scope of this study is limited to considering optimization performed off line, although the dynamic optimization part bears a resemblance to real-time optimization (RTO) or, particularly, to dynamic real-time optimization (DRTO), see e.g. the study by Tosukhowong, Lee, Lee, and Lu (2004), and references therein. On-line control, e.g. nonlinear model predictive control (NMPC, NLMPC), studied e.g. by BenAmor, Doyle, and McFarlane (2004) is not considered in more detail. However, the transitions obtained from the dynamic optimization could be used as reference trajectories in an (N)MPC application. Furthermore, supervisory systems, such as the one described by Srinivasan, Viswanathan, Vedam, and Nochur (2005), could be used for switching between different control applications for the transition stages and the steady-state stages. The approach is demonstrated on a case with three polymerization reactors operated in parallel, for which a short campaign of producing nine products is calculated. The illustrated case is of relatively small scale and of short time horizon. In an industrial application, the aim might be to solve problems involving tens to hundreds of orders over time horizons of weeks or months. The approach can be adapted to large problems by resorting to practical measures, such as increasing the use of stored (a priori) information, splitting up the dynamic optimizations into even smaller subproblems, and heuristics for making the MILP scheduling problem tractable. From this viewpoint, the ability to decouple the problem and use state-of-the art software for the dynamic optimization and scheduling subproblems is an advantage.

8. Conclusions and outlook

Acknowledgements

A method for production optimization for continuously operated processes has been presented. The method is applicable to cases where a multitude of similar units are operated in parallel. The whole production optimization problem is decoupled into a master problem, representing the scheduling problem, and a primal problem, representing the lower-level dynamic optimization for determining optimal operating points and transition trajectories. The master problem belongs to the class of MILP. The primal problem consists of a number of nonlinear DO subproblems. A solution to the overall problem is obtained through iteration between the master problem and the primal problem, as in the conventional methods for solving MINLP problems. The presented method utilizes the knowledge of the structure of the problem to provide a compact master problem MILP using estimated key parameters from the primal problem. The binary variables from the master problem do not need to enter the DO subproblems of the primal problem. Instead, the appropriate constraints are programmed by the DO preprocessor. The method illustrates a possible path for future planning applications, where the interaction between advanced planning and scheduling and dynamic optimization grows stronger, evolving into dynamic model-centric planning and scheduling. This type of approach supports swift reaction to changes in the production

The European Union is gratefully acknowledged for the support of R. Nystr¨om under the Marie Curie Fellowship, contract no. HPMI-CT-2001-00138. We thank Professor Dr. Marquardt and Dr. Jan Oldenburg at LFPT Aachen for valuable discussions, Dr. R¨udiger Franke for his support with the dynamic optimization and HQP, Mr. Adrian Prata at ABB Corporate Research for his work on the process model and his alternative multistage formulation in another environment, Dr. Alexander Frick, and Dr. Pousga Kabor´e at ABB Corporate Research for general discussion and help with the manuscript. Finally, we thank the two anonymous reviewers, whose feedback greatly helped to improve the manuscript. References Allgor, R. J., & Barton, P. I. (1999). Mixed-integer dynamic optimization. I. Problem formulation. Computers and Chemical Engineering, 23 (4/5), 567– 584. Backx, T., Bosgra, O., & Marquardt, W. (1998). Towards intentional dynamics in supply chain conscious process operations. Proceedings of FOCAPO Snowbird, Utah. Bansal, V., Sakizlis, V., Ross, R., Perkins, J. D., & Pistikopoulos, E. N. (2003). New algorithms for mixed-integer dynamic optimization. Computers and Chemical Engineering, 27 (5), 647–668.

406

R.H. Nystr¨om et al. / Computers and Chemical Engineering 30 (2006) 392–406

Bemporad, A., & Morari, M. (1999). Control of systems integrating logic, dynamics and constraints. Automatica, 35, 407–427. BenAmor, S., Doyle III, F. J., & McFarlane, R. (2004). Polymer grade transition control using advanced real-time optimization software. Journal of Process Control, 14, 349–364. Bhatia, T., & Biegler, L. T. (1996). Dynamic optimization in the design and scheduling of multiproduct batch plants. Industrial and Engineering Chemistry Research, 35, 2234–2246. Brendel, M., Oldenburg, J., Schlegel, M., & Stockmann, K. (2003). DyOS User’s Guide Version 2.1. LFPT, Prof. Dr. -Ing. Wolfgang Marquardt: RWTHAachen University. Cervantes, A. M., Tonelli, S., Brandolin, A., Bandoni, J. A., & Biegler, L. T. (2002). Large-scale dynamic optimization for grade transitions in a low density polyethylene plant. Computers and Chemical Engineering, 26, 227– 237. Chatzidoukas, C., Perkins, J. D., Pistikopoulos, E. N., & Kiparissides, C. (2003a). Optimal grade transition and selection of closed-loop controllers in a gas-phase olefin polymerization fluidized bed reactor. Chemical Engineering Science, 58, 3643–3658. Chatzidoukas, C., Kiparissides, C., Perkins, J. D., & Pistikopoulos, E. N. (2003b). Optimal grade transition campaign scheduling in a gas-phase polyolefin FBR using mixed-integer dynamic optimization. In B. Chen, & A. W. Westerberg (Eds.), Process Systems Engineering (pp.744–747). Elsevier. Emet, S., & Westerlund, T. (2004). Comparisons of solving a chromatographic separation problem using MINLP methods. Computers and Chemical Engineering, 28, 673–682. Feather, D., Harrell, D., Lieberman, R., & Doyle III, F. J. (2004). Hybrid approach to polymer grade transition control. AIChE Journal, 50 (10), 2502– 2513. Floudas, C. A. (1995). Nonlinear and Mixed-Integer Optimization. New York: Oxford University Press. Franke, R., & Arnold, E. (1997). Applying new numerical algorithms to the solution of discrete-time optimal control problems. In K. Warwick, & M. Karny (Eds.), Computer-Intensive Methods in Control and Signal Processing: The Curse of Dimensionality (pp.105–118). Basel: Birkh¨auser Verlag See also http://hqp.sourceforge.net/for download of HQP. Gallestey, E., Stothert, A., Castagnoli, D., Ferrari-Trecate, G., & Morari, M. (2003). Using model predictive control and hybrid systems for optimal scheduling of industrial processes. Automatisierungstechnik, 51 (6), 285– 293. Grossmann, I. E. (2002). Review of nonlinear mixed-integer and disjunctive programming techniques. Optimization and Engineering, 3, 227–252. ILOG, Inc. (2002). ILOG CPLEX 8. 0 User’s Manual. Mahadevan, R., Doyle III, F. J., & Allcock, A. C. (2002). Control-relevant scheduling of polymer grade transitions. AIChE Journal, 48 (8), 1754–1764. The MathWorks, Inc., (2002). Simulink® Model-Based and System-Based Design: Writing S-Functions. McAuley, K. B., & MacGregor, J. F. (1992). Optimal grade transitions in a gas phase polyethylene reactor. AIChE Journal, 38 (10), 1564–1576.

Mohideen, M. J., Perkins, J. D., & Pistikopoulos, E. N. (1996). Optimal design of dynamic systems under uncertainty. AIChE Journal, 42 (8), 2251– 2272. Nystr¨om, R. H., Franke, R., Harjunkoski, I., & Kroll, A. (2005). Production campaign planning including grade transition sequencing and dynamic optimization. Computers and Chemical Engineering, 29 (10), 2163–2179. Oldenburg, J., Marquardt, W., Heinz, D., & Leineweber, D. B. (2003). Mixedlogic dynamic optimization applied to batch distillation process design. AIChE journal, 11, 2900–2917. Pekny, J. F., & Miller, D. L. (1991). Exact solution of the no-wait flowshop scheduling problem with a comparison to heuristic methods. Computers and Chemical Engineering, 15 (11), 741–748. Pinto, J. M., & Grossmann, I. E. (1998). Assignment and sequencing models for the scheduling of process systems. Annals of Operations Research, 81, 433–466. Prata, A., Oldenburg, J., Marquardt, W., Nystr¨om, R., & Kroll, A. (2005). Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor, submitted for publication. Schulz, C., Engell, S., & Rudolf, R. (1998). Scheduling of a multi-product polymer batch plant. In: J. F. Pekny, G. E. Blau, B. Carnahan (Eds.),Proceedings of FOCAPO Snowbird, Utah (pp. 75–90). Michigan: CACHE Publications. Sinclair, K. B. (1987). Grade change flexibility—defined, determined, compared. In Proceedings of Fifth International Conference on Polyolefins (pp. 2–19). Houston, Texas. Smania, P., & Pinto, J. M. (2003). Mixed-integer nonlinear programming techniques for the short term scheduling of oil refineries. In B. Chen, & A. W. Westerberg (Eds.), Process Systems Engineering (pp.744–747). Elsevier. Srinivasan, R., Viswanathan, P., Vedam, H., & Nochur, A. (2005). A framework for managing transitions in chemical plants. Computers and Chemical Engineering, 29, 305–322. Tosukhowong, T., Lee, J. M., Lee, J. H., & Lu, J. (2004). An introduction to a dynamic plant-wide optimization strategy for an integrated plant. Computers and Chemical Engineering, 29, 199–208. Tousain, R. L., (2002). Dynamic Optimization in Business-wide Process Control. Ph.D. Thesis, Technical University Delft. Van Brempt, W., Van Overschee, P., Backx, T., Ludlage, J., Hayot, P., Oostvogels, L., & Rahman, S. (2003). Grade-change control using INCA model predictive controller: application on a DOW polystyrene process model. Proceedings of the American Control Conference (pp. 5411–5416). Denver, Colorado. Vassiliadis, V. S., Sargent, R. W. H., & Pantelides, C. C. (1994). Solution of a class of multistage dynamic optimization problems 1. Problems without path constraints. Industrial and Engineering Chemistry Research, 33, 2111– 2122. Viswanathan, J., & Grossmann, I. E. (1990). A combined penalty function and outer-approximation method for MINLP optimization. Computers and Chemical Engineering, 14, 769–782. Winston, W. L. (1994). Operations Research—Applications and Algorithms (3rd ed.). Belmont, California: Duxbury Press.