Suboptimal Control Strategies for Finite-Time Nonlinear Processes with Input Constraints

Abstract

Novel techniques for the optimization and control of finite-time processes in real-time are pursued. These are developed in the framework of the Hamiltonian optimal control. Two methods are designed. The first one constructs the reference control trajectory as an approximation of the optimal control via the Riccati equations in an adaptive fashion based on the solutions of a set of partial differential equations called the 𝛼 and 𝛽 matrices. These allow calculating the Riccati gain for a range of the duration of the process 𝑇 and the final penalization 𝑆. The second method introduces input constraints to the general optimization formulation. The notions of linear matrix inequalities allow us to recuperate the Riccati gain as in the first method, but using an infinite horizon optimization method. Finally, the performance of the proposed strategies is illustrated through numerical simulations applied to a batch reactor and a penicillin fed-batch reactor.

1. Introduction

Batch processes have received much attention during the past two decades due to developing chemical and pharmaceutical products, new polymers, and recent biotechnological applications. Usually, the engineers state that a batch process has three operative stages: start-up, batch run, and shutdown. While these three stages are widely studied by the engineers for each particular batch process, it is important to note that in a wide number of cases, the batch run stage is far from an optimal operation, supported only with the experience of operators and engineers.

Control techniques have been made inroad in the industry in order to improve the performance of the process. In batch context, interesting ideas such as controllability, observability, and stability have been introduced despite the duration of batch processes. Numerous control techniques have been adapted to the peculiarities of these kinds of systems. One example is the model predictive control (MPC). This method is essentially numerical, usually implemented on-line. Most successful industrial applications of MPC reported so far are in refining and petrochemical plants (see for instance [1, 2]), where processes are run near optimal steady-states and model linearizations are reliable approximations. In spite of its many advantages, some undesired disadvantages (such as the high computational cost) have caused a major problem: only a few of the available commercial software packages are cautiously suggested for truly nonlinear batch processes as it is reported in [3].

Here, the optimal control paradigm is tackled, where the Hamiltonian formalism plays a major role. If the optimal control problem is regular, that is, a unique control 𝑢(𝑡) minimizes the cost function of the problem, then there is an explicit form for the control law. For linear systems, it has been shown that this control strategy is in fact a proportional feedback law 𝑢(𝑡)=𝐾𝑥(𝑡) (see [4] for details). For nonlinear systems, the resulting control strategy is not a strict feedback control law depending on its state. Indeed, it is an open-loop strategy depending on a nonmeasurable co-state of Hamiltonian system associated with the process as it is illustrated in [4, 5].

This problem cannot generally be solved in real-time due to the missing initial condition for the costate. Several numerical approaches have been designed to cope with this problem, and many of them were based on the shooting methods (see [6]), recursive methods with generalized functions (see [7]), and techniques coming from artificial intelligence as the results presented by [8] and their references. The main disadvantage of these techniques, in the context of on-line applications, is the computational effort required. Therefore, these techniques remain as off-line procedures.

When input constraints are imposed to the problem, even in the linear control problem, the control law becomes nonlinear and phenomena like saturation appear. Many methods deal with this problem, for instance, control parametrization introduced by [7], neural networks (see [9] and their references), and other numerical approaches described by [10]. However, few of them are designed for on-line applications. Another disadvantage presented in these methods can be the saturation of the actuators. For many batch processes, this situation is not desired.

MPC strategies and other techniques based on nonlinear programming schemes are proposed in the literature to solve both optimal problems in real-time (see for instance [11, 12]). However, the model process is converted to discrete time to pose the nonlinear programming problem and a large computation is required at each sampling time. In this paper, the continuous differential equations are maintained, analytical PDEs are used to calculate the control law, and lower computation requirement is needed in comparison with these methods. Recently, an analytical method to solve these problems was proposed by [13]. This approach solves efficiently the optimal control problem but taking into account just the optimal objectives without considering the dynamic behavior of the solution. Here, the controllers obtained have a simple form (similar to PID controllers), and its dynamics can be adjusted to obtain smooth behaviors with a good performance. An interesting survey paper written by [14] can be seen for obtaining more details about the methods to solve both problems described before either in an off-line or in a real-time framework.

The contributions of this paper are as follows. (i) A method to generate control strategies that can be used to obtain reference trajectories (or batch recipes) by means of a linear approximation of the nonlinear optimal solution using auxiliary matrices 𝛼(𝑇,𝑆), 𝛽(𝑇,𝑆) solutions of a set of partial differential equations (PDEs). The time-varying Riccati gain is calculated in parallel with the process based on 𝛼(𝑇−𝑡,𝑆), 𝛽(𝑇−𝑡,𝑆). (ii) A method to solve the constrained version of the optimization problem when restrictions in the manipulated variable are taken into account, exploiting the ideas of LMIs. A new LMI is proposed which can be designed to avoid saturation by tuning a parameter called 𝛾. Although the control strategy is completely calculated and adapted for a finite-time horizon problem, the optimization problem is solved in an infinite horizon fashion. (iii) The controllers obtained in both cases have a simple form and can be implemented in real-time. Numerical procedures to implement the strategies are provided.

The rest of the paper has the following structure. In Section 2, a brief description of the Hamiltonian formalism is elaborated. Then the first method is proposed based on a set of PDEs and their relationships with the Riccati equations. Next, the restricted optimization problem is analyzed through the lines of linear matrix inequalities and optimal control. In Section 3, examples of batch processes are discussed. The strategies are reviewed by means of numerical simulations. Finally, in Section 4 the conclusions are presented.

2. Model Formulation and Solution Methodology

The dynamics of a fed-batch reactor usually is a nonlinear control system, which can be written in the affine form as follows: 𝑥̇(𝑡)=𝑓(𝑥(𝑡))+𝑔(𝑥(𝑡))𝑢(𝑡);𝑥(0)=𝑥₀, (1) where 𝑥∈ℝ^𝑛 is the vector of state variables with initial conditions 𝑥₀, 𝑢∈𝒰∈ℝ^𝑚 is the control vector, and 𝑓(𝑥)∈ℝ^𝑛 and 𝑔(𝑥)∈ℝ^𝑛×𝑚 are smooth vector fields. The general dynamic optimization of a bioprocess needs to find the control vector 𝑢, which minimizes a performance index as 𝒥𝑇,0,𝑥₀,𝑢(⋅)=^𝑇₀𝐿(𝑥(𝜏),𝑢(𝜏))𝑑𝜏+𝒦(𝑥(𝑇)) (2) subject to the dynamics (1). The function 𝐿 is known as the Lagrangian, 𝒦 is a final penalization on the states, and 𝑇 is the duration of the batch run. Although 𝐿 and 𝒦 can be general functions, here only quadratics forms will be considered, namely, 𝐿(𝑥,𝑢)=𝑥^′𝑄𝑥+𝑢^′𝑅𝑢,𝒦=𝑥(𝑇)^′𝑆𝑥(𝑇), (3) where 𝑄 and 𝑆 are 𝑛×𝑛 symmetric and positive semidefinite matrices and 𝑅 is a 𝑛×𝑚 positive-definite matrix. Under these conditions the Hamiltonian ℋ of such problem defined as ℋ(𝑥,𝜆,𝑢)≜𝐿(𝑥,𝑢)+𝜆^′𝑓(𝑥)+𝑔(𝑥)𝑢 (4) is known to be regular; that is, there is a unique 𝐻-optimal control 𝑢⁰ that minimizes the Hamiltonian ℋ, 𝑢⁰(𝑥,𝜆)≜argmin_𝑢ℋ(𝑥,𝜆,𝑢), (5) and the derivative of ℋ with respect to 𝑢 vanishes at (𝑥,𝜆,𝑢⁰(𝑥,𝜆)). For control affine systems (1), the optimal control takes the following form: 𝑢⁰()1𝑥,𝜆=−2𝑅⁻¹𝑔(𝑥)^′𝜆. (6) The control-Hamiltonian ℋ⁰(𝑥,𝜆)≜ℋ𝑥,𝜆,𝑢⁰(,ℋ𝑥,𝜆)₀(𝑥,𝜆)=𝑥^′𝑄𝑥+𝜆^′𝑓(𝑥)−14𝜆^′𝑔(𝑥)𝑅⁻¹𝑔(𝑥)^′𝜆, (7) gives rise to the Hamiltonian canonical equations (HCEs) [4]. Consider 𝑥̇(𝑡)=𝜕ℋ⁰𝜕𝜆^′;𝑥(0)=𝑥₀,𝜆̇(𝑡)=−𝜕ℋ⁰𝜕𝑥^′;𝜆(𝑇)=2𝑆𝑥(𝑇). (8)

2.1. Optimal Control for Linear Systems

In [4, 15] it is shown that the Hamiltonian canonical equations (HCEs) present several inconveniences when they are integrated. For instance, the adequate initial conditions of the costates are not known, and even if there are methods to obtain them, the HCEs are an unstable system. Besides, to construct the optimal trajectory with (6) for systems with many states can be difficult, and the optimal control problem becomes an unsolved problem. In contrast, to calculate the optimal control 𝑢⁰ for a linear system is easier than the nonlinear case. The main idea is then to use the tools from linear systems to approximate the feed-forward control equation (6) to a linear feedback law as follows.(i)Only when the control action is calculated, the system considered will be a linear approximation of (1) around a small neighborhood of (𝑥, 𝑢), that is, 𝑥̇(𝑡)≈𝐴𝑥(𝑡)+𝐵𝑢(𝑡),𝑥(0)=𝑥₀,𝐴=𝜕𝑓(𝑥)+𝜕𝑥𝜕𝑔(𝑥)𝑢𝜕𝑥_𝑥,𝑢,𝐵=𝑔(𝑥)_𝑥,𝑢, (9) where (𝑥, 𝑢) is a fixed point of the trajectory of the system (1) at time 𝑡≥0.(ii)The costates 𝜆(𝑡) can be linearized in a small neighborhood of (𝑥,𝜆); that is, 𝜆(𝑡)≈𝜆(𝑡)=2𝑃(𝑡)𝑥(𝑡). As a result, the optimal control law equation (6) becomes 𝑢≈𝑢=−𝑅⁻¹𝐵^′𝑃(𝑡)𝑥(𝑡), (10) where 𝑃(𝑡) is the solution to the differential Riccati equation (DRE): 𝑃̇(𝑡)=𝑃(𝑡)𝑊𝑃(𝑡)−𝑃(𝑡)𝐴−𝐴^′𝑃(𝑡)𝑃(𝑇)−𝑄,=𝑆,𝑊=𝐵𝑅⁻¹𝐵^′. (11)Notice that, in this paper, 𝑢 is the manipulated variable corresponding to nonlinear problem and 𝑢 is the approximated linear control law.

The main challenge here is either to solve the HCEs (8) corresponding to the linear problem (equivalently to find the missing conditions 𝑥(𝑇) and 𝜆(0)) or to calculate the solution of differential Riccati equation 𝑃(𝑡).

2.2. Method 1: Generating of Suboptimal Control Reference Trajectory

In order to find the missing boundary values 𝑥(𝑇) and 𝜆(0), the classical approach presented in [4, pages 370-371]) is revised. This method transforms the original boundary-value problem (8) into an initial-value one, by introducing the following auxiliary objects:(i)the Hamiltonian matrix 𝐇, (which comes from (8) applied to the system (9)) 𝐇=^def1𝐴−2𝑊−2𝑄−𝐴^′, (12)(ii)the augmented Hamiltonian system (a linear 2𝑛-dimensional matrix ODE with a final condition) defined for the combination of two 𝑛×𝑛 matrices 𝑋(𝑡),Λ(𝑡),𝑡∈[0,𝑇] through 𝑋̇Λ̇𝑋Λ;𝑋(𝑇)Λ(𝑇)=𝐼=𝐇2𝑆. (13)The solution to the augmented (linear) problem (13) being unique must verify 𝑋(0)Λ(0)=𝑒^−𝐇𝑇𝐼2𝑆, (14) and since in this linear case equation (8) also reads as follows: 𝑥̇𝜆̇𝑥𝜆=𝐇, (15) then the missing boundary conditions can be explicitly found, namely, 𝑥₀𝜆(0)=𝑒^−𝐇𝑇𝑥(𝑇)(𝑇)=𝑋(0)Λ(0)𝑥(𝑇)𝑥(𝑇)2𝑆𝑥,=𝑋⁻¹(0)𝑥₀,𝜆(0)=Λ(0)𝑋⁻¹(0)𝑥₀, (16)(17) (see [4] for a proof of the invertibility of 𝑋 and other details). Actually, the whole solution to the DRE can be recovered from the solution to problem (13), explicitly, as follows: 𝑃(𝑡)=Λ(𝑡)𝑋⁻¹(𝑡),𝑡∈[0,𝑇]. (18) As it will become clear in the following sections, it is desirable to have access to the missing boundary values for different values of the parameters (𝑇,𝑆) without solving either the DRE or the HCEs. To solve the whole (𝑇,𝑆)-family of LQR problems (with common 𝐴,𝐵,𝑄,𝑅,𝑥₀ values) consider each particular (𝑇,𝑆)-problem: 𝜌(𝑇,𝑆)=^def𝑥^𝑇,𝑆(𝑇);𝜎(𝑇,𝑆)=^def𝜆^𝑇,𝑆(0);𝑈(𝑇)=^def𝑒^−𝐇𝑇,𝛼()𝛽()𝑇,𝑆𝑇,𝑆^def=𝑋^𝑇,𝑆(0)Λ^𝑇,𝑆(0), (19)(20)(21) where the superscript (^𝑇,𝑆) refers to the individual problem of duration 𝑇 and final penalty matrix 𝑆. Then (14) can be rewritten in the form 𝛼()𝛽()(𝑇)⋅𝑇,𝑆𝐼𝑇,𝑆=𝑈2𝑆. (22) Since the subjacent Hamiltonian system is linear, its solutions depend smoothly on parameters and initial conditions, and then the (partial) derivative of (22) with respect to 𝑇 (denoted by the subscript (_𝑇)) results in 𝛼_𝑇𝛽_𝑇𝐼=−𝐇𝑈2𝑆, (23) where the variables (𝑇,𝑆) have been dropped for clarity here and whenever convenient and clear through the rest of the paper. At this point, it is convenient to start by adopting the simplest expressions for matrix 𝑆 to make sense of the partial derivatives (𝜕/𝜕𝑆) with respect to the second parameter. Consequently, we will consider for simplicity that the form of the final penalty matrix will be scalar; that is, 𝑆=𝑠𝐼, 𝑠∈ℝ⁺₀. The partial derivatives of (22) with respect to the real variable 𝑠 are then 𝛼_𝑠𝛽_𝑠0=𝑈2𝐼. (24) Now, by partitioning 𝑈 into 𝑛×𝑛 submatrices 𝑈_𝑖,𝑖=1,…4 in the form 𝑈𝑈=¹𝑈₂𝑈₃𝑈₄, (25) then (24) reads as follows: 12𝛼_𝑠=𝑈₂,12𝛽_𝑠=𝑈₄, (26) which combined with (22) gives 𝑈₁=𝛼−𝑠𝛼_𝑠,𝑈₃=𝛽−𝑠𝛽_𝑠, (27) and by inserting these results in (23), the following relations are obtained: 𝛼_𝑇−𝛼_𝑠𝛽𝑀=−𝛼𝑁,_𝑇−𝛽_𝑠𝑀=−𝛽𝑁, (28) where 𝑀=^def𝐴^′𝑆+𝑆𝐴+𝑄−𝑆𝑊𝑆,𝑁=^def𝐴−𝑊𝑆. (29) These are first-order quasilinear PDEs for the matrices 𝛼,𝛽. Boundary conditions for a process of zero horizon are imposed in view of (21) and (17), that is, 𝛼(0,𝑠)=𝐼,𝛽(0,𝑠)=2𝑆. (30) Now, the missing values, the final state and initial costate, for any (𝑇,𝑠)-problem can be recuperated immediately as 𝜌=𝛼⁻¹𝑥₀,𝜎=𝛽𝛼⁻¹𝑥₀=𝛽𝜌. (31) Concerning the existence and uniqueness of solutions to the quasilinear PDEs equation (28), there are only local results (see [16, page 51]), although the field of vector and matrix PDE integration is in active development (see for instance [17]). The solutions 𝛼 and 𝛽 illustrated before can be alternatively obtained from the procedure detailed by [18, 19] for the one- and 𝑛-dimensional cases, respectively. The advantage of the procedure developed here is that it is not necessary to compute the PDEs for the final state 𝜌 and the initial coestate 𝜎 as it is required in [18, 19].

2.2.1. Relations amongst (𝛼, 𝛽) Solutions and Riccati Equations

For any (𝑇,𝑆)-LQR problem the optimal initial costate has the form given by (10) 𝜎=𝜆^𝑇,𝑆(0)=2𝑃(0)𝑥₀, (32) where 𝑃(⋅) is the numerical solution of the corresponding DRE, that is, the final value matrix ODE in (11). Therefore, from (31) and (32), for each (𝑇,𝑆)-problem the Riccati matrix 𝑃(⋅) should also verify 𝑃(0)=12𝛽(𝑇,𝑆)[𝛼(𝑇,𝑆)]⁻¹, (33) allowing to solve the DRE as an initial-value problem (eventually on-line with the dynamics). But actually, as it is shown below (see (37)) the (𝛼, 𝛽) matrices for missing boundary conditions avoid solving the DRE for each particular (𝑇,𝑆)-problem. They also avoid storing, necessarily as an approximation, the Riccati matrix 𝑃(𝑡) for the values of 𝑡∈[0,𝑇] chosen by the numerical integrator, possibly different from the time instants for which the control 𝑢(𝑡) is constructed. Instead, the HCEs (8) can be integrated with initial conditions 𝑥(0)=𝑥₀,𝜆(0)=𝜎(𝑇,𝑆), (34) rendering the optimal trajectories 𝑥^∗(𝑡),𝜆^∗(𝑡) for 0≤𝑡≤𝑇, which allow to generate the optimal control at each sampling time 𝑢^∗(𝑡)=𝑢⁰𝑥^∗(𝑡),𝜆^∗(𝑡1)=−2𝑅⁻¹𝐵^′𝜆^∗(𝑡). (35) Even better, the feedback form for the control becomes directly available due to the linear dependence of (31) on initial conditions, namely 𝑢^∗(𝑡)1=−2𝑅⁻¹𝐵^′𝛽(𝑇−𝑡,𝑆)[𝛼(𝑇−𝑡,𝑆)]⁻¹𝑥^∗(𝑡), (36) since 𝜆^∗(𝑡) is also the optimal initial costate for the problem with parameters (𝑇−𝑡,𝑆) starting at 𝑥(0)=𝑥^∗(𝑡). Then, as a side-product, an alternative formula for the Riccati matrix results: 𝑃(𝑡)=12𝛽(𝑇−𝑡,𝑆)[𝛼(𝑇−𝑡,𝑆)]⁻¹,∀𝑡∈[0,𝑇]. (37)

2.2.2. Numerical Method for Nonlinear Systems

Now, a suboptimal control strategy for nonlinear systems can be developed based on the linearization (9) around a chosen pair (𝑥,𝑢). After that, the DRE is solved by the matrices 𝛼, 𝛽 according to (37) for a range of 𝑇 and 𝑆, and its solution is kept in memory.

Normally, this linearization is performed around an equilibrium point (or an steady-state point). And, this will work well near that point. However, the initial condition is far from the equilibrium and in fed-batch rectors equilibrium (or steady state) does not even exist as in continuous reactors. Here, the linearization will be calculated around the initial condition, so updating it in the running time is recommended.

In the following procedure, a basic adaptive scheme to update the controller gain as well as the matrices 𝐴, 𝐵 at arbitrary time instants 𝜏_𝑖, 𝑖=0,1,…,𝑛 is given:In: 𝑥(0), 𝑄, 𝑅, 𝑇, 𝑆, {𝜏₀,𝜏₁,…,𝜏_𝑛},Out: 𝑢(𝑡), for all 𝑡≥0. Step 1: Take 𝑖=0. At time 𝑡=𝜏₀=0, the first pair (𝐴,𝐵) is computed from (9) taking 𝑥=𝑥(0) and 𝑢=0 or as solution of 𝑓(𝑥(0))+𝑔(𝑥(0))𝑢=0. Step 2: Compute (𝛼,𝛽) using (22) or as a solution of (28) by mean of the two matrices (𝐴,𝐵), and Step 3: run in parallel the process sending the control signal 𝑢(𝑡) = (−1/2)𝑅⁻¹𝐵^′𝛽(𝑇−𝑡,𝑆)[𝛼(𝑇−𝑡,𝑆)]⁻¹𝑥(𝑡) during the interval 𝜏_𝑖≤𝑡≤𝜏_𝑖+1. Step 4: at time instant 𝑡=𝜏_𝑖+1, update the matrices 𝐴, 𝐵, 𝛼, and 𝛽 of the control signal of Step 3 taking 𝑥=𝑥(𝜏_𝑖+1) and 𝑢=𝑢(𝜏_𝑖+1). Update 𝑖=𝑖+1. Step 5: return to Step 3. This procedure continues until completing the whole process horizontally, that is, until to reach 𝑡=𝑇=𝜏_𝑛. Step 6: if 𝑡=𝑇, then the process is over.

Notice that 𝑃(𝑡)=(1/2)𝛽(𝑇−𝑡,𝑆)[𝛼(𝑇−𝑡,𝑆)]⁻¹ describes a time-variant feedback control law, where the matrices 𝐴, 𝐵, 𝛼, and 𝛽 are updated at each instant 𝜏_𝑖. In the scheme above, the time instants (𝜏_𝑖) are chosen arbitrarily. More sophisticated algorithms to detect the time of updating the linearization can be found in the specialized literature of adaptive control, but that is out of scope of this paper.

2.3. Method 2: Handling Input Constraints

The goal in this section is to develop a suboptimal feedback control for the nonlinear systems (1) but taking into account restrictions in the manipulated variable. This last objective cannot be included in Method 1 due to the lack of efficient tools for on-line optimization. However, in the last years, numerical approaches (as neural networks (see [9])) have been introduced as alternative tools to cope with this problem in the off-line context.

Here, an on-line procedure to generate suboptimal feedback control laws is developed. To incorporate restrictions, the manner to solve the optimization problem was changed, basically modifying its horizon to infinite (see [12]). So, the optimization problem to solve reads as follows: min_𝑢∈𝒰𝒥=^∞₀𝑥^′𝑄𝑥+𝑢^′𝑅𝑢𝑑𝜏, (38) subject to dynamics (9) and ‖𝑢(𝑡)‖⩽𝑢_max. (39)

Note that the optimization horizon time was changed to infinity only to evaluate the cost function and to generate the control law, but the problem remains as a finite-time process. The finite-time horizon with input constraints case and on-line applications is still an open problem. Some first results in this context are illustrated in [20].

The solution to the optimal problem control described above is not easy to find, even when the system is linear. For these kinds of systems, the traditional Ricatti equations are no longer valid, and then a number of numerical approaches attempt to solve this problem like the MPC formulation [1] or the methods described in [7].

Linear matrix inequalities (LMIs) are an alternative method to treat the optimization problem posed above with a low computational effort, at least in the infinite horizon set-up and for linear time-invariant systems.

A linear matrix inequality or LMI is a matrix inequality of the form: ℱ(𝜋)≜ℱ₀+^𝑚_𝑖=1𝜋_𝑖ℱ_𝑖>0, (40) where 𝜋∈ℝ^𝑚 and 𝜋_𝑖 are the optimization variables and the positive and symmetric matrices ℱ_𝑖=ℱ^′_𝑖∈ℝ^𝑚×𝑚,  𝑖=1,…,𝑚, are given. The LMI (40) is a convex constraint on 𝜋; that is, the set {𝜋∣ℱ(𝜋)>0} is convex (for more details, it is suggested to see two interesting contributions [21, 22]).

In particular, the following control problems are handled as a single LMI-based optimization problem in this paper.

(i) If 𝑇 is considered as an infinite horizon, it is known that when the control weighting matrix 𝑅 in the cost functional is positive definite, the state weighting matrix 𝑄 is nonnegative definite and the pair (𝐴,𝐵) is controllable; then the LQR problem is well posed and can be solved via the classic algebraic Riccati equation (ARE) (see [4] for more details). This algebraic equation can be converted to a constraint (inequality) [21] as follows: Π𝐴+𝐴^′Π+𝑄−Π𝐵𝑅⁻¹𝐵^′Π≥0, (41) where Π>0 is the solution of the last inequality and Π∈ℝ^𝑛×𝑛 is a positive symmetric matrix. In this section, the Riccati gain 𝑃 used to generate the control action is renamed as Π because it is the solution of the ARE equation.

A common procedure to convert convex nonlinear (matrix norm, quadratic, etc.) inequalities to a LMI form is through the well-known Schur complements; that is, a LMI of the form 𝑍(𝜋)𝐷(𝜋)𝐷(𝜋)^′𝑊(𝜋)>0, (42) where 𝑍(𝜋)=𝑍(𝜋)^′, 𝑊(𝜋)=𝑊(𝜋)^′, and 𝐷(𝜋) depend on 𝜋, is equivalent to the matrix inequality 𝑊(𝜋)>0,𝑍(𝜋)−𝐷(𝜋)𝑊(𝜋)⁻¹𝐷(𝜋)^′>0, (43) or equivalently one has 𝑍(𝜋)>0,𝑊(𝜋)−𝐷(𝜋)^′𝑍(𝜋)⁻¹𝐷(𝜋)>0. (44)

Since 𝑅>0, and by using (44), the algebraic Riccati inequality (41) can be written as a LMI: ℱ¹(Π)=Π𝐴+𝐴^′𝐵Π+𝑄Π𝐵_′Π𝑅>0. (45)

Note that to minimize the optimal control problem without constraints, it is equal to solve the ARE equation. And in [21] it is proved that to find the ARE’s solution is equivalent to minimize the LMI equation (45).

(ii) Physical limitations inherent in process equipment invariably impose hard constraints on the manipulated variable 𝑢(𝑡) as in (39). These constraints can be seen as a LMI.

Considering that 𝑢(𝑡) is a stabilizing control law, then one has the following.

Lemma 1. Assuming that a 𝛾 parameter exists, the following LMI ℱ²(Π)=𝑢²_max𝛾𝐼𝑅⁻¹𝐵^′Π𝑅⁻¹𝐵^′Π^′𝐼≥0 (46) is equivalent to ‖𝑢(𝑡)‖≤𝑢_max,𝑡≥0. (47)

Proof. Consider that ‖𝑢(𝑡)‖²=𝑅⁻¹𝐵^′Π𝑥(𝑡)²=𝑥(𝑡)^′𝑅⁻¹𝐵^′Π^′𝑅⁻¹𝐵^′(𝑡)(𝑡)Π𝑥=𝑥^′𝑁^′𝑁𝑥(𝑡)=𝑥(𝑡)^′𝐷𝑥(𝑡). (48) Since 𝐷 is a symmetric and square constant matrix, it is possible to decompose it into its singular values, and then the following inequality holds: 𝜈_min𝐼≤𝑀𝐷𝑀⁻¹≤𝜈_max𝐼, (49) where 𝜈 is a singular value; then 𝑥(𝑡)^′𝑁^′𝑁𝑥(𝑡)≤𝜈_max𝐼𝑥(𝑡)^′𝑥(𝑡)≤𝜈_max𝜐𝐼≤𝑢²_max𝜈𝐼,_max𝑅⁻¹𝐵^′Π^′𝑅⁻¹𝐵^′Π𝜐≤𝑢²_max𝑅𝐼,₋₁𝐵^′Π^′𝑅⁻¹𝐵^′Π𝜐≤𝑢²_max𝐼. (50) Using (44), ‖𝑢(𝑡)‖≤𝑢_max results equivalently to 𝑢²_max𝛾𝐼𝑅⁻¹𝐵^′Π𝑅⁻¹𝐵^′Π^′𝐼≥0, (51) where 𝛾=1/𝜐 is a real positive constant that satisfies a stabilization condition ‖𝑥(𝑡)‖²≤𝜐 for any 𝑡≥0.

Note that to find 𝜐, it is necessary to involve another procedure to design an invariant ellipsoid to guarantee either ‖𝑥(𝑡)‖²≤𝜐 or 𝑥(𝑡)^′𝑃₁𝑥(𝑡)≤𝜐, where 𝑃₁ is a positive definite matrix. But it requires a higher computational effort. However, a simple criterion can be used to approximate the 𝜐 value; that is, since normally the linear system involved is a dissipative system to open loop, then 𝜐=‖𝑥(0)‖² as it is suggested in [23].

One important advantage of the LMI is that the problem with multiple constrains can be expressed with multiple LMIs; that is, ℱ¹(𝜋),…,ℱ^𝑛(𝜋)>0 as a single LMI as follows: diagℱ¹(𝜋),…,ℱ^𝑛(𝜋)>0, (52) which makes us introduce the input constraints to the optimal problem as ℱ(Π)=Π𝐴+𝐴^′𝐵Π+𝑄Π𝐵00_′Π𝑅0000𝑢²_max𝛾𝐼𝑅⁻¹𝐵^′Π00𝑅⁻¹𝐵^′Π^′𝐼>0. (53)

To summarize, the optimal control problem with input constraints expressed by (38) and (39) is equivalent to the minimizing of the LMI (53) with respect to Π.

The control action that minimizes this problem can be expressed as 𝑢(𝑡)=−𝑅⁻¹𝐵^′Π𝑥(𝑡), (54) where Π is the solution of minimizing of the LMI (53).

The procedure to implement these results in nonlinear systems is as follows.In: 𝑥(0), 𝑄, 𝑅, 𝑇,  𝑢_max, {𝜏₀,𝜏₁,…,𝜏_𝑛}.Out: 𝑢(𝑡), for all 𝑡≥0. Step 1: Take 𝑖=0. At time 𝑡=𝜏₀=0, the first pair (𝐴,𝐵) is obtained either from (9) taking 𝑥=𝑥(0) and 𝑢=0 or as solution of 𝑓(𝑥(0))+𝑔(𝑥(0))𝑢=0. Step 2: Compute the LMI (53) by mean of the two matrices (𝐴,𝐵), Step 3: run in parallel the process sending the control signal 𝑢(𝑡)=−𝑅⁻¹𝐵^′Π𝑥(𝑡) during the interval 𝜏_𝑖≤𝑡≤𝜏_𝑖+1. Step 4: at time instant 𝑡=𝜏_𝑖+1, update the matrices 𝐴, 𝐵, and Π of the control signal of Step 3 taking 𝑥=𝑥(𝜏_𝑖+1) and 𝑢=𝑢(𝜏_𝑖+1). Take 𝑖=𝑖+1. Step 5: return to Step 3. This procedure continues until completing the whole optimization horizon, that is, until to reach 𝑡=𝑇=𝜏_𝑛. Step 6: if 𝑡=𝑇, then the optimization process is over.

The time instants 𝜏_𝑖 when the matrices are updated are taken in the same fashion as Method 1. The numerical methods to solve LMI are well-known and they are already available in a standard software like MATLAB.

3. Numerical Examples

3.1. Nonlinear Batch Reactor

Consider a batch reactor with a nonlinear dynamic where an exothermic and irreversible second order chemical reaction 𝐴→𝐵 takes place. It is assumed that the reactor has a cooling jacket whose temperature (𝑢) can be directly manipulated. The goal is to regulate the reactor temperature (𝑥₂) by means of inlet coolant temperature (𝑢). Furthermore, the manipulated variable has minimum and maximum constraints −10≤𝑢≤10. The concentration is denoted by the state 𝑥₁.

The reactor dynamics is modeled by the following equations: 𝑥̇₁(𝑡)=−𝑘₀𝑒^{𝐸𝑅/𝑥2(𝑡)}𝑥²₁(𝑡),𝑥̇₂(𝑡)=−𝑈𝐴𝑀𝑐_𝑝𝑥₂(𝑡)−()−𝑢+298.16Δ𝐻𝑀𝑐_𝑝𝑉𝑘₀𝑒^{𝐸𝑅/𝑥2(𝑡)}𝑥²₁(𝑡), (55) subject to initial conditions 𝑥₁(0)=0.9 and 𝑥₂(0)=307.16 and the nominal values of the batch reactor are the same as taken by [24, 25].

Initially, the batch reactor is on the nominal values, and the goal is to control the reaction temperature without restrictions in 𝑢. The desired value for the temperature of the reactor is 𝑥₂(𝑇)=303.16. For this objective, Method 1 is applied to the process. Two scenarios were simulated. In the first one, the linearization was not updated, and the control system works with the linearization around the initial condition. In the second one, it was updated at each 𝜏_𝑖=1. The optimal solution was computed by methods described in [19]. The final penalization was chosen as 𝑠=22408. The results are illustrated in Figures 1, 2, and 3.

Figure 1 

State trajectories for Method 1 and a comparison against the optimal solution with 𝑠=22408.

Figure 2 

Control trajectories for Method 1 and a comparison against the optimal solution with 𝑠=22408.

Figure 3 

Difference between optimal cost and suboptimal cost trajectories generated by Method 1.

Dynamically, the optimal solution has a worse behavior than that of the suboptimal trajectories because its state and control have more variability than the state and control of the suboptimal trajectories. However, the final error between the final desired state and the reached state is lower for the optimal solution as well as its total cost. Note that the solution with 𝜏_𝑖=1 is nearly arrived at the final desired state (𝑥₂(𝑇)=303.16).

The relative error between the optimal cost (𝐽_opt) and the suboptimal cost (𝐽_meth1) was calculated as %error=100(|𝐽_opt−𝐽_meth1|/𝐽_opt). When the linearization is updated with more frequency (for instance 𝜏_𝑖=1), the error decreases as it is illustrated in Figure 3. For both updating schemes of linearization, the error is smaller than 3%. That means that Method 1 approximates appropriately the optimal cost of the nonlinear system.

When the input constraints are active, two control strategies are proposed. A traditional predictive controller (MPC) according to [12, 26] and a suboptimal controller following Method 2 are proposed here. The suboptimal controller can incorporate the possibility of satisfy restrictions in the manipulated variable handling saturation using the tuning of the parameter 𝛾. Two 𝛾-values are selected (𝛾=0.10and0.15). Furthermore, the linearization is updated at each 𝜏_𝑖=1 and the time operation for the reactor is 𝑇=50.

Notice that neither the MPC nor the suboptimal controller with 𝛾=0.10 can completely reject that the initial disturbance due to the saturation of the manipulated variable. This phenomenon is shown in Figure 4. When the reaction begins, both control schemes try to correct the increase of temperature in the reactor, but quickly the manipulated variable is saturated. However, when 𝛾=0.15 the saturation of the manipulated variable is avoided.

Figure 4 

State trajectories for Method 2 with 𝛾=0.1 and 𝛾=0.15 and comparison to MPC.

Note also that the saturation time with the suboptimal controller is less than that with the MPC, and this could be modified by setting the parameter 𝛾. This is certainly an advantage in systems where it may result beneficially in an initial saturation.

Figure 5 shows the dynamic responses of the reactor temperature for the three cases here studied. Even though the performance is similar in all three cases, the functional objective 𝐽 was evaluated and a worse performance than the MPC was obtained. Anyway, the improvement in the functional is approximately 2% with the method proposed here.

Figure 5 

Control trajectories for Method 2 with 𝛾=0.1 and 𝛾=0.15 and comparison to MPC.

3.2. Case Study: Penicillin Fed-Batch Reactor

In this section, a fed-batch reactor for penicillin production is considered (see [27, 28]). The nonlinear process model is given by the following ordinary differential equations (ODEs): 𝑥̇₁(𝑡)=ℎ₁𝑥₁𝑥−𝑢¹500𝑥₄,𝑥₁(0)𝑥̇=1.5,₂(𝑡)=ℎ₂𝑥₁−0.01𝑥₂𝑥−𝑢²500𝑥₄,𝑥₂(0)𝑥̇=0,₃(𝑡)ℎ=−¹𝑥₁0.47−𝑥₁0.029𝑥₃0.0001+𝑥₃−⋯ℎ₂𝑥₁+𝑢1.2𝑥₄𝑥1−³500,𝑥₃(0)𝑥̇=0,₄(𝑡)=𝑢500,𝑥₄(0)ℎ=7.₁𝑥=0.11³0.006𝑥₁+𝑥₃,ℎ₂𝑥=0.0055³0.0001+𝑥₃1+10𝑥₃. (3.2)

In this model 𝑥₁,𝑥₂, 𝑥₃, and 𝑥₄ are the biomass, penicillin concentration, substrate concentration (g/L), and the reactor volume (L), respectively. The control variable 𝑢 is the feed rate of substrate.

The goals for the control of the penicillin reactor were extracted from [28], and they are as follows.(i)To reach in 132 hours (𝑇=132 is the time horizon) a final penicillin concentration of 8 g/L; that is, the desired output is 𝑥₂=8.(ii)The batch process includes an input constraint, which must be satisfied during the operation process |𝑢|⩽15. (57)The values of the parameters of the cost function are 𝑄=𝐼_4×4,𝑅=3, and 𝑆=20𝐼_4×4.

3.2.1. Results of Method 1

The strategy starts by computing the linearization of the system given by (1); that is, 𝐴 and 𝐵 are initially calculated at 𝑡=0, namely, ,01𝐴=0018.300−0.0182.5000−542.700000𝐵=−0.000471.500 (58)

The two equations in (28) were calculated for this approximation and the Riccati gain is now available. Four time instants were selected: 𝜏₁=30, 𝜏₂=60, 𝜏₃=90, and, 𝜏₄=120. In these time instants the controller gain was updated.

In Figure 6, the evolution of the states for the finite horizon optimization problem and the control reference trajectory are depicted. The batch process has a typical response (according to [28]), but for this set-up there is an error between the desired state 𝑥₂=8 and the real state reached, 𝑥₂=6.47, at the end of the operation time of the batch reactor. So, the parameter 𝑆 is increased and the results are shown in Figure 7. Clearly, the goals are now satisfied. Note that these strategies accept higher control values than 15. In other words, the constraint imposed in the manipulated variable is not satisfied yet.

Figure 6 

States and control trajectories for finite horizon set-up. 𝑇=132 and 𝑆=20.

Figure 7 

States and control trajectories for finite horizon set-up. 𝑇=132 but 𝑆=60.

3.2.2. Results of Method 2

In this subsection, the results of Method 2 are described. The optimization problem via LMI is solved by traditional numerical methods (see [21]) and implemented with standard software like MATLAB, as well as the first strategy.

The controller gain was updated four times according to the time instants 𝜏_𝑘 from the first method solving the LMI ℱ(Π) (53) at each time.

The parameter 𝛾 in (46) was taken as [23] suggests and that ensures the validity of the assumption of Lemma 1. The strategy for the second method was simulated, and in Figure 8 the results are depicted. The reader can notice that the objectives fixed for the problem are satisfied; that is, the output 𝑥₂ is near the desired value and the manipulated variable satisfies the constraint avoiding saturation.

Figure 8 

States and control trajectories for the second method via LMIs in an infinite horizon set-up.

4. Conclusions

In this paper, complete control schemes have been developed and illustrated. The innovation in this paper is the novel way to obtain the control reference trajectory (or recipe), which is necessary in the context of batch reactors. These kinds of trajectories are constructed by using a set of new algebraic equations which are described above for the linear case. The other contribution is to handle constraints in the manipulated variable when the process is controlled and optimized. All tools developed to design come from modern optimal control theory based on the Hamiltonian formalism.

It is important to remark that the entire proposed schemes are designed to be implemented in real-time, with the controller parameters being updated on-line. This is because the resulting controller is based on the simplifications made of costate on the system, and, in consequence, a simple controller can be online-tuned. Thus, the two resulting strategies for tuning simple controllers are able to generate suboptimal trajectories such that, in the second case, different constraints can be satisfied.

Also, it is significant to comment that the first strategy was developed using finite horizon optimal control theory with the restrictions that support the definition of the objective function. Since the control action written in costate terms and assuming linearity in the costates, the resulting control law has time-varying gains. On the other hand, the second strategy is presented under infinite horizon optimal control framework and allows the inclusion of new constraints by means of LMI as for example maximum and minimum constraints in control variable. Maintaining the linearity in the co-states, it is possible to find a suboptimal control law, but with constant feedback gains. At this point, notice that the simplifications may lead to a suboptimal solution to the optimization problem. The counterpart of this is that the reduction in the complexity of the mathematical problem treatment leads to less computational effort.

Finally, the numerical results show a good performance in the controlled variables, with the suboptimal strategies, and, especially, in the second strategy, when the control is bounded. In general, all variables have a good performance inside the operation time fixed for the batch reactor achieving the requirements imposed upon the system.