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

Authors: Pablo S. Rivadeneira; Eduardo J. Adam
Journal: Journal of Nonlinear Dynamics
 (2013)

---

## 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 T and the final penalization S. 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.


---

## Body

## 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 controlu(t) 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 u(t)=Kx(t) (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α(T,S), β(T,S) solutions of a set of partial differential equations (PDEs). The time-varying Riccati gain is calculated in parallel with the process based on α(T-t,S), β(T-t,S). (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 Section2, 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:(1)x˙(t)=f(x(t))+g(x(t))u(t);x(0)=x0,
where x∈ℝn is the vector of state variables with initial conditions x0, u∈𝒰∈ℝm is the control vector, and f(x)∈ℝn and g(x)∈ℝn×m are smooth vector fields. The general dynamic optimization of a bioprocess needs to find the control vector u, which minimizes a performance index as
(2)𝒥(T,0,x0,u(·))=∫0TL(x(τ),u(τ))dτ+𝒦(x(T))
subject to the dynamics (1). The function L is known as the Lagrangian, 𝒦 is a final penalization on the states, and T is the duration of the batch run. Although L and 𝒦 can be general functions, here only quadratics forms will be considered, namely,
(3)L(x,u)=x′Qx+u′Ru,𝒦=x(T)′Sx(T),
where Q and S are n×n symmetric and positive semidefinite matrices and R is a n×m positive-definite matrix. Under these conditions the Hamiltonian ℋ of such problem defined as
(4)ℋ(x,λ,u)≜L(x,u)+λ′(f(x)+g(x)u)
is known to be regular; that is, there is a unique H-optimal control u0 that minimizes the Hamiltonian ℋ,
(5)u0(x,λ)≜argminuℋ(x,λ,u),
and the derivative of ℋ with respect to u vanishes at (x,λ,u0(x,λ)). For control affine systems (1), the optimal control takes the following form:
(6)u0(x,λ)=-12R-1g(x)′λ.
The control-Hamiltonian
(7)ℋ0(x,λ)≜ℋ(x,λ,u0(x,λ)),ℋ0(x,λ)=x′Qx+λ′f(x)-14λ′g(x)R-1g(x)′λ,
gives rise to the Hamiltonian canonical equations (HCEs) [4]. Consider
(8)x˙(t)=(∂ℋ0∂λ)′;x(0)=x0,λ˙(t)=-(∂ℋ0∂x)′;λ(T)=2Sx(T).


### 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 u0 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 (x-, u-), that is,
(9)x˙(t)≈Ax(t)+Bu(t),x(0)=x0,A=(∂f(x)∂x+∂g(x)∂xu)|x-,u-,B=g(x)|x-,u-,
where (x-, u-) is a fixed point of the trajectory of the system (1) at time t≥0.(ii)
The costatesλ(t) can be linearized in a small neighborhood of (x,λ); that is, λ(t)≈λ~(t)=2P(t)x(t). As a result, the optimal control law equation (6) becomes
(10)u≈u~=-R-1B′P(t)x(t),
where P(t) is the solution to the differential Riccati equation (DRE):
(11)P˙(t)=P(t)WP(t)-P(t)A-A′P(t)-Q,P˙(t)=P(t)WP(t)-cciP(t)A-P(T)=S,W=BR-1B′.Notice that, in this paper, u is the manipulated variable corresponding to nonlinear problem and u~ 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 x(T) and λ(0)) or to calculate the solution of differential Riccati equation P(t).


### 2.2. Method 1: Generating of Suboptimal Control Reference Trajectory

In order to find the missing boundary valuesx(T) 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 matrixH, (which comes from (8) applied to the system (9))
(12)H=def(A-12W-2Q-A′),(ii)
the augmented Hamiltonian system (a linear2n-dimensional matrix ODE with a final condition) defined for the combination of two n×n matrices X(t),Λ(t),t∈[0,T] through
(13)(X˙Λ˙)=H(XΛ);(X(T)Λ(T))=(I2S).The solution to the augmented (linear) problem (13) being unique must verify
(14)(X(0)Λ(0))=e-HT(I2S),
and since in this linear case equation (8) also reads as follows:
(15)(x˙λ˙)=H(xλ),
then the missing boundary conditions can be explicitly found, namely,
(16)(x0λ(0))=e-HT(x(T)2Sx(T))=(X(0)Λ(0))x(T),(17)x(T)=X-1(0)x0,λ(0)=Λ(0)X-1(0)x0,
(see [4] for a proof of the invertibility of X and other details). Actually, the whole solution to the DRE can be recovered from the solution to problem (13), explicitly, as follows:
(18)P(t)=Λ(t)X-1(t),t∈[0,T].
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 (T,S) without solving either the DRE or the HCEs. To solve the whole (T,S)-family of LQR problems (with common A,B,Q,R,x0 values) consider each particular (T,S)-problem:
(19)ρ(T,S)=defxT,S(T);σ(T,S)=defλT,S(0);(20)U(T)=defe-HT,(21)(α(T,S)β(T,S))=def(XT,S(0)ΛT,S(0)),
where the superscript (T,S) refers to the individual problem of duration T and final penalty matrix S. Then (14) can be rewritten in the form
(22)(α(T,S)β(T,S))=U(T)·(I2S).
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 T (denoted by the subscript (T)) results in
(23)(αTβT)=-HU(I2S),
where the variables (T,S) 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 S to make sense of the partial derivatives (∂/∂S) 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, S=sI, s∈ℝ0+.  The partial derivatives of (22) with respect to the real variable s are then
(24)(αsβs)=U(02I).
Now, by partitioning U into n×n submatrices Ui,i=1,…4 in the form
(25)U=(U1U2U3U4),
then (24) reads as follows:
(26)12αs=U2,12βs=U4,
which combined with (22) gives
(27)U1=α-sαs,U3=β-sβs,
and by inserting these results in (23), the following relations are obtained:
(28)αT-αsM=-αN,βT-βsM=-βN,
where
(29)M=defA′S+SA+Q-SWS,N=defA-WS.
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,
(30)α(0,s)=I,β(0,s)=2S.
Now, the missing values, the final state and initial costate, for any (T,s)-problem can be recuperated immediately as
(31)ρ=α-1x0,σ=βα-1x0=βρ.
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 n-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(T,S)-LQR problem the optimal initial costate has the form given by (10)
(32)σ=λT,S(0)=2P(0)x0,
where P(·) is the numerical solution of the corresponding DRE, that is, the final value matrix ODE in (11). Therefore, from (31) and (32), for each (T,S)-problem the Riccati matrix P(·) should also verify
(33)P(0)=12β(T,S)[α(T,S)]-1,
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 (T,S)-problem. They also avoid storing, necessarily as an approximation, the Riccati matrix P(t) for the values of t∈[0,T] chosen by the numerical integrator, possibly different from the time instants for which the control u(t) is constructed. Instead, the HCEs (8) can be integrated with initial conditions
(34)x(0)=x0,λ(0)=σ(T,S),
rendering the optimal trajectories x*(t),λ*(t) for 0≤t≤T, which allow to generate the optimal control at each sampling time
(35)u*(t)=u0(x*(t),λ*(t))=-12R-1B′λ*(t).
Even better, the feedback form for the control becomes directly available due to the linear dependence of (31) on initial conditions, namely
(36)u*(t)=-12R-1B′β(T-t,S)[α(T-t,S)]-1x*(t),
since λ*(t) is also the optimal initial costate for the problem with parameters (T-t,S) starting at x(0)=x*(t). Then, as a side-product, an alternative formula for the Riccati matrix results:
(37)P(t)=12β(T-t,S)[α(T-t,S)]-1,∀t∈[0,T].


#### 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 (x-,u-). After that, the DRE is solved by the matrices α, β according to (37) for a range of T and S, 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 matricesA, B at arbitrary time instants τi, i=0,1,…,n is given:In:
x(0), Q, R, T, S, {τ0,τ1,…,τn},Out:
u(t), for all t≥0.Step 1: Take i=0. At time t=τ0=0, the first pair (A,B) is computed from (9) taking x-=x(0) and u-=0 or as solution of f(x(0))+g(x(0))u=0.Step 2: Compute (α,β) using (22) or as a solution of (28) by mean of the two matrices (A,B), andStep 3: run in parallel the process sending the control signal u(t) = (-1/2)R-1B′β(T-t,S)[α(T-t,S)]-1x(t) during the interval τi≤t≤τi+1.Step 4: at time instant t=τi+1, update the matrices A, B, α, and β of the control signal of Step 3 taking x-=x(τi+1) and u-=u(τi+1). Update i=i+1.Step 5: return to Step 3. This procedure continues until completing the whole process horizontally, that is, until to reach t=T=τn.Step 6: if t=T, then the process is over.Notice thatP(t)=(1/2)β(T-t,S)[α(T-t,S)]-1 describes a time-variant feedback control law, where the matrices A, B, α, and β are updated at each instant τi. In the scheme above, the time instants (τi) 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:
(38)minu∈𝒰𝒥=∫0∞(x′Qx+u′Ru)dτ,
subject to dynamics (9) and
(39)∥u(t)∥⩽umax.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:(40)ℱ(π)≜ℱ0+∑‍mi=1πiℱi>0,
where π∈ℝm and πi are the optimization variables and the positive and symmetric matrices ℱi=ℱi′∈ℝm×m,  i=1,…,m, 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) IfT is considered as an infinite horizon, it is known that when the control weighting matrix R in the cost functional is positive definite, the state weighting matrix Q is nonnegative definite and the pair (A,B) 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:
(41)ΠA+A′Π+Q-ΠBR-1B′Π≥0,
where Π>0 is the solution of the last inequality and Π∈ℝn×n is a positive symmetric matrix. In this section, the Riccati gain P 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(42)(Z(π)D(π)D(π)′W(π))>0,
where Z(π)=Z(π)′, W(π)=W(π)′, and D(π) depend on π, is equivalent to the matrix inequality
(43)W(π)>0,Z(π)-D(π)W(π)-1D(π)′>0,
or equivalently one has
(44)Z(π)>0,W(π)-D(π)′Z(π)-1D(π)>0.SinceR>0, and by using (44), the algebraic Riccati inequality (41) can be written as a LMI:
(45)ℱ1(Π)=(ΠA+A′Π+QΠBB′ΠR)>0.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 variableu(t) as in (39). These constraints can be seen as a LMI.Considering thatu(t) is a stabilizing control law, then one has the following.Lemma 1.
Assuming that aγ parameter exists, the following LMI
(46)ℱ2(Π)=(umax2γIR-1B′Π(R-1B′Π)′I)≥0
is equivalent to
(47)∥u(t)∥≤umax,t≥0.Proof.
Consider that(48)∥u(t)∥2=∥R-1B^′Πx(t)∥2=x(t)′(R-1B^′Π)′R-1B^′Πx(t)=x(t)′N′Nx(t)=x(t)′Dx(t).
Since D is a symmetric and square constant matrix, it is possible to decompose it into its singular values, and then the following inequality holds:
(49)νminI≤MDM-1≤νmaxI,
where ν is a singular value; then
(50)x(t)′N′Nx(t)≤νmaxIx(t)′x(t)≤νmaxυI≤umax2I,νmax{(R-1B′Π)′R-1B′Π}υ≤umax2I,((R-1B′Π)′R-1B′Π)υ≤umax2I.
Using (44), ∥u(t)∥≤umax results equivalently to
(51)(umax2γIR-1B′Π(R-1B′Π)′I)≥0,
where γ=1/υ is a real positive constant that satisfies a stabilization condition ∥x(t)∥2≤υ for any t≥0.Note that to findυ, it is necessary to involve another procedure to design an invariant ellipsoid to guarantee either ∥x(t)∥2≤υ or x(t)′P1x(t)≤υ, where P1 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 υ=∥x(0)∥2 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,ℱ1(π),…,ℱn(π)>0 as a single LMI as follows:
(52)diag(ℱ1(π),…,ℱn(π))>0,
which makes us introduce the input constraints to the optimal problem as
(53)ℱ(Π)=(ΠA+A′Π+QΠB00B′ΠR0000umax2γIR-1B′Π00(R-1B′Π)′I)>0.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)u(t)=-R-1B′Πx(t),
where Π is the solution of minimizing of the LMI (53).The procedure to implement these results in nonlinear systems is as follows.In:
x(0), Q, R, T,  umax, {τ0,τ1,…,τn}.Out:
u(t), for all t≥0.Step 1: Take i=0. At time t=τ0=0, the first pair (A,B) is obtained either from (9) taking x-=x(0) and u-=0 or as solution of f(x(0))+g(x(0))u=0.Step 2: Compute the LMI (53) by mean of the two matrices (A,B),Step 3: run in parallel the process sending the control signal u(t)=-R-1B′Πx(t) during the interval τi≤t≤τi+1.Step 4: at time instant t=τi+1, update the matrices A, B, and Π of the control signal of Step 3 taking x-=x(τi+1) and u-=u(τi+1). Take i=i+1.Step 5: return to Step 3. This procedure continues until completing the whole optimization horizon, that is, until to reach t=T=τn.Step 6: if t=T, then the optimization process is over.The time instantsτi 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.


## 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 u0 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 (x-, u-), that is,
(9)x˙(t)≈Ax(t)+Bu(t),x(0)=x0,A=(∂f(x)∂x+∂g(x)∂xu)|x-,u-,B=g(x)|x-,u-,
where (x-, u-) is a fixed point of the trajectory of the system (1) at time t≥0.(ii)
The costatesλ(t) can be linearized in a small neighborhood of (x,λ); that is, λ(t)≈λ~(t)=2P(t)x(t). As a result, the optimal control law equation (6) becomes
(10)u≈u~=-R-1B′P(t)x(t),
where P(t) is the solution to the differential Riccati equation (DRE):
(11)P˙(t)=P(t)WP(t)-P(t)A-A′P(t)-Q,P˙(t)=P(t)WP(t)-cciP(t)A-P(T)=S,W=BR-1B′.Notice that, in this paper, u is the manipulated variable corresponding to nonlinear problem and u~ 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 x(T) and λ(0)) or to calculate the solution of differential Riccati equation P(t).


## 2.2. Method 1: Generating of Suboptimal Control Reference Trajectory

In order to find the missing boundary valuesx(T) 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 matrixH, (which comes from (8) applied to the system (9))
(12)H=def(A-12W-2Q-A′),(ii)
the augmented Hamiltonian system (a linear2n-dimensional matrix ODE with a final condition) defined for the combination of two n×n matrices X(t),Λ(t),t∈[0,T] through
(13)(X˙Λ˙)=H(XΛ);(X(T)Λ(T))=(I2S).The solution to the augmented (linear) problem (13) being unique must verify
(14)(X(0)Λ(0))=e-HT(I2S),
and since in this linear case equation (8) also reads as follows:
(15)(x˙λ˙)=H(xλ),
then the missing boundary conditions can be explicitly found, namely,
(16)(x0λ(0))=e-HT(x(T)2Sx(T))=(X(0)Λ(0))x(T),(17)x(T)=X-1(0)x0,λ(0)=Λ(0)X-1(0)x0,
(see [4] for a proof of the invertibility of X and other details). Actually, the whole solution to the DRE can be recovered from the solution to problem (13), explicitly, as follows:
(18)P(t)=Λ(t)X-1(t),t∈[0,T].
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 (T,S) without solving either the DRE or the HCEs. To solve the whole (T,S)-family of LQR problems (with common A,B,Q,R,x0 values) consider each particular (T,S)-problem:
(19)ρ(T,S)=defxT,S(T);σ(T,S)=defλT,S(0);(20)U(T)=defe-HT,(21)(α(T,S)β(T,S))=def(XT,S(0)ΛT,S(0)),
where the superscript (T,S) refers to the individual problem of duration T and final penalty matrix S. Then (14) can be rewritten in the form
(22)(α(T,S)β(T,S))=U(T)·(I2S).
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 T (denoted by the subscript (T)) results in
(23)(αTβT)=-HU(I2S),
where the variables (T,S) 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 S to make sense of the partial derivatives (∂/∂S) 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, S=sI, s∈ℝ0+.  The partial derivatives of (22) with respect to the real variable s are then
(24)(αsβs)=U(02I).
Now, by partitioning U into n×n submatrices Ui,i=1,…4 in the form
(25)U=(U1U2U3U4),
then (24) reads as follows:
(26)12αs=U2,12βs=U4,
which combined with (22) gives
(27)U1=α-sαs,U3=β-sβs,
and by inserting these results in (23), the following relations are obtained:
(28)αT-αsM=-αN,βT-βsM=-βN,
where
(29)M=defA′S+SA+Q-SWS,N=defA-WS.
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,
(30)α(0,s)=I,β(0,s)=2S.
Now, the missing values, the final state and initial costate, for any (T,s)-problem can be recuperated immediately as
(31)ρ=α-1x0,σ=βα-1x0=βρ.
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 n-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(T,S)-LQR problem the optimal initial costate has the form given by (10)
(32)σ=λT,S(0)=2P(0)x0,
where P(·) is the numerical solution of the corresponding DRE, that is, the final value matrix ODE in (11). Therefore, from (31) and (32), for each (T,S)-problem the Riccati matrix P(·) should also verify
(33)P(0)=12β(T,S)[α(T,S)]-1,
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 (T,S)-problem. They also avoid storing, necessarily as an approximation, the Riccati matrix P(t) for the values of t∈[0,T] chosen by the numerical integrator, possibly different from the time instants for which the control u(t) is constructed. Instead, the HCEs (8) can be integrated with initial conditions
(34)x(0)=x0,λ(0)=σ(T,S),
rendering the optimal trajectories x*(t),λ*(t) for 0≤t≤T, which allow to generate the optimal control at each sampling time
(35)u*(t)=u0(x*(t),λ*(t))=-12R-1B′λ*(t).
Even better, the feedback form for the control becomes directly available due to the linear dependence of (31) on initial conditions, namely
(36)u*(t)=-12R-1B′β(T-t,S)[α(T-t,S)]-1x*(t),
since λ*(t) is also the optimal initial costate for the problem with parameters (T-t,S) starting at x(0)=x*(t). Then, as a side-product, an alternative formula for the Riccati matrix results:
(37)P(t)=12β(T-t,S)[α(T-t,S)]-1,∀t∈[0,T].


### 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 (x-,u-). After that, the DRE is solved by the matrices α, β according to (37) for a range of T and S, 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 matricesA, B at arbitrary time instants τi, i=0,1,…,n is given:In:
x(0), Q, R, T, S, {τ0,τ1,…,τn},Out:
u(t), for all t≥0.Step 1: Take i=0. At time t=τ0=0, the first pair (A,B) is computed from (9) taking x-=x(0) and u-=0 or as solution of f(x(0))+g(x(0))u=0.Step 2: Compute (α,β) using (22) or as a solution of (28) by mean of the two matrices (A,B), andStep 3: run in parallel the process sending the control signal u(t) = (-1/2)R-1B′β(T-t,S)[α(T-t,S)]-1x(t) during the interval τi≤t≤τi+1.Step 4: at time instant t=τi+1, update the matrices A, B, α, and β of the control signal of Step 3 taking x-=x(τi+1) and u-=u(τi+1). Update i=i+1.Step 5: return to Step 3. This procedure continues until completing the whole process horizontally, that is, until to reach t=T=τn.Step 6: if t=T, then the process is over.Notice thatP(t)=(1/2)β(T-t,S)[α(T-t,S)]-1 describes a time-variant feedback control law, where the matrices A, B, α, and β are updated at each instant τi. In the scheme above, the time instants (τi) 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.2.1. Relations amongst (α, β) Solutions and Riccati Equations

For any(T,S)-LQR problem the optimal initial costate has the form given by (10)
(32)σ=λT,S(0)=2P(0)x0,
where P(·) is the numerical solution of the corresponding DRE, that is, the final value matrix ODE in (11). Therefore, from (31) and (32), for each (T,S)-problem the Riccati matrix P(·) should also verify
(33)P(0)=12β(T,S)[α(T,S)]-1,
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 (T,S)-problem. They also avoid storing, necessarily as an approximation, the Riccati matrix P(t) for the values of t∈[0,T] chosen by the numerical integrator, possibly different from the time instants for which the control u(t) is constructed. Instead, the HCEs (8) can be integrated with initial conditions
(34)x(0)=x0,λ(0)=σ(T,S),
rendering the optimal trajectories x*(t),λ*(t) for 0≤t≤T, which allow to generate the optimal control at each sampling time
(35)u*(t)=u0(x*(t),λ*(t))=-12R-1B′λ*(t).
Even better, the feedback form for the control becomes directly available due to the linear dependence of (31) on initial conditions, namely
(36)u*(t)=-12R-1B′β(T-t,S)[α(T-t,S)]-1x*(t),
since λ*(t) is also the optimal initial costate for the problem with parameters (T-t,S) starting at x(0)=x*(t). Then, as a side-product, an alternative formula for the Riccati matrix results:
(37)P(t)=12β(T-t,S)[α(T-t,S)]-1,∀t∈[0,T].


## 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 (x-,u-). After that, the DRE is solved by the matrices α, β according to (37) for a range of T and S, 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 matricesA, B at arbitrary time instants τi, i=0,1,…,n is given:In:
x(0), Q, R, T, S, {τ0,τ1,…,τn},Out:
u(t), for all t≥0.Step 1: Take i=0. At time t=τ0=0, the first pair (A,B) is computed from (9) taking x-=x(0) and u-=0 or as solution of f(x(0))+g(x(0))u=0.Step 2: Compute (α,β) using (22) or as a solution of (28) by mean of the two matrices (A,B), andStep 3: run in parallel the process sending the control signal u(t) = (-1/2)R-1B′β(T-t,S)[α(T-t,S)]-1x(t) during the interval τi≤t≤τi+1.Step 4: at time instant t=τi+1, update the matrices A, B, α, and β of the control signal of Step 3 taking x-=x(τi+1) and u-=u(τi+1). Update i=i+1.Step 5: return to Step 3. This procedure continues until completing the whole process horizontally, that is, until to reach t=T=τn.Step 6: if t=T, then the process is over.Notice thatP(t)=(1/2)β(T-t,S)[α(T-t,S)]-1 describes a time-variant feedback control law, where the matrices A, B, α, and β are updated at each instant τi. In the scheme above, the time instants (τi) 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:
(38)minu∈𝒰𝒥=∫0∞(x′Qx+u′Ru)dτ,
subject to dynamics (9) and
(39)∥u(t)∥⩽umax.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:(40)ℱ(π)≜ℱ0+∑‍mi=1πiℱi>0,
where π∈ℝm and πi are the optimization variables and the positive and symmetric matrices ℱi=ℱi′∈ℝm×m,  i=1,…,m, 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) IfT is considered as an infinite horizon, it is known that when the control weighting matrix R in the cost functional is positive definite, the state weighting matrix Q is nonnegative definite and the pair (A,B) 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:
(41)ΠA+A′Π+Q-ΠBR-1B′Π≥0,
where Π>0 is the solution of the last inequality and Π∈ℝn×n is a positive symmetric matrix. In this section, the Riccati gain P 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(42)(Z(π)D(π)D(π)′W(π))>0,
where Z(π)=Z(π)′, W(π)=W(π)′, and D(π) depend on π, is equivalent to the matrix inequality
(43)W(π)>0,Z(π)-D(π)W(π)-1D(π)′>0,
or equivalently one has
(44)Z(π)>0,W(π)-D(π)′Z(π)-1D(π)>0.SinceR>0, and by using (44), the algebraic Riccati inequality (41) can be written as a LMI:
(45)ℱ1(Π)=(ΠA+A′Π+QΠBB′ΠR)>0.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 variableu(t) as in (39). These constraints can be seen as a LMI.Considering thatu(t) is a stabilizing control law, then one has the following.Lemma 1.
Assuming that aγ parameter exists, the following LMI
(46)ℱ2(Π)=(umax2γIR-1B′Π(R-1B′Π)′I)≥0
is equivalent to
(47)∥u(t)∥≤umax,t≥0.Proof.
Consider that(48)∥u(t)∥2=∥R-1B^′Πx(t)∥2=x(t)′(R-1B^′Π)′R-1B^′Πx(t)=x(t)′N′Nx(t)=x(t)′Dx(t).
Since D is a symmetric and square constant matrix, it is possible to decompose it into its singular values, and then the following inequality holds:
(49)νminI≤MDM-1≤νmaxI,
where ν is a singular value; then
(50)x(t)′N′Nx(t)≤νmaxIx(t)′x(t)≤νmaxυI≤umax2I,νmax{(R-1B′Π)′R-1B′Π}υ≤umax2I,((R-1B′Π)′R-1B′Π)υ≤umax2I.
Using (44), ∥u(t)∥≤umax results equivalently to
(51)(umax2γIR-1B′Π(R-1B′Π)′I)≥0,
where γ=1/υ is a real positive constant that satisfies a stabilization condition ∥x(t)∥2≤υ for any t≥0.Note that to findυ, it is necessary to involve another procedure to design an invariant ellipsoid to guarantee either ∥x(t)∥2≤υ or x(t)′P1x(t)≤υ, where P1 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 υ=∥x(0)∥2 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,ℱ1(π),…,ℱn(π)>0 as a single LMI as follows:
(52)diag(ℱ1(π),…,ℱn(π))>0,
which makes us introduce the input constraints to the optimal problem as
(53)ℱ(Π)=(ΠA+A′Π+QΠB00B′ΠR0000umax2γIR-1B′Π00(R-1B′Π)′I)>0.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)u(t)=-R-1B′Πx(t),
where Π is the solution of minimizing of the LMI (53).The procedure to implement these results in nonlinear systems is as follows.In:
x(0), Q, R, T,  umax, {τ0,τ1,…,τn}.Out:
u(t), for all t≥0.Step 1: Take i=0. At time t=τ0=0, the first pair (A,B) is obtained either from (9) taking x-=x(0) and u-=0 or as solution of f(x(0))+g(x(0))u=0.Step 2: Compute the LMI (53) by mean of the two matrices (A,B),Step 3: run in parallel the process sending the control signal u(t)=-R-1B′Πx(t) during the interval τi≤t≤τi+1.Step 4: at time instant t=τi+1, update the matrices A, B, and Π of the control signal of Step 3 taking x-=x(τi+1) and u-=u(τi+1). Take i=i+1.Step 5: return to Step 3. This procedure continues until completing the whole optimization horizon, that is, until to reach t=T=τn.Step 6: if t=T, then the optimization process is over.The time instantsτi 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 reactionA→B takes place. It is assumed that the reactor has a cooling jacket whose temperature (u) can be directly manipulated. The goal is to regulate the reactor temperature (x2) by means of inlet coolant temperature (u). Furthermore, the manipulated variable has minimum and maximum constraints -10≤u≤10. The concentration is denoted by the state x1.The reactor dynamics is modeled by the following equations:(55)x˙1(t)=-k0eER/x2(t)x12(t),x˙2(t)=-UAMcp(x2(t)-(u+298.16))-ΔHMcpVk0eER/x2(t)x12(t),
subject to initial conditions x1(0)=0.9 and x2(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 inu. The desired value for the temperature of the reactor is x2(T)=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 τi=1. The optimal solution was computed by methods described in [19]. The final penalization was chosen as s=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 withs=22408.Figure 2
Control trajectories for Method 1 and a comparison against the optimal solution withs=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τi=1 is nearly arrived at the final desired state (x2(T)=303.16).The relative error between the optimal cost (Jopt) and the suboptimal cost (Jmeth1) was calculated as %error=100(|Jopt-Jmeth1|/Jopt). When the linearization is updated with more frequency (for instance τi=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 τi=1 and the time operation for the reactor is T=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.Figure5 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 J 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):
(3.2)x˙1(t)=h1x1-u(x1500x4),x1(0)=1.5,x˙2(t)=h2x1-0.01x2-u(x2500x4),x2(0)=0,x˙3(t)=-h1x10.47-x1(0.029x30.0001+x3)-⋯h2x11.2+ux4(1-x3500),x3(0)=0,x˙4(t)=u500,x4(0)=7.h1=0.11(x30.006x1+x3),h2=0.0055(x30.0001+x3(1+10x3)).In this modelx1,x2, x3, and x4 are the biomass, penicillin concentration, substrate concentration (g/L), and the reactor volume (L), respectively. The control variable u 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 in132 hours (T=132 is the time horizon) a final penicillin concentration of 8 g/L; that is, the desired output is x-2=8.(ii)
The batch process includes an input constraint, which must be satisfied during the operation process(57)|u|⩽15.The values of the parameters of the cost function are Q=I4×4,R=3, and S=20I4×4.


#### 3.2.1. Results of Method 1

The strategy starts by computing the linearization of the system given by (1); that is, A and B are initially calculated at t=0, namely,
(58)A=(0018.300-0.0182.5000-542.700000),B=(-0.00040171500).The two equations in (28) were calculated for this approximation and the Riccati gain is now available. Four time instants were selected: τ1=30, τ2=60, τ3=90, and, τ4=120. In these time instants the controller gain was updated.In Figure6, 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 x-2=8 and the real state reached, x2=6.47, at the end of the operation time of the batch reactor. So, the parameter S 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.T=132 and S=20.Figure 7
States and control trajectories for finite horizon set-up.T=132 but S=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τk 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 x2 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.


## 3.1. Nonlinear Batch Reactor

Consider a batch reactor with a nonlinear dynamic where an exothermic and irreversible second order chemical reactionA→B takes place. It is assumed that the reactor has a cooling jacket whose temperature (u) can be directly manipulated. The goal is to regulate the reactor temperature (x2) by means of inlet coolant temperature (u). Furthermore, the manipulated variable has minimum and maximum constraints -10≤u≤10. The concentration is denoted by the state x1.The reactor dynamics is modeled by the following equations:(55)x˙1(t)=-k0eER/x2(t)x12(t),x˙2(t)=-UAMcp(x2(t)-(u+298.16))-ΔHMcpVk0eER/x2(t)x12(t),
subject to initial conditions x1(0)=0.9 and x2(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 inu. The desired value for the temperature of the reactor is x2(T)=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 τi=1. The optimal solution was computed by methods described in [19]. The final penalization was chosen as s=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 withs=22408.Figure 2
Control trajectories for Method 1 and a comparison against the optimal solution withs=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τi=1 is nearly arrived at the final desired state (x2(T)=303.16).The relative error between the optimal cost (Jopt) and the suboptimal cost (Jmeth1) was calculated as %error=100(|Jopt-Jmeth1|/Jopt). When the linearization is updated with more frequency (for instance τi=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 τi=1 and the time operation for the reactor is T=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.Figure5 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 J 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):
(3.2)x˙1(t)=h1x1-u(x1500x4),x1(0)=1.5,x˙2(t)=h2x1-0.01x2-u(x2500x4),x2(0)=0,x˙3(t)=-h1x10.47-x1(0.029x30.0001+x3)-⋯h2x11.2+ux4(1-x3500),x3(0)=0,x˙4(t)=u500,x4(0)=7.h1=0.11(x30.006x1+x3),h2=0.0055(x30.0001+x3(1+10x3)).In this modelx1,x2, x3, and x4 are the biomass, penicillin concentration, substrate concentration (g/L), and the reactor volume (L), respectively. The control variable u 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 in132 hours (T=132 is the time horizon) a final penicillin concentration of 8 g/L; that is, the desired output is x-2=8.(ii)
The batch process includes an input constraint, which must be satisfied during the operation process(57)|u|⩽15.The values of the parameters of the cost function are Q=I4×4,R=3, and S=20I4×4.


### 3.2.1. Results of Method 1

The strategy starts by computing the linearization of the system given by (1); that is, A and B are initially calculated at t=0, namely,
(58)A=(0018.300-0.0182.5000-542.700000),B=(-0.00040171500).The two equations in (28) were calculated for this approximation and the Riccati gain is now available. Four time instants were selected: τ1=30, τ2=60, τ3=90, and, τ4=120. In these time instants the controller gain was updated.In Figure6, 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 x-2=8 and the real state reached, x2=6.47, at the end of the operation time of the batch reactor. So, the parameter S 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.T=132 and S=20.Figure 7
States and control trajectories for finite horizon set-up.T=132 but S=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τk 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 x2 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.


## 3.2.1. Results of Method 1

The strategy starts by computing the linearization of the system given by (1); that is, A and B are initially calculated at t=0, namely,
(58)A=(0018.300-0.0182.5000-542.700000),B=(-0.00040171500).The two equations in (28) were calculated for this approximation and the Riccati gain is now available. Four time instants were selected: τ1=30, τ2=60, τ3=90, and, τ4=120. In these time instants the controller gain was updated.In Figure6, 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 x-2=8 and the real state reached, x2=6.47, at the end of the operation time of the batch reactor. So, the parameter S 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.T=132 and S=20.Figure 7
States and control trajectories for finite horizon set-up.T=132 but S=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τk 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 x2 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.



---

*Source: 182034-2013-11-17.xml*