Efficient Linear Programming Formulations of Model Predictive Control

2134

Ind. Eng. Chem. Res. 2005, 44, 2134-2144

Efficient Linear Programming Formulations of Model Predictive Control for Cross-Directional Control of a Paper Machine Daniel R. Saffer, II† Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716

Francis J. Doyle, III* Department of Chemical Engineering, University of California, Santa Barbara, Santa Barbara, California 93106

Efficient formulations of linear programming based model predictive control are considered for cross direction (CD) control of a dual headbox paper machine. A number of techniques including coincidence points, efficient mathematics, implementation strategies, and the solution of progressively more complicated optimization problems are considered for both their performance and efficiency. 1. Overview The paper machine is the final unit operation involved in the manufacturing of paper from wood pulp. In a typical machine, a mixture of wood pulp, water, and chemical additives are delivered to the headbox of the machine. The pulp is delivered to a fourdrinier table by way of an opening in the headbox and from there to a series of heated cylinders where the paper sheet is dried and collected in rolls at the end of the machine opposite of the headbox. The control of this process is of significant importance since the properties of end product directly depend on the regulation of these properties within the paper machine. The paper machine control problem is typically divided into two directional domains running parallel with and perpendicular to the direction of flow of the sheet through the machine, otherwise known as machine direction (MD) and cross direction (CD), respectively. One of the properties of interest in paper production is basis weight, which is the weight per area of the final product. MD control of the average basis weight of the sheet is performed by controlling the bulk flow of pulp to the headbox and as such is a single-input-single-output control problem. The focus of this study is the much more complicated problem of CD control of basis weight. CD control uses of a series of actuators to manipulate a slice lip at the headbox of the paper machine to control the local flow of pulp leaving the headbox. In this paper, the focus is on the constrained model predictive control (MPC) for the CD control of a dual headbox paper machine that uses full array wet end sensors to measure basis weight. The large scale nature, constrained control space, and fast sampling rate makes this problem challenging for MPC implementation. Recent efforts by a number of groups have made strides in achieving the goal of on-line implementation of MPC to the large-scale CD control problem. Doyle, * To whom correspondence should be addressed. Tel. (805) 893-8133. Fax: (805) 893-4731. E-mail: doyle@ engineering.ucsb.edu. † Current address: Weyerhaeuser, Red River, Highway 480, P.O. Box 377, Campti, LA 71411.

Pekny, and co-workers explore the nontraditional use of a linear programming based MPC (LP-MPC) CD control scheme1 and later refined the LP algorithm for fast on-line implementation through a four-stage calculation.2 More recently, Biegler and co-workers3 explore a computationally efficient quadratic programming based MPC (QP-MPC) algorithm based on a dual active set method deployed with efficient object-oriented programming. To reduce the number of constraints to search over, their QP solver uses soft input constraints rather than hard input constraints. The studies of both the Doyle and Biegler research groups are developed for simulated single headbox paper machines with scanning sensor technologies. One of the first documented implementations of QP-MPC to a dual headbox CD control problem is attributed to Backstro¨m et al.4 The MPC algorithm is applied to CD control of basis weight and moisture profiles using the actuators of one headbox of a dual-headbox machine. In their approach, the effects of the secondary headbox are introduced to this machine as input disturbances. While the previously noted research in MPC for CD control aimed to solve the full-scale CD control problem efficiently, other groups have explored techniques to reduce the size of the MPC programming problem through approximations to either the input-output parameter space or the constraint set. Rigopoulos et al.5 explore the use of adaptive principle components analysis to reduce the number of inputs and outputs that are optimized in the CD control problem. Van Antwerp and Braatz6 approximate the full-scale constraint space by an ellipsoid in order to compute the optimal control solution by a bisection technique. While these examples provide efficient solutions to the CD control problem, they both do so at the expense of suboptimal control calculations as compared to techniques that solve the full-scale MPC problem for CD control. The problem of computationally efficient linear programming solutions for LP-MPC and QP-MPC is addressed throughout this paper. Although industrial standards for CD quality are standard deviation values, σCD and σMD, based on a 2-norm objective, MPC based on linear objectives, either the sum of absolute values

10.1021/ie049303i CCC: $30.25 © 2005 American Chemical Society Published on Web 02/23/2005

Ind. Eng. Chem. Res., Vol. 44, No. 7, 2005 2135

Figure 1. Paper machine with two headboxes consisting of primary (up) and secondary slice lip actuators (us), as well as wet end (yp and ys) full array sensors. In the top view (a), the dashed line represents an actuator-sensor lane, indicating the center of the response of the actuator at the sensor location. The times in the side view (b) represent approximate time delays at various points in the model.

(1-norm) or the minimization of the maximum value (1norm), is explored here. These σ values, originally defined in a TAPPI standard,7 are

σMD )

σCD )

x∑ x∑ N

(xji - xj)2

i)1

N-1

M

(xjj - xj)2

j)1

M-1

(1)

(2)

where the overall average of the individual data points N M (xij) is xj ) ∑i)1 ∑j)1 xij/NM, the sample average of the ith N sample is xji ) ∑i)1 xij/M, the lane average of the jth M lane is xjj ) ∑j)1xij/N, and i and j are indices for the sample time and CD position, respectively. These LPMPC algorithms are compared to standard QP-MPC for their ability to efficiently reduce product variance in the face of sustained disturbances. To address the issue of computational efficiency, the four-stage algorithm that was originally developed by Dave, Doyle, and Pekny2 is extended to the dualheadbox paper machine benchmark problem with full array sensors. The technique is also compared to other algorithms, including a simpler algorithm that utilizes

a computationally efficient commercial LP solver, much like the work of Backstro¨m et al.4 The use of coincidence points is explored to reduce the computational burden of the four-stage and the simpler algorithm just mentioned by eliminating a number of constraints and variables within the linear program as an alternative to hard constraints in an effort to efficiently solve the LP-MPC problem for CD control. Finally, the LDMC approach of Morshedi et al.8 is explored for its ability to either minimize the CD quality variables by itself or to warm start the solution to a quadratic program in QP-MPC. 1.1. Benchmark CD Control Problem. The benchmark problem studied in this paper is a proprietary Matlab and Simulink model of dual-headbox paper machine, which was first introduced in Saffer et al.9 Like many other sheet and film processes, the simulation model assumes first-order-plus-time-delay dynamics between the slice lip actuator and full array sensor locations with a Toeplitz symmetric input-output interaction matrix. Measurement of cross-direction basis weight variations is performed with full-array wet-end sensors at the end of each of the Fourdrinier tables (Figure 1). These sensors have a sample time of 1 s. The objective in each case is to control the CD variation of basis weight in the presence of sustained disturbances. While robustness measures are not addressed in this work, it is well-known that CD processes are ill-

2136

Ind. Eng. Chem. Res., Vol. 44, No. 7, 2005

Figure 2. Profiles of the input disturbance at the primary headbox, dp, and secondary headbox, ds, that are employed in the case studies throughout this paper. A representative closed-loop results using the four stage technique of Section 2 is shown at the primary (up) and secondary (us) slice lip actuators as well as the primary (yp) and secondary (ys) full array wet-end sensors.

conditioned due to the Toeplitz symmetric gain, with a large number of singular values of the plant model close to zero.10 As such, any component of the error profile that aligns with the direction of these singular values will not be attenuated effectively. With a properly identified model, these algorithms have been found to be stable to many disturbance profiles, although zero steady-state error is rarely, if ever, attainable owing to the nonsquare nature of the problem. The primary headbox has 45 slice lip actuators while the secondary headbox has 36 actuators, both of which are used to control basis weight at the primary and secondary wet-end full array sensor locations. Each sensor bank contains 360 databoxes, yielding an 81 input by 720 output CD control problem. Each closedloop trial is subjected to the input disturbances shown

in the top two plots of Figure 2 at the primary headbox, dp, and secondary headbox, ds. The disturbances at both headboxes were generated by a summation of sine waves of various amplitudes and spatial frequencies. In addition, the primary headbox had a dominantsfive times the maximum amplitude of the summation of wavesssine wave that spanned the entire array much like a headbox imbalance disturbance. 1.2. Optimization and Model Predictive Control. MPC has become a viable alternative to classical control techniques for multivariable processes that require optimal control in the presence of process constraints. As opposed to other available control techniques, MPC provides the only methodology for handling process constraints in a manner that is consistent with the design and implementation of the controller.11 Other


advantages of using MPC include the scalability to a nonsquare multivariable framework, the ease of handling difficult process dynamics, and the ability to compensate for measured and unmeasured disturbances.12 MPC uses a dynamic model process to solve an on-line optimization within the controller at each time step. The control law that is generated from MPC is open-loop optimal with respect to the objective function that is used within the controller. At the core of any MPC problem is an optimization algorithm. MPC uses a receding horizon approach that samples the output (y˜ k) at time k to predict the process response over a prediction horizon (p). Control values (uk) are optimized over a move horizon (m e p) to minimize the difference between the predicted outputs (yk+i|k) and the reference trajectory (r). A penalty for the control moves and/or control values is usually added to the objective function so that the controller is not overly aggressive or sensitive to noise. The future control moves after time k + m - 1 are set to zero. Once an optimal solution is found, the first control move is implemented and the controller is re-initialized at time k + 1. Originally studied in Dave et al.’s work,2 the LP to be solved within LP-MPC for the CD control problem throughout this work is a mixed ∞/∞-norm penalty on basis weight deviations from setpoint and 1/1-norm penalty on control move velocity: m

min uk|k

fobjective(r,∆u,u) ) ξ +

R∆u[∆ui+ + ∆ui-] ∑ i)1

(3)

s.t. R+ + R- + s1 - Ξ ) 0 -P∆u+ + P∆u- - R+ + R- ) y˜ - r ˜ u - ∆u+ + ∆u- ) u A ˆ u + s2 ) δmax {R+, R-, s1} g 0 0 e {∆u+, ∆u-} e ∆umax 2

0 e s e 2δ u

min

max max

eueu

In the above LP, A ˆ ) toeplitz [1, 2, 1] is a constraint matrix corresponding to penalties on adjacent slice lip actuators on each headbox so that the slice lip is not permanently damaged. The variable ξ is the scalar value representing the maximum deviation from setpoint, Ξ is a one-column vector of this quantity. The values R+ and R- represent the auxiliary variables to represent the absolute value of the predicted error term while the same is true for ∆ui+ and ∆u- for the absolute value of the control move ∆u. The variables s1 and s2 are slack variables associated with the definition of ξ and the bending moment constraints. For this work, the model, P, is a step-response model based directly on the simulation model that relates the predicted inputs to the predicted outputs of the process. This assumption is ideal, and the model would need to be determined through means of identification such as those described elsewhere.9 The step-response model of the dual-headbox benchmark CD control model uses a

move horizon of m ) 1 and a prediction horizon of p ) 7. This model deviates from that in the four-stage algorithm in Dave et al.2 in that P must contain the transient response for a dual headbox machine with full array sensors with a sample time of Ts ) 1 s. Since the previous work focused on slower scanning sensor technology, P is chosen to be the steady-state input-output matrix with m ) p ) 1. Although the newer technology provides more data to the controller, the algorithms to be used within the controller must now be able to compute solutions for larger models, in this case a 14fold increase in the size of P. For the benchmark example in this paper, the current formulation contains 2pNy + 2mNu - 4m ) 10238 constraints, 3pNy + 4mNu - 4m + 1 ) 15441 variables. This is compared to the problem in ref 2 that has 2Ny + 2Nu - 4 ) 1598 constraints and 3Ny + 4Nu - 4 + 1 ) 2481 variables for the same problem using a steady-state model. With this optimization problem comes two sequential operations that must be performed at each time step. The first step is to set up the optimization problem to be solved while the second step is to use the appropriate solver to determine the optimal control move. Each of these two steps can potentially be time-consuming, depending on the values that change from sample time to sample time. One advantage to LP-MPC is the minimal number of changes in the constraint coefficients that need to be made from sample to sample within the LP problem. When using LP-MPC, the right-hand side constraint vector is the only component that changes in LP-MPC formulations at every sample time. The LP problems to be solved throughout this paper make use this property to quickly set up the LP problem at each sample time. In contrast to this, QP-MPC formulations require changes to multiple components that make up the QP problem. The corresponding QP-MPC problem addressed in this work is

min fobjective(r, uk+n|k) ) uk

T uk+n|k (PTP + RuTRu)uk+n|k + (r - y˜ )Puk+n|k (4)

s.t. A ˆ u + s1 ) δmax s1 g 0 -∆umax e ∆u e ∆umax 0 e s1 e 2δmax umin e u e umax In QP-MPC, the objective function matrixes as well as the right-hand side constraint vector all require modification from sample time to sample time. This work develops different solution methodologies to the linear programming problem within LP-MPC. Each LP within the tested algorithms utilizes the simplex routines of the commercial linear programming solver CPLEX version 7.113 that is accessed through a communication routine with Matlab 6.0 without its Java-enabled features. CPLEX is chosen as the solver in these examples since its algorithms have been optimized for efficient computation of large-scale sparse linear programs, which is the case in LP-MPC applications to CD control. The disadvantage to using CPLEX,

2138


however, is that only one LP problem can remain within the CPLEX workspace at any given instant. The advantages and disadvantages that this might bring are discussed throughout this paper. All of the closed-loop simulations are run on an AMD Athlon Thunderbird processor running Linux 6.0 at 800 MHz with 255 MB of RAM and 500 MB of swap space. During the simulations, Matlab and CPLEX were the only user processes running on this machine. 2. Four Stage Approach to Efficient LP-MPC As mentioned in section 1, Doyle, Pekny, and coworkers1,2 were the first to propose that the CD control problem can be controlled through an LP-MPC approach. The four-stage approach of ref 2 attempts to solve the LP-MPC problem by starting with an easily solved LP problem that progressively adds complexity of the full-scale problem along the stages, using the solution from the prior stage to warm-start the next, until the full scale LP is solved. The work in ref 2 demonstrates that the computers of the time (an HP9000/ 770) solve the four-stage algorithm within 0.58-0.73 s for a 79 × 79 input to output CD control problem as compared to 0.9-1.31 s for a standard LP approach to the same size problem. The efficiency of the four-stage algorithm is more evident for larger problems. For example, a 399 × 399 problem is solved in 4.12-43.43 s using the four-stage approach as compared to 76.89130.58 s for the standard LP. The following explains the four-stage technique and the modifications in this work for the dual headbox CD control problem. Stage one of the four-stage algorithm solves a reducedscale version of the LP problem in 3, yielding an intermittently optimal, or reduced-scale solution (ur), to the full-scale problem. The number of slice lip actuators on each slice lip and the number of sensor locations in the full array sensor are reduced by evenly removing rows and columns of the model, P, along with the corresponding constraints, based on a scale factor. The second stage of the four-stage algorithm is designed to obtain a full-scale feasible solution, u ˆ , to the LP in problem 3 using ur: p

min

u ˆ {u ˆ k|k,...,u ˆ k+m-1|k}

fobjective(u ˆ) )

[di+ + di-] ∑ i)1

(5)

s.t. ui - di+ + di- ) uri A û ˆ + s ) δumax {d+,d-} g 0 0 e s e 2δmax The variables d+ and d- are the auxiliary variables that correspond to the absolute value of the error between the reduced control move ur and the corresponding values from the feasible solution vector u ˆ . This problem ensures that none of the bending moment constraints are violated using the bending moment constraints. A û ˆ + s ) δumax. The third stage obtains a vertex solution uv using uˆ . Already accounting for the slice lip constraints in the second stage, the matrix A ˆ is partitioned into its active rows, B ˆ , and inactive rows, C ˆ , based on the solution of the second stage. The matrix V ˆ is basis for null (B ˆ ) and

is used to account for any additional rows of the bending moment constraints that may become active within the third stage. The vertex solution is computed as follows:

min

t{tk|k,...,tk+m-1|k}

fobjective(t) ) ξv

(6)

s.t. Rv+ + Rv- + sv1 - Ξ ) 0 Rv+ - Rv- + Mt ) b Nt + sv2 ) d {Rv+, Rv-, sv1} g 0 0 e sv2 e 2δmax b ) r - y˜ + Pu ˆ - Pu ˆ M ) PV ˜ N)C ˜V ˜ ˜u ˆ d ) δmax - C Once solved, uv can be found by uv ) V ˆt + u ˆ . As before, the variables Rv+ and Rv- are the auxiliary variables that correspond to absolute value of the error of the predicted error term based on the vertex solution. The variables sv1 and sv2 correspond to the slack variables for the bounds on the maximum predicted error and the bending moment constraints. Finally, t is a directionality vector that brings u ˆ to the vertex solution, uv. Finally, stage 4 solves the full-scale LP problem 3, using the vertex solution as a starting basis for the optimal solution. A summary of this algorithm is demonstrated in the flowchart in Figure 3. This flowchart shows the number of steps that the controller must go through to reach the ultimate goal of finding an optimal solution to the LP problem 3. Since CPLEX is designed to only solve one LP at a time, the code first initializes all of the objective function coefficients and left-hand side constraint matrixes to the Matlab workspace. At each controller calculation, the right-hand side constraint vector of the reduced LP is built and all of the arguments of the reduced LP are sent to CPLEX to be solved. This pattern continues for each of the other three stages until the LP problem 3 is solved. It should be noted that additional mathematics are necessary to calculate the basis for the null set of B ˆ . Although there exists a communication overhead in the CPLEX calls, this technique is still three to four times faster than using Matlab’s linear programming solver. Closed-loop simulation results at both sensor locationssyp and yssand the corresponding slice lip profiless up and ussare shown in the lower four plots of Figure 2. These are the typical closed-loop profiles for all of the results in this paper. The results are interpreted for three criteria throughout this paper: calculation time per control move and the two-directional variance measuressσMD and σCDsas defined in eqs 1 and 2. The controller calculation time and directional variance results corresponding to the lower two plots in Figure 2 are shown in Figure 4. The top plot of Figure 4 shows the calculation time of the four stages as well as the total calculation time while the lower plots show the σCD and σMD values, respectively. The LP-MPC algo-


Figure 3. Flowchart for the four-stage algorithm.

Figure 4. Controller calculation time and closed-loop performance results using the four-stage technique of section 2.

rithm is able to bring the basis weight to its steadystate values within 70 s of the disturbance measurement. The controller reduces the σCD value to the same levels as the QP-MPC result within 120 s of simulation. The controller does not, however, provide a similar level of σMD. This can be addressed by suitable design of an MD-controller. The profiles demonstrate that there are times that the secondary headbox behaves somewhat

erratically, attempting to minimize a variation within the profile that causes high-frequency ringing in the final basis weight CD profiles. These issues may be tuning related and are not as striking an issue as the slow calculation timesson average of 100 s per calculationsshown in Figure 4. 2.1. Modifications to the Four-Stage Approach. An analysis of Figure 4 shows that the vertex and fullsize LP problems are the most time-consuming steps of

2140


Figure 5. Controller calculation time and closed-loop performance results using variants of the four stage technique of section 2.

the four stage algorithm applied to this benchmark problem. Since this is also demonstrated in the original work on the four-stage algorithm,2 two solutions are now proposed to bypass the vertex LP problem 6. The first proposed solution to speed up the four stage algorithm is to instead use a three-stage approachs eliminating the vertex LP problem 6 and using the feasible LP problem 5 to warm start the full LP problem 3. The other proposed solution is to perform a two-stage algorithm using the feasible LP problem 5 directly as a control law. This assumes that the reduced LP problem 3 is near the optimal value of the full scale LP problem 3 prior to further optimization. In Figure 5 comparisons of the computation times and performance for the solutions proposed above are compared to the results for the original four-stage example. The three-stage algorithm yields the best tradeoff between acceptable performance and computational efficiency. This technique, however, is far from the 1 s computation times that are required for industrial implementation. The two-stage technique, however, is found to perform rather poorly although its computation time is quite fast. 3. Direct Solution of LP-MPC For problems of a large scale, the four-stage approach may not necessarily yield a more efficient solution to the LP-MPC CD control problem owing to communication bottlenecks between CPLEX and Matlab in the proposed architecture. The algorithm proposed in this section solves the full scale LP problem 3 directly at each sample time. The algorithm proposed here takes advantage of the efficiencies within the simplex algorithm of CPLEX version 7.1.13 The initialization of the controller code builds the matrixes and vectors describing the objective function, the left-hand side matrix of the constraint equalities and the variable bounds since

these do not change from sample to sample unless the model is revised. Figure 6 summarizes the implementation of this approach. This technique exploits the advantages of the CPLEX solver in ways that the four-stage algorithm cannot. First, the matrixes that are built within the initialization of the controller remain within the CPLEX memory space since there is no need to alter the matrixes from one sample time to the next. This is in contrast to the four-stage algorithm which stores the model matrixes in the Matlab workspace and communicates all four LP problems to CPLEX at each sample time. Also, the direct solution algorithm proposed here uses the solution from the previous sample time to warm start the full scale LP problem 3. Computation times and directional variance measures for this controller on the closed-loop example are shown in Figure 7. Since the objective of both the four-stage algorithm and the direct LP-MPC algorithm are the same, the performance measures of the two algorithms are equal. The direct LP-MPC algorithm, however, is at least 1-2 orders of magnitude faster in its calculation of the LP problem 3 at each sample time. The computation times for this algorithm are still 1 order of magnitude slower than those necessary for industrial implementation within the given sample time of 1 s. 4. Coincidence Points and LP-MPC The original LP-MPC research by Dave et al.2 uses a steady-state model. This poses a question with respect to the current work: Even with a faster sensor, is it necessary to minimize the error over the transient response? This section introduces the technique of coincidence points to answer this question. A coincidence point refers to the technique of setting the weights in the optimization problem such that only one or a few key times along the prediction horizon are weighted


Figure 6. Flowsheet for the direct LP-MPC algorithm.

state coincidence point, from now on referred to as only the coincidence point. In essence, this optimization reduces the number of rows and columns within the LP problem to be solved at each sample time. The LP problem to be solved at each sample time is a modification of the LP problem 3: m

min uk|k

fobjective(rP,∆u,u) ) ξ +

R∆u[∆ui+ + ∆ui-] ∑ i)1

(7)

s.t. Rp+ + Rp- + s1 - Ξ ) 0 -Pp∆u+ + Pp∆u- - Rp+ + Rp- ) y˜ - rp Figure 7. Controller calculation time and performance using only the direct LP-MPC algorithm of Section 3.

within the objective function of the optimization problem to be solved at each sample time. The first published uses of coincidence points are in the descriptions of the industrial MPC algorithms14 by both SetPoint’s (now owned by AspenTech) Multivariable Identification and Command (IDCOM-M) routines14,15 and Adersa Inc.’s Predictive Functional Controller (PFC),14 both of which cite Richalet’s Model Predictive Heuristic Control (MPHC) and IDCOM16 as a basis for their design. Gupta uses a technique in which only the predicted steadystate values are considered in the objective function.17 Bamberger and Isermann18 develop an on-line steadystate optimization algorithm with the static part of a nonlinear model and a higher order gradient technique. Marchenko and Rybashov19 seek optimal inputs corresponding to a steady-state solution of a dynamic linear model. Yousfi and Tournier20 adapt a MPC algorithm to solve dynamic and steady-state objectives within the same objective function. Finally, Muske21 uses linear programming to perform an economic optimization of steady-state target values then sends these values to QP-MPC for dynamic optimization. Although more than one coincidence point in time can be chosen for various reasons, Gupta’s technique explores the extreme case of only minimizing the steady

u - ∆u+ + ∆u- ) u ˜ A ˆ u + s2 ) δmax {R+, R-, s1} g 0 0 e {∆u+, ∆u-} e ∆umax 0 e s2 e 2δmax umin e u e umax The above LP is modified only in that the prediction model, Pp, and the corresponding Rp+ and Rp- are only of the size Ny corresponding to the coincidence point value instead of pNy of the original LP problem 3 based on minimizing values across the entire prediction horizon. The steady-state gain matrix, Pp, for the LPMPC CD control problems in this section is the part of the original gain matrix, P, corresponding to the value at k ) p ) 7. This coincidence point prediction model is equivalent to the model used in the original research regarding the four stage algorithm2 although the sample time of the true process is still 1s as compared to the much slower sample rates for scanning sensors, which range from 20 to 40 s per scan or longer.

2142


Figure 8. Controller calculation time and performance comparison of the coincidence point technique with the previous technique of optimizing over the full prediction horizon.

The coincidence point technique is tested here on the CD control problem using the four stage LP solution method described in section 2 and the direct LP-MPC method described in Section 3. The resulting computation times and performance measures for these formulations are compared to those using a full prediction horizon in Figure 8. Using the coincidence point technique, each algorithm yields a 2 orders of magnitude decrease in computation times without any degradation in the performance of the control action. In fact, using an objective function that weights the steady-state profile much more than the dynamic response yields better σMD values. This is likely due to the fact that the objective function minimizes one profile rather than CD errors over many profiles over the dynamic response. The direct LP-MPC method combined with the coincidence point technique is the first to yield both an acceptable variance reduction and a fast (0.1 s) computation time. 5. LDMC and QP-MPC Despite the efforts made above, these techniques all focus on nontraditional linear programming implementations of the MPC problem. In this section, a QP-MPC formulation and a QP motivated LP-MPC formulation are developed and tested on the benchmark problem. With the exception of the research by Biegler,3 and earlier efforts by Backstro¨m,4 QP-MPC applications to the CD control problem were not implemented due to lack in solver efficiency. The CD control in their research still use scanning sensor technologies that require less stringent sample times. Also, the example that Biegler and co-workers explore is a single headbox machine.3 While Baskstro¨m and co-workers4 address a multiple headbox problem, they only use one headbox for control purposes.4 As is evident from the above sections, optimal control utilizing all of the inputs and outputs from full array measurements provide a challenge that is not addressed on in these works.

The technique developed in this section addresses the issue of QP-MPC implementation as defined in QP problem 4. The idea here is to use a modified LP-MPC objective function that minimizes an objective that is close in form to the MPC optimization problem based on a quadratic criterion to warm start a standard QPMPC problem. The LP-MPC problem uses the same modified objective function used in Saffer et al.:9 m

min fobjective(r, uk+n|k) ) uk

{Rj+ + Rj-} ∑ j)1

(8)

s.t. Rj+ - Rj- - (PTP + RuTRu)uk+n|k ) PT(r - y˜ ) A ˆ u + s1 ) δmax {R+, R-, s1} g 0 -∆umax e ∆u e ∆umax 0 e s1 e 2δmax umin e u e umax This objective function is first used by Morshedi et al. in their proposed LDMC routine.8 Their motivation is to use a move suppression factor, PT, to reduce the computational demands of the optimization routine. This proposal is motivated from the idea that the objective function of the LP problem 8 minimizes the matrix derivative of the QP problem 4 with respect to the manipulated variables subject to the process constraints. While QP theory states that the constrained minimal slope obtained by eq 8 of the quadratic objective function eq 4 is not necessarily the minimum of said quadratic objective, eq 4, the convexity of the QP-MPC problem suggests that the solution to the LP Problem 8 might be close to the QP problem 4 minimum.


Figure 9. Comparison of controller calculation times and performance for the controllers testing QP-MPC with and without the warmstart of a LDMC algorithm as well as the LDMC algorithm itself.

The closed-loop examples for this section demonstrate three alternatives. First, QP-MPC using the QP problem 4 is performed to provide a basis for comparison. Second, LP-MPC using the proposed modified LP problem 8 is tested alone using an implementation that is the same as that in Figure 6 with the LP problem 8 replacing LP problem 3. Finally, the proposed technique of using the LP problem 8 to warm start the QP problem 4. Just as in section 3, the modified LP problem 8 in this section takes advantage of the efficiencies within the simplex algorithm of CPLEX version 7.1.13 The initialization of the controller code again builds the matrixes and vectors describing the objective function, the left-hand side matrix of the constraint equalities and the variable bounds since these do not change from sample to sample unless the model is re-identified. Figure 9 shows the computation time, σCD, and σMD performance values for the three controllers. While the LDMC-warm-started-QP is not computed within the sample time of 1 s, it does reduce the computation time to half of the QP-MPC without warm starting. Although it is not a rigorous QP-MPC formulation, the LDMC algorithm performs as well as the QP formulations for this closed-loop example while providing a solution that has the fastest efficiency of any previously tested methods. Rigorous testing of the modified LP problem 8 would be necessary, however, prior to implementation of this in an industrial setting. 6. Summary A number of computationally efficient linear programming solutions for LP-MPC and QP-MPC are addressed in this paper in the context of the CD control of a dualheadbox benchmark paper machine simulation using fast full array sensors. While not all algorithms provide a sufficiently fast solution or a desirable performance, two of the combinations of techniques that are exploreds the technique that combines the coincidence point

approach with the direct solution approach as well as the LDMC algorithmsdemonstrate both efficient computation times and acceptable performance. The extension of the four-stage algorithm developed by Dave et al.2 in its original form with a dynamic prediction model is found to provide solutions that are acceptable in performance, but unacceptable in efficiency. Bypassing the most time-consuming stage of the four-stage approach still yields unacceptable computational times and/or performance. A simpler algorithm that utilizes a computationally efficient commercial LP solver, CPLEX, is found to perform similarly to the four-stage approach while producing the results at times that are1-2 orders of magnitude more efficient. The approach does not always compute the solution within the 1 s sample rate, but does so more than 50% of the time. The use of coincidence points is explored to reduce the computational burden of the four-stage algorithm and the direct LP-MPC algorithm. While the problem size is significantly reduced, the performance of the resulting controllers does not degrade to any significant degree. The LDMC approach of Morshedi et al.8 is also found to yield an efficient and high performance solution to the CD control problem. This algorithm actually performs more comparable to the QP-MPC approach than the algorithms that strictly solve for a 1/∞-norm, either with or without using coincidence points. Acknowledgment We acknowledge the financial support of the University of Delaware and the technical and financial support of Weyerhaeuser. Literature Cited (1) Dave, P.; Willig, D.; Kudva, G.; Pekny, J. F.; Doyle, F. J., III. LP methods in MPC of large scale systems - Application to paper machine CD control. AIChE J. 1997, 43 (4), 1016-1031.

2144


(2) Dave, P.; Doyle, F. J., III; Pekny, J. F. Customization strategies for the solution of linear programming problems arising from large scale model predictive control of a paper machine. J. Proc. Cont. 1999, 9 (5), 385-396. (3) Bartlett, R. A.; Biegler, L. T.; Backstrom, J.; Gopal, V. Quadratic programming algorithims for large-scale model predictive control. J. Proc. Cont. 2002, 12 (7), 775-795. (4) Backstro¨m, J. U.; Gheorghe, C.; Stewart, G. E.; Vyse, R. N. Constrained model predictive control for cross directional multiarray processes. Pulp Paper Canada 2001, 102 (5), 33-36. (5) Rigopoulos, A.; Arkun, Y.; Kayihan, F. Model predictive control of CD profiles in sheet forming processes using full profile disturbance models identified by adaptive PCA. Proc. American Control Conf. Albuquerque, NM, 1997, pp 1468-1472. (6) Van Antwerp, J. G.; Braatz, R. D. Model predictive control of large scale processes. J. Proc. Cont. 2000, 10 (1), 1-9. (7) TAPPI. Calculation and partitioning of variance using paper machine scanning sensor measurements; Technical Report TIP 1101-01; Process Control Committee of the Process Control, Electical and Information Division, TAPPI: Atlanta, GA, 1997. (8) Morshedi, A.; Cutler, C.; Skrovanek, T. Optimal solution of dynamic matrix control with linear programming techniques (LDMC). In Proc. American Control Conf. Boston, MA, 1985, pp 199-208. (9) Saffer, D. R., II; Doyle, F. J., III; Rigopoulos, A.; Wisnewski, P. MPC study for a dual headbox CD control problem. Pulp Paper Canada 2001, 102 (12), 97-101. (10) Featherstone, A. P.; Van Antwerp, J. G.; Braatz, R. D. Identification and Control of Sheet and Film Processes; Springer: London, 2000. (11) Garcıá, C. E.; Prett, D. M.; Morari, M. Model predictive control: Theory and practice - a survey. Automatica 1989, 25 (3), 335-348.

(12) Ogunnaike, B. A.; Ray, W. Process Dynamics, Modeling, and Control; Oxford University Press: New York, 1994. (13) CPLEX. ILOG CPLEX 7.1 User’s Manual, 2001. (14) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control Eng. Prac. 2003, 11 (7), 733764. (15) Grosdidier, P.; Froisy, B.; Hammann, M. The IDCOM-M controller. Proc. IFAC Workshop Model Based Process Control Atlanta, GA, 1988, pp 31-36. (16) Richalet, J.; Rault, A.; Testud, J. L.; Papon, J. Model predictive heuristic control: Applications to industrial processes. Automatica 1978, 14 (5), 413-428. (17) Gupta, Y. P. A simplified predictive control approach for handling constraints through linear programming. Computers Ind. 1993, 21 (3), 255-265. (18) Bamberger, W.; Isermann, R. Adaptive on-line optimization of the steady-state behavior of slow dynamic processes with process computers. Digital Computer Applications to Process Control, North Holland, 1977. (19) Marchenko, M. I.; Rybashov, M. V. Optimization of the quasistationary mode in linear systems. Autom. Remote Control 1988, 48 (12), 1605-1614. (20) Yousfi, C.; Tournier, R. Steady state optimization indise model predictive control. Proc. Am. Control Conf. Boston, MA, 1991, pp 1866-1870. (21) Muske, K. R. Steady-state target optimization in linear model predictive control. In Proc. Am. Control Conf. Albuquerque, NM, 1997, pp 3597-3601.

Received for review August 3, 2004 Revised manuscript received December 30, 2004 Accepted January 14, 2005 IE049303I

Efficient Linear Programming Formulations of Model Predictive Control

Recommend Documents