Offset-Free Model Predictive Control with Explicit ... - ACS Publications

Matt Wallace, Steven Spielberg Pon Kumar, and Prashant Mhaskar. Department of Chemical Engineering, McMaster University, Hamilton, Ontario L8S 4L7, ...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/IECR

Offset-Free Model Predictive Control with Explicit Performance Specification Matt Wallace, Steven Spielberg Pon Kumar, and Prashant Mhaskar* Department of Chemical Engineering, McMaster University, Hamilton, Ontario L8S 4L7, Canada ABSTRACT: This work considers the problem of designing a control law to yield a desired closed-loop response subject to input constraints and plant-model mismatch. To this end, a two-tiered model predictive controller is designed that first computes the best response of a certain desired form, and then computes the control action to achieve the prescribed closed−loop response. The key idea is first illustrated in the absence of plant−model mismatch for a prescribed first- and second-order response of the controlled variables, and shown using a simulation example. Subsequently, a formulation is presented to handle plant−model mismatch. Finally, the proposed formulation is applied to a chemical process example.

1. INTRODUCTION The operation of chemical processes must combat many complex challenges. The challenges arise due to a myriad of issues, including nonlinearity, disturbances, and multivariable interactions, as well as constraints that must be satisfied to safely and optimally meet product specifications. One control method well-suited to handling these challenges is model predictive control (MPC), where the control action is computed by solving an optimization problem that minimizes a given cost function (often used to represent the regulation objective), while accounting for the process dynamics (through a prediction model), and accommodating state and input constraints. The industrial implementation of MPC has contributed to research toward and benefited from a better characterization of the rigorous stability and performance properties of the MPC. In one direction, this has included the development of Lyapunov-based MPC formulations1 that enable an explicit characterization of the region from which stability can be guaranteed in the presence of input2 and state3 constraints. These approaches also allow uncertainty handling through the robust control approach4 and via accounting for the stochastic nature of the disturbances.5 In contrast to the robust control approach, where the control action is computed to counteract the worst-case effect of the uncertainty, the offset-free MPC approach6 has focused on the inclusion of so-called “augmented states” (akin to the integral mode in classical proportional integral controllers) to estimate and counteract the uncertainty. The offset-free MPC approaches have been implemented using both linear7,8 and nonlinear models, as well as in conjunction with robust control approaches.9 In the aforementioned MPC approaches, the cost function is used as a means to specify and enforce the regulation objective. Thus, even though the control action is computed to be optimal, because the weights in the objective function are initially selected to satisfy the regulation objective, they must be further tuned to try and match a desired closed-loop behavior. The recognition of process economics as the overriding objective has fostered the development of economic MPC formulations, with the key feature that the MPC internally determines the appropriate setpoint at which the process © 2016 American Chemical Society

variables should be kept. Results in this direction have focused on ensuring stability in an appropriate sense for linear (see, e.g., ref 10) and nonlinear (see, e.g., refs 11 and 12) systems. Notwithstanding the plethora of results on stability analysis for standard or economic MPC formulations, relatively little attention has been given to incorporating desired closed-loop behavior explicitly in the MPC formulation, especially for multiple-input, multiple-output (MIMO) systems. The important issue of explicit specification of the control performance has been addressed for single-input, single-output (SISO) classical control structures. In particular, the notion of loop shaping13 enables determining the control law to achieve a certain desired closed-loop transfer function (via closed-loop system analysis in the frequency domain). Generally, there is no guarantee that the resultant control law will be consistent with a proportional-integral-derivative (PID) controller structure. Furthermore, the approach does not generalize readily to MIMO systems with constraints. In one effort, MIMO systems are handled through the so-called “funnel control approach”, where the controller is designed (e.g., see ref 14) to bound the output error of specific control variables within a “funneling” function, which dictate the controller gain. However, the results do not explicitly account for the presence of input constraints and plant−model mismatch. Motivated by the above considerations, this work considers the problem of explicit specification of the performance in the control design for linear time-invariant systems subject to input constraints and plant−model mismatch. To this end, we first present the system description in section 2, followed by a brief review of offset-free MPC. We then present in section 3 a twotiered MPC framework that enables determining the best achievable performance subject to constraints, and then implements it in the closed loop. The design is illustrated using an example in section 3.4. We next address the problem of handling plant−model mismatch within the proposed Received: Revised: Accepted: Published: 995

October 8, 2015 January 6, 2016 January 7, 2016 January 7, 2016 DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research

rate of change of input moves versus the state deviation through the relative values of the Q and R matrices. While these characteristics go significantly further in terms of the ability to specify the performance objectives, the closed-loop behavior is often judged (and controllers fine-tuned) against criteria such as “a smooth but fast first-order response”, etc., and in the existing MPC formulations, there is no mechanism to directly incorporate these considerations in the control calculations. One way of handling plant−model mismatch in MPC formulations with linear models is the offset-free approach,6,7,17 where additional “integrating” states, θ (assumed constant over the prediction horizon), are incorporated in the model. In particular, the model is augmented as follows:

framework in section 4, followed by an illustrative example in section 4.2. In section 5, the MPC is implemented on a chemical process example. We finally present conclusions in section 6.

2. PRELIMINARIES In this section, we present the system description, followed by a review of linear offset-free MPC approaches. 2.1. System Description. We consider linear systems of the form (1) ẋ = Ax + Bu n m where x ∈ R , u ∈ R , and A and B are constant matrices. We assume that the system is controllable, and the input is subject to the constraints, umin ≤ u ≤ umax, where umin, umax ∈ Rm are the bounds on the manipulated inputs. In keeping with the discrete implementation of MPC with a sampling time Δt, u is piecewise constant and defined over a sampling instance l as u (t ) = u (l )

x k = A*x k − 1 + B*uk − 1 + Gθθk − 1 θk = θk − 1

The augmented model can be represented as ⎡ x̂l ⎤ x̃l = ⎢ ⎥ = AΣx̃l − 1 + BΣul − 1 ⎢⎣ θl̂ ⎥⎦

lΔt ≤ t < (l + 1)Δt

2.2. Linear Offset-Free MPC. In this section, a representative linear offset-free MPC formulation is described, where, for the sake of simplicity, and without loss of generality, x = 0 is the prescribed setpoint. In particular, the control action at a sampling instance l is computed by solving the following optimization problem: P

2

min ∑ xk̂ û

k=1

(3)

(4)

where ⎡ A* Gθ ⎤ AΣ = ⎢ ⎥ ⎣0 I⎦

2

+

and

Δuk̂ − 1

Q

R

⎡ B* ⎤ BΣ = ⎢ ⎥ ⎣0⎦

(2)

subject to

The augmented states act like an integrator and, under certain assumptions, can eliminate offset at steady state. The current value of the augmented states can be estimated using a variety of approaches, such as a moving-horizon-estimator (MHE) or a Luenberger observer. As an illustration, when using a Luenberger observer, the state estimation equations take the following form:

u min ≤ uk̂ − 1 ≤ u max x̂k = A*xk̂ − 1 + B*uk̂ − 1 with x̂ 0 = x l Stability Constraint

where at sampling instance l, the implemented input to the process is given as ul = û0 + unom, where unom is the nominal input, x̂ denotes the prediction vector for all states, i = 1, ..., n, Δu(k) = û(k) − û(k − 1), k > 1 denotes the rate of change of successive inputs for a candidate input profile denoted by û. To capture the difference between the first value of the candidate input and current input, Δu(0) ≡ û(0) + unom − u(l − 1). Q and R are corresponding penalty matrices on state deviation and rate of change of the input, where the notation ∥·∥2Q is the weighted norm, defined by ∥x∥2Q = xTQx. P is the prediction horizon, and A* and B* are the discretized versions of the matrices A and B, respectively, with a sampling time Δt. The stability constraint is representative of various formulations that guarantee stability by incorporating some form of a stability constraint (see, e.g., refs 15 and 16). Remark 1. As stated earlier, the problem of stability is handled either by explicitly incorporating a stability constraint, or by augmenting a terminal penalty in the objective function to ensure that the optimal solution is also stabilizing. In either case, the closed-loop implementation is optimal, primarily in the sense of ensuring closed-loop stability. The other mechanisms through which performance considerations are incorporated are, for instance, the ability to specify known constraints on the state and input variables to satisfy safety considerations, or through specifying the relative importance of

⎡ x̂l ⎤ ⎡ x̂l − 1 ⎤ ⎡L ⎤ ⎥ + BΣul − 1 + ⎢ x ⎥[x l − 1 − x̂l − 1] x̃l = ⎢ ⎥ = AΣ⎢ ⎢⎣ θl̂ ⎥⎦ ⎢⎣ θl̂ − 1⎥⎦ ⎣ Lθ ⎦ (5)

⎡L ⎤ where ⎢ L x ⎥ is the observer gain matrix. The choice of Gθ and ⎣ θ⎦ the poles of the observer gain matrix are tunable parameters. The augmented state estimates are used to initialize the following MPC optimization problem: 2

P

min û

∑ k=1

xk̂

2

+

Δuk̂ − 1

Q

R

(6)

subject to u min ≤ uk̂ − 1 ≤ u max ∼x̂ = A ∼x̂ Σ k − 1 + BΣ uk̂ − 1 k

with

∼x̂ = x̃ 0 l where 996

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research ∼x̂ = [x ̂ , x ̂ , ... x ̂ , θ ̂ , θ ̂ , ... θ ̂ ]′ k 1k 2k nk 1k 2k nk Stability Constraint

The above formulation(s) considers the case when the desired state setpoint(s) is defined at the origin (thus, the steady-state/nominal value of the input also is zero). In particular, the states, and the inputs (as well as the constraints) are defined in terms of deviation variables. For the purpose of implementation, for any given setpoint, the corresponding steady-state value of the input must be computed to implement the above formulation. In particular, for a given setpoint xSP, the corresponding steady-state input (for the case where B* ∈ Rn×n) is computed as follows: uss, xSP = B*−1[xSP − A*xSP]

Figure 2. Proposed MPC approach after initialization.

3.1. MPC Formulation: Tier 2. In the proposed MPC formulation, the control action is computed by solving the following optimization problem: P

min ∑ x̂k − x̅ SP û

(7)

2 Q

(8)

k=1

where uss,xSP = [uss,x1,SP, ..., uss,xn,SP] and xSP = [x1,SP, ..., xn, SP]. If B* is not square and the number of inputs is more than the number of states, that simply allows more degrees of freedom in achieving a desired behavior. On the other hand, with fewer inputs than states, one must recognize that not all the states can be made to evolve in a prescribed manner. Thus, either the prescribed behavior must be described cognizant of this fact, or recognized that the prescribed behavior can only be achieved in some best fit sense. The details for this scenario remain outside the scope of the present work.

subject to u min ≤ û k − 1 ≤ u max

where xk̂ = A*x̂k − 1 + B*û k − 1

with x̂0 = xi, where x̅ i,SPk = Desired behavior

where, at sampling instant l, the implemented input to the process is ul = û0 + unom, x̂i denotes the prediction for state i, xl is the measured states used to initialize the state predictions, and xS̅ P is the desired behavior of the controlled variables. The desired behavior is directly incorporated into the MPC objective function by only penalizing the deviation of the predicted states in the objective function. Remark 2. Note that the objective function in the proposed formulation does not include a rate of change of input penalty term, nor does it require a “stability constraint”. These requirements are taken care of when the best achievable desired performance is computed. For instance, for a desired first-order response, the maximum allowable input movement between two successive control calculations can be incorporated as a hard constraint when determining the fastest firstorder response achievable. Another key point to recognize is that the proposed MPC is “driven” by the desired behavior, and to specifically enable the implementation of readily acceptable, and intuitive, closed-loop behavior (such as a stable first-order response, a stable second-order response, etc.). Thus, the stability of the closed-loop system is derived, not from stability constraints as in the traditional MPC formulations, but from the nature of the desired behavior being asked of the MPC (see section 3.2 for an illustration). 3.2. Achieving a First-Order Response. We first consider the case where each controlled state is required to approach the setpoint in a first-order fashion:

3. MODEL PREDICTIVE CONTROL (MPC) FORMULATION WITH PERFORMANCE SPECIFICATION In this section, we first consider the case of a perfect model to illustrate the key idea of explicit incorporation of the control performance specification, where the desired performance is described in terms of a desired behavior of the state evolution. The proposed MPC formulation recognizes that a desired behavior is often expressed in a general functional form (such as first-order response/second-order response with minimum overshoot ratio, first order with time delay, etc.). However, the system dynamics and presence of constraints limit the best (e.g., fastest first-order response/fastest second-order response with minimum overshoot ratio, etc.) achievable performance. Therefore, a two-tiered framework is constructed, where, in the first tier (referred to as the performance optimizer block), the best achievable performance subject to constraints is computed at controller initialization, and then the second tier implements that specified control performance. The two-tiered control structure is schematically presented in Figures 1 and 2. The

⎡ ⎛ −Δt ⎞⎤ ⎛ −Δt ⎞ xi̅ ,SPk = ⎢1 − exp⎜ ⎟⎥xi ,SPk + ⎜exp ⎟xi ,SP ⎢⎣ τi ⎠ ̅ k −1 ⎝ τi ⎠⎥⎦ ⎝

Figure 1. Proposed MPC approach at initialization.

(9)

where xi,SP is the desired setpoint for each state i, with τi denoting the time constant. The key idea in the proposed framework is to require the MPC to force the setpoints to be tracked in this fashion, with the performance dictated by how quickly the setpoints can be achieved, which, in turn, are impacted by the values of τi. For a given linear system subject to constraints, the “best achievable” values of τi can be computed.

MPC described by eq 8 is the second tier in the proposed structure, where a specific desired closed-loop behavior computed in Tier 1 is implemented. Specific instances of the Tier 1 formulations for first-order and second-order desired behavior are presented next. 997

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research In particular, the control action required to achieve a first-order response is given by uk − 1 = B*−1[ΓFO − A*]x k − 1

function can be made to decay by using the available control action. By virtue of fulfilling this requirement (i.e., finding a feasible solution to this optimization problem), stability of the closed-loop system follows. In other words, consider a Lyapunov function of the form V(x) = x12 + x22 + ... xn2. By virtue of how τ is computed in the performance optimizer (satisfying input constraints), it is known that a set of input moves exist for which the value of the objective function is exactly zero, or that the state exactly tracks the prescribed firstorder response, with a finite τ. Thus, upon implementing that input, each state will approach their respective desired setpoints in a first-order fashion. This leads to successive decay in the value of the Lyapunov function, leading to stabilization. Remark 5. The τ values are computed by minimizing the norm of the τ vector. If the τ values were recomputed at all subsequent times, it is possible that the subsequently computed τ values give a lower norm than at the previous time(s), with some of the τi values potentially increasing. To guarantee successive nonincreases of each τi, a constraint of the following form could be added to the optimization problem of eq 11:

(10)

where ⎡ e−Δt / τ1 0 ⎢ −Δt / τ2 ⎢ 0 e ΓFO = ⎢ ⋮ ⎢ ⋮ ⎢⎣ 0 0

⎤ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ e−Δt / τn ⎦ ⋯

0

Thus, the smallest achievable time constants can be determined by solving the following optimization problem at controller initialization: min τ τ

(11)



subject to û k − 1 = B*−1[ΓFO − A*]x̂k − 1 u min ≤ û k − 1 ≤ u max x̂k = A*x̂k − 1 + B*û k − 1

∀ i = 1 , ..., n

τik ≤ τik −1,min

(12)

Remark 6. After the optimal τ values are computed, the second-tier MPC optimization problem must be implemented with a sample time of at least some fraction of the fastest optimal τ value. This is important to ensure that the controller can provide the control action that will satisfy the optimal performance specifications for each state. As a result, a guideline of Δt = 1/3 min(τopt,i) (i = 1, ..., n) is recommended (and utilized for the simulations). 3.3. Explicit Tuning Approach for an Underdamped Second-Order Specification. As another example of an acceptable response, an underdamped second-order response is considered. In particular, we consider the scenario where each state is desired to evolve in accordance with the following differential equation:

∀ k = 1 , ..., P2

x̂ 0 = x l 0≤τ

where Qτ is the matrix that specifies the relative importance of one state versus the other (an alternative approach could be, for instance, to minimize the largest of the resultant τ values) and P2 is the prediction horizon. These best values of the time constants are then used to compute the specific profile that is implemented by the second tier of the controller described by eq 8. Note that, in the absence of plant−model mismatch, it is sufficient to solve this optimization problem only once. The prediction horizon must be of an adequate length to ensure that the slowest desired performance specification state is able to converge to the desired setpoint. Note that if the horizon is not sufficiently long, then it will not be guaranteed that it is possible to drive the states to the desired setpoint with the computed τ. One way to ensure this is to use a multiple of the slowest time constant, for example, P2 = 5 × τmax, where τmax corresponds to the slowest open-loop time constant. The idea being to exploit the fact that in a closed loop, it should be possible to drive all the states to the setpoint within a time that is smaller than the largest open-loop time constant. Remark 3. Note that the key objective of the first tier is to determine if the desired behavior is achievable and then determine solutions with improvements. Given that the optimization problem defined in eq 10 is not convex, it is understood that the solution can, at best, be guaranteed to be a local minimum. However, the key is to recognize that the solution would be a feasible one for the closed-loop system, and thus is able to provide the desired closed-loop behavior. Furthermore, advances in the field of nonconvex optimization (see, e.g., refs 18−21) can be utilized to compute improved solutions to the optimization problem. Remark 4. Note that the manner in which the closed-loop performance is prescribed (first-order response), the first tier can be considered to be an optimization problem that searches for a quadratic Lyapunov function such that the Lyapunov

τi 2

∂ 2xi ∂t

2

+ 2ζiτi

∂xi + xi = 0 ∂t

(13)

where τi is the time constant and ζi represents the damping factors (note that 0 < ζi < 1). Using the variable transformation, zi = dxi/dt, eq 13 can be transformed to the following firstorder representation: ⎡ dxi ⎤ ⎢ ⎥ x ⎢ dt ⎥ = W ⎡ i ⎤ i⎢ z ⎥ ⎢ d zi ⎥ ⎣ i⎦ ⎢ ⎥ ⎣ dt ⎦

(14)

where ⎡ 0 1 ⎤ ⎥ ⎢ Wi = ⎢ 1 2ζτ ⎥ ⎢⎣− τ 2 − τ 2 ⎥⎦

and discretizing eq 14 yields ⎡ xik ⎤ ⎡ xik −1 ⎤ ⎢ ⎥ = W *i ⎢ ⎥ ⎣ z ik ⎦ ⎣ z ik −1 ⎦

(15)

where 998

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research ⎡ w* w* ⎤ 1,1 1,2 ⎥ * Wi = ⎢ ⎢ w* w* ⎥ ⎣ 2,1 2,2 ⎦

order specification, the second-tier MPC optimization problem is implemented with a sample time of Δt = 1/3 min(τopt,i) (i = 1, ..., n). Remark 7. The main difference between the proposed formulation and that of a traditional MPC is that the performance specification is explicitly included in the objective function. The performance objective can be that the controlled variable approach the setpoint in a first-order or second-order fashion (as in sections 3.2 and 3.3), or any other form. The key is that, given a particular form, the proposed formulation allows one to determine the “best” performance that is achievable, subject to the constraints, and then implements that. In the present formulation, the performance specification is described in terms of the controlled output behavior. If it is desired that the inputs also behave in a specific fashion (beside the satisfaction of constraints, which have already been taken into account), then these could also be directly incorporated in the objective function. Remark 8. Note that the proposed MPC is different, both in concept and implementation from trajectory tracking or other MPC approaches where a filtered version of the setpoint is provided to the MPC. In particular, in trajectory tracking implementations, the MPC is required to track a time-varying setpoint, instead of a constant one, but using the same control design. On the other hand, in other approaches, an arbitrarily filtered trajectory (to the setpoint) is provided to the MPC to track, which it again tracks using the traditional control structures. In contrast, in the proposed framework, the control design is given a setpoint and a desired closed-loop behavior, and the formulation determines the best closed-loop implementation and executes it. 3.4. Illustrative Simulation Results. To illustrate the proposed approach, a three-input, three-output linear dynamical system was used (see Table 1 for the system matrices). For the purpose of comparison, two standard linear MPC (eq 2) implementations are also shown (with different penalties on the rate of change of input movement) in Figure 3. This is done to illustrate the fact that while the speed of the response can be altered by tuning parameters such as Q and R in standard MPC

For a particular xi0 and a choice of z10 (in the simulation example, we choose z10 = 0; note that this value becomes nonzero over the duration of the specification of the desired behavior), xik and zik can be computed. For simplification, we define the following matrix (note that the matrix entries are dependent on the choice of τi, ζi, xi0, and zi0): ⎡ x1 0 ⎢ k ⎢ 0 x 2k ΓSO = ⎢ ⎢⋮ ⋮ ⎢0 0 ⎣

0 ⎤ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⋯ xnk ⎥⎦ ⋯

(16)

Then, by equating eq 16 with the discrete representation of eq 1, it can be shown that the control action required to achieve a second-order response is given by uk − 1 = B*−1[ΓSO − A*x k − 1]

(17)

As stated earlier, a given second-order response is achievable only if the input, as computed by the above equation, satisfies the constraints. In addition, to ensure that the controller provides a minimal oscillatory response, terms were included in the objective function, as well as constraints. In particular, a term penalizing the difference between each ζi value and that of their theoretical maximum (unity) was introduced in the objective function and peak overshoot constraints for each state’s response (POi) were taken into account. With these considerations, the smallest achievable time constants and within the allowable damping factor can be determined by solving the following optimization problem at controller initialization: n

min ∑ τi τi , ζi

i=1

+

(1 − ζi)

Q

Z

Table 1. System and Filter Characteristics

subject to û k − 1 = B*−1[ΓSO − A*x̂k − 1]

parameter

u min ≤ û k − 1 ≤ u max

A

x̂k = A*x̂k − 1 + B*û k − 1

∀ k = 1 , ..., P2

x̂ 0 = x l

B

POmin ≤ PO ≤ 1

x(0) τ1st,opt (s)

∀ i = 1 , ..., n

0 ≤ τi 0 ≤ ζi < 1

∀ i = 1 , ..., n

POi =

umin umax POmin

[1.51 × 10−1 1.68 × 10−1 1.18 × 10−1] [ 5.33 × 10−2 5.96 × 10−2 2.42 × 10−2 ] [8.26 × 10−1 8.26 × 10−1 9.97 × 10−1] [− 100 − 110 − 120] [150 140 130] [1.00 × 10−2 1.00 × 10−2 1.00 × 10−2 ]

POmax Δt (s) Rlo Δu Rhi Δu τmax (s)

[1.00 1.00 1.00] 5.00 × 10−4 diag{1, 1 × 10−1, 1} diag{100, 10, 100} 4

τ2nd,opt (s) ζ2nd,opt

where −ζiπ 1 − ζi 2

∀ i = 1 , ..., n

PO = [ PO1 ⋯ POn ]

and Qτ and Z are weighting matrices of appropriate dimensions, and P2 is the prediction horizon, of length 5 × τmax (where τmax is the dominant open-loop time constant). Identical to the first999

matrix/value

⎡− 6 − 1 1⎤ ⎥ ⎢ ⎢ 3 −5 3⎥ ⎢⎣ 4 1 − 2 ⎥⎦ ⎡− 2 − 1 − 5 ⎤ ⎢ ⎥ ⎢ 4 2 7⎥ ⎢⎣− 2 − 3 3⎥⎦ [− 1.73 3.84 3.46]

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research

Figure 3. Closed-loop input and state profiles.

4. OFFSET-FREE MPC FORMULATION Having illustrated the key idea in the proposed approach, we now present a formulation designed to handle plant−model mismatch. To this end, the model is augmented with states as in the standard offset-free approach described in section 2.2. At a sampling instance, l, an observer is first used to generate estimates of the states (for the purpose of illustration, we use a simple Luenberger design as below):

formulations, a rigorous procedure to select them, to provide a desired nature of the response, simply does not exist. The color vector defined by = [MPCτ, SP, MPClo Δu, MPChi Δu, MPCτ,ζ], represents the closed-loop profiles of the first-order performance specification-based MPC, the state setpoints, standard MPC (eq 2) with low Δu objective function weighting, standard MPC with high Δu objective function weighting, and, finally, the second-order performance specification-based MPC, respectively. As can be seen in Figure 3, the proposed MPC provides a first-order response for each of the states. In comparison, the standard linear MPC approaches are able to have certain states achieve faster convergence (x2), albeit with undesired behavior of x1. This behavior is present through a range of tuning settings for the standard linear MPC designs (both high and low Δu weights in the objective function) and underscores the inherent inability to specify a desired behavior in existing MPC approaches. Thus, while it is possible to achieve fast and slow response by changing the relative values of the tuning parameters in the standard MPC formulations, there does not exist a way to ensure that the behavior is acceptable for all the states. In contrast, the proposed MPC formulation implements a control action that ensures that the response is smooth and, yet, the fastest-possible first-order response. The other key point is that if a second-order response is desirable/allowable, the proposed approach allows readily accommodating such a behavior. Following the design procedure outlined in section 3.3, the proposed approach is thus able to achieve just as fast a response as the standard MPC, albeit with the overshoot within the acceptable limit specified in the controller design (see Figure 4).

⎡ x̂l ⎤ ⎡ x̂l − 1 ⎤ ⎥ + B Σ ul − 1 + x̃l = ⎢ ⎥ = AΣ⎢ ⎢⎣ θl̂ ⎥⎦ ⎢⎣ θl̂ − 1⎥⎦

⎡ Lx ⎤ ⎢ ⎥[x l − 1 − x̂l − 1] ⎣ Lθ ⎦ (18)

The state estimates are then used as initial conditions, and the control action computed by solving the following MPC optimization problem: P

2

min ∑ xî k − xi̅ ,SPk û k=1

q

subject to u min ≤ û k − 1 ≤ u max ∼x̂ = A ∼x̂ Σ k − 1 + BΣ û k − 1 k

with x∼0̂ = xl̃ ∼x̂ = [x ̂ , x ̂ , ..., x ̂ , θ ̂ , θ ̂ , ..., θ ̂ ] k 1k 2k nk 1k 2k nk ⎡ ⎛ Δt ⎞⎤ ⎛ Δt ⎞ xi̅ ,SPk = ⎢1 − exp⎜ − ⎟⎥xi ,SPk + ⎜exp − ⎟xi ,SP ⎢⎣ τi ⎠⎥⎦ τi ⎠ ̅ k −1 ⎝ ⎝

where, at sampling instant l, the implemented input to the process is ul = û0 + unom, x̂i denotes the prediction for state i, xl̃ is the state estimate (obtained using eq 18), and xi̅ ,SP is the desired trajectory of the control variables to reach the setpoint xi,SP (for the sake of illustration, we consider a first-order specification) with a prediction horizon P. For each of i = 1, ..., n states, individual time constants (τi) dictate the prescribed speed of the response. 4.1. Specifying a Closed-Loop Trajectory. Following a similar idea as that in section 3.2, the control action required to achieve a particular response can be computed (we illustrate

Figure 4. Magnified x1 closed-loop profiles. 1000

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research using a desired first-order response). However, the key difference is that, in computing the required control action, the current best estimate of the plant dynamics is utilized, which is described by x̂l = A*x̂l − 1 + B*ul − 1 + Gθ θl̂ − 1

aggressive, and minimizing the negative impact of plant−model mismatch. Remark 10. The presentation of the offset-free MPC formulation is a way of illustrating how the main idea behind the approach (that of providing the best prescribed desired behavior) can be integrated with existing MPC approaches. Future work will consider the problem of generalizing the present formulation to explicitly address other important issues, such as handling nonlinearity and the lack of availability of state measurements. 4.2. Illustrative Simulation Results. To illustrate the offset-free approach, we again consider a three-input, threeoutput linear dynamical system, where the model matrices are different from the plant (see Table 2). In designing the offset-

(19)

The control action that would result in a first-order response is given by ul − 1 = B*−1([ΓFO − A*]x̂l − 1 − Gθ θl̂ − 1)

(20)

Thus, the best achievable response is computed in Tier I by solving the following optimization problem: min τ

(21)

τ

Table 2. Plant and Model Matrices

subject to uk̂ − 1 = B*−1([ΓFO − A*]xk̂ − 1 − Gθ θk̂ − 1)

parameter

u min ≤ uk̂ − 1 ≤ u max

Bvary, model

xk̂ = A*xk̂ − 1 + B*uk̂ − 1 + Gθ θk̂ − 1

∀ k = 1 , ..., P2

x0̂ = x l

Avary, model

θk̂ − 1 = θ0̂

∀ k = 1 , ..., P2

Gθ Rstd Pstd

0≤τ

As improved estimates of θ become available, the achievable τ values are updated by solving the first-tier optimization repeatedly. Using a guideline similar to that used in section 3.3, P2 was chosen as 5 × τmax, where τmax is the dominant openloop time constant. When dealing with plant−model mismatch, care must be taken to not implement control action that is too aggressive. This would be particularly important initially (when the plant− model mismatch has not been adequately estimated). As the most natural form of detuning, during periods when significant plant−model mismatch exist, the computed τ values obtained from the τ-based optimization problem are relaxed. The model predictions are used as a means of determining the current plant−model mismatch, and if the mismatch is significant, the control action is detuned, by using larger values of the prescribed time constant as follows: if

l−1

l−2

then τi = τi ,comp else

τi = ki , ττi ,comp

Poles([ x1 x 2 x3 θ1 θ2 θ3 ])

Bvary, model diag{30, 50,30} 3 [0.66 0.72 0.55 0.6 0.54 0.72]

umin umax Δt (s) τmax (s)

[− 100 − 110 − 120] [150 140 130] 1.25 × 10−4 4

free mechanism, Gθ is chosen as B. The most effective observer Poles are presented in Table 2. The closed-loop measured and augmented state, as well as input profiles, are presented in = [MPCτ, SP, Figure 5. The color vector defined by MPCstd] represents the closed-loop profiles of the first-order performance specification-based MPC, the state setpoints, and that of a standard offset-free MPC approach (eq 6), respectively. With the first-order specification in place, the τ-based MPC is seen for many prescribed setpoints to lessen the suddenness in which actuator movement occurs. By choosing poles of an aggressive nature, the state estimates converge quickly to their corresponding measured values (and before the desired setpoint is reached), allowing for the specified designations to be seen in the closed-loop responses of the τ-based MPC. For certain setpoints, the standard offset-free formulation is able to bring the some of the states to the setpoint faster, albeit with significant overshoot in the other states. In contrast, the proposed τ-based formulation enables quick and yet smooth transients to the prescribed setpoints. Remark 11. In the above illustrative example with full state feedback, the primary purpose of the state estimator is to provide the estimates of the augmented states, which are then used by the MPC in the prediction model. Generally, the estimation scheme could perform multiple tasks including estimation of unmeasured state and noise handling, along with the estimation of the augmented states.

|γi | ≤ Toli and |γi | ≤ |γi | l−1

value

⎡ −1 −1 −2 ⎤ ⎢ ⎥ ⎢ 4 2 6⎥ ⎢⎣− 2 − 4 5⎥⎦ ⎡− 6.4 − 2 − 2.8 ⎤ ⎢ ⎥ ⎢ 2 − 4 3.3 ⎥ ⎢⎣− 3.7 1.5 − 2.1⎥⎦

(22)

where γil−1 is the latest estimate of the plant−model mismatch for state i, defined as γil−1 = x̂il − xil−1. The term Toli is a predefined tolerance that acts as a tuning parameter for each state i, τi,comp is the τ resulting from the τ-based optimization (eq 21), and ki,τ is a tuning parameter for each state i. Remark 9. In existing MPC formulations, one of the means of achieving detuning is to either decrease the penalty on state deviation or increase the penalty on input moves. In the present formulation, for a given first-order response, smaller values of the first-order time constants result in the prescribed control variable trajectory being a “fast” trajectory to the setpoint, resulting, in turn, in aggressive control action. Therefore, relaxing (or making larger) the values of the first-order time constants is a natural means of making the control action less 1001

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research

temperature, and CA, which is the concentration of A in the reactor. For the control design, a linearized model of the system (around the stable steady state, [Q, CA0, CA, TR] = [−6.55 × 104 kJ/s, 4.25 kmol/m3, 4.12 kmol/m3, 358 K]) was used. The closed-loop profiles are presented in Figure 6, with a magnified state profile in Figure 7.

Figure 6. State, input, and fictitious state closed-loop profiles for the nonlinear CSTR example.

Figure 5. Offset-free closed-loop state; input and fictitious state profiles.

5. APPLICATION TO A CHEMICAL PROCESS EXAMPLE Lastly, the proposed MPC structure was implemented on a r1 jacketed CSTR example where a reaction of the form A → B occurs, and a mole and energy balance results in the following model: CȦ =

⎛ F E ⎞ (CA0 − CA ) − k 0 ⎜exp − ⎟CATṘ V RTR ⎠ ⎝

Figure 7. Magnified TR and CA closed-loop profiles.

As evident in the closed-loop responses (see, e.g., the CA profile in Figure 7), the proposed MPC formulation enables smoother, and yet reasonably fast tracking of the set-point profile, in comparison to traditional offset-free MPC.

( −ΔH ) ⎛ Q E ⎞ F k 0⎜exp − = (TA0 − TR ) + ⎟CA + ρcp ρcpV V RTR ⎠ ⎝ (23)

6. CONCLUSION In this work, an MPC framework has been developed, explicitly incorporating performance specification in the control design. The approach is first illustrated using examples of prescribed closed-loop response in the absence of plant−model mismatch. Subsequently, a formulation is presented to account for plant−

where the manipulated inputs are the heat addition/removal due to the jacket (Q) and the inlet concentration of A in the feed (CA0), subject to the constraints −9.05 × 104 kJ/s ≤ Q ≤ −4.05 × 104 kJ/s, and 0.25 kmol/m3 ≤ CA0 ≤ 8.25 kmol/m3. The control objective is that of tracking a given set-point profile for the two state variables, TR, which is the reactor 1002

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003

Article

Industrial & Engineering Chemistry Research

Lyapunov-Based Predictive Control. Syst. & Contr. Lett. 2006, 55, 650−659. (17) Pannocchia, G.; Rawlings, J. Disturbance models for offset-free model predictive control. AIChE J. 2003, 49, 426−437. (18) Flores-Tlacuahuac, A.; Moreno, S. T.; Biegler, L. T. Global optimization of highly nonlinear dynamic systems. Ind. Eng. Chem. Res. 2008, 47, 2643−2655. (19) Wachter, A.; Biegler, L. T. Line search filter methods for nonlinear programming: Motivation and global convergence. SIAM J. Optim. 2005, 16, 1−31. (20) Li, X.; Chen, Y.; Barton, P. I. Nonconvex generalized benders decomposition with piecewise convex relaxations for globaloptimization of integrated process design and operation problems. Ind. Eng. Chem. Res. 2012, 51, 7287−7299. (21) Singer, A. B.; Barton, P. I. Global optimization with nonlinear ordinary differential equations. J. Glob. Optim. 2006, 34, 159−190.

model mismatch using the offset-free approach. The proposed approach is compared against traditional offset-free MPC and shown to provide superior (less oscillatory and often faster steady-state convergence) performance via application to a chemical process example.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Financial support from the Natural Sciences and Engineering Research Council of Canada through the Collaborative Research and Development Program (in collaboration with Johnson Controls, Inc.) is gratefully acknowledged. Also, insightful discussions with researchers John House and Tim Salsbury (Johnson Controls, Inc.) are gratefully acknowledged.



REFERENCES

(1) Mahmood, M.; Mhaskar, P. Enhanced Stability Regions for Model Predictive Control of Nonlinear Process Systems. AIChE J. 2008, 54, 1487−1498. (2) Mhaskar, P.; El-Farra, N. H.; Christofides, P. D. Predictive control of switched nonlinear systems with scheduled mode transitions. IEEE Trans. Autom. Control 2005, 50, 1670−1680. (3) El-Farra, N. H.; Mhaskar, P.; Christofides, P. D. Output feedback control of switched nonlinear systems using multiple Lyapunov functions. Syst. Control Lett. 2005, 54, 1163−1182. (4) Mhaskar, P. Robust model predictive control design for faulttolerant control of process systems. Ind. Eng. Chem. Res. 2006, 45, 8565−8574. (5) Mahmood, M.; Mhaskar, P. Lyapunov-based model predictive control of stochastic nonlinear systems. Automatica 2012, 48, 2271− 2276. (6) Muske, K. R.; Badgwell, T. Disturbance Modeling for Offset-free Linear Model Predictive Control. J. Process Control 2002, 12, 617−632. (7) Muske, K. R.; Rawlings, J. B. Model predictive control with linearmodels. AIChE J. 1993, 39, 262−287. (8) Wallace, M.; Das, B.; Mhaskar, P.; House, J.; Salsbury, T. Offsetfree model predictive control of a vapor compression cycle. J. Process Control 2012, 22, 1374−1386. (9) Das, B.; Mhaskar, P. Lyapunov-based offset-free model predictive control of nonlinear process systems. Can. J. Chem. Eng. 2015, 93, 471−478. (10) Angeli, D.; Amrit, R.; Rawlings, J. B. On average performance and stability of economic model predictive control. IEEE Trans. Autom. Control 2012, 57, 1615−1626. (11) Ellis, M.; Durand, H.; Christofides, P. D. A tutorial review of economic model predictive control methods. J. Process Control 2014, 24, 1156−1178. (12) Ellis, M.; Christofides, P. D. Selection of control configurations for economic model predictive control systems. AIChE J. 2014, 60, 3230−3242. (13) Skogestad, S., Postlethwaite, I. Multivariable Feedback Control; Wiley: New York, 2005. (14) Hopfe, N.; Ilchmann, A.; Ryan, E. Funnel Control With Saturation: Linear MIMO Systems. IEEE Trans. Autom. Control 2010, 55, 532−538. (15) Mayne, D. Q.; Rawlings, J. B.; Rao, C.; Scokaert, P. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789−814. (16) Mhaskar, P.; El-Farra, N. H.; Christofides, P. D. Stabilization of Nonlinear Systems with State and Control Constraints Using 1003

DOI: 10.1021/acs.iecr.5b03772 Ind. Eng. Chem. Res. 2016, 55, 995−1003