# Data-Driven Approximated Optimal Control for Chemical Processes with State and Input Constraints

Authors: Mingfang He
Journal: Complexity
 (2019)

---

## Abstract

Input and state constraints widely exist in chemical processes. The optimal control of chemical processes under the coexistence of inequality constraints on input and state is challenging, especially when the process model is only partially known. The objective of this paper is to design an applicable optimal control for chemical processes with known model structure and unknown model parameters. To eliminate the barriers caused by the hybrid constraints and unknown model parameters, the inequality state constraints are first transformed into equality state constraints by using the slack function method. Then, adaptive dynamic programming (ADP) with nonquadratic performance integrand is adopted to handle the augmented system with input constraints. The proposed approach requires only partial knowledge of the system, i.e., the model structure. The value information of the model parameters is not required. The feasibility and performance of the proposed approach are tested using two nonlinear cases including a continuous stirred-tank reactor (CSTR) example.


---

## Body

## 1. Introduction

Constraints on input and state commonly exist in chemical processes due to finite capability of the actuators [1], safety limits [2], requirements on product quality, and environmental regulations [3]. For example, in practice, it is required to keep the outlet species concentrations of a chemical reactor in reasonable and stable ranges. In addition, the value of manipulated variables is constrained in a certain range defined by the operating instructions. Therefore, the ability to handle constraints is an essential concern in the control design and synthesis of real chemical processes [4, 5].The study of constrained optimal control has undergone different stages and can be classified into different categories, e.g., problem transformation methods, Lyapunov function-based methods, state-dependent Riccati equation (SDRE) methods, model predictive control (MPC), and machine learning-based methods, to name a few. Since the 1960s, numerous approaches have been proposed to handle the input and/or state constraints in the optimal control of linear/nonlinear systems. Leitmann [6] and Bryson and Denham [7] proposed penalty-function technique for the optimal control of systems with state-variable inequality constraints; the solution of the penalized optimal control problem converges to the constrained optimal control problem if the penalty multiplier approaches infinity [8]. Sakawa [9] transformed the optimal control of linear systems with input constraints into an infinite-dimensional nonlinear programming problem by integrating the linear differential equation of the system states. Jacobson and Lele [10] used slack variable method to eliminate the scalar inequality constraints on the state variables. Hager [11] applied the Ritz–Trefftz method in the optimal control of system with both state and control constraints. Vlassenbroeck [12] transformed the state-variable inequality constrained optimal control problem into a parameter optimization problem using Chebyshev series expansion. Lim et al. [13] generalized the separation theorem to the constrained linear-quadratic (LQ) and linear-quadratic-Gaussian (LQG) optimal control problem. Manousiouthakis and Chmielewski [14] approximated the optimal control of nonlinear systems subject to pointwise-in-time inequality constraints in the SDRE framework. El-Farra and Christofides [3] developed a unified framework for constrained optimal control of nonlinear uncertain systems based on a general state-space Lyapunov approach. Balestrino et al. [15] designed a control Lyapunov R-functions (CLRF) solution for constrained linear systems with input and state constraints. Kiefer et al. [16] first transformed the input and state constraints into the constraints on the outputs and their higher-order derivatives, and then, special saturation functions were utilized to incorporate the constraints into overall control design. Stathopoulos et al. [17] studied the linear-quadratic-regulation (LQR) with input and state constraints. By using proximal algorithms and duality, they decomposed the corresponding quadratic program (QP) into two subproblems including an infinite-dimensional least squares problem and a simple clipping of an infinite sequence to the nonpositive orthant.MPC is a standard tool to handle input and state constraints within an optimal control setting [15], especially in industry applications [18], e.g., paper and pulp [19], minerals processing [20], chemical engineering [21–23], renewable energy [24], mechatronics engineering [25], and urban water supply [26]. One problem in the MPC implementation is the high computational load for large-scale or fast-sampling systems. Scokaert and Rawlings [27] solved the constrained LQR problem in a finite-dimensional MPC setting with optimality and stability. Wang and Wan [28] developed a structured neural network to solve the QP problem in a constrained MPC. Bemporad [29] moved the MPC computations offline and proposed a technique to obtain the piecewise explicit optimal control law for both MPC and constrained linear-quadratic-regulation (CLQR). By offline approximation, Pin et al. [30] designed a MPC scheme for constrained nonlinear discrete-time systems which allows coping with discontinuous control laws. Mhaskar [31] combined control Lyapunov function (CLF) and MPC to form a Lyapunov-based MPC approach that guarantees stability and constraint satisfaction from an explicitly characterized set of initial conditions. More detailed review of constrained optimal control could be found in [32–34] and the references therein.Application of the aforementioned methods relies on the full knowledge of system dynamics. With the increased unavailability of quality raw materials, it is imperative that raw materials of low grade with large variations should be employed in the production in order to maximize the use of resources. In this context, the operation of some chemical processes exhibit complexity in terms of variable dynamic characteristics, strong nonlinearities, heavy coupling, unclear mechanism, and mathematically unmodelable parts. Therefore, a common situation is that a chemical process has various working conditions. The structure of the process can be derived by applying conservative laws and its essential physicochemical mechanisms. However, the model parameters are unknown under some working conditions, e.g., insufficient data samples for model identification and undetermined reaction mechanism. This constitutes challenges to the existing control theory and technology.In recent years, the development of machine learning supported and gave rise to the emergence of data-driven constrained optimal control methods and the integration of machine learning and MPC [35]. Chakrabarty et al. [36] used a support-vector machine (SVM) to learn feasible region boundaries for an explicit nonlinear model predictive control (ENMPC). Lin and Zheng [37] devised a reinforcement learning agent to obtain the optimal control strategy for nonlinear systems with inequality constraints on input and state via cycle-by-cycle finite-time optimization. Full knowledge of the system model is still required in this approach. Abu-Khalaf and Lewis [38], Modares et al. [39], Luo et al. [40], Yang et al. [41], Zhang et al. [42], and Zhu et al. [43] used adaptive dynamic programming (ADP) [44–46] to approximate the optimal state-feedback controller for input-constrained nonlinear systems. Fan and Yang [47] applied ADP to approximate the optimal control law for state-constrained nonlinear systems. These ADP-based approaches are model-free or do not need partial knowledge of the system model; however, they provide the optimal control solution for either input-constrained or state-constrained systems. Chi et al. [48] studied the constrained data-driven iterative learning control (ILC) for point-to-point optimal control problem of discrete-time nonlinear systems with input and output constraints.The aim of this study is to design an approximated optimal control algorithm for continuous-time nonlinear chemical processes with both state and input constraints. The proposed algorithm requires only the structural knowledge of the system model, the value of the model parameters is not needed, which is a common situation in industrial applications. The input and state constraints were handled in a sequential manner. The state constraints were first eliminated by introducing slack functions to form an augmented system without inequality state constraints. Then, a nonquadratic performance integrand was adopted in the ADP framework to account for the input constraints in the augmented system.The rest of this paper is organized as follows. In Section2, the constrained optimal control is formulated, some preliminaries are introduced. The constraint handling approach and the approximated optimal control design are described in Section 3. Simulation results are presented and discussed in Section 4, followed by concluding remarks in Section 5.


## 2. Problem Formulation and Preliminaries

### 2.1. Formulation of the Nonlinear Constrained Optimal Control Problem

Consider a chemical process described by following continuous-time constrained system:(1)x˙=fx+gxu,where x∈X is the system state and u∈U is the manipulated input; fx∈ℝn and gx∈ℝn×m are differentiable and Lipschitz continuous functions. X and U are compact and contain the equilibrium point, i.e., the origin in their interiors. More specifically, X=x=x1,…,xn∈ℝn:Cmin⩽Cx⩽Cmax and U=u=u1,…,um∈ℝm:umin⩽u⩽umax denote the hard constraints on the states and inputs, respectively. Cx is a Lipschitz continuous function of x, which demonstrates the physical or technical constraints on the system states. Cmax, Cmin, umax, and umin are the maximum and minimum value vectors of Cx and input variables.For chemical processes described by (1) with system state and input constrained in X and U, the aim of optimal control is to find a control policy which minimizes the following performance index:(2)Vx0=∫0tfWxτ+Muτdτ,where Vx0 is the value function and W⋅ and M⋅ are positive definite weighting functions of the immediate states and inputs, respectively. Vx0 is an overall index of the control performance. A lower value of Vx0 indicates better control performance. Wx evaluates the deviation of the system state with the origin, e.g., the difference between the actual outlet concentration of a chemical reactor with its setting value. Mu accounts for the control effort to be placed on the process, e.g., the dosage of additive for the reactor to change the outlet concentration. 0,tf is the time interval of interest.


### 2.2. Preliminaries on Optimal Control of Nonlinear Systems

To start with, the preliminaries on optimal control of nonlinear systems are introduced without considering the system and input constraints. When the system is unconstrained,Wx and Mu are usually selected as Wx=xTQx and Mu=uTRu, with Q⩾0 and R>0, the symmetric weighting matrices. The solution of the unconstrained optimal control problem can be obtained by solving a Hamilton–Jacobi–Bellman (HJB) equation:(3)VxTfx+xTQx−14VxTgxR−1gxTVx=0,with the boundary condition V0=0. Vx is the partial derivative of the cost function V regarding system state x. If the HJB equation has an optimal solution V∗, then the optimal control is(4)u∗=−12gxTVx∗.The HJB equation involves a nonlinear differential equation whose analytic solution is difficult to obtain. However, for the HJB equation, if a control policyu has been improved using Vx to generate a better policy u′, then Vx′ corresponding to the improved policy u′ can be used to yield an even better policy u″. This approach of solving the optimal control policy is called policy iteration, which can gradually converge to an optimal policy and optimal value function.Therefore, successive approximation methods are developed to iteratively improve the control policy [49]. For system (1), starting from an initial admissible control u0x, two sequences of Vixi=0∞ and ui+1xi=0∞ could be generated via the policy iteration (PI) approach defined in [49]:(1)
Fori⩾0, solve the Lyapunov equation (LE) which is linear in V:(5)VixTfx+gxui+xTQx+uiTRui=0.(2)
Use the solutionVix to update the control policy:(6)ui∗=−12gxTVix∗,and(1)
0⩽Vi+1x⩽Vix(2)
uix is admissible(3)
ifV∗ and u∗ exist, Vix⟶V∗ and uix⟶u∗ [50]The above is the general iterative framework for solving the optimal control problem of nonlinear processes. Different from the unconstrained idea case, for practical chemical processes, the evolution trajectory of the system must obey the constraints. So, the feasible region of state and input is shrunk. These constraints must be taken into account in the optimal control design in order to avoid violating physical or technical limits.


## 2.1. Formulation of the Nonlinear Constrained Optimal Control Problem

Consider a chemical process described by following continuous-time constrained system:(1)x˙=fx+gxu,where x∈X is the system state and u∈U is the manipulated input; fx∈ℝn and gx∈ℝn×m are differentiable and Lipschitz continuous functions. X and U are compact and contain the equilibrium point, i.e., the origin in their interiors. More specifically, X=x=x1,…,xn∈ℝn:Cmin⩽Cx⩽Cmax and U=u=u1,…,um∈ℝm:umin⩽u⩽umax denote the hard constraints on the states and inputs, respectively. Cx is a Lipschitz continuous function of x, which demonstrates the physical or technical constraints on the system states. Cmax, Cmin, umax, and umin are the maximum and minimum value vectors of Cx and input variables.For chemical processes described by (1) with system state and input constrained in X and U, the aim of optimal control is to find a control policy which minimizes the following performance index:(2)Vx0=∫0tfWxτ+Muτdτ,where Vx0 is the value function and W⋅ and M⋅ are positive definite weighting functions of the immediate states and inputs, respectively. Vx0 is an overall index of the control performance. A lower value of Vx0 indicates better control performance. Wx evaluates the deviation of the system state with the origin, e.g., the difference between the actual outlet concentration of a chemical reactor with its setting value. Mu accounts for the control effort to be placed on the process, e.g., the dosage of additive for the reactor to change the outlet concentration. 0,tf is the time interval of interest.


## 2.2. Preliminaries on Optimal Control of Nonlinear Systems

To start with, the preliminaries on optimal control of nonlinear systems are introduced without considering the system and input constraints. When the system is unconstrained,Wx and Mu are usually selected as Wx=xTQx and Mu=uTRu, with Q⩾0 and R>0, the symmetric weighting matrices. The solution of the unconstrained optimal control problem can be obtained by solving a Hamilton–Jacobi–Bellman (HJB) equation:(3)VxTfx+xTQx−14VxTgxR−1gxTVx=0,with the boundary condition V0=0. Vx is the partial derivative of the cost function V regarding system state x. If the HJB equation has an optimal solution V∗, then the optimal control is(4)u∗=−12gxTVx∗.The HJB equation involves a nonlinear differential equation whose analytic solution is difficult to obtain. However, for the HJB equation, if a control policyu has been improved using Vx to generate a better policy u′, then Vx′ corresponding to the improved policy u′ can be used to yield an even better policy u″. This approach of solving the optimal control policy is called policy iteration, which can gradually converge to an optimal policy and optimal value function.Therefore, successive approximation methods are developed to iteratively improve the control policy [49]. For system (1), starting from an initial admissible control u0x, two sequences of Vixi=0∞ and ui+1xi=0∞ could be generated via the policy iteration (PI) approach defined in [49]:(1)
Fori⩾0, solve the Lyapunov equation (LE) which is linear in V:(5)VixTfx+gxui+xTQx+uiTRui=0.(2)
Use the solutionVix to update the control policy:(6)ui∗=−12gxTVix∗,and(1)
0⩽Vi+1x⩽Vix(2)
uix is admissible(3)
ifV∗ and u∗ exist, Vix⟶V∗ and uix⟶u∗ [50]The above is the general iterative framework for solving the optimal control problem of nonlinear processes. Different from the unconstrained idea case, for practical chemical processes, the evolution trajectory of the system must obey the constraints. So, the feasible region of state and input is shrunk. These constraints must be taken into account in the optimal control design in order to avoid violating physical or technical limits.


## 3. Approximated Constrained Optimal Controller Design with Partially Known System Dynamics

In this section, an optimal controller design approach without knowing the full system dynamics is proposed for constrained nonlinear systems. The state constraints are eliminated first by introducing slack functions. Then, the input constraints are managed by using ADP in the optimal control of the resulted augmented system. The iterative solution procedure introduced in Section2 relies on the precise knowledge of the system dynamics fx and gx. In order to eliminate this dependence, in this section, an approximated optimal control is studied and designed.


### 3.1. Handling State Constraints Using Slack Functions

Consider the state constraintCmin⩽Cx⩽Cmax; it contains the following inequality constraints:(7)Cmin−C⩽0,C−Cmax⩽0.There exists slack functionliαit ([10, 47]), which satisfy, for each state constraint, (8)ηix+liαit=0,where ηix could be Cix−Cimax, Cimin−Cix, or Cix−CimaxCimin−Cix, and according to the structure of ηix and the value of Cimin and Cimax, the slack function liαi⩾0 could be selected as γieαit or 1/2αi2, γi>0.By using (8), the inequality state constraints could be transformed into equality state constraints. For example, consider a single inequality constraint:(9)Cimin−Cix⩽0.If(10)liαi=γieαit,then γi is chosen that(11)γi=−ηi0.In addition, the relationship between the time derivative ofαit and the input could be obtained by differentiating (8) pipi⩾1 times until the appearance of input u:(12)∂li∂αiα˙i+ηi1x=0,∂2li∂2αiα˙i2+∂li∂αiα¨i+ηi2x=0,⋮∂li∂αiαipi+lipiαi−∂li∂αiαipi+ηipix,u=0,where ηij is the jth derivative of ηix1⩽j⩽pi and lipiαi is the pith derivative of liαi.Denote thejth order derivative of αi as αij, and decompose ηipix,u as ηipix,u=Γi1x+Γi2xu. If there exist k (k⩾1) state constraints and Ns(Ns=∑i=1kpi) slack variables, then combining (12) with system (1) could form the following augmented system:(13)x˙s=fsxs+gsxsu,where xs=x,α1,α11,α12,…,α1p1,…,αk,αk1,αk2,…,αkpkT,(14)fsxs=fxα11α12⋮−∂l1∂α1−1Γ11x+l1p1α1−∂l1∂α1α1p1⋮αk1αk2⋮−∂lk∂αk−1Γk1x+lkpkαk−∂lk∂αkαkpk,(15)gsxs=gx00⋮−∂l1∂α1−1Γ12x⋮00⋮−∂lk∂αk−1Γk2x.From (8), for each i, if αit exists and is bounded, then the state constraints are not violated. Therefore, regarding the following performance index, if there exists an optimal controller u∗ for system (13) without state constraints, then u∗ is also optimal for the state-constrained system (1):(16)Vx0=∫0tfWsxsτ+Muτdτ,where Wsxs=xsTQsxs and Qs=diagQ,0 with a dimension n+Ns.


### 3.2. Handling Input Constraints by Using Nonquadratic Performance Integrand

In order to incorporate the control constraints in the optimal control design, the nonquadratic sufficiently differentiable performance integrand could be applied [51]:(17)Mu=2∫0uH−TvRdv,where H⋅=h1⋅,h2⋅,…,hm⋅T are vectors of bounded, monotonic Cϵ functions (ζ⩾1 and ϵ⩾1) to approximate the hard constraints(18)umin⩽H⋅⩽umax.The corresponding Hamiltonian and HJB equations using the nonquadratic performance integrand are(19)Vxs=xsTQsxs+∫0uH−TvRdv+λTfsxs+gsxsu,(20)VxsTfxs−gxsH12R−1gTxsVxs+xsTQsxs+2∫0−H12R−1gTxsVxsH−1vTRdv=0,where λt∈ℝn is the Lagrange multiplier and V0=0.If the solutionV∗x of the HJB equation exists, the optimal control law is(21)u∗x=−H12R−1gTxsVxs∗,which can be approximated by using a computational approach based on PI [38]:(1)
Choose an initial stabilizing admissible controlleru0xs.(2)
Fori⩾0, solve the following LE for Vi:(22)VxsTfsxs+gsxsuixs+xsTQsxs+2∫0uixsH−1vTRdv=0.(3)
Update the control law:(23)ui+1xs=−H12R−1gTxsVxs.


### 3.3. Approximated Optimal Control Design

In Section 3.1 and 3.2, the problems caused by the constraints on the states and inputs are eliminated, respectively. However, the PI denoted by (22) and (23) relies on the full knowledge of the system dynamics. To eliminate this requirement, a computational approach [50] is applied and tailored. For denotation simplicity, in this subsection, xs, fsxs, gsxs, and Wsxs are denoted as x, fx, gx, and Wx, respectively.To start with, decomposeu as(24)u=ui+vi,then(25)x˙s=fsxs+gsxsui+gsxsvi,where u is the control forced actually on the system and ui is the control policy to be improved iteratively.Consider a cost functionVi satisfying (22):(26)V˙i=VxsiTfsxs+gsxsui+gsxsvi=−Wsxs−Mui−2H−Tui+1Rvi.IntegrateV˙i on a time interval t,t+Δt:(27)Vit+Δt−Vit=∫tt+Δt−Wsxs−Mui−2H−Tui+1Rvidτ.Assumeϕjxsj=1∞ and ψjxsj=1∞ are two infinite sequences of linearly independent smooth basis functions on X, where ϕj0=0 and ψj0=0 for all j=1,2,…. Approximate the cost function and the control policy as(28)V^i=∑j=1N1c^jiϕjxs=ΦxsTc^i,where c^i=c^1ic^2i⋯c^N1iT, Φx=ϕ1xϕ2x⋯ϕN1xT, and(29)u^i=−HR−Tϒ^ixs,where(30)ϒ^ixs=k^iΨxs,with(31)k^i=k^11ik^12i⋯k^1N2ik^21ik^22i⋯k^2N2i⋮k^m1ik^m2i⋯k^mN2i,(32)Ψxs=ψ1xsψ2xs⋯ψN2xsT.In (28)–(32), N1>0 and N2>0 are two sufficiently large integers. Substitute equations (28)–(32) into (27), and denote t and t+Δt as tk and tk+1:(33)∑j=1N1c^jiϕjxstk+1−ϕjxstk−∫tktk+1k^i+1Ψxsujxs−ujixsdt=−∫tktk+1Wsxs+Muidt+εi,k.In (33), there are N1+N2 variables to be determined in c^i and k^ji+1j=1,2,…,m, considering a sufficient long time sequence tkk=0L with L⩾N1+N2; then, according to (33),(34)IΦ−IΨu+IΨuic^iveck^i+1T=−IW−IMui+Δi,where(35)IΦ=Φxs1T−Φxs0TΦxs2T−Φxs1T⋮ΦxsLT−ΦxsL−1TL×N1,(36)IΨu=∫t0t1uTIm⊗ΨxsTdt∫t1t2uTIm⊗ΨxsTdt⋮∫tL−1tLuTIm⊗ΨxsTdtL×mN2,(37)IΨui=∫t0t1H12R−1k^iΨxsTIm⊗ΨxsTdt∫t1t2H12R−1k^iΨxsTIm⊗ΨxsTdt⋮∫tL−1tLH12R−1k^iΨxsTIm⊗ΨxsTdtL×mN2,(38)IW=∫t0t1Wxsdt∫t1t2Wxsdt⋮∫tL−1tLWxsdtL×1,(39)IMui=t0∫t1Muidt∫t1t2Muidt⋮∫tL−1tLMuidtL×1,and Δi=εi,1εi,1⋯εi,LT∈ℝL is the approximation error.From (34), if(40)rankΘk^i=N1+mN2,where Θk^i=IΦ−IΨu+IΨui, then c^i and k^i+1 can be obtained as(41)c^iveck^i+1T=Θk^iTΘk^i−1Θk^iTΠk^i,where Πk^i=−IW−IMui.Equation (41) along with (28)–(32) serves as a computational approach to approximate Vi and generate an improved controller ui+1 when an admissible control policy ui is given. By using this approach iteratively, starting from an admissible bounded control u0, a control sequence uixi=0∞ could be obtained, and it satisfies the following theorem.Theorem 1.
Consider system (1) with both state and input constraints, and augment it as (13) using an appropriate number of slack functions; if there exists an optimal controller and Assumption 1-2 hold, then(42)ui+1=−HR−Tk^iΨxs,which is generated by the PI ((41), (28), and (29)), is an admissible control sequence for (1) on X. The cost function Vi=Jt0u=ui satisfies the following Lyapunov equation (LE)(43)VxsiTfsxs+gsxsui+Wsxs+Mui=0,and then V∗xs⩽Vi+1xs⩽Vixs for ∀xs∈X. In addition, Vi⟶V∗ and ui⟶u∗ uniformly on X.Proof.
The proof of this theorem follows the same lines of reasoning as in the proof of Lemma 1 and Theorem 1 in [38] and is omitted here for brevity.Assumption 1.
There existL0>0 and δ>0, such that for all L>L0,(44)1L∑k=1Lθi,kTθi,k⩾IN1+N2,where θi,k is the kth row of Θk^i.Assumption 2.
The closed-loop system is composed of (18), and(45)u=u0+eis ISS (input-to-state stable) [52] when the exploration noise e is considered as the input.


## 3.1. Handling State Constraints Using Slack Functions

Consider the state constraintCmin⩽Cx⩽Cmax; it contains the following inequality constraints:(7)Cmin−C⩽0,C−Cmax⩽0.There exists slack functionliαit ([10, 47]), which satisfy, for each state constraint, (8)ηix+liαit=0,where ηix could be Cix−Cimax, Cimin−Cix, or Cix−CimaxCimin−Cix, and according to the structure of ηix and the value of Cimin and Cimax, the slack function liαi⩾0 could be selected as γieαit or 1/2αi2, γi>0.By using (8), the inequality state constraints could be transformed into equality state constraints. For example, consider a single inequality constraint:(9)Cimin−Cix⩽0.If(10)liαi=γieαit,then γi is chosen that(11)γi=−ηi0.In addition, the relationship between the time derivative ofαit and the input could be obtained by differentiating (8) pipi⩾1 times until the appearance of input u:(12)∂li∂αiα˙i+ηi1x=0,∂2li∂2αiα˙i2+∂li∂αiα¨i+ηi2x=0,⋮∂li∂αiαipi+lipiαi−∂li∂αiαipi+ηipix,u=0,where ηij is the jth derivative of ηix1⩽j⩽pi and lipiαi is the pith derivative of liαi.Denote thejth order derivative of αi as αij, and decompose ηipix,u as ηipix,u=Γi1x+Γi2xu. If there exist k (k⩾1) state constraints and Ns(Ns=∑i=1kpi) slack variables, then combining (12) with system (1) could form the following augmented system:(13)x˙s=fsxs+gsxsu,where xs=x,α1,α11,α12,…,α1p1,…,αk,αk1,αk2,…,αkpkT,(14)fsxs=fxα11α12⋮−∂l1∂α1−1Γ11x+l1p1α1−∂l1∂α1α1p1⋮αk1αk2⋮−∂lk∂αk−1Γk1x+lkpkαk−∂lk∂αkαkpk,(15)gsxs=gx00⋮−∂l1∂α1−1Γ12x⋮00⋮−∂lk∂αk−1Γk2x.From (8), for each i, if αit exists and is bounded, then the state constraints are not violated. Therefore, regarding the following performance index, if there exists an optimal controller u∗ for system (13) without state constraints, then u∗ is also optimal for the state-constrained system (1):(16)Vx0=∫0tfWsxsτ+Muτdτ,where Wsxs=xsTQsxs and Qs=diagQ,0 with a dimension n+Ns.


## 3.2. Handling Input Constraints by Using Nonquadratic Performance Integrand

In order to incorporate the control constraints in the optimal control design, the nonquadratic sufficiently differentiable performance integrand could be applied [51]:(17)Mu=2∫0uH−TvRdv,where H⋅=h1⋅,h2⋅,…,hm⋅T are vectors of bounded, monotonic Cϵ functions (ζ⩾1 and ϵ⩾1) to approximate the hard constraints(18)umin⩽H⋅⩽umax.The corresponding Hamiltonian and HJB equations using the nonquadratic performance integrand are(19)Vxs=xsTQsxs+∫0uH−TvRdv+λTfsxs+gsxsu,(20)VxsTfxs−gxsH12R−1gTxsVxs+xsTQsxs+2∫0−H12R−1gTxsVxsH−1vTRdv=0,where λt∈ℝn is the Lagrange multiplier and V0=0.If the solutionV∗x of the HJB equation exists, the optimal control law is(21)u∗x=−H12R−1gTxsVxs∗,which can be approximated by using a computational approach based on PI [38]:(1)
Choose an initial stabilizing admissible controlleru0xs.(2)
Fori⩾0, solve the following LE for Vi:(22)VxsTfsxs+gsxsuixs+xsTQsxs+2∫0uixsH−1vTRdv=0.(3)
Update the control law:(23)ui+1xs=−H12R−1gTxsVxs.


## 3.3. Approximated Optimal Control Design

In Section 3.1 and 3.2, the problems caused by the constraints on the states and inputs are eliminated, respectively. However, the PI denoted by (22) and (23) relies on the full knowledge of the system dynamics. To eliminate this requirement, a computational approach [50] is applied and tailored. For denotation simplicity, in this subsection, xs, fsxs, gsxs, and Wsxs are denoted as x, fx, gx, and Wx, respectively.To start with, decomposeu as(24)u=ui+vi,then(25)x˙s=fsxs+gsxsui+gsxsvi,where u is the control forced actually on the system and ui is the control policy to be improved iteratively.Consider a cost functionVi satisfying (22):(26)V˙i=VxsiTfsxs+gsxsui+gsxsvi=−Wsxs−Mui−2H−Tui+1Rvi.IntegrateV˙i on a time interval t,t+Δt:(27)Vit+Δt−Vit=∫tt+Δt−Wsxs−Mui−2H−Tui+1Rvidτ.Assumeϕjxsj=1∞ and ψjxsj=1∞ are two infinite sequences of linearly independent smooth basis functions on X, where ϕj0=0 and ψj0=0 for all j=1,2,…. Approximate the cost function and the control policy as(28)V^i=∑j=1N1c^jiϕjxs=ΦxsTc^i,where c^i=c^1ic^2i⋯c^N1iT, Φx=ϕ1xϕ2x⋯ϕN1xT, and(29)u^i=−HR−Tϒ^ixs,where(30)ϒ^ixs=k^iΨxs,with(31)k^i=k^11ik^12i⋯k^1N2ik^21ik^22i⋯k^2N2i⋮k^m1ik^m2i⋯k^mN2i,(32)Ψxs=ψ1xsψ2xs⋯ψN2xsT.In (28)–(32), N1>0 and N2>0 are two sufficiently large integers. Substitute equations (28)–(32) into (27), and denote t and t+Δt as tk and tk+1:(33)∑j=1N1c^jiϕjxstk+1−ϕjxstk−∫tktk+1k^i+1Ψxsujxs−ujixsdt=−∫tktk+1Wsxs+Muidt+εi,k.In (33), there are N1+N2 variables to be determined in c^i and k^ji+1j=1,2,…,m, considering a sufficient long time sequence tkk=0L with L⩾N1+N2; then, according to (33),(34)IΦ−IΨu+IΨuic^iveck^i+1T=−IW−IMui+Δi,where(35)IΦ=Φxs1T−Φxs0TΦxs2T−Φxs1T⋮ΦxsLT−ΦxsL−1TL×N1,(36)IΨu=∫t0t1uTIm⊗ΨxsTdt∫t1t2uTIm⊗ΨxsTdt⋮∫tL−1tLuTIm⊗ΨxsTdtL×mN2,(37)IΨui=∫t0t1H12R−1k^iΨxsTIm⊗ΨxsTdt∫t1t2H12R−1k^iΨxsTIm⊗ΨxsTdt⋮∫tL−1tLH12R−1k^iΨxsTIm⊗ΨxsTdtL×mN2,(38)IW=∫t0t1Wxsdt∫t1t2Wxsdt⋮∫tL−1tLWxsdtL×1,(39)IMui=t0∫t1Muidt∫t1t2Muidt⋮∫tL−1tLMuidtL×1,and Δi=εi,1εi,1⋯εi,LT∈ℝL is the approximation error.From (34), if(40)rankΘk^i=N1+mN2,where Θk^i=IΦ−IΨu+IΨui, then c^i and k^i+1 can be obtained as(41)c^iveck^i+1T=Θk^iTΘk^i−1Θk^iTΠk^i,where Πk^i=−IW−IMui.Equation (41) along with (28)–(32) serves as a computational approach to approximate Vi and generate an improved controller ui+1 when an admissible control policy ui is given. By using this approach iteratively, starting from an admissible bounded control u0, a control sequence uixi=0∞ could be obtained, and it satisfies the following theorem.Theorem 1.
Consider system (1) with both state and input constraints, and augment it as (13) using an appropriate number of slack functions; if there exists an optimal controller and Assumption 1-2 hold, then(42)ui+1=−HR−Tk^iΨxs,which is generated by the PI ((41), (28), and (29)), is an admissible control sequence for (1) on X. The cost function Vi=Jt0u=ui satisfies the following Lyapunov equation (LE)(43)VxsiTfsxs+gsxsui+Wsxs+Mui=0,and then V∗xs⩽Vi+1xs⩽Vixs for ∀xs∈X. In addition, Vi⟶V∗ and ui⟶u∗ uniformly on X.Proof.
The proof of this theorem follows the same lines of reasoning as in the proof of Lemma 1 and Theorem 1 in [38] and is omitted here for brevity.Assumption 1.
There existL0>0 and δ>0, such that for all L>L0,(44)1L∑k=1Lθi,kTθi,k⩾IN1+N2,where θi,k is the kth row of Θk^i.Assumption 2.
The closed-loop system is composed of (18), and(45)u=u0+eis ISS (input-to-state stable) [52] when the exploration noise e is considered as the input.


## 4. Case Study

In this section, two cases, including a CSTR (continuous stirred-tank reactor) model [14], are selected to test the feasibility and performance of the proposed control design.


### 4.1. A Nonlinear System

Consider the following nonlinear system:(46)x˙1=x2,x˙2=1−x12x2−x1+u,with state and control constraints x2⩾−0.4 and −1⩽u⩽1, and initial condition x1,x2=1,0.To start with, one slack function was introduced:(47)−x2−0.4+γeαt=0.According to (11) and (12), γ=0.4 and α0=0. Denote x3=α, and combine x˙3 with the original system:(48)x˙1=x2,x˙2=1−x12x2−x1+u,x˙3=10.4ex31−x12x2−x1+u.For simulation purpose, we setQs=diag1,1,0, so the augmented state x3 is not taken into account in the value function, R=0.5; so, the value function mainly concerns the deviation between the actual state and its setting point, Hv=tanhv. The basis functions for approximation of value function and control input were selected as polynomials of system states:(49)Φ=x12x22x1x2x14x24x12x22x1x23x13x2,Ψ=x1x2x1x2x1x22x12x2.The initial weightk^0=0.8−1.520−1. The initial weight was determined using the structural information of system model and the basis functions. The actual input was set such that u=−H0.5R−1k^0Ψ+0.008e with the exploration noise(50)e=rand⋅rand⋅sin0.1t+rand⋅sint+rand⋅sin5t+rand⋅sin10t.The ADP iteration started after the learning finished att=20s. After 20 iterations, the difference between k^20 and k^19 was less than 10−7, and the ADP algorithm stopped. The resultant approximated optimal controller was k^20=0.45863.391418.310719.56157.1479.The state trajectories driven by the initial admissible controller and the approximated optimal controller were given in Figure1. The approximated cost difference was illustrated in Figure 2. The control input using the two controllers were shown in Figure 3. From the figures, both the initial admissible controller and the approximated optimal controller can force the system to the equilibrium point x1=0 and x2=0. However, the state trajectories of x1 and x2 driven by the approximated optimal control has less fluctuations and less approximated control cost compared with the initial admissible controller.Figure 1
State trajectories driven by initial admissible control and approximated optimal control.Figure 2
Cost difference between initial admissible control and approximated optimal control.Figure 3
Control inputs using initial admissible controller and approximated optimal controller.


### 4.2. A CSTR Model

Consider a reversible reaction2A⇌B taken place in a CSTR, and the forward and reverse kinetic constants with respect to species B are ka and kb, respectively. The dynamics of the process is(51)x˙1=−2kax12−c1x1+2kbx2+FVu,x˙2=kax12+c2x1−FV+kbx2,where x1=CA−C¯A, x2=CB−C¯B, and u=CAin−C¯Ain are the deviation of the outlet concentrations of species A and B and the inlet concentration of species A, respectively. c1=4kaC¯A+F/V and c2=2kaC¯A are some constant variables. V and F are the volume and inlet flow-rate of the reactor, respectively. The constraints, initial conditions, and the values of model parameters were shown in Table 1.Table 1
Model parameters, initial conditions, and constraints.
Variable
Value
k
a
0.05
k
b
0.01
c
1
0.338
c
2
0.159
F
/
V
0.02
C
¯
A
1.59
C
¯
B
4.21
C
A
in
10
x
1
0
0
x
2
0
−1.71
x
1
t
[−0.09, 0.16]
x
2
t
⩾
−
4.21
u
t
[−10, 10]According to Section3.1, two slack functions were introduced first to transform the inequality constraints into equality constraints:(52)x1+0.09x1−0.16+γ1eα1t=0,(53)−x2−4.21+γ2eα2t=0.The values ofγ1 and γ2 were calculated by setting both x1, α1 and x2, α2 as 0: γ1=0.0144 and γ2=4.21. The initial value of α1 and α2 were calculated using the initial value of x1 and x2: α10=0 and α20=−0.5212. By differentiating (51) and (52), the derivatives of α1 and α2 were obtained:(54)α˙1=−2x1−0.070.0144eα1−2kax12−c1x1+2kbx2+FVu,(55)α˙2=14.21eα2kax12+c2x1−FV+kbx2.Then, (54) was differentiated until the appearance of input u which is present in x˙1:(56)α¨2=−α˙22+14.21eα20.1x1+0.159x˙1−0.03x˙2,and α˙20=0.0205.Denoteα1, α2, and α˙2 as x3, x4, and x5, and combining the original system with the derived equality constraints, the following augmented system was obtained:(57)x˙1=−2kax12−c1x1+2kbx2+FVu,x˙2=kax12+c2x1−FV+kbx2,x˙3=−2x1−0.070.0144ex3−2kax12−c1x1+2kbx2+FVu,x˙4=x5,x˙5=−x42+14.21ex40.1x1+0.159−2kax12−c1x1+2kbx2+FVu−0.03kax12+c2x1−FV+kbx2.In the simulation,Qs=diag5,5,0,0,0; so, the augmented states were not considered when evaluating the control performance. Compared with the evolution of outlet concentrations to their set points in time, the selection of input variables is not that important, so R=0.5 and Hv=10tanhv/10. The basis functions in Φ include x12, x22, x1x2, x14, x24, x12x22, x1x23, x13x2, x16, x26, x1x25, x12x24, x13x23, x14x22, x15x2, and Ψ=x1x2x12x22. The initial weight k^0=−1−1−1−1. The exploration noise was set such that u=−H0.5R−1k^0Ψ+0.1e with e=rand⋅rand⋅sin0.1t+rand⋅sint+rand⋅sin5t+rand⋅sin10t.The ADP iteration started after the learning finished att=40s. After 20 iterations, the difference between k^20 and k^19 was less than 10−7, and the ADP algorithm stopped. The resultant approximated optimal controller was k^20=−6.13732.7789−3.20451.4604.The state trajectories ofx1 and x2 under the initial admissible controller and the approximated optimal controller were given in Figures 4 and 5, respectively. The approximated cost difference was illustrated in Figure 6. The control input using the two controllers were shown in Figure 7. Compared with the initial admissible controller, the approximated optimal controller could force the CSTR to the desired working point CA=1.59 and CB=4.21 with less steady-state error and better performance regarding to the performance index.Figure 4
Trajectory ofCA driven by initial admissible control and approximated optimal control.Figure 5
Trajectory ofCB driven by initial admissible control and approximated optimal control.Figure 6
Cost difference between initial admissible control and approximated optimal control.Figure 7
Control inputs using initial admissible controller and approximated optimal controller.In the simulation, structural knowledge of the system model was required to derive the explicit formulations of (47) and (56). However, the numerical value of the slack variables can be obtained without knowing the full system dynamics; only the model structure is required. Thus, via learning from online input and state information, the performance of the controller was improved.


## 4.1. A Nonlinear System

Consider the following nonlinear system:(46)x˙1=x2,x˙2=1−x12x2−x1+u,with state and control constraints x2⩾−0.4 and −1⩽u⩽1, and initial condition x1,x2=1,0.To start with, one slack function was introduced:(47)−x2−0.4+γeαt=0.According to (11) and (12), γ=0.4 and α0=0. Denote x3=α, and combine x˙3 with the original system:(48)x˙1=x2,x˙2=1−x12x2−x1+u,x˙3=10.4ex31−x12x2−x1+u.For simulation purpose, we setQs=diag1,1,0, so the augmented state x3 is not taken into account in the value function, R=0.5; so, the value function mainly concerns the deviation between the actual state and its setting point, Hv=tanhv. The basis functions for approximation of value function and control input were selected as polynomials of system states:(49)Φ=x12x22x1x2x14x24x12x22x1x23x13x2,Ψ=x1x2x1x2x1x22x12x2.The initial weightk^0=0.8−1.520−1. The initial weight was determined using the structural information of system model and the basis functions. The actual input was set such that u=−H0.5R−1k^0Ψ+0.008e with the exploration noise(50)e=rand⋅rand⋅sin0.1t+rand⋅sint+rand⋅sin5t+rand⋅sin10t.The ADP iteration started after the learning finished att=20s. After 20 iterations, the difference between k^20 and k^19 was less than 10−7, and the ADP algorithm stopped. The resultant approximated optimal controller was k^20=0.45863.391418.310719.56157.1479.The state trajectories driven by the initial admissible controller and the approximated optimal controller were given in Figure1. The approximated cost difference was illustrated in Figure 2. The control input using the two controllers were shown in Figure 3. From the figures, both the initial admissible controller and the approximated optimal controller can force the system to the equilibrium point x1=0 and x2=0. However, the state trajectories of x1 and x2 driven by the approximated optimal control has less fluctuations and less approximated control cost compared with the initial admissible controller.Figure 1
State trajectories driven by initial admissible control and approximated optimal control.Figure 2
Cost difference between initial admissible control and approximated optimal control.Figure 3
Control inputs using initial admissible controller and approximated optimal controller.


## 4.2. A CSTR Model

Consider a reversible reaction2A⇌B taken place in a CSTR, and the forward and reverse kinetic constants with respect to species B are ka and kb, respectively. The dynamics of the process is(51)x˙1=−2kax12−c1x1+2kbx2+FVu,x˙2=kax12+c2x1−FV+kbx2,where x1=CA−C¯A, x2=CB−C¯B, and u=CAin−C¯Ain are the deviation of the outlet concentrations of species A and B and the inlet concentration of species A, respectively. c1=4kaC¯A+F/V and c2=2kaC¯A are some constant variables. V and F are the volume and inlet flow-rate of the reactor, respectively. The constraints, initial conditions, and the values of model parameters were shown in Table 1.Table 1
Model parameters, initial conditions, and constraints.
Variable
Value
k
a
0.05
k
b
0.01
c
1
0.338
c
2
0.159
F
/
V
0.02
C
¯
A
1.59
C
¯
B
4.21
C
A
in
10
x
1
0
0
x
2
0
−1.71
x
1
t
[−0.09, 0.16]
x
2
t
⩾
−
4.21
u
t
[−10, 10]According to Section3.1, two slack functions were introduced first to transform the inequality constraints into equality constraints:(52)x1+0.09x1−0.16+γ1eα1t=0,(53)−x2−4.21+γ2eα2t=0.The values ofγ1 and γ2 were calculated by setting both x1, α1 and x2, α2 as 0: γ1=0.0144 and γ2=4.21. The initial value of α1 and α2 were calculated using the initial value of x1 and x2: α10=0 and α20=−0.5212. By differentiating (51) and (52), the derivatives of α1 and α2 were obtained:(54)α˙1=−2x1−0.070.0144eα1−2kax12−c1x1+2kbx2+FVu,(55)α˙2=14.21eα2kax12+c2x1−FV+kbx2.Then, (54) was differentiated until the appearance of input u which is present in x˙1:(56)α¨2=−α˙22+14.21eα20.1x1+0.159x˙1−0.03x˙2,and α˙20=0.0205.Denoteα1, α2, and α˙2 as x3, x4, and x5, and combining the original system with the derived equality constraints, the following augmented system was obtained:(57)x˙1=−2kax12−c1x1+2kbx2+FVu,x˙2=kax12+c2x1−FV+kbx2,x˙3=−2x1−0.070.0144ex3−2kax12−c1x1+2kbx2+FVu,x˙4=x5,x˙5=−x42+14.21ex40.1x1+0.159−2kax12−c1x1+2kbx2+FVu−0.03kax12+c2x1−FV+kbx2.In the simulation,Qs=diag5,5,0,0,0; so, the augmented states were not considered when evaluating the control performance. Compared with the evolution of outlet concentrations to their set points in time, the selection of input variables is not that important, so R=0.5 and Hv=10tanhv/10. The basis functions in Φ include x12, x22, x1x2, x14, x24, x12x22, x1x23, x13x2, x16, x26, x1x25, x12x24, x13x23, x14x22, x15x2, and Ψ=x1x2x12x22. The initial weight k^0=−1−1−1−1. The exploration noise was set such that u=−H0.5R−1k^0Ψ+0.1e with e=rand⋅rand⋅sin0.1t+rand⋅sint+rand⋅sin5t+rand⋅sin10t.The ADP iteration started after the learning finished att=40s. After 20 iterations, the difference between k^20 and k^19 was less than 10−7, and the ADP algorithm stopped. The resultant approximated optimal controller was k^20=−6.13732.7789−3.20451.4604.The state trajectories ofx1 and x2 under the initial admissible controller and the approximated optimal controller were given in Figures 4 and 5, respectively. The approximated cost difference was illustrated in Figure 6. The control input using the two controllers were shown in Figure 7. Compared with the initial admissible controller, the approximated optimal controller could force the CSTR to the desired working point CA=1.59 and CB=4.21 with less steady-state error and better performance regarding to the performance index.Figure 4
Trajectory ofCA driven by initial admissible control and approximated optimal control.Figure 5
Trajectory ofCB driven by initial admissible control and approximated optimal control.Figure 6
Cost difference between initial admissible control and approximated optimal control.Figure 7
Control inputs using initial admissible controller and approximated optimal controller.In the simulation, structural knowledge of the system model was required to derive the explicit formulations of (47) and (56). However, the numerical value of the slack variables can be obtained without knowing the full system dynamics; only the model structure is required. Thus, via learning from online input and state information, the performance of the controller was improved.


## 5. Conclusions

An approximated optimal control approach for chemical processes with both input and state constraints was proposed. Its feasibility and performance were tested via two nonlinear examples. The proposed approach requires only the structure knowledge of the system model. The value information of the model parameters is not needed. This indicates the proposed approach can be applied to practical systems or a single unit (like a CSTR in this study) with determined model structure but has unknown parameters under some working conditions. However, the proposed approach is applicable only when Assumptions 1 and 2 hold. The global stability issue still needs future study.



---

*Source: 1396913-2019-10-28.xml*