This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.
Article pubs.acs.org/IECR
Real-Time Nonlinear Model Predictive Control of a Transport− Reaction System Andreas Steinboeck,*,† Martin Guay,‡ and Andreas Kugi† †
Automation and Control Institute, Vienna University of Technology, 1040 Vienna, Austria Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6
‡
ABSTRACT: Two real-time nonlinear model predictive control (NMPC) algorithms for a transport−reaction system are designed. The system is modeled by a hyperbolic partial differential equation and discretized by means of a two-timelevel semi-implicit semi-Lagrangian scheme. For the resulting lumped-parameter system, a constrained optimal control problem is formulated and state constraints are implemented in the form of barrier functions. The NMPC algorithms perform a single step or several steps of an iterative solution routine of the optimal control problem at every sampling point. With this suboptimal solution strategy, a fixed maximum evaluation time and execution in real time are guaranteed. An analysis of the nominal stability is provided for one NMPC scheme. The robustness of the controllers is evaluated for an example problem, where a nonisothermal plug-flow reactor with irreversible exothermic reactions is considered. The control objectives are to limit the maximum reactor temperature (avoid hot spots) and to maximize the process output.
1. INTRODUCTION Applications in the chemical and in the process industry have been major drivers of the development of model predictive control (MPC) since its beginnings in the 1970s.1−4 In chemical engineering applications, current issues of MPC include the development of nonlinear model predictive control (NMPC), guaranteed closed-loop stability, efficient computation, fast optimization methods, and constraints.5 While linear MPC is widely accepted by control engineers, the theoretical complexity and the computational demands associated with NMPC hamper its widespread application.6 A common challenge in NMPC is the typically high computational load associated with the model-based prediction of the state trajectory and the solution of the underlying optimal control problem. Specifically, the numerical or analytical computation of gradients of the objective function with respect to the control inputs is a demanding task.7 These challenges become more severe if the system considered has distributed parameters, i.e., if it is infinite dimensional. In this case, a common solution is to use lumpedparameter approximations of the system in the MPC formulation.8 To achieve the desired degree of accuracy, these models are typically high dimensional, which again entails a high computational load. For transport−reaction systems, an alternative integration of the partial differential equations utilizing the method of characteristics may help to keep the model order low. Direct application of the method of characteristics with Lagrangian coordinates is not recommended since the state vector may have a variable dimension and meaning, i.e., states may continuously enter and leave the system. This makes stability proofs difficult. The problem can © XXXX American Chemical Society
be avoided if an upwind discretization scheme with an equidistant spatial grid, a fixed time step, and a Courant number9 equal to 1 is used.10,11 In recent times, economic MPC12 moved into scientific focus for the control of chemical processes. Economic MPC has been used in combination with a backup stabilizing state-feedback controller for the control of concentrations and temperatures in continuous stirred tank reactors.13 Economic MPC has also been applied to distributed-parameter chemical processes, e. g., to diffusion-convection-reaction systems using a high-order Galerkin discretization.14 The use of high-order lumpedparameter models is necessary if state constraints are to be satisfied. In many chemical engineering applications, controllers must ensure that temperature constraints are satisfied. In plug-flow reactors, so-called hot spots may typically occur at unknown positions that depend on the boundary conditions and the system inputs. The control of hot spots is of ongoing scientific interest and can be done by PI control,15 nonlinear feedback controllers,16 input−output linearization,17 generalized predictive control,18 and economic MPC.14 If the heat transfer system is appropriately designed, hot spots in plug-flow reactors are best controlled by means of the cooling rate.15 Given this state of the art, the current paper aims at the development of real-time NMPC and its application to transport−reaction systems. Here, the primary control objective Received: February 12, 2016 Revised: June 20, 2016 Accepted: June 23, 2016
A
DOI: 10.1021/acs.iecr.6b00592 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
simply connected set < . Clearly, we assume x s ∈ ? and us ∈ < . System 1 may be an early lumping, discrete-time approximation of a distributed-parameter system, e.g., a transport−reaction system. 2.2. Receding Horizon. The receding prediction horizon used in MPC ranges from tk to tk+Nk, i.e., it covers Nk sampling intervals. The integer Nk may depend on k. Figure 1 shows an
is the avoidance of hot spots and the secondary control objective is the maximization of the process output. The following questions are to be answered in this work: • Can the control objectives be satisfied by means of realtime NMPC? • Is (nominal) closed-loop stability of the NMPCcontrolled system ensured? • Is the design robust against disturbances like model− plant mismatches or changing operating conditions? The paper is organized as follows: A real-time MPC algorithm with guaranteed nominal closed-loop stability and an ad-hoc real-time MPC algorithm are developed in section 2 for a general discrete-time dynamical model. The distributedparameter model of a transport−reaction system and its discretization are introduced in section 3. An example problem of a nonisothermal plug-flow reactor with hot spots is analyzed in section 4. The two control designs and an open-loop controller are compared, and the robustness of these control strategies against model−plant mismatches and variations of boundary conditions is analyzed in five different simulation scenarios. Conclusions are summarized in section 5.
Figure 1. Time grids during the prediction horizon, shown for Nk = 5, M = 2, and K = 3.
2. REAL-TIME MPC MPC typically involves the solution or suboptimal solution of a dynamic optimization problem. Especially for fast dynamical systems and distributed-parameter systems, this is still a major challenge for MPC. To ensure that the MPC can be executed in real time, we describe a discrete-time system formulation in section 2.1, define a receding horizon in section 2.2, formulate an optimal control problem in section 2.3, and present a realtime MPC algorithm based on a suboptimal nonlinear programming method in section 2.4. In this paper, we describe two very similar suboptimal MPC formulations. Sections 2.1−2.4 are equally applicable to both formulations. Sections 2.5 and 2.6 describe the differences between the two formulations. The first formulation, which will be explained in section 2.5, is inspired by DeHaan and Guay,19 who proposed a continuous-time MPC method for a lumpedparameter system. This approach facilitates a closed-loop stability proof but requires a specific formulation of the objective functions, which can limit the performance of the control system. The second approach, which will be explained in section 2.6, involves an ad-hoc formulation of the objective functions and is inspired by economic MPC strategies.12,20 2.1. System. The control system uses the time grid tk = kΔt with k = 0, 1, 2, ..., and the sampling time Δt, which is constant in this paper. The prediction model is integrated on this time grid and the MPC is executed at the sampling points tk. Consider a discrete-time model of the form 0 = fk(x k + 1, x k , uk)
example with Nk = 5. The planned or predicted state and input trajectories during the horizon are summarized in the vectors xk̅ = [xTk,0, ..., xTk,Nk]T and uk̅ = [uTk,0, ..., uTk,Nk−1]T, respectively. Let xk,0 = xk be the initial state at the beginning of the horizon. Throughout this paper, we assume that this initial state xk is known with sufficient accuracy and can thus be used for feedback. In practice, xk has to be measured or estimated, e. g., by Kalman filter techniques for nonlinear systems,21 particle filters,22 or moving horizon estimation.23 In many MPC formulations, xk comes from a prediction made before the time tk to allow for the non-negligible time span required to execute the MPC control law that defines the input uk,0.24,25 This strategy ensures a systematic consideration of the computational delay. Consider that the input during the horizon [tk, tk+Nk] is defined by the input parameter vector Uk = [UTk,0, ..., UTk,M−1]T. Its subvectors Uk,i with i = 0, ..., M − 1, the dimension m, and the constant M define a piecewise constant input on a coarse local time grid Tk,0, ..., Tk,M. That is, the parameters Uk,i define the constant input values during the interval [Tk,i, Tk,i+1) with i = 0, ..., M − 1. The coarse time grid is called local because it is only defined during the respective prediction horizon and it is generally updated for each prediction horizon. As indicated in Figure 1 for M = 2, the coarse local time grid does not need to be equidistant. The coarse local time grids are chosen so that each sampling point Tk,i of a coarse grid is also a sampling point of the fine time grid and that {Tk,1, ..., Tk,M−1} ⊂ {Tk+1,0, ..., T k+1,M−1 } is satisfied. Hence, there exists a matrix Mk ∈ mNk × mM with Mk,ij ∈ {0, 1} and ∑mM j=1 Mk,ij = 1 so that u̅k = MkUk maps the input parameter vector Uk to the input vector u̅k. In MPC, it can be useful to constrain the predicted state xk,Nk at the end of the receding horizon. For instance, this can be a requirement for closed-loop stability. Therefore, we stipulate x k , Nk ∈ ?f , where ?f ⊆ ? is referred to as terminal constraint set. We assume that x s ∈ ?f holds for the target steady state xs. 2.3. Optimal Control Problem. Consider the optimal control problem
(1)
n
where x k ∈ ? ⊆ is the state at the time tk, x 0 ∈ ? is the initial state, and uk ∈ < ⊆ m defines the input during the interval [tk, tk+1). Let fk : n × n × m → n be sufficiently smooth in its arguments and consider that xs is a steady state associated with the constant input us, i.e., 0 = fk(xs,xs,us). The primary control objective is that the system converges to a suitable steady state. The selection of a target steady state may be influenced by economic considerations. The states are restricted to the closed and simply connected set ? . Similarly, the inputs are restricted to the closed and B
DOI: 10.1021/acs.iecr.6b00592 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research min J(x k , Uk)
Uk ∈ < M
3. Repeat the following step of an iterative optimization routine for a fixed number of times: • Compute the gradient gk = dJ/dUk. • Compute the step size
(2a)
s.t. 0 = fk + j(x k , j + 1, x k , j, uk , j)
∀ j = 0, ..., Nk − 1
(2b)
x k ,0 = x k
(2c)
u̅ k = MkUk
(2d)
αk = arg min J(x k , Proj(Uk − α Hgk)) α∈ ≥ 0
• Set Uk ← Proj(Uk − αkHgk). 4. Apply the input uk,0 according to (2d) to the plant. Here, H is a user-defined gain matrix, which is typically an approximation of the inverse Hessian. Moreover, κ: ?f → < is a user-defined discrete-time feedback law and ΔT is the corresponding constant sampling time. To satisfy the conditions described in section 2.2, ΔT must be an integer multiple of Δt, i.e., ΔT= KΔt with the constant K ∈ . If K = 3 and M = 2 as in Figure 1, Nk takes the values 6, 5, 4, 6, etc. During the computation of gk and αk in step 3 of the MPC algorithm, the equality constraints 2b−2d have to be considered. Because adjoint variables are used in this paper and because 2c and 2d are eliminated by substitution, the analytical computation of gk is straightforward.28 The operator Proj: mM → < M is defined in the form
with Nk − 1
J(x k , Uk) = V (x k , Nk) +
∑ j=0
(3)
vj(x k , j, uk , j) (2e)
V: n → ≥ 0 and vj : n × < → ≥ 0 are user-defined objective functions. Generally, formulation 2 requires additional measures to ensure that x k ∈ ? and x k , Nk ∈ ?f . In this paper, x k , j ∈ ? is ensured by requiring J(xk,Uk) < ∞ and by choosing vj to be a barrier function for ? , i.e., vj → ∞ as xk,j approaches ∂? . Similarly, V is formulated as a barrier function for the constraint x k , Nk ∈ ?f . When solving (2), we consider equality constraints 2b by using Lagrange multipliers λk,j (adjoint variables) with j = 0, ..., Nk − 1. Then the Lagrangian follows in the form L(xk,Uk) = Nk−1 T J(xk,Uk) + ∑j=0 λk,j fk+j(xk,j+1, xk,j, uk,j). Alternatively, if constraints 2b can be rewritten in an explicit form, they may also be eliminated by recursive substitution of xk,j. However, in section 3, we use an implicit method for time integration implying that (2b) cannot be written in an explicit form. The equality constraints 2c and 2d are eliminated by straightforward substitution. To keep the computational load as low as possible, we refrain from an exact solution of (2). Alternatively, we use a suboptimal solution strategy, which significantly reduces the computational load and ensures that the controller can be executed in real-time. In fact, we perform a single step or several steps of an iterative optimization routine. This reduces J, which is the main property required to achieve the desired control performance. We assume that an initial solution U0 ∈ < M satisfying all state constraints as well as a corresponding coarse local time grid T0,0, ..., T0,M exist and are known at the time t0 = T0,0. This assumption is commonly used in suboptimal MPC.26,27 2.4. Real-Time MPC Algorithm. At each sampling point tk, the following suboptimal MPC algorithm is executed. Unknown parameters and functions used in this algorithm are explained afterward. 1. If k = 0, use • the initial solution U0 and • the corresponding coarse local time grid T0,0, ..., T0,M else if tk = Tk−1,1, set • Tk,j = Tk−1,j+1 ∀ j = 0, ..., M − 1 • Tk,M = Tk,M−1 + ΔT • Uk,j = Uk−1,j+1 ∀ j = 0, ..., M − 2 • Uk,M−1 = κ(xk−1,Nk−1) else set • Tk,0 = tk • Tk,j = Tk−1,j ∀ j = 1, ..., M • Uk,j = Uk−1,j ∀ j = 0, ..., M − 1. 2. Measure or estimate the current state xk.
Proj(Uk) = arg minM ∥U − Uk∥1 U∈
0 ∀ xk,Nk ≠ xs, v(xs, uk,j) = 0, and v(xk,j, uk,j) > 0 ∀xk,j ≠ xs for arbitrary uk , j ∈ < . Clearly, these properties are equivalent to positive definiteness of V and v in terms of the deviations xk,Nk − xs and xk,j − xs, respectively. V and v have to be chosen in accordance with the feedback controller κ from section 2.5.1. In fact, it is a sufficient requirement for closedloop stability that for any x k ∈ ?f
Δuk = κ(x k) − us = K(x k − x s) = KΔx k
and the positive definite weighting matrix P so that Ṽ (Δxk) = ΔxTk PΔxk is a Lyapunov function candidate. For the linear closed-loop system 6 with 7, Ṽ (Δxk) represents the infinitehorizon cost function and the identity Ṽ (Δx k + 1) − Ṽ (Δx k) ≤ −Δx kT(Q + KTRK)Δx k
∑ v(xk+ j, κ(xk)) − V (xk) ≤ 0
⎛ Ṽ (Δx) ⎞2 ⎛ Ṽ (Δx) ⎞ ⎟⎟ ⎟⎟ ⎜⎜3 − 2 b(Δx) = ⎜⎜ β1 ⎠ ⎝ β1 ⎠ ⎝
(4)
j=0
with 0 = fk + j(x k + j + 1, x k + j , κ(x k))
it follows from (8) and the monotonicity of b with respect to Ṽ that
∀ j = 0, ..., K − 1
and ΔT = KΔt holds. From these requirements, it follows that J(xk, Uk) is indeed a suitable Lyapunov function candidate and an upper bound on the objective function value of a fictitious infinite prediction horizon. 2.5.3. Stability Analysis. In this section, we outline the proof of stability of the nominal closed-loop system.33 The objective function J from (2e) with the functions v and V according to section 2.5.2 is used as a discrete-time Lyapunov function candidate and the Lyapunov difference is shown to be negative definite. First, we consider the case tk = Tk−1,1 of the MPC algorithm from section 2.4: Then, the interval [tk−1, tk) is truncated from the receding horizon and the interval (Tk−1,M, Tk,M] is appended. The latter does not increase the objective value J because of condition 4. Moreover, steps 3−4 of the algorithm given in section 2.4 do not increase the value of J. Consequently, we have the inequality J(x k , Uk) − J(x k − 1, Uk − 1) ≤ −v(x k , uk)
⎛ ⎞ β2 ⎟⎟ b(Δx k + 1)⎜⎜Ṽ (Δx k + 1) + β1 − Ṽ (Δx k + 1) ⎠ ⎝ ⎞ ⎛ β2 ⎟⎟ − b(Δx k)⎜⎜Ṽ (Δx k) + β1 − Ṽ (Δx k) ⎠ ⎝ ≤ −b(Δx k)Δx kT(Q + KTRK)Δx k
holds for all Δx k + x s ∈ ?f and any constant β2 > 0. This motivates the choice ⎛ ⎞ β2 ⎟⎟ V̅ (Δx) = b(Δx)⎜⎜Ṽ (Δx) + β1 − Ṽ (Δx) ⎠ ⎝ v ̅ (Δx , Δu) = b(Δx)(Δx TQΔx + Δu TRΔu)
and demonstrates that V̅ (Δxk) is a Lyapunov function for the linear closed-loop system 6 with 7, as long as Δx k + x s ∈ ?f . However, it has to be checked whether V̅ (Δxk) is also a Lyapunov function for the nonlinear closed-loop system 1 with 7, in the whole domain Δx k + x s ∈ ?f . If this is not the case, a different choice of the parameters β1,β2, and ρ may help. If this is the case and if κ: ?f → < , the requirements of sections 2.5.1 and 2.5.2 are satisfied with the cost functions V(xk) = V̅ (xk − xs) and v(xk, uk) = v(xk − xs, uk − us). Clearly, V(xk,Nk) is a barrier function for the terminal constraint x k , Nk ∈ ?f . In the same way, v(xk,j, uk,j) can be extended by a barrier function term for the constraint x k , j ∈ ? if ? ≠ n. It is recommended to design this term to be 0 in the domain ?f and strictly positive in the domain ?\?f . Otherwise, it has to be checked whether the requirements of sections 2.5.1 and 2.5.2, and especially condition 4, are satisfied. In this section, the choice K = 1 was made for the sake of clarity. It allowed us to use the discrete-time LQR synthesis with the sampling time ΔT = Δt. For K > 1, we can use an
(5)
Recall that both J and v are positive definite in terms of the control error xk − xs. Inequality 5 shows that the Lyapunov difference is negative definite. In the second case tk < Tk−1,1, we follow the same line of reasoning. The only difference compared to the first case is that Tk−1,M = Tk,M, i.e., just the interval [tk−1, tk) is truncated from the receding horizon. This shows that inequality 5 is also applicable in the second case, which concludes the proof. 2.5.4. Design Example of Feedback Controller and Objective Functions. This section describes a strategy for designing κ, V, v, and ?f if ? = n and K = 1. Linearization of (1) at the point xs with the corresponding input us yields Δx k + 1 = AΔx k + BΔuk
(8)
proves closed-loop asymptotic stability. For ρ = 1, (8) holds with identity. L e t ?f b e d e fi n e d a s t h e l e v e l s e t ?f = {x ∈ ?|Ṽ (x − x s) ≤ β1} with an appropriate constant β1 > 0. If we introduce the function
K−1
V (x k + K ) +
(7)
(6)
with Δxk = xk − xs and Δuk = uk − us. Using the stage cost ΔxTk QΔxk + ΔuTk RΔuk with positive definite weighting matrices Q and R, a discrete-time LQR steady-state optimal controller34 can be designed for system 6 and used as a feedback controller D
DOI: 10.1021/acs.iecr.6b00592 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research emulation design method,34 i.e., a discrete-time LQR synthesis based on an equivalent cost function evaluated with the sampling time ΔT = KΔt. 2.5.5. Discussion. The MPC formulation proposed in this section features provable nominal closed-loop stability, which is clearly a desirable property. However, it entails also some disadvantages. Finding a feedback law κ that satisfies all stipulated requirements, especially (4) and the forward invariance property of ?f , can be a challenge. In fact, the choice of ?f , V , and v depends on the design of κ. In section 2.5.4, ?f is chosen to be a level set of a Lyapunov function used in the stability proof of κ. The converse Lyapunov theorem35 guarantees the existence of such a Lyapunov function. These dependencies may interfere with the required barrier shapes of V and v and their user-defined design to pursue specific control objectives. Using a small region ?f , which is desirable in terms of both control performance and satisfaction of (4), can render it a challenging task to find a feasible initial solution when iteratively solving (2). In the following section, we circumvent these problems by formulating an ad-hoc MPC algorithm. A rigorous proof of closed-loop stability of the proposed ad-hoc approach is beyond the scope of this study. 2.6. Ad-hoc MPC Formulation. Inspired by economic MPC strategies,12,20 we outline a real-time MPC formulation with additional freedom in the design of ?, ?f ∈ ?, < , V , and vj, e.g., there is no need for ?f to be a level set. This additional freedom should be utilized for better satisfaction of (economic) control objectives. In many cases, this entails the disadvantage that there is no rigorous proof of closed-loop stability.12 The proposed real-time MPC formulation uses the general algorithm given in section 2.4. The domains ?, ?f ∈ ? , and < can be freely chosen depending on the respective control problem and economic objectives. Note that there is generally no need for ?f to be a level set. The objective functions V and vj can also be freely chosen, except that they have to be barrier functions for the constraints x k , j ∈ ? and x k , Nk ∈ ?f . In cases where it is known that some steady state xs is the optimal operational strategy, V and vj can be designed to have their minimum at xs. However, in economic MPC, the objective is often to maximize some output quantity without a-priori knowledge of xs. Recall that the control law κ in the MPC algorithm from section 2.4 is not directly used for feedback but just to generate a feasible initial guess for Uk,M−1. There are various possibilities for the design of κ and ΔT: First, we could use κ = 0. Second, we could use κ = us. Third, we could use κ = Uk,M−2. As a fourth alternative, we could use some feedback control law depending on the prediction xk−1,Nk−1. This feedback controller could be designed based on linearized system 6 using any standard method from linear control theory.34 Clearly, for κ to yield a feasible guess for Uk,M−1, it is important that ?f is a forward invariant set of system 1 with the feedback κ. A proof of closed-loop stability of this ad-hoc MPC formulation has not been found so far. Therefore, it has to be individually decided whether it is suitable to operate a system with this control approach. An appropriate choice of ?, ?f , and < can at least ensure safe operation of the closedloop system.
3. TRANSPORT−REACTION SYSTEM We consider a transport−reaction system modeled by the hyperbolic partial differential equation ∂ξ(z , t ) ∂ξ(z , t ) + w(t ) − r(ξ(z , t ), u(t )) = 0 ∂t ∂z
(9)
with the state ξ(z, t), the spatial coordinate z ∈ [0, L], the time t ≥ 0, a single characteristic curve (defined by dz/dt = w(t)), and a generally nonlinear reaction term r. We assume that the velocity w(t) satisfies w(t) > 0 at all times. Equation 9 is supplemented by the initial condition ξ(z , 0) = ξ0(z)
∀ z ∈ [0, L]
(10a)
and the boundary condition ξ(0, t ) = ξb(t )
∀t≥0
(10b)
Here, ξb(t) can be considered as a system input, which may or may not be controllable. Consider an equidistant spatial grid zl = lΔz with l = 0, ..., Nz and the step size Δz = L/Nz. To integrate (9), we use the modif ied method of characteristics36 combined with the trapezoidal rule37 for time integration. The approach is also known as two-time-level semi-implicit semi-Lagrangian scheme.38 This approach features high accuracy, low numerical diffusion, no grid distortion, implicit time integration, and benign numerical stability properties. The characteristic curves of (9) are defined by ∫ w(t) dt, which can be approximated in the form ∫ ttk+1 w(t) dt ≈ (w(tk) + k w(tk+1))Δt/2. Based on this formulation, the Courant number9 follows in the form Ck =
w(tk) + w(tk + 1) Δt 2 Δz
Let ξlk be an abbreviation for ξ(zl, tk)with l = 0, 1, ..., Nz and k = 0, 1, 2, ... From the initial and boundary conditions, we have ξl0 = ξ0(zl) and ξ0k = ξb(tk), respectively. Given that ξlk ∀ l = 1, ..., Nz is known from previous calculations, we are interested in ξlk+1 ∀ l = 1, ..., Nz. By back-tracing the characteristic curve that starts at zl and tk+1, we find the so-called departure point zl ̅ = (1 − Ck)zl + Ckzl−1. For simplicity, we assume 0 ≤ Ck ≤ 1, which implies zl ̅ ∈ [zl−1, zl]. Let ξlk̅ = ξ(zl,̅ tk) be the state at the departure point zl,̅ which is normally not a grid point. This is why we apply linear interpolation to obtain the approximation ξkl ̅ = (1 − Ck)ξkl + Ckξkl − 1
based on known grid-point values ξlk and ξl−1 k . Application of the trapezoidal rule 37 (Crank−Nicolson method) for time integration of (9) along characteristic curves and consideration of a piecewise constant input u(t) = uk ∀ t ∈ [tk, tk+1) give 0=
ξkl+ 1 − ξkl ̅ Δt
−
∀ l = 1, ..., Nz
1 (r(ξkl+ 1 , uk) + r(ξkl ̅ , uk)) 2 (11)
Equation 11 is successively evaluated for k = 0, 1, 2, ... N Assembling the states in the vector xk = [(ξ1k )T, ..., (ξk z)T]T and arranging (11) in a system of equations, we obtain formulation 1. The sparse structure of the resulting system can be exploited in computer implementations to improve the computational efficiency. E
DOI: 10.1021/acs.iecr.6b00592 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research ξb = [4 mol/L, 0 mol/L, 320 K]T
4. EXAMPLE PROBLEM 4.1. Control of an Exothermic Plug-Flow Reactor. In this section, we apply the proposed real-time MPC algorithms to a nonisothermal plug-flow reactor with two irreversible exothermic first-order reactions.39 These reactions may cause undesired local temperature elevations, so-called hot spots.40,41 Therefore, an important control objective or constraint is that the local temperature T(z, t) in the reactor never exceeds a certain critical value T̅ , i.e., T (z , t ) ≤ T ̅
(14)
i.e., the incoming fluid does not contain the species B. The value of ξb is constant and cannot be used for control. The reactor operates with a constant flow velocity w, which implies that the Courant number Ck = C is also constant. The nominal model parameters are summarized in Table 1. Table 1. Nominal Parameters of Plug-Flow Reactor
(12)
Karafyllis and Daoutidis16 addressed a similar control task but their method requires regularity assumptions on the boundary condition (cf. (10b)) and may violate (12) during transient operating situations. However, undesirable temperature elevations should be avoided at all times, i.e., also during transient periods. In fact, the control of transitions between different steady-state operating points is a frequent task in continuous production processes. Performing this task in an optimal way is important for the consumption of resources, the production costs, the amount of off-spec products, etc. Bendjaouahdou and Bendjaouahdou18 used generalized predictive control to control hot spots in a fixed bed reactor modeled as a diffusion−convection−reaction system. The model is discretized by finite differences and the Crank− Nicolson method. The formulated optimization problem implements a set-point tracking objective. Constraints are not explicitly considered, i.e., hot spots are indirectly avoided. Lao et al.14 used the Galerkin method to obtain a lumpedparameter model of a tubular reactor and controlled the system by economic MPC with a boundary input and a state constraint to avoid hot spots, i.e., hot spots are directly avoided. This state constraint is formulated as an inequality constraint of the optimization problem. This typically requires an exact solution of the considered optimization problem and can thus be incompatible with the idea of real-time suboptimal MPC. To circumvent this problem, barrier functions are used in the current paper to implement temperature constraint 12. The first considered reaction converts the supplied reactant species A into the desired product B. The second reaction converts the species B into an undesired product. Both reactions are irreversible, exothermic, and first-order. The objective is to produce the desired product B. The system can be modeled in forms 9 and 10 with the state vector ξ = [CA, CB, T]T. It contains the molar concentrations CA and CB of the species A and B, respectively, and the local temperature T of the fluid in the reactor. The reaction term reads as
parameter
variable
value
flow velocity length of reactor reaction constant reaction constant activation energy activation energy enthalpy of reaction enthalpy of reaction universal gas constant volume-specific heat capacity volume-specific heat transfer coefficient maximum reactor temperature
w L k1 k2 E1 E2 H1 H2 Ru ρc h T̅
1 m/min 1m 5 × 1012 L/min 5 × 103 L/min 20000 kcal/kmol 25000 kcal/kmol −0.274 kcal/mol −0.0986 kcal/mol 1.987 kcal/(kmol K) 0.006 kcal/(L K) 0.17 kcal/(min L K) 370 K
In this paper, we consider a transition from an original steady state x0 with u = u0 = 293.15 K and CB(L, t) = 0.03 mol/L (off state of the reactor) to a desired final steady state xd with maximum output concentration CB(L, t). At the desired final steady state xd, the temperature constraint 12, holds with equality. To find x0, the steady-state formulation of the prediction model 11, i.e., 0=w
ξ0l − ξ0l − 1 Δz
+ Cξ0l − 1 , u0))
1 (r(ξ0l , u0) + r((1 − C)ξ0l 2
−
∀ l = 1, ..., Nz
is solved numerically. The desired final steady state xd = xs = [(ξ1s )T, ..., (ξNs z)T]T and the corresponding input ud = us can be easily computed by solving the static optimization problem
max[010]ξsNz
(15a)
us ∈
