Reference trajectory optimization under constrained predictive control.

Lam, David K. ; Baker, Rhoda ; Swartz, Christopher L.E. 等

INTRODUCTION

Chemical process systems often undergo transitions from one steady-state operating point to another. This could be prompted, for example, by a shift in the steady-state economic optimum due to changing process or economic parameters, or a response to a market demand for different product specifications. The latter is particularly common in polymer plants, where there has been a shift from single products toward the production of multiple polymer grades from the same process. Chatzidoukas et al. (2003) report as many as 30-40 polymer grades being produced in a polyolefin plant. This motivates careful consideration of the cost of transitions, and the development and application of operating practices that minimize this cost.

Several studies have applied dynamic optimization to grade transitions (McAuley and MacGregor, 1992; Takeda and Ray, 1999; Wang et al., 2000; Cervantes et al., 2002; Chatzidoukas et al., 2003; Flores-Tlacuahuac et al., 2006; Asteasuain et al., 2006). A set of process variables is computed that minimizes a measure of the cost of the transition, subject to constraints on the inputs and possibly other specification and/or operational constraints. The model states are related to the inputs through a dynamic model. The decision space in the above-cited studies includes the open-loop trajectories of certain inputs.

McAuley and MacGregor (1992) show that plant/model mismatch could result in deviation of product quality variables from their desired values, and advocate the use of feedback control in the implementation of computed optimal transitions. In a subsequent paper (McAuley and MacGregor, 1993), they develop a non-linear model-based controller for a polymerization process, and apply it to track the profiles of output variables determined from open-loop dynamic optimization. Wang et al. (2000) propose the use of a feedforward-feedback control scheme to implement optimal open-loop trajectories computed via dynamic optimization in a polymer grade transition application. Chatzidoukas et al. (2003) formulate grade transition as a mixed-integer dynamic optimization problem, where the determination of a multi-loop PI control structure is included within the optimization framework. Four PI control loops are considered in an application study, with two input trajectories computed for feed-forward control of the polymer density and melt index. Flores-Tlacuahuac et al. (2006) first compute economically optimal steady-state operating conditions and design parameters. Optimal transitions are computed in a subsequent dynamic optimization problem. Open-loop trajectories are computed, and the economics of the transition are not directly considered. Asteasuain et al. (2006) consider grade transition within an optimization-based design and control framework, similar to that of Chatzidoukas et al. (2003). Steady-state operating points are included as optimization decision variables, and an [epsilon]-constraint multi-objective optimization approach is followed in which steady-state economic and transition performance objectives are considered. Kadam et al. (2007) propose a grade transition approach that seeks to satisfy the necessary conditions for optimality associated with an optimal control problem that minimizes the transition time. The resulting policy includes utilization of PI-type controllers over different operating regimes.

In this paper, we consider optimal transitions of a process regulated via constrained predictive control. This is motivated by the necessity for feedback control to achieve effective set-point tracking in the face of disturbances and plant/model mismatch, and the widespread adoption of model predictive control (MPC) as the advanced control strategy of choice within the chemical process and several other industries (Qin and Badgwell, 2003). A key feature is that the reference trajectories of controlled variables are computed, rather than the process inputs themselves (which are determined by the controller). Moreover, the closed-loop dynamics are taken into account in the reference trajectory optimization. Since the process inputs under constrained MPC are determined from the solution of a quadratic programming (QP) problem at every sampling period, the overall reference trajectory optimization problem is multi-level in nature. In the next sections, we describe the formulation of the reference trajectory optimization problem that we consider, discuss an effective solution strategy, and illustrate the performance of the method through application to two case studies.

Our approach follows a similar structure to the dynamic real-time optimization strategy proposed in Kadam et al. (2003). In both approaches, economically optimal set-point trajectories are determined at an upper plant optimization level, and are passed to a model predictive controller. However, Kadam et al. (2003) do not include the dynamics of the MPC system in the calculation of the reference trajectories, although they do propose a methodology for re-calculation of the trajectories during the course of the transition.

We outline, in the remainder of this section, the key difference between the method proposed here, and reference management or command governors that have appeared in relatively recent control literature (Bemporad et al., 1997; Angeli et al., 1998; Bemporad and Mosca, 1998; Sugie and Yamamoto, 2001). In reference management, the primal control system is an unconstrained, typically linear, control system, and the reference signal is manipulated in order to handle constraints on the closed-loop response. In the present application, constrained MPC is applied as the regulatory controller. The key objective of the reference optimization is to effect a required transition in an optimal (typically economic based) manner.

FORMULATION

Reference Trajectory Optimization

Our objective is to compute an optimal reference (set-point) trajectory that is tracked by a constrained MPC controller. The objective function would typically reflect the cost of the transition and would, in general, be a function of the plant inputs, outputs, and states over the optimization horizon considered. Since the plant is assumed to be controlled, it is the closed-loop response that will be considered during the set-point trajectory optimization. The control structure is illustrated in Figure 1. The reference trajectory optimization problem takes the following form:

Minimize { Cost of transition } Subject to

* Bounds on process outputs

* Bounds on process inputs

* Dynamic process model relating inputs to outputs

* Controller equations relating set-point trajectory and measured plant outputs to plant inputs

where the optimization decision variables correspond to the set-point trajectory. For a discrete-time dynamic system, this may be stated mathematically as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (1)

where x [member of] [R.sup.nx.N] nx.N is a vector of plant states over the reference trajectory optimization horizon, N; y [member] [R.sup.ny.N] is a corresponding vector of plant outputs; u [member of] [R.sup.nu.N] is a vector of plant inputs; [y.sub.sp] [member of] [R.sup.ny.N] is a vector of set-point trajectories; [x.sub.mpc] [member of] [R.sup.nx*P.N] is a vector of MPC model states over the prediction horizon, P, for each time point in the outer-level reference trajectory optimization; [y.sub.mpc] [member of] [R.sup.ny.P.M] is a corresponding vector of MPC model outputs; [u.sub.mpc] [member of] [R.sup.ny.P.N] is a vector of MPC inputs over the input move horizon, M, for each time point in the outer-level optimization; and r [member of] [R.sup.ny.P.N] is a vector of set-point trajectories utilized at the MPC control level. The trajectories, r, are directly related to ysp, but shifted in time to account for the moving horizon of the MPC controller.

h includes the MPC controller equations, which for constrained MPC cannot be expressed as an explicit, continuous function. The control moves at every time step are determined through the solution of a quadratic programming problem. Since the dynamic optimization problem consists of a primary objective, [PHI], as well as MPC optimization subproblems, the result is a multi-level optimization problem. This is illustrated in Figure 2, where the constraints involving the MPC variables would be included within the corresponding MPC subproblems. The following sections describe the MPC algorithm, and a solution approach for this multi-level problem.

[FIGURE 1 OMITTED]

It is also possible to impose further constraints on the reference trajectory. In this paper we explore the following variations:

* Number of allowed changes

Instead of including the full reference trajectory in the design variables, only a limited number of set-point changes are allowed. This may be implemented through the constraints

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (2)

where [y.sub.sp](k) defines the set-point trajectory over the optimization horizon and NAC is the number of allowed set-point changes.

* Set-point hold

In this variation, the reference trajectory is held at a particular value for a specified number of sampling periods, SPH, before being allowed to change again. This may be stated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

where K = {k| k mod SPH [not equal to] 0}. The mod operator in the definition of K gives the remainder upon dividing k by SPH.

* First-order filter

The first-order low pass exponential reference filter is similar to the structure discussed within traditional reference management literature, and is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (4)

where [y.sub.tgt] is the desired set-point target, and the optimal closedloop filter time constant, [f.sub.i], is determined to shape the closedloop response. The single tuning parameter is very appealing, because it offers simplicity of design and tuning, and is easily

[FIGURE 2 OMITTED]

tuned on-line. However, the arbitrary structure of the first-order filter may limit performance. The corresponding difference equation in the discrete-time domain is given by:

[y.sub.sp] (k) = (1 - [f.sub.i]) [y.sub.tgt] + [f.sub.i][y.sub.sp] (k - 1) (5)

Model Predictive Control

Model predictive control (MPC) utilizes an internal dynamic model to predict future process outputs over a prediction horizon, P, in response to future input changes over a control move horizon, M. An optimization problem is formulated, typically to minimize a scalar measure of the deviation of the predicted outputs from a desired set-point and the severity of the input action, and solved to give an optimal set of input moves. The calculated inputs corresponding to the first sampling period are applied to the plant. At the end of the sampling period, the process is repeated, with the most recently measured outputs used to adjust the predicted outputs. This results in a receding horizon control strategy with feedback to compensate for disturbances and model uncertainty. A comprehensive treatment of MPC is given in Maciejowski (2002).

We consider here the following state-space MPC formulation, where the optimization problem to be solved at each sampling period takes the form,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (6)

where x [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is a vector of predicted states, u [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is a vector of predicted inputs, y [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is a vector of the predicted outputs, and r [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is a specified reference trajectory. A [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and C [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] are linear(ized), discrete-time, statespace matrices. The norms in the objective function are defined as

[[parallel]x[parallel].sub.Q] = [x.sup.T]Qx

and y(k + i|k) represents the predicted value of the outputs at time step k + i, based on information available at time step k. A similar definition applies to the state and input vectors, x and u, respectively.

d(k + i|k) represents a disturbance estimate, which in the original dynamic matrix control (DMC) formulation (Cutler and Ramaker, 1979), is taken to be constant over the prediction horizon and computed as the difference between the measured outputs and predicted outputs, using information available at the previous time step. Using the present notation, this would correspond to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

While hard constraints on the outputs may be included, this is often avoided in practice as it may result in closed-loop instability or infeasible QP problems (Zafiriou, 1990; Muske and Rawlings, 1993; Qin and Badgwell, 2003); we therefore impose these instead at the outer set-point optimization level.

We note also that referring to the notation of the previous section, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

that is, a composite vector comprising all the MPC outputs for

k = 1, ... ,N. [x.sub.mpc] and [u.sub.mpc] are similarly defined.

Solution Approach

As discussed earlier, the set-point trajectory optimization problem we consider is multi-level in nature, due to the MPC quadratic programming subproblems required for generation of the closed-loop response. One method for solving this is to use a sequential approach in which the set-point trajectory optimization and closed-loop simulation are carried out iteratively. A potential drawback with this approach is the presence of derivative discontinuities induced by the hard constraints within the MPC controller.

Baker and Swartz (2005) present a simultaneous optimization formulation that accounts for the closed-loop behaviour of a constrained model predictive controller by rewriting the constrained MPC quadratic programming problem at each time step in terms of the Karush-Kuhn-Tucker (KKT) optimality conditions. Consider the general form of a quadratic programming problem,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where x [member of] [R.sup.n]. The KKT conditions for this QP problem can be written as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (7)

where [lambda] and v are the associated Lagrange multipliers.

Reformulation of the inner optimization problems by their KKT conditions results in a single-level optimization problem with complementarity constraints ([x.sub.i][v.sub.i] = 0). Mathematical programs with complementarity constraints (MPCCs) generally cannot be solved directly using standard non-linear algorithms and may require reformulation of the NLP or alternative algorithmic strategies.

Baker and Swartz (2005) applied an interior-point algorithm as implemented in the software package, IPOPT-C (Raghunathan and Biegler, 2003), that is tailored to handle the complementarity constraints in the primal problem; this is the method that is used in the present study. The KKT conditions of the MPCC contain the complementarity constraints present in the primal problem, as well as those arising from the primal problem inequality constraints. The approach relaxes the complementarity constraints as

[x.sub.i][v.sub.i] [less than or equal to] i = 1, ... ,n,

where [micro] is a small positive parameter, and solves a series of subproblems during which the value of ??is gradually reduced toward zero. While this is similar to the manner in which complementarity constraints are dealt with in standard interior point algorithms (Wright, 1997), a key feature is that the complementarity constraints in the primal problem are handled differently from the general, non-linear equality constraints.

Alternatively, Ralph and Wright (2004) showed that, under similar assumptions of regularity and MPCC linear independence as used in Baker and Swartz (2005), the solution to the original MPCC can be obtained if the complementarity constraints are moved to the objective function and are multiplied by a sufficiently large penalty parameter. However, the size of this parameter is not known a priori, and choosing too large a parameter can potentially lead to scaling issues.

Baker and Swartz (2005) used this single-level formulation to solve constraint back-off problems in which the operating point is treated as a design variable and "backed-off" from the steady-state economic optimum in order that output constraints are not violated in the event of a disturbance. Thus, the principle design variable was the controlled variable set-point, which is kept constant over the prediction horizon. In this study, the model predictive controller reference trajectory is included as part of the decision space. The extension, moreover, does not add any further non-linearity to the single-level problem.

The case studies will demonstrate how the single-level formulation can be used to determine optimal reference trajectories for problems with set-point transitions. All case studies are solved with AMPL as the modelling language on a 3.0GHz Pentium IV with 1GB RAM running Debian Linux.

Reference Management with Dynamic Economics

In this section we formulate an objective function that is based directly on the economics of a grade transition. Consider the scenario illustrated in Figure 3, where the transition from product A to product B takes place. The band around [y.sup.A.sub.1] represents the limits within which quality variable y1 must lie in order to meet the quality specifications for product A; a similar specification region is shown for product B. Acceptable product further requires variable [y.sub.2] to lie within the region shown in the lower diagram in Figure 3.

[FIGURE 3 OMITTED]

In the formulation below, we use continuous approximations to capture the discontinuous switching of contributions to the objective function as variables move into and out of the specification regions. Referring to the top diagram in Figure 3, we define the continuous function, [R.sub.1]([y.sub.1]), that is approximately zero for [y.sub.1] < [y.sup.A.sub.1] - [delta][y.sup.A.sub.y], and unity for [y.sub.1] > [y.sup.A.sub.1] - [delta]py.sup.A.sub.1]. This may be achieved by the hyperbolic tangent function,

[R.sub.1] = 1/2 tanh [[gamma]([y.sub.1] - [y.sup.A.sub.1] + [delta][y.sup.A.sub.1])] + 1/2 (8)

where the parameter, [gamma], determines the steepness of the switching function. The shape of [R.sub.1] is shown in Figure 4.

We similarly define [R.sub.2]([y.sub.1]) such that [R.sub.2] [approximately equal to] 1 for [y.sub.1] < [y.sup.A.sub.1] < [y.sup.A.sub.1] + [delta]py.sup.A.sub.1], and [approximately equal to] 0 for [y.sub.1] > [y.sup.A.sub/1] + [delta][y.sup.A.sub.1]:

[R.sub.1] = 1/2 tanh [[gamma]([y.sup.Asub.1] + [delta][y.sup.A.sub.1] - [y.sub.1])] + 1/2 (9)

Thus, [R.sub.1][R.sub.2] [approximately equal to] 1 if [y.sub.1] lies within the specification region for product A and is approximately zero otherwise. We similarly define [R.sub.3] and [R.sub.4] such that [R.sub.3][R.sub.4] [approximately equal to] 1 if [y.sub.1] lies within the specification region for product B; and [R.sub.5] and [R.sub.6] such that [R.sub.5][R.sub.6] [approximately equal to] 1 if [y.sub.2] lies within its specification bounds, which here are common for products A and B. The revenue may thus be formulated as

R(t) = F[P.sub.A][R.sub.1][R.sub.2][R.sub.5][R.sub.6] + F[P.sub.B][R.sub.3][R.sub.4][R.sub.5][R.sub.6] (10)

where F is the product flow rate, and [P.sub.A] and [P.sub.B] are the prices ($ per unit of F) of products A and B, respectively. We wish to maximize the profit over the transition,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

where C(t) represents the costs associated with the input streams and utilities.

[FIGURE 4 OMITTED]

We use an economic objective function of this form (with a discrete approximation of the integral) in the polymerization grade transition case study to follow. A similar construct for handling the discontinuities in the grade transition objective function is proposed in Tousain (2002), who also investigates a mixed-integer formulation. The methodology readily extends to more than two specification variables.

Two-Tier Optimization

Depending on the objective function and application, the optimal reference trajectory optimization formulation may yield set-point trajectories with relatively high variation and, in some cases, a non-unique solution. We describe here a hierarchical two-tier approach, in which the economic objective is first maximized, followed by a subsequent optimization problem in which the setpoint variation is minimized, subject to the economic performance meeting a threshold determined from the optimal objective function value obtained in the first optimization.

Consider the set-point trajectory optimization problem (1) with an economic objective function, such as given by (11). Denote the optimal objective function value as ?*. We then pose the second-tier optimization problem as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (12)

This is similar to the ?-constraint approach in multi-objective optimization (Luyben and Floudas, 1994; Asteasuain et al., 2006). By solving the first problem, we obtain an upper bound on economic performance and, in the subsequent problem, specify by how much we wish to sacrifice economic performance to reduce set-point variability. We explore this approach in the second case study.

CASE STUDIES

Case Study 1

We consider here the single-input, single-output, non-minimum phase system,

[G.sub.p] (s) -1.2s + 0.1601/[s.sup.2] + 0.40s + 0.16 (13

The system is controlled using input-constrained model predictive control, with output constraints imposed at the outer, setpoint trajectory optimization level. A prediction horizon of 30 and an input horizon of 10 are chosen, with a sampling time of 1 min. The model predictive control tuning uses an output to input weighting ratio of 1:4, with the closed-loop output required to be within the range [-0.50, 1.0], while the inputs are constrained to lie within the range [0, 1.0]. The system is initialized with the output at 0, and the target output value is 1.0. The objective function was taken as the sum of the squared errors (SSE) between the actual output trajectory and the target value over the set-point optimization time horizon:

[PHI] = [N.summation over (k=1)] [(y (k) - [y.sub.tgt).sup.2] (14)

Single step change

In the first scenario, the reference trajectory is constrained to be a single step change from zero to the target value. This corresponds to a simulation; there are no bounds imposed on the output trajectory at the outer level since there are insufficient degrees of freedom to modify the solution if an output constraint is violated. The input and output trajectories are plotted in Figure 5, from which it can be seen that the lower bound is violated at t = 3 min. This implies that the controller tunings would have to be adjusted in order to ensure that the bound is not violated if a single step reference trajectory is used. While it is possible to detune the controller to prevent the constraint violation, this may be undesirable in practice, resulting in increased transition times. By optimizing the reference trajectory, there is an opportunity to maintain aggressive performance while satisfying output constraints.

Full reference trajectory optimization

Next, we examine the effect of optimizing the full reference trajectory. Additional constraints were imposed on the system, namely, the requirement that the final value of the reference trajectory be equal to the target value, and that the reference trajectory not exceed the upper and lower bounds of the output. The optimization was solved using IPOPT-C in 1.32 s, and returned an objective value of 16.186.

The input and output trajectories are plotted in Figure 6. Because the system exhibits inverse response behaviour, it is not surprising to see that the optimal reference trajectory initially moves in the direction opposite to that of the target value. From these graphs we also see that the optimal reference trajectory is adjusted until the end of the simulation horizon, instead of immediately settling to the target value. Since this example has the most degrees of freedom, it is expected that the objective value obtained will act as an upper bound on the achievable performance .

Limited number of allowed changes

Here, we examine the effect of limiting the number of changes that the reference trajectory is allowed to make. The optimization was solved using IPOPT-C in 1.00 s for the case with 10 allowed set-point moves, for which the input and output trajectories are plotted in Figure 7. Once again an initial reference trajectory move is made in the opposite direction to the final target value.

[FIGURE 5 OMITTED]

The effect of varying the number of allowed changes was also investigated, with the results illustrated Figure 8. It should be noted that case studies with fewer than 7 allowed changes were infeasible. From Figure 8 it is apparent that the more degrees of freedom afforded to the reference trajectory optimization, the better the performance. However, this improvement in performance eventually reaches a point of diminishing returns. It should be mentioned again that if the number of allowed changes is too few, feasible output and input trajectories might not exist.

Set-point hold

In this study, we examine the effect of maintaining the reference trajectory at a particular value for a specified number of steps before allowing the trajectory to change. The optimization was solved using IPOPT-C in 1.00 s for the case with a set-point hold of 3 sampling periods, with the corresponding input and output trajectories plotted in Figure 9. In this case, unlike the previous two examples, it is no longer optimal for the reference trajectory to move initially in the direction opposite to that of the target value. However, with a set-point hold of 2 sampling periods, the optimal reference trajectory again exhibits an inverse response characteristic.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

Further case studies were run in order to investigate the effect of varying set-point hold duration on the performance; the results are reported in Table 1. The general trend indicates that the longer a particular value is held, the worse the performance. However, the performance for a set-point duration of 4 is worse than the performance for a duration of 5. This could possibly be explained by noting that the locations of the set-point changes have an effect on the performance as well as the lengths for which the set-point values are maintained constant.

Discrete reference filter

In this example, the discrete reference filter is applied to the reference trajectory optimization problem, with the filter parameter [f.sub.i] as the optimization variable, rather than the individual elements of the set-point trajectory. The system was constrained so that the initial value of the reference trajectory be equal to zero, but it was not required that the final value of the reference trajectory be equal to the target value since this value would be determined by the discrete filter. The optimization was solved using IPOPT-C in 2.45 s and returned a value of 0.9097 for [f.sub.i] with an SSE of 16.417. The input and output trajectories are plotted in Figure 10 and demonstrate a gradual increase toward the target value.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

Discussion

SSE values for various models of the reference trajectory are given in Table 2. It is not valid to compare the SSE for the case with a single step in the reference trajectory since this example violates the output constraints.

From Table 2 we see that optimization of the full reference trajectory profile results in the lowest sum of squared errors and acts as an upper bound on the achievable performance. However, the set-point variation displayed in Figure 6 is potentially undesirable, although the input and output trajectories do not appear to be unduly affected by this.

In contrast to the full reference trajectory optimization case study, the first-order filter method has the highest sum of squared errors. Since the closed-loop system exhibits inverse response behaviour rather than a first-order response, this could have had a negative impact on the performance. However, the reference trajectory profile is smooth in comparison to the other cases.

[FIGURE 10 OMITTED]

Case Study 2

We consider here a styrene polymerization reaction process based on Maner et al. (1996), to which an approximate linear model was fitted,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (15)

where the outputs, [y.sub.1] and [y.sub.2], represent the number average molecular weight (NAMW) and temperature, respectively, controlled by inputs [u.sub.1] and [u.sub.2], representing the initiator flow rate, [Q.sub.i], and coolant flow rate, [Q.sub.c].

The following constraints (in deviation form) are applied:

-80 [less than or equal to] [u.sub.1] [less than or equal to] 42 -471.6 [less than or equal to] [u.sub.2] [less than or equal to] 28.4 -8.5 [less than or equal to] [y.sub.1] [less than or equal to] 21.5 -2 [less than or equal to] [y.sub.2] [less than or equal to] 2

The following constraints were imposed on the reference trajectory ,

-20 [less than or equal to] [y.sub.sp,1] [less than or equal to] 20 -2 [less than or equal to] [y.sub.sp,2] [less than or equal to] 2

First, we consider the optimization of the set-point trajectory when the objective function is based on minimizing the deviation of the outputs from their target values, then consider the effect of using an economic based objective function.

Objective based on deviation from target

The system was controlled using input-constrained MPC, with a prediction horizon, control move horizon and sampling period of P = 20, M = 5, and [DELTA]t = 1 h. The desired target specifications are [y.sub.1] = 20 and [y.sub.2] = -1. The output weights are Q = diag(10, 100), with input move suppression weights, R = I. The performance objective, [PHI], was weighted consistent with the output weights used in the MPC controller.

Single step change

We first constrain the reference trajectory to be a single step change from zero to the target value. The reference trajectory is therefore not optimized, and no bounds are imposed on the output trajectory. The input and output trajectories are plotted in Figure 11. The optimal objective function value is 18 048.

Reference trajectory optimization

Figure 12 shows the results of reference trajectory optimization with a set-point hold of 5 sampling periods. The trajectory optimization results in an improved objective function value of 17 772.

The input and output trajectories for the single step and optimized reference trajectory cases are plotted together in Figure 13. In both cases, u1 touches lower bound, illustrating the use of the constraint-handling feature of the MPC controller.

The primary benefit of the reference trajectory optimization is that the temperature profile returns to the target value approximately 10 h sooner than the single step case. It does so by temporarily reducing the flow rate of cooling water below the final steady state value, whereas the profile for the single setpoint change is more gradual. This comes at the expense of slightly oscillatory behaviour of the NAMW around the target value; however, this difference is hardly noticeable when the two profiles are compared.

Next we consider the same polymerization example using an economic-based objective function.

Objective based on economics

In this case study, we use the economic-based objective function (11) to take into account the loss in revenue for producing off specification product outside product quality tolerances. The following cost function, based on the cost of monomer, initiator and cooling water, is used:

C (t) = [C.sub.m][Q.sub.m] + [C.sub.i][Q.sub.i] (t) + [C.sub.c][Q.sub.c] (t) (16)

where [Q.sub.m], [Q.sub.i] and [Q.sub.c] are the monomer, initiator and cooling water flow rates, respectively, with the coefficients preceding them representing the corresponding costs. [Q.sub.m] has a constant value of 378 L/h for this example, and the initial steady-state value of the initiator is 108 L/h. The revenue is calculated using (10), with the product flow rate, F, calculated as

F (t) = [Q.sub.i] (t) + [Q.sub.m] (17)

Solvent that is added (Maner et al., 1996) is assumed to be recovered and recycled, and its effect on the transition cost neglected. The feed, product and utility cost information is reported in Table 3.

[FIGURE 11 OMITTED]

Here, the end-point condition on the desired grade is enforced via a constraint on the reference trajectory. The revenue from the initial steady state and final desired grade is taken into account, and the value of off-specification product created during the transition is assumed to be negligible.

Single step change

The input and output trajectories resulting from a single setpoint change and shown in Figure 11 result in an objective value of $3963. This provides a basis of comparison against which the subsequent reference trajectory optimization results may be compared.

[FIGURE 12 OMITTED]

Reference trajectory optimization

In this case, the reference trajectory is optimized to allow 30 changes before it is fixed at the target value, with input and output constraints as described previously. The input and output trajectories are plotted in Figure 14, and give a significantly improved objective function value of $4723.

In this example, u1 touches the lower bound at t = 3 h, illustrating that the constraint-handling functionality of the model predictive controller is utilized. The optimal reference trajectory for the NAMW appears to be a delayed single step change, most likely because valuable product is being produced at the initial steady state as well as the final steady state. In contrast, the reference trajectory for the temperature is highly oscillatory and frequently lies at the upper and lower bounds of reference trajectory profile. We notice that the temperature is slightly higher than in the base case. This may be explained as follows. Any temperature profile that lies between the maximum and minimum specification bounds is considered acceptable in terms of its impact on the objective function through the product revenue. However, by operating at an increased temperature, a reduction in the cooling water requirement and corresponding cost may be achieved.

[FIGURE 13 OMITTED]

Two-tier optimization

We next explore the application of the two-tier optimization approach described earlier to mitigate the temperature set-point variation. The secondary objective is to minimize the the changes between subsequent set-points while maintaining the profit within 5% of the previously determined value. The resulting trajectories are plotted in Figure 15, which yield a profit of $4486.

We observe that the changes in the reference trajectory for the temperature profile are significantly reduced, and that the set-point trajectory does not hit the constraints. More importantly, this approach allows users to quantify exactly how much profit they are willing to sacrifice for a more conservative set-point trajectory profile. Despite the decrease in profit, this optimal trajectory still offers an improvement over the single step change.

[FIGURE 14 OMITTED]

We recognize that there is a degree of error associated with the hyperbolic tangent smoothing function. Large values of ? give a tighter approximation, but at the expense of a higher degree of non-linearity and potential numerical difficulties. Alternative formulation strategies are being explored.

CONCLUSIONS AND FUTURE WORK

Many chemical processes undergo frequent transitions, motivating the application of economic optimization to the transition problem. In this paper, a strategy has been proposed for reference trajectory optimization that includes the closed-loop dynamics of a constrained MPC controller used to regulate the underlying process. A simultaneous solution framework has been described, whereby the resulting multi-level optimization problem is transformed into a single-level mathematical program with complementarity constraints (MPCC). Application of an interior point approach to the solution of the MPCC was found to be reliable and efficient. The approach was applied to two case studies, the second of which considered the economics of a grade transition.

[FIGURE 15 OMITTED]

While the applications in this paper were to linear dynamic systems, the method is readily extended to non-linear processes, with the non-linear dynamic model discretized using a technique such as orthogonal collocation on finite elements. This is planned as the subject of a future communication.

Manuscript received January 31, 2007; revised manuscript received May 7, 2007; accepted for publication June 8, 2007.

REFERENCES

Angeli, D., A. Casavola and E. Mosca, "Enforcing Constraints in Industrial Control Systems: Command Governor Tools," Proc. IEEE International Conference on Control Applications, 746-750 (1998).

Asteasuain, M., A. Bandoni, C. Sarmoria and A. Brandolin, "Simultaneous Process and Control System Design for Grade Transition in Styrene Polymerization," Chem. Eng. Sci. 61, 3362-3378 (2006).

Baker, R. and C. L. E. Swartz, "Interior Point Solution of Multilevel QP Problems Arising in Embedded MPC Formulations," Paper 12e, AIChE Annual Meeting, Cincinnati (2005).

Bemporad, A., A. Casavola and E. Mosca, "Nonlinear Control of Constrained Linear Systems via Predictive Reference Management," IEEE Trans. Auto. Control 42(3), 340-349 (1997).

Bemporad, A. and E. Mosca, "Fulfilling Hard Constraints in Uncertain Linear Systems by Reference Managing," Automatica 34(4), 451-461 (1998).

Cervantes, A. M., S. Tonelli, A. Brandolin, J. A. Bandoni and L. T. Biegler, "Large-Scale Dynamic Optimization for Grade Transitions in a Low Density Polyethylene Plant," Comp. Chem. Eng. 26, 227-237 (2002).

Chatzidoukas, C., J. D. Perkins, E. N. Pistikopoulos and C. Kiparissides, "Optimal Grade Transition and Selection of Closed-Loop Controllers in a Gas-Phase Olefin Polymerization Fluidized Bed Reactor," Chem. Eng. Sci. 58, 3643-3658 (2003).

Cutler, C. R. and B. L. Ramaker, "Dynamic Matrix Control--A Computer Control Algorithm," AIChE National Meeting, Houston, Texas (1979).

Flores-Tlacuahuac, A., L. T. Biegler and E. Saldivar-Guerra, "Optimal Grade Transitions in the High-Impact Polystyrene Polymerization Process," Ind. Eng. Chem. Res. 45, 6175- 6189 (2006).

Kadam, J. V., W. Marquardt, M. Schlegel, O. H. Bosgra, T. Backx, P.-J. Brouwer, G. Dunnebier, D. van Hessem, A. Tiagounov and S. de Wolf, "Towards Integrated Dynamic Real-Time Optimization and Control of Industrial Proceses," I. E. Grossmann and C. M. McDonald, Eds., "FOCAPO 2003," 593-596 (2003).

Kadam, J. V., W. Marquardt, B. Srinivasan and D. Bonvin, "Optimal Grade Transition in Industrial Polymerization Processes via NCO Tracking," AIChE J. 53(3), 627-639 (2007).

Luyben, M. L. and C. A. Floudas, "Analyzing the Interaction of Design and Control - 1. A Multiobjective Framework and Application to Binary Distillation Synthesis," Comp. Chem. Eng. 18(10), 933-969 (1994).

Maciejowski, J. M., "Predictive Control with Constraints," Prentice Hall (2002).

Maner, B. R., F. J. Doyle III, B. A. Ogunnaike and R. K. Pearson, "Nonlinear Model Predictive Control of a Simulated Multivariable Polymerization Reactor using Second-Order Volterra Models," Automatica 32(9), 1285-1301 (1996).

McAuley, K. B. and J. F. MacGregor, "Optimal Grade Transitions in a Gas Phase Polyethylene Reactor," AIChE J. 38(10), 1564-1576 (1992).

McAuley, K. B. and J. F. MacGregor, "Nonlinear Product Property Control in Industrial Gas-Phase Polyethylene Reactors," AIChE J. 39(5), 855-866 (1993).

Muske, K. R. and J. B. Rawlings, "Model Predictive Control with Linear Models," AIChE J. 39(2), 262-287 (1993).

Qin, S. J. and T. A. Badgwell, "A Survey of Industrial Model Predictive Control Technology," Control Eng. Prac. 11, 733-764 (2003).

Raghunathan, A. U. and L. T. Biegler, "Mathematical Programs with Equlibrium Constraints (MPECs) in Process Engineering," Comp. Chem. Eng. 27, 1381-1392 (2003).

Ralph, D. and S. J. Wright, "Some Properties of Regularization and Penalization Schemes for MPECs," Optimization Methods and Software 19(5), 527-556 (2004).

Sugie, T. and H. Yamamoto, "Reference Management for Closed Loop Systems with State and Control Constraints," Proc. ACC, 1426-1431 (2001).

Takeda, M. and W. H. Ray, "Optimal-Grade Transition Strategies for Multistage Polyolefin Reactors," AIChE J. 45(8), 1776-1793 (1999).

Tousain, R. L., "Dynamic Optimization in Business-Wide Process Control," PhD Thesis, Delft University of Technology, Delft University Press, Delft, The Netherlands (2002).

Wang, Y., H. Seki, S. Ohyama, K. Akamatsu, M. Ogawa and M. Ohshima, "Optimal Grade Transition Control for Polymerization Reactors," Comp. Chem. Eng. 24, 1555-1561 (2000).

Wright, S. J., "Primal-Dual Interior Point Methods," SIAM, Philadelphia, PA (1997).

Zafiriou, E., "Robust Model Predictive Control of Processes with Hard Constraints," Comp. Chem. Eng. 14, 359-371 (1990).

David K. Lam ([dagger]), Rhoda Baker and Christopher L. E. Swartz * Department of Chemical Engineering, McMaster University, 1280 Main Street West, Hamilton, ON, Canada L8S 4L7

* Author to whom correspondence may be addressed.

E-mail address: swartzc@mcmaster.ca

Table 1. SSE with increasing duration of the Set-Point Hold (SPH)

SPH SSE

1 16.186
2 16.187
3 16.294
4 16.394
5 16.385

Table 2. SSE for different reference trajectory scenarios for Example 1

Case SSE

Full 16.186
NAC 10 16.250
SPH 3 16.294
Filter 16.417

Table 3. Pricing information for polymerization example

Parameter Value

[C.sub.m] Monomer cost 0.98 $/L
[C.sub.i] Initiator cost 0.07 $/L
[C.sub.c] Cooling water cost 0.05 $/[m.sup.3]
[P.sub.A] Price of products A and B 1.76 $/L