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

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. Twomethods 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 processT 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.


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 costate 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.

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: where  ∈ R  is the vector of state variables with initial conditions  0 ,  ∈ U ∈ R  is the control vector, and () ∈ R  and () ∈ R × are smooth vector fields.The general dynamic optimization of a bioprocess needs to find the control vector , which minimizes a performance index as subject to the dynamics (1).The function  is known as the Lagrangian, K is a final penalization on the states, and  is the duration of the batch run.Although  and K can be general functions, here only quadratics forms will be considered, namely, where  and  are  ×  symmetric and positive semidefinite matrices and  is a  ×  positive-definite matrix.Under these conditions the Hamiltonian H of such problem defined as is known to be regular; that is, there is a unique -optimal control  0 that minimizes the Hamiltonian H, and the derivative of H with respect to  vanishes at (, ,  0 (, )).For control affine systems (1), the optimal control takes the following form: The control-Hamiltonian H 0 (, ) ≜ H (, ,  0 (, )) , gives rise to the Hamiltonian canonical equations (HCEs) [4].Consider 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  0 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, 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 where () is the solution to the differential Riccati equation (DRE): 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 ().

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 H, (which comes from (8) applied to the system ( 9)) (ii) the augmented Hamiltonian system (a linear 2dimensional matrix ODE with a final condition) defined for the combination of two  ×  matrices (), Λ(),  ∈ [0, ] through The solution to the augmented (linear) problem ( 13) being unique must verify and since in this linear case equation ( 8) also reads as follows: then the missing boundary conditions can be explicitly found, namely, (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: 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 , , , ,  0 values) consider each particular (, )-problem: where the superscript ( , ) refers to the individual problem of duration  and final penalty matrix .Then ( 14) can be rewritten in the form 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 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,  = ,  ∈ R + 0 .The partial derivatives of ( 22) with respect to the real variable  are then Now, by partitioning  into  ×  submatrices   ,  = 1, . . . 4 in the form then (24) reads as follows: which combined with (22) gives and by inserting these results in (23), the following relations are obtained: where 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, Now, the missing values, the final state and initial costate, for any (, )-problem can be recuperated immediately as 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].

Relations amongst (𝛼, 𝛽) Solutions and Riccati Equations.
For any (, )-LQR problem the optimal initial costate has the form given by ( 10) 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 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 rendering the optimal trajectories  * (),  * () for 0 ≤  ≤ , which allow to generate the optimal control at each sampling time Even better, the feedback form for the control becomes directly available due to the linear dependence of (31) on initial conditions, namely 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: 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.
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)(−, )[( − , )] −1 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.

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: subject to dynamics (9) and 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 online 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 setup and for linear time-invariant systems.
A linear matrix inequality or LMI is a matrix inequality of the form: where  ∈ R  and   are the optimization variables and the positive and symmetric matrices F  = F   ∈ R × ,  = 1, . . ., , are given.The LMI (40) is a convex constraint on ; that is, the set { | F() > 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: where Π > 0 is the solution of the last inequality and Π ∈ R × 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 where () = ()  , () = ()  , and () depend on , is equivalent to the matrix inequality or equivalently one has  () > 0,  () − ()  () −1  () > 0. (44) Since  > 0, and by using (44), the algebraic Riccati inequality (41) can be written as a LMI: 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
Proof.Consider that Since  is a symmetric and square constant matrix, it is possible to decompose it into its singular values, and then the following inequality holds: where ] is a singular value; then Using (44), ‖()‖ ≤  max results equivalently to where  = 1/ is a real positive constant that satisfies a stabilization condition ‖()‖ 2 ≤  for any  ≥ 0.
Note that to find , it is necessary to involve another procedure to design an invariant ellipsoid to guarantee either ‖()‖ 2 ≤  or ()   1 () ≤ , where  1 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)‖ 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, F 1 (), . . ., F  () > 0 as a single LMI as follows: which makes us introduce the input constraints to the optimal problem as 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 where Π is the solution of minimizing of the LMI (53).
The procedure to implement these results in nonlinear systems is as follows.
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.

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 ( 2 ) by means of inlet coolant temperature ().Furthermore, the manipulated variable has minimum and maximum constraints −10 ≤  ≤ 10.The concentration is denoted by the state  1 .
The reactor dynamics is modeled by the following equations: subject to initial conditions  1 (0) = 0.9 and  2 (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  2 () = 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.
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 ( 2 () = 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.10 and 0.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.
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.Error between optimal solution and Method 1 Error between optimal solution and Method 1 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.

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): ) ,  1 (0) = 1.5, In this model  1 , 2 ,  3 , and  4 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  2 = 8.

(58)
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 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  2 = 8 and the real state reached,  2 = 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.

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.Final error ≈ 0.04 The input constraint is not satisfied (u ≤ 15) x 2 (132) ≈ 7.96 The controller gain was updated four times according to the time instants   from the first method solving the LMI F(Π) (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  2 is near the desired value and the manipulated variable satisfies the constraint avoiding saturation.

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.

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

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

Figure 3 :Figure 4 :
Figure 3: Difference between optimal cost and suboptimal cost trajectories generated by Method 1.

Figure 5 :
Figure 5: Control trajectories for Method 2 with  = 0.1 and  = 0.15 and comparison to MPC.