Fast Offset-Free Nonlinear Model Predictive Control Based on Moving

Mar 10, 2010 - One strategy for a constrained state estimation is the moving horizon estimation (MHE). MHE is based on solving a dynamic optimization ...
6 downloads 0 Views 2MB Size
Download PDF
7882

Ind. Eng. Chem. Res. 2010, 49, 7882–7890

Fast Offset-Free Nonlinear Model Predictive Control Based on Moving Horizon Estimation Rui Huang,†,‡ Lorenz T. Biegler,*,†,‡ and Sachin C. Patwardhan§ Collaboratory for Process and Dynamic Systems Research, National Energy Technology Laboratory, P.O. Box 880, Morgantown, West Virginia 26507-0880, Department of Chemical Engineering, Carnegie Mellon UniVersity, Pittsburgh, PennsylVania 15213, and Department of Chemical Engineering, Indian Institute of Technology, Bombay, Mumbai, Maharashra, 400076, India

To deal with plant-model mismatches in control practice, this paper proposes two variations of an offsetfree framework which integrates nonlinear model predictive control (NMPC) and moving horizon estimation (MHE). We prove that the proposed method achieves offset-free regulatory behavior, even in the presence of plant-model mismatches. If the plant uncertainty structure is known, the MHE can be tuned to estimate uncertainty parameters, to remove the plant-model mismatch online. In addition, we incorporate the advanced step NMPC (as-NMPC) and the advanced step MHE (as-MHE) strategies into the proposed method to reduce online computational delay. Finally, the proposed method is applied on a large scale air separation unit, and the steady state offset-free behavior is observed. Introduction Electricity is not readily stored and must be used or wasted after production. Therefore, power plants must be able to ramp up and down frequently to meet the fluctuating power demand. The combination of nonlinear model predictive control (NMPC) has emerged as a promising approach to meet this control requirement. Recently, with advanced algorithm and computational power, many successful NMPC implementations of industrial systems have been reported.1-3 However, as larger and more complicated process models are considered, computational delay associated with solving the NMPC problem online becomes an issue, resulting in degrated NMPC performance.4,5 To overcome this difficulty, Zavala and Biegler6 proposed the so-called adVanced step NMPC (as-NMPC) based on NLP sensitivity. It separates the computational burden into background and inexpensive online tasks. The only computational delay is the result of the online task. In recent work,7,8 we have illustrated that the computational delay of as-NMPC can be reduced by 2 orders of magnitude as compared to ideal NMPC, which solves the full NMPC problem online. Huang et al.8 reported a successful application of as-NMPC based on a rigorous dynamic model for an air separation unit (ASU) in an integrated gasification combined cycle (IGCC) power plant. It shows that the as-NMPC strategy is able to steer the ASU during the ramping processes and demonstrates superior performance over linear MPC. Moreover, although as-NMPC yields similar performance as ideal NMPC, the online computational delay is reduced from 4 min to less than 1 s. The previous study is based on the assumption that all the states of the ASU process are measured, and the control action is calculated on the basis of the perfect model. In practical applications, however, state information is usually not completely available for measurement, and disturbances and modeling errors are usually present and often not predictable. In addition, nonvanishing plant-model mismatch can develop over time, especially with large changes in operating conditions. Thus, it is critical to ensure controller’s * Address correspondence to [email protected]. † National Energy Technology Laboratory. ‡ Carnegie Mellon University. § Indian Institute of Technology-Bombay.

performance in the presence of plant-model mismatches, with a state estimator running in parallel to reconstruct the state information from the plant output. Offset-free behavior, that is, maintaining the controlled outputs at their desired steady-state set points, is an important requirement in control applications. To achieve offset-free behavior for linear systems, Muske and Badgwell9 presented a general disturbance model that accommodates unmeasured disturbances by state augmentation. Pannocchia and Bemporad10 proposed an integrated design strategy for disturbance model and dynamic observer. Moreover, Pannocchia and Kerrigan11 proposed a target setting strategy based on state augmentation to guarantee offset-free behavior for linear systems. For nonlinear systems, Meadows and Rawlings12 summarized a conventional offset-free (N)MPC method which adds an output disturbance into the objective function for the entire predictive horizon. The output disturbance is generated by comparing the measured plant output to the model prediction at the current time step. Huang et al.13 proposed an offset-free NMPC formulation that integrates both the state and output disturbances from the extended Kalman filter (EKF). It has the advantage that it is able to be applied to open-loop unstable systems. However, the design of EKF for a complex nonlinear system can be a challenging problem, especially when inequality constraints are enforced on state variables. One strategy for a constrained state estimation is the moving horizon estimation (MHE). MHE is based on solving a dynamic optimization problem online, and state constraints can be easily introduced in the optimization problem.14 On the basis of NLP sensitivity, Zavala et al.15 proposed advanced step MHE (as-MHE), to reduce the online computational delay for solving the optimization problem. In this work, we propose a general offset-free NMPC framework using MHE as the state estimator. In addition to the formulation with state and output disturbances, MHE can estimate the plant state and uncertainty parameter simultaneously, if the uncertainty structure is known. Therefore, the plant-model mismatch is adaptively removed. This paper is organized as follows: in the next section, we propose two variations of offset-free NMPC formulation with MHE as the estimator. Then, the as-NMPC and as-MHE are integrated into the framework to reduce the online computational time. Finally,

10.1021/ie901945y  2010 American Chemical Society Published on Web 03/10/2010

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

the proposed method is applied to the large scale air separation unit in the IGCC power plant.8 Offset-Free NMPC with MHE This work considers a general nonlinear discrete-time system with possible plant-model mismatches: xk+1 ) f(xk, uk, θk) yk ) h(xk),

(1a)

kg0

(1b)

where xk ∈ Rnx, uk ∈ Rnu, yk ∈ Rny and θk ∈ Ωθ ⊂ Rnθ are the plant states, controls, outputs, and uncertainty parameters, respectively, defined at time steps k g 0, and Ωθ is a compact set. Here if θk ) 0, (1) represents the nominal plant. Without losing generality, we assume that the given plant (1) has an equilibrium point at the origin, that is f(0, 0, 0) ) 0. Moreover, no trajectory of this system exhibits finite escape time. Depending on the prior knowledge of the plant-model mismatch, two variations of the offset-free NMPC with MHE are proposed in this section. In general, the uncertainty structure of a process is unknown; this could be due to sensor failures and process disturbances. In this case, we introduce state and output disturbances in both the NMPC and MHE problems to compensate for plant-model mismatch and achieve offset-free control behavior. On the contrary, some uncertainty structure of a process can be determined beforehand, for example, a certain parameter drifts away from its nominal value along a chemical reaction. Then we propose an MHE formulation to estimate the uncertainty parameter and plant state simultaneously. This formulation yields offset-free control behavior as well, due to the removal of the plant-model mismatch. Formulation with State and Output Disturbances. In the first scenario, we assume that the uncertainty structure is unknown. Hence, the nominal zero value of the uncertainty parameter is used in the predictive model for both MHE (2b) and NMPC (3b). To compensate for the plant-model mismatch, state and output disturbances are introduced. At time step k, with the measured output sequence yk-Ne, yk-Ne+1, ..., yk, the MHE problem based on the nominal uncertainty value is formulated as Ne

min

∑ (ζ

T k-Ne+jΠyζk-Ne+j

j)0

T + ξk-N Πxξk-Ne+j) + e+j

(xˆk-Ne - jxk-Ne)TΠ0(xˆk-Ne - jxk-Ne) (2a) subject to xˆk-Ne+j+1 ) f(xˆk-Ne+j, uk-Ne+j, 0) + ξk-Ne+j

(2b)

yˆk-Ne+j ) h(xˆk-Ne+j)

(2c)

ζk-Ne+j ) yk-Ne+j - yˆk-Ne+j

(2d)

xˆk-Ne+j ∈ X,

ζk-Ne+j ∈ Ωζ, j ) 0, ..., Ne

ξk-Ne+j ∈ Ωξ

(2e) (2f)

where Ne is the estimation horizon length, Πy, Πx, and Π0 are symmetric positive definite tuning matrices, xˆk and yˆk are the estimated state and output values, and ξk and ζk are the state and output disturbances which are assumed to be bounded in compact sets Ωζ and Ωξ, respectively. In addition, jxk-Ne is the most likely prior value of xk-Ne. After the MHE problem is

7883

solved at time step k, we choose xˆk-Ne+1 as the prior value jxk-Ne+1 for the arrival cost in (2a) at time step k + 1. Note that though we consider a noise-free plant (1) for notational simplicity, the proposed MHE and NMPC can easily incorporate the state and output noises, due to the introduction of the state and output disturbances. In addition, state noise can be considered as a special form of the uncertainty parameter θ. After the MHE problem (2b) is solved, it yields the optimal estimated state (xˆk), as well as the state and output disturbances (ξk, ζk) at the current time step k. Therefore, NMPC with state and output disturbances is formulated as Nc-1

Np

min

∑ (l

k+j

- yr)TΓy(lk+j - yr) +

j)0

∑ ∆V

T k+iΓu∆Vk+i

i)0

(3a) subject to zk+j+1 ) f(zk+j, Vk+j, 0) + ξk, zk ) xˆk, lk+j ) h(zk+j) + ζk,

j ) 0, ..., Np - 1

zk+j ∈ X ∆Vk+i ) Vk+i+1 - Vk+i

(3b) (3c) (3d)

Vk+j ) Vk+ti for ti e j < ti+1, Vk+i ∈ U

(3e)

t0 ) 0 e t1 e t2 e , ..., e Np - 1

(3f)

where Np and Nc are the prediction and control horizon length, respectively; yr is the set point for the output; zk, lk, and Vk are the predicted state, output, and control movement, respectively. In addition, the predicted state and control are required to satisfy the inequality constraints zk+j∈X, and Vk+I∈U. NMPC is initialized with the estimated state xˆk. In typical NMPC applications, fewer degrees of freedom are available for the control movement. The last two constraints indicate that the control action is the input blocking form, ensuring that the available degrees of freedom spread over the entire prediction horizon. After the NMPC problem is solved, the first manipulated variable Vk is injected into the plant, that is, uk ) Vk. Note that unlike conventional NMPC formulations, the predictive model is perturbed by the state and output disturbances (ζk and ξk) to compensate for the unknown uncertainty parameter θk. It is worth mentioning that ζk and ξk, which are the calculated value at the end of the estimation horizon in MHE, are fixed parameters in the NMPC problem. We now show that the proposed method yields zero steadystate offset. The analysis is similar to that in Meadows and Rawlings.12 However, it does not depend on a target setting optimization problem. Theorem 1. If (a) the set point yr is reachable for the perturbed predictive model zk+j+1 ) f(zk+j, Vk+j, 0) + ξk and lk+j ) h(zk+j) + ζk, (b) the NMPC controller (3) is asymptotically stabilizing for the perturbed predictive model, (c) the closedloop system goes to a steady state, (d) the perturbed predictive model is observable at the steady state, then the system controlled by the MHE (2) and the NMPC (3) has zero steadystate offset. Proof: In the following analysis, the superscript ss denotes the steady-state value. Since yr is reachable for the perturbed predictive model and the NMPC control law is asymptotically stable, the stage cost in NMPC (3) is zero at steady state, that is, lss ) yr

(4)

7884

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

Moreover, at the steady state, the predictive state remains constant, zss ) xˆss. The control action is also a constant, that is, uss ) Vss. The state and output disturbances (ξss and ζss) can be estimated from the assumption that the perturbed predictive model is observable. Thus the following equations hold true: xˆss ) f(xˆss, uss, 0) + ξss

(5a)

lss ) h(xˆss) + ζss

(5b)

In addition, at the steady state the MHE evolves according to

(2b), ζk-Ne+j is the output disturbance, bounded in the compact set Ωζ. However, there is no state disturbance in this formulation. MHE (9) yields the optimal uncertainty parameter θˆ k that minimizes the difference between the estimated output yˆ, and the measured output y over the entire estimation horizon. As a result, the estimated state xˆ is also smoothed over the entire estimation horizon. With the optimal estimated state xˆk, output disturbance ζk and the optimal uncertainty parameter θˆ k at the current time step k, the NMPC with parameter estimation is formulated as the following: Nc-1

Np

xˆss ) f(xˆss, uss, 0) + ξss

(6a)

yˆss ) h(xˆss)

(6b)

ζss ) yss - yˆss

(6c)

∑

min

(lk+j - yr)TΓy(lk+j - yr) +

j)0

Since the predictive model in the MHE (6a) is exactly the same as that in the NMPC (5a), by combining (5b), (6b), and (6c), we see lss ) yss

(7)

Then by virtue of (4) and (7), the following equation can be derived: yss ) yr

(8)

It indicates that the plant output is equal to the set point at the steady state. 9 Remark 1. If the set point yr is not feasible, then the proposed approach minimizes the steady-state output difference, that is, (yss - yr)TΓy(yss - yr). Remark 2. Note that in this analysis, the observer is in a general formulation (6). Thus, the Theorem also applies to the NMPC incorporated with nonlinear recursive observers, such as extended Kalman filter and extended Luenberger Observer. Formulation with State and Parameter Estimation. In the second scenario, we assume that the uncertainty parameter structure is known. Then instead of compensating for the uncertainty by state and output disturbances, MHE can estimate both the state and uncertainty simultaneously. Thus, the model for the MHE and the NMPC is modified adaptively online. At time step k, with the measured output sequence yk-Ne, yk-Ne+1, ..., yk, the MHE is formulated as

∑ ∆V

T k+jΓu∆Vk+j

j)0

(10a) subject to zk+j+1 ) f(zk+j, Vk+j, θˆ k), zk ) xˆk, lk+j ) h(zk+j) + ζk,

j ) 0, ..., Np - 1

zk+j ∈ X

(10c)

∆Vk+j ) Vk+j+1 - Vk+j

Vk+j ) Vk+ti for ti e j < ti+1,

(10b)

Vk+j ∈ U

t0 ) 0 e t1 e t2 e , ..., e Np - 1

(10d) (10e) (10f)

Similar to the NMPC formulation with state and output disturbances (3b), this NMPC problem is initialized with the estimated state variable xˆk and the planned control movement is chosen as the input-blocking form. However, this NMPC formulation does not perturb the predictive model with the state disturbance. Instead, it adaptively updates the uncertainty parameter with the optimal value from the MHE. Similar to the analysis in Theorem 1, we can show that the NMPC and MHE with parameter estimation (MHE (9) and NMPC (10)) is able to provide offset-free control behavior. The proof follows in the same way as that in Theorem 1, and is omitted here. Illustrative Simulation Example. In this subsection, we illustrate the proposed two variations of offset-free NMPC with MHE through a simulation example. Here, a simulated NMPC scenario with a nonlinear continuous stirred tank reactor (CSTR) is considered. The CSTR is represented by the following differential equations:

Ne

min

∑ (ζ j)0

T k-Ne+jΠyζk-Ne+j)

+ θˆ Tk Πθθˆ k +

(xˆk-Ne - jxk-Ne) Π0(xˆk-Ne - jxk-Ne) (9a) T

subject to xˆk-Ne+j+1 ) f(xˆk-Ne+j, uk-Ne+j, θˆ k)

(9b)

yˆk-Ne+j ) h(xˆk-Ne+j)

(9c)

ζk-Ne+j ) yk-Ne+j - yˆk-Ne+j

(9d)

xˆk-Ne+j ∈ X,

ζk-Ne+j ∈ Ωζ, j ) 0, ..., Ne

θˆ k ∈ Ωθ

(9e) (9f)

where θˆ k is the estimated uncertainty parameter which is bounded in a compact set Ωθ. Similar to the MHE formulation

dzc ) (zc - 1)/u2 + k0zc exp(-Ea /zT) dt

(11a)

dzT ) (zT - zfT)/u2 + k0zc exp(-Ea /zT) + νu1(zT - zcw T ) dt (11b) This system involves two states z )[zc, zT] corresponding to dimensionless concentration and temperature, and two manipulated inputs, corresponding to the cooling water flow rate u1 and the inverse of dilution rate u2. The model parameters are zcw T ) 0.38, zfT ) 0.395, and Ea ) 5, ν ) 1.95 × 104, and k0 is an uncertainty parameter in the plant with nominal value k0 ) 300. The system is operated at steady state zss ) [0.12466, 0.74068] corresponding to uss )[378, 20]. In this simulation, it is assumed that both the states are measured, that is, the output mapping function in (1b) is chosen to be h( · ) ) diag[1, 1] × [zc, zT]T. This model is regulated by the proposed two variations of offset-free NMPC with MHE. The horizon of the MHE, Ne, is

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

Figure 1. State profile in scenario 1.

Figure 3. Error profile in scenario 1.

Figure 2. Control profile in scenario 1.

Figure 4. State profile in scenario 2.

chosen as 6 time units with sampling time being equal to 1 time unit. Let Q ) diag[1, 1] and R ) diag[1, 1], and define Ass ) ∂f/∂z|zss,uss, Bss ) ∂f/∂u|zss,uss, Css ) ∂h/∂z|zss and Vss ) R T + CssQCss . The weighting matrices are chosen to be inverse of the covariance information which is calculated similarly to the T T extended Kalman filter, that is, Π0-1 ) AssQAss + BssQBss T T AssQCss Vss-1CssQAss , Πx-1 ) Q, Πy-1 ) R and Πθ ) 0. The NMPC is tuned with prediction horizon Np chosen to be 10 time units and control horizon Nc chosen to be 5 time units and equally distributed over the entire prediction horizon. The tuning matrices are Γy ) diag[1 × 106, 1 × 106], Γu ) 0. In the first simulation, the plant is controlled by the formulation with state and output disturbances (MHE (2) and NMPC (3)). The plant starts from the nominal steady-state value. At time step 10, the uncertainty parameter k0 is reduced to 70% of its nominal value as shown at the bottom of Figure 3. The resulting closed-loop responses are shown in Figure 1 and Figure 2. It is observed that although both states are measured, the state estimates are biased after the plant-model mismatch is introduced. However, these differences are utilized in the NMPC formulation to remove the steady-state offset of the plant outputs. As shown in Figure 1, the proposed method is able to regulate the plant outputs at the desired set points. The second simulation scenario is the same as the first one, except that the plant is regulated by the formulation with state and parameter estimation (MHE (9) and NMPC (10)). As shown

7885

Figure 5. Control and uncertainty profile in scenario 2.

at the bottom of Figure 5, the estimated uncertainty gradually converges to the plant value after the plant-model mismatch is introduced at time step 10. Thus, after time step 25, the uncertainty parameter in the model is equal to that in the plant, removing the plant-model mismatch. Figure 4 shows that the proposed method quickly rejects the disturbance and yields

7886

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

offset-free control behavior. In addition, the estimated states converge to the measured plant states without any error after the plant-model mismatch is eliminated. It is interesting to compare Figure 1 and Figure 4 to see that the formulation with parameter estimation (9)-(10) rejects the plant-model mismatch faster than the formulation with state and output disturbances (2)-(3). Offset-Free as-MHE-NMPC In the previous section, the CSTR is controlled by the hypothetical ideal-MHE-NMPC, which is based on the assumption that the MHE and NMPC problems can be solved instantaneously, that is, no computational delay. However, in practice, the computational delay resulting from solving the optimization problems associated with MHE and NMPC may degrade the controller’s performance,4,5 especially for large scale processes. As-NMPC6 and as-MHE,15 which could reduce the computational delay by 2 orders of magnitude, are alternatives to approximate the solution of hypothetical ideal-MHE-NMPC. In this section, we first recall the as-MHE15 and as-NMPC16 algorithm, which is based on the IPOPT16 algorithm and NLP sensitivity.17 The as-MHE and as-NMPC are incorporated with the proposed offset-free framework. Then the illustrative simulation example in the previous section is revisited. Offset-Free as-MHE-NMPC Algorithm. The MHE and NMPC formulations can be represented as a general parametric nonlinear programming (NLP) problem of the form min f(w(p), p) w(p)

(12a)

subject to c(w(p), p) ) 0,

w(p) g 0

(12b)

with a fixed parameter vector p. Note the implicit dependence of the problem variables on the particular value of the parameter. For a given nominal parameter vector p ) p0, IPOPT solves the NLP problem through Newton step K*(p0)∆s ) -φ(s*(p0), p0)

(13)

where s*(p0) is the optimal solution vector while φ( · , · ) and K*( · ) are the KKT conditions and KKT matrix evaluated at this solution, respectively. Here ∆s ) 0. Assume that the objective and constraints in NLP (12) are at least twice continuously differentiable with respect to p and that a given nominal solution s*(p0) satisfies the linear independence constraint qualification (LICQ), the sufficient second order conditions (SSOC), and strict complementary slackness. This assumption implies that the corresponding system for NMPC and MHE is locally controllable and observable, respectively. Moreover, if the parameter vector p0 enters linearly in the constraints, it can be shown17 that replacing the right hand side in (13) by -φ(s*(p0),p) leads to the first order step ∆s ) s˜(p) - s*(p0). Here, s˜(p) is an approximation to the true optimal solution s*(p). This result has important practical implications since the factorization of the KKT matrix in (13) is already available as a natural outcome of the NLP solver. Consequently, the approximated solution s˜(p) can be obtained through a single backsolve, which can be performed over 2 orders of magnitude faster than the solution of the full NLP problem. Hence, the as-MHE and as-NMPC can be incorporated into the previously proposed framework, yielding offset-free asMHE-NMPC.

Figure 6. State profile in the revisiting of scenario 1.

In Background, Between Time Step k and k + 1. (i) asMHE: having xˆk and uk, compute the disturbance-free extrapolation of the state x˜k+1 and the corresponding output y˜k+1. Solve an extended MHE problem (2) or (9) with horizon Ne + 1 and output sequence yk-Ne, yk-Ne + 1, ..., yk, y˜k+1. Let p0 ) y˜k+1 and hold the KKT matrix at the solution. (ii) as-NMPC: having xˆk and uk, predict the future state of the system zk+1 using the predictive model (3b) or (10b). Solve the NMPC (3) or (10) with p0 ) zk+1. Hold the KKT matrix at the solution. Online Update, at Time Step k + 1. (i) as-MHE: obtain the true measurement yk+1. Set p ) yk+1 and compute the fast approximation solution using (13) for MHE. Extract the estimated state xˆk+1. (ii) as-NMPC: obtain the state estimate xˆk+1. Set p ) xˆk+1 and use (13) for NMPC to get the fast updated solution. Extract the control action uk+1 from the approximate solution vector and inject to the plant. It is observed that the online computational delay of this asMHE-NMPC scheme is only two single backsolves associated to the MHE and NMPC problems, which is significantly faster compared to solving the MHE and NMPC problems online. In addition, the background solution of as-MHE-NMPC can be further parallelized on different computing units to ensure that it can be finished within one sampling time. Revisiting the Illustrative Simulation Example. In this subsection, we revisit the two simulation scenarios of the illustrative CSTR example (11) with offset-free as-MHE-NMPC. The model parameters are the same, and k0 is still considered as the uncertainty parameter. In addition, the tuning parameters in the two simulations are chosen exactly the same as before. In scenario 1, we use as-MHE-NMPC formulation with state and output disturbances. The as-MHE-NMPC is not able to handle the drastic change of the uncertainty parameter k0 by 30% at time step 10. It loses the robust stability due to the large perturbation. Consequently, in this simulation, the uncertainty parameter k0 is reduced by 20% at time step 10 as shown at the bottom of Figure 8. The closed-loop plant response is shown in Figure 6 and Figure 7. Here, the profiles are the plant responses controlled by ideal-MHE-NMPC and as-MHENMPC; the profile of state estimation is omitted. It is observed that both as-MHE-NMPC and ideal-MHE-NMPC are able to regulate the plant outputs at their desired set points without any offset. In addition, note in Figure 7 that the control action from as-MHE-NMPC differs slightly from ideal-MHE-NMPC between time steps 10 and 20 and is otherwise the same. This is due to the sensitivity error introduced by as-MHE and as-NMPC.

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

Figure 7. Control profile in the revisiting of scenario 1.

7887

Figure 9. State profile in the revisiting of scenario 2.

Figure 8. Error profile in the revisiting of scenario 1. Figure 10. Control and uncertainty profile in the revisiting of scenario 2.

In scenario 2, the as-MHE-NMPC formulation with parameter estimation is utilized. Since the as-MHE (9) is able to update the uncertainty parameter and remove the plant-model mismatch online, the as-MHE-NMPC is stable even when the uncertainty parameter k0 is reduced drastically by 30%, as shown at the bottom of Figure 10. The closed-loop plant response is shown in Figure 9 and Figure 10. Both as-MHE-NMPC and ideal-MHE-NMPC are able to regulate the plant outputs at their desired set points without any offset. Moreover, as-MHE-NMPC has almost identical performance with ideal-MHE-NMPC because the plant-model mismatch is removed. Hence, in practical applications, it is recommended to estimate parameters with the state simultaneously when possible. ASU Simulation In this study, we consider an industrial-size air separation unit in IGCC power plants.8 The unit contains two integrated cryogenic distillation columns as shown in Figure 11. The high pressure column (bottom) has 40 trays and operates at 5-6 bar, while the low pressure column (top) operates at 1-1.5 bar and also has 40 trays. An air feed flow is split into two substreams. The high pressure air (MA) enters the bottom of the high pressure column and the expanded air (EA) enters the 20th tray of the low pressure column. Crude nitrogen gas (GN) from the main heat exchanger is also added to the 25th tray of the high

pressure column. The reboiler of the low pressure column is integrated with the condenser of the high pressure column. The main products of the high pressure column are pure nitrogen (PNI) (>99.99%) and crude liquid oxygen (∼50%). The crude oxygen stream is fed into the 19th tray of the low pressure column. In addition, an intermediate side stream from the 15th tray of the high pressure column (LN) is fed into the top of the low pressure column. A high purity separation is achieved in the low pressure column, leading to a nitrogen product with ∼99% purity and an oxygen product (POX) with ∼97% purity. The ASU model is represented by tray-by-tray equations consisting of mass balances (overall and component balances of nitrogen, oxygen, and argon), energy balances, phase equilibrium, hydraulic, and summation equations. Overall Mass Balance. dMi ) Li-1 + Vi+1 - Li - Vi + Fi dt

(14)

where i is the index of each tray, starting from the top of the column. Mi is the liquid mole holdup ([mol]) on tray i, Li and Vi are liquid and vapor molar flow rates, respectively, and Fi is the molar feed ([mol/min]). If there is no feed to tray i, then Fi ) 0. In this case, the only nonzero values of Fi are those corresponding to expanded air (EA, U1), pure air (MA, U2),

7888

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

Summation Equation. 1)

∑

di,j

(19)

j∈COMP

Hydraulic Equation. Li ) kdMi

(20)

where kd is considered as an uncertainty parameter with kd ) 0.5 min-1 as the nominal value. Vapor-Liquid Equilibrium. di,jpi ) γi,jci,jpsat i,j

Figure 11. Air separation unit flow sheet.

liquid nitrogen (LN, U3), crude nitrogen gas (GN, U4), and crude oxygen stream as shown in Figure 11. Component Balances. Mi

dci,j ) Li-1(ci-1,j - ci,j) + Vi+1(di+1,j - ci,j) dt Vi(di,j - ci,j) + Fi(cfi,j - ci,j) (15)

where j ∈ COMP is the index of each component, ci, j and di, j are component mole fractions in the liquid and vapor phases, and ci,f j is the mole fraction of the feed. Energy Balance. Mi

(

∂hLi j + T ∂Ti i

)

∂hLi L jc ) Li-1(hi-1 - hLi ) + ∂ci,j i,j j∈COMP

∑

V Vi+1(hi+1

-

hLi )

- Vi(hVi - hLi ) + Fi(hfi - hli) (16)

where hiL ) f hl(Ti, Pi) and hiV ) f hv(Ti, Pi) are liquid and vapor enthalpies in [kJ/mol], pi is the total pressure on tray i, hif is the feed enthalpy, and Ti is the tray temperature. Expressions and L data to compute hV i and hi can be found in a number of standard references. Moreover jci,j :)

dci,j ) (Li-1(ci-1,j - ci,j) + Vi+1(di+1,j - di,j) dt Vi(di,j - ci,j) + Fi(cfi,j - ci,j))/Mi (17)

and

j i :) T

∑

dTi j∈COMP )dt

[

ci,j

∑

k∈COMP

∑

(

)

∂Ki,j jc + Ki,jjci,k ∂ci,k i,k

]

ci,j∂Ki,j /∂Ti

j∈COMP

(18) Here we define Ki, j ) γi, jpi,satj /pi, where pi,satj ) f p(Ti) is the saturation pressure of pure component j on tray i, and γi, j denotes the liquid activity coefficient describing the nonideal vapor-liquid equilibrium calculated as in Huang et al.8

(21)

After modeling each tray in the distillation columns by (14)-(21), the ASU model is composed of 320 differential equations and 1200 algebraic equations. The control structure is set the same as that reported in Huang et al.8 We choose the molar flow rate of product oxygen (POXY1), product nitrogen (PNI-Y2), the temperature at 30th tray in the low pressure column (Tl30-Y3), and temperature at the 15th tray in the high pressure column (Th15-Y4) as output variables. Four stream flow rates are considered as manipulated variables, including the expanded air feed (EA-U1), main air feed (MAU2), reflux liquid nitrogen (LN-U3), and crude nitrogen gas (GNU4). The objective is to regulate the outputs at their set points in the presence of the plant-model mismatch. Since there are many uncertainties in the ASU process, for example, thermodynamic properties, tray efficiencies, etc., it is not trivial to determine the uncertainty structure of the ASU model in practice. Therefore, we choose to use the formulation with state and output disturbances (MHE (2) and NMPC (3)). The prediction (Np) horizon and control (Nc) horizon in the NMPC formulation are chosen to be 20 time units, each with a 5 min sampling time. Γy is chosen to be a diagonal matrix with 3 × 10-2 corresponding to the Tl30, Th15, and 1 × 10-4 corresponding to the POX, PNI, while Γu is set to be a null matrix. For the MHE formulation, the sampling time is still 5 min, but the estimation horizon (Ne) is chosen to be 5 time units. The MHE is tuned with Πy as a diagonal matrix with 1 × 10-2 corresponding to Tl30, Th15, and 1 × 10-4 corresponding to POX, PNI; Πx is chosen as a diagonal matrix with 1 corresponding to the ci, j,∀i, j and 1 × 10-5 corresponding to Mi,∀i; Π0 is a diagonal matrix with 5 corresponding to the ci, j,∀i, j and 5 × 10-5 corresponding to Mi,∀i. It is worth mentioning that in this problem, beyond reasonable guidelines, the as-MHENMPC performance is not very sensitive to the tuning parameters. The simulation starts from the nominal steady state with kd ) 0.5. At 25 min, the kd in the plant is increased to 0.6, while kd ) 0.5 in the model, introducing a plant-model mismatch. The closed-loop plant output is presented in Figure 12, and the control action is shown in Figure 13. It is observed that both as-MHE-NMPC and ideal-MHE-NMPC in the proposed offsetfree framework reject the disturbance and regulate the plant outputs at their set points without any steady-state offset. Moreover, the as-MHE-NMPC yields comparable performance as the ideal-MHE-NMPC. In addition, from the output of IPOPT, we see that the solutions of both MHE and NMPC problems satisfy LICQ and SSOC at each time step, meaning that this system is locally observable and controllable at each time step. This performance arises because the ASU is openloop stable and the MHE and NMPC problems at each time step are well initialized from their previous solutions.

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

Figure 12. Output profile of the ASU.

7889

Figure 14. Product purity profile of the ASU.

NMPC, the online computational delay is reduced by 150 times, from around 5 min to less than 2 s. Summary

Figure 13. Control profile of the ASU.

In this control structure setting, there are only four inputs, which means that the offset-free control behavior can be achieved for up to four outputs. We choose to use the same four measurements as in the previous study.8 The purities in the product streams POX and PNI are not directly measured. Nevertheless, it is interesting to note in Figure 14 that the oxygen and nitrogen purity in the product streams satisfy the requirement (oxygen purity g96%). In this work, the ASU model which consists of 1520 differential algebraic equations (DAE), is discretized using a finite element discretization at Radau collocation points.18 After discretization, the resulting nonlinear programming (NLP) problem corresponding to the NMPC problem (20 finite elements and 3 collocations) has 116 900 constraints and 117 140 variables; while the NLP corresponding to the MHE problem (five finite elements and three collocations) contains 29 285 constraints and 30 885 variables. Both the NLPs are solved using AMPL and IPOPT on an Intel DuoCore 2.4 GHz personal computer. The ideal-MHE problems take up to 15 iterations and 90 CPU seconds to solve, while ideal-NMPC problems take up to 6 IPOPT iterations and 200 CPU seconds to solve. On the contrary, the online computational time is less than 1 CPU second for as-MHE and around 1 CPU second for as-NMPC, despite that as-MHE-NMPC has similar performance with ideal-MHE-NMPC. It indicates that by using as-MHE-

This work addresses the practical issue of achieving offsetfree behavior in NMPC applications when plant-model mismatch is present and a state estimator is integrated to reconstruct plant state from output. Using MHE as the general state estimator, we propose an offset-free NMPC framework based on state and output disturbances. It can be shown that the proposed method guarantees offset-free behavior if the set point is feasible for the observable perturbed predictive model and the control law asymptotically converges to the steady state. Moreover, another variation of offset-free NMPC framework based on integrated state and parameter estimation is proposed. Here, MHE is tuned to estimate the plant state and uncertainty parameter simultaneously, in order to remove the plant-model mismatch online. Provided that the uncertainty structure is known, we observe that this framework has better performance than the formulation based on state and output disturbances. In addition, as-NMPC and as-MHE is incorporated into the proposed offset-free frameworks to reduce online computational delay associated with solving the MHE and NMPC problems. We see that the framework with parameter estimation can improve the robust stability of as-MHE-NMPC. Finally, the proposed as-MHE-NMPC based on state and output disturbances is successfully implemented on a large scale air separation unit. The online computational delay is reduced by at least 150 times by using as-MHE-NMPC. Acknowledgment Partial support of this work was provided by the National Energy Technology Laboratory’s ongoing research in Process and Dynamic Systems Research under the RDS Contract No. DE-AC26-04NT41817. Literature Cited (1) Bartusiak, R. NLMPC: A platform for optimal control of feed- or product-flexible manufacturing. In Assessment and Future Directions of Nonlinear Model PredictiVe Control; Findeisen, R., Allgo¨wer, F., Biegler, L., Eds.; Springer: New York, 2007; pp 1-16. (2) Nagy, Z.; Mahn, B.; Franke, R.; Allgo¨wer, F. Real-time implementation of nonlinear model predictive control of batch process in an industrial framework. In Assessment and Future Directions of Nonlinear Model

7890

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

PredictiVe Control; Findeisen, R., Allgo¨wer, F., Biegler, L., Eds.; Springer: New York, 2007; pp 465-472. (3) Diehl, M.; Uslu, I.; Findeisen, R.; Schwarzkopf, S.; Allgo¨wer, F.; Bock, H. G.; Burner, T.; Gilles, E. D.; Kienle, A.; Schloder, J. P. Realtime optimization for large scale processes: nonlinear model predictive control of a high purity distillation column. J. Process Control 2002, 12, 577–585. (4) Findeisen, R.; Allgo¨wer, F. Computational delay in nonlinear model predictive control. Proceedings of the International Symposium Advanced Control of Chemical Processes, Hong Kong, 2004. (5) Santos, L.; Afonso, P.; Castro, J.; Oliveira, N.; Biegler, L. On-line implementation of nonlinear MPC: An experimental case study. Control Eng. Practice 2001, 9, 847–857. (6) Zavala, V.; Biegler, L. The advanced step NMPC Controller: Optimality, stability and robustness. Automatica 2009, 45, 86–93. (7) Zavala, V.; Laird, C.; Biegler, L. Fast implementations and rigorous models: Can both be accommodated in NMPC? Int. J. Robust Nonlinear Control 2007, 18, 800–815. (8) Huang, R.; Zavala, V.; Biegler, L. Advanced step nonlinear model predictive control for air separation units. J. Process Control 2009, 19, 678– 685. (9) Muske, K.; Badgwell, T. Disturbance modeling for offset-free linear model predictive control. J. Process Control 2002, 12, 617–632. (10) Pannocchia, G.; Bemporad, A. Combined design of disturbance model and observer for offset-free model predictive control. IEEE Trans. Automat. Control 2007, 52, 1048–1053.

(11) Pannocchia, G.; Kerrigan, E. Offset-free receding horizon control of constrained linear systems. AIChE J. 2005, 51, 3134–3146. (12) Meadows, E.; Rawlings, J. Model predictive control. In Nonlinear Process Control; Henson, M., Seborg, D., Eds.; Prentice Hall: Upper Saddle River, NJ, 1997; pp 233-310. (13) Huang, R.; Patwardhan, S.; Biegler, L. Robust extended Kalman filter based nonlinear model predictive control formulation. Proceedings of the 48th IEEE Conference on Decision and Control, Shanghai, China, 2009. (14) Rao, C.; Rawlings, J.; Mayne, D. Constrained state estimation for nonlinear-discrete-time systems: stability and moving horizon approximations. IEEE Trans. Automat. Control 2003, 48, 246–258. (15) Zavala, V.; Laird, C.; Biegler, L. A fast moving horizon estimation algorithm based on nonlinear programming sensitivity. J. Process Control 2008, 18, 876–884. (16) Wa¨chter, A.; Biegler, L. On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Program 2006, 106, 25–27. (17) Fiacco, A. Introduction to SensitiVity and Stability Analysis in Nonlinear Programming; Academic Press: New York, 1983. (18) Biegler, L.; Cervantes, A.; Wa¨chter, A. Advances in simultaneous strategies for dynamic process optimization. Chem. Eng. Sci. 2002, 57, 575– 593.

ReceiVed for reView December 8, 2009 ReVised manuscript receiVed February 19, 2010 Accepted February 22, 2010 IE901945Y