# Hybrid Functions Direct Approach and State Feedback Optimal Solutions for a Class of Nonlinear Polynomial Time Delay Systems

Authors: Mohamed Karim Bouafoura; Naceur Benhadj Braiek
Journal: Complexity
 (2019)

---

## Abstract

The aim of this paper is to determine the optimal open loop solution and a nonlinear delay-dependent state feedback suboptimal control for a class of nonlinear polynomial time delay systems. The proposed method uses a hybrid of block pulse functions and Legendre polynomials as an orthogonal base for system’s states and input expansion. Hence, the complex dynamic optimization problem is then reduced, with the help of operational properties of the hybrid basis and Kronecker tensor product lemmas, to a nonlinear programming problem that could be solved with available NLP solvers. A practical nonlinear feedback controller gains are deduced with respect to a least square formalism based on the optimal open loop control results. Simulation results show efficiency of the proposed numerical optimal approach.


---

## Body

## 1. Introduction

Time delays affect systems dynamics in many engineering applications like chemical control systems, biology, and medicine [1, 2]. Delays are also encountered in communication and information technologies like high-speed communication networks [3]. It should be noted that time delay may be, in some applications like communication lines, a source of instability and performance degradation [4]. Time delay system is therefore a very important class of processes whose stabilization [5] and optimization [6, 7] have been of interest to many researchers.Particularly, many attempts have been made in literature to solve optimal control problems for many classes of linear [8–11] and nonlinear [6, 12, 13] time delay systems. Among them, we recall the application of Pontryagins maximum principle to the optimization of control systems with time delays which was firstly proposed by [14]. It had been shown that it results in a system of coupled two-point boundary-value (TPBV) problem involving both delay and advance terms whose exact solution, except in very special cases, is very difficult to determine (see [15]). Perhaps one of the most effective techniques is dynamic programming approaches (see [2]) for overcoming the complexity of the nonlinear time delay systems in optimal control problems. Of course, application of dynamic programming methods has some difficulties due to the need to provide an appropriate level model and also to define recursive relationships for each case problem. Also a computational algorithm considering a linear approximation of the original system which is defined about a nominal trajectory is offered by [16]. Clearly, using the linear approximation is not reliable and may lead to large errors. Reference [7] proposed an approach based on discretization techniques and necessary conditions to obtain approximate optimal control and the state for optimal control problems with nonlinear delay systems. Despite the good performance of this method, achieving the necessary conditions in some problems and the implementation of approach may be faced with difficulties. So different numerical methods have been proposed to avoid the problems arising from the applications of analytical methods. It is then straightforward that many of the numerical methods dedicated to solving classical optimal control problems have been extended to handle optimal control problems governed by time delay systems.Typically, direct methods are based on converting the dynamical optimal control problem into static optimization problem. Among direct methods, parametrization technique [17–19] is known to minimize decision variables compared to the discretization of the problem [7]. It is worth noting that parametrization relies basically on orthogonal functions or wavelets [20–22]; however that tool have been used to solve various other problems of dynamic systems like identification (see [8]), tracking control (see [23]), observer based control (see [24]), or minimum time control (see [25]). The main characteristic of this pseudo-spectral technique is that it allows transforming complex dynamic optimization problems to solving a set of algebraic equations in the least square sense in the linear systems case [26, 27] or permits formulating an equivalent nonlinear static programming problem for problems related to nonlinear systems [13, 28, 29].In recent years, a growing interest has been appeared toward the application of hybrid functions, which is a combination of block pulse and an orthogonal polynomials basis [26]. In the nonlinear time delay optimal control problems context, an approach using hybrid functions which consist of block pulse functions and orthonormal Taylor series (see [15, 29]) had been proposed, where authors propose to solve the necessary and sufficient condition equations for stationary emanating from the Hamiltonian based on state and control coefficients over the basis. Similarly, [28] propose a direct approach based on a hybrid of block pulse functions and Lagrange interpolating polynomials in order to convert the original optimal problem containing multiple delay into a mathematical programming one, where the resulting optimization problem is solved numerically by the Lagrange multipliers method. Reference [27] proposed similar approach based on hybrid functions of block pulse and Bernouilli polynomials, while [30] uses biorthogonal cubic Hermite spline multiwavelets in addition to block pulse functions to constitute the hybrid basis. Although above contributions treat some nonlinear delayed optimal control problem, they do not propose any general nonlinear programming problem that could handle all examples depicted in their works. In fact, for each considered nonlinear system, a nonlinear optimization problem is formulated and then solved with an NLP solver. Furthermore, only open loop control solutions are investigated therein, which is not of great interest in practice.In the present paper, we introduce a direct method to solve forwardly the finite time quadratic optimal control problem of polynomial systems with delayed state, by the use of hybrid functions of block pulse and Legendre polynomials. The operational matrices of delay and Kronecker product specific to that basis are recalled. At first, the open loop solution of the nonlinear time delay optimal control problem is investigated. Secondly, a suboptimal nonlinear state feedback is determined based on the first part results. Hence, the main contributions in this work could be summarized as follows:(a) expressing the constraint of the formulated NLP problem properly for the class of polynomial systems; thus the proposed formulation could handle a wide range of nonlinear analytic nonlinear systems. Then, a unified development is carried for that class of systems,(b) deriving a nonlinear polynomial delay-dependent nonlinear suboptimal state feedback that reproduce the optimal state trajectories determined in the open loop framework,(c) using an hybrid basis with reduced number of elementary functions, which makes open loop synthesis faster, with a good enough accuracy compared to other approaches, and closed loop solution within a simpler formulation and resolution.The remainder of the paper is organised as follows. In the second section, hybrid functions and their properties are introduced. In the third section, the open loop numerical solution of the nonlinear time delay optimal control problem is detailed. The suboptimal closed loop framework is presented in the fourth section. In the fifth section, computational results are depicted. Finally, concluding remarks and future works are presented.


## 2. Hybrid Functions

Hybrid functionshij(t), i=1,2,⋯N, j=0,1⋯M-1, have three arguments; i and j are the order of block pulse functions and Legendre polynomials, respectively, and t is the normalized time. They are defined on the interval t=[0,tf] as [26](1)hijt=Lj2Ntft-2i+1,t∈i-1Ntf,iNtf,0otherwise.Here,Lj(t) are the well-known Legendre polynomials of order M which constitute the base L(t) and satisfy the following recursive formula:(2)L0t=1,L1t=t,Lj+1t=2j+1j+1tLjt-jj+1Lj-1t,j=1,2,3,⋯We defineΦ(t) the vector of N block pulse functions ϕi(t), i=0,1,⋯N-1, as follows:(3)ϕit=1,∀t∈i-1Ntf,iNtf,0,otherwise.Sincehij(t) is the combination of Legendre polynomials and block pulse functions which are both complete and orthogonal, then the set of hybrid functions is a complete orthogonal system.In the rest of the paper we notatew=NM the dimension of the hybrid basis.


### 2.1. Operational Matrix of Integration

The integration ofh(t) can be approximated by [26](4)∫0thνdν≅PhtwhereP is the integration operational matrix of order w×w(5)Pw×w=THH⋯H0TH⋯H00T⋯H⋮⋮⋮⋱⋮000⋯Twhere(6)HM×M=tfN100⋯0000⋯0000⋯0⋮⋮⋮⋮⋮000⋯0and(7)TM×M=tf2N.1100⋯000-130130⋯0000-15015⋯000⋮⋮⋮⋮⋮⋮⋮⋮0000⋯-12N-3012N-30000⋯0-12N-10


### 2.2. Delay Modeling with Hybrid Functions

A vector functiong(t) of r dimensional components which are square integrable in [0,tf] can be approximated by a block pulse series as(8)gt≅∑i=1Ngiϕit=GTΦtwhereG=[g1,g2,…,gN]T.For anr component delay vector variable g(t-τ) with(9)gt=ζt∀t∈-τ,0the block pulse series approximation ofg(t-τ) may be defined as [31](10)gt-τ≅∑i=1Ng∗τϕit=G∗τΦtwhere(11)gi∗τ=NT∫itf/Ni+1tf/Ngt-τdt=ζiτfori<μgi-μfori≥μwith(12)ζiτ=NT∫itf/Ni+1tf/Nζt-τdtfori<μandμ is the number of block pulse functions considered over 0≤t≤τ, and G∗(τ)=[g1∗(τ),g2∗(τ),…,gN∗(τ)]T.Let(13)ζμτ=ζ1∗τ,ζ2∗τ,…,ζμ-1∗τTThen, it comes [31](14)vecG∗τ=Er,μvecζμτ+Dr,μvecGE ( r , μ ) and D(r,μ) are called the shift operational matrices, given by(15)Er,μ=Irμ×rμ⋯0rN-μ×rμand(16)Dr,μ=0rμ×rN-μ⋮0rμ×rμ⋯⋮⋯IrN-μ×rN-μ⋮0rN-μ×rμIt is worth noticing that the Shift operational matrices of hybrid functions could be derived forwardly from those of block pulse functions by(17)Eh=IM⊗Er,μDh=IM⊗Dr,μwhere⊗ stands for the Kronecker product.However, it should be noticed that block pulse functions are fundamental for delay modeling. The choice ofN depends on τ and tf, which issue had been addressed in [26]. In this framework, we propose to choose N as follows:(18)N=a.ninttfτ,a∈N∗wherea is a nonnegative integer, to be chosen bigger than one if possible in order to improve approximation and nint(.) denotes the nearest integer function [22] (implemented by round routine in MATLAB).


### 2.3. The Integration of the Cross Product

The integration of the cross product of two hybrid functions vectorsh(t) can be obtained as [26](19)C=∫0tfhthTtdt=L0M⋯0M0ML⋯0M⋮⋮⋱⋮0M0M⋯LwhereC is an w×w matrix. 0M stands for the zero M×M matrix and L is an M×M diagonal matrix that is given by(20)L=tfN10⋯0013⋯0⋮⋮⋱⋮00⋯12M-1


### 2.4. Kronecker Product Operational Matrix

It would be interesting, for bilinear systems, as it will be proven later, to investigate the Kronecker product operational matrix for hybrid functions. This particular matrix operator derivation, as it is the case for the integration, cross product, and delay operators, is highly inspired of the Kronecker operational matrix of both Legendre polynomials and block pulse functions.For the block pulse functions, we can state(21)Φt⊗Φt=ϕ1tΦtϕ2tΦt⋮ϕNtΦt=KΦΦtwith(22)KΦ=E1,1N×NE2,2N×N⋮EN,NN×NwhereKΦ∈RN2×N is the Kronecker product operational matrix of block pulse functions and the matrix Ei,jN×N is defined in Appendix.On the other hand, the product of two Shifted Legendre PolynomialsLi(t) and Lj(t) can be expressed by(23)LitLjt≅∑k=0M-1ψijkLktwith(24)ψijk=2k+1tf∫0tfLitLjtLktdtA practical implementation of the latter scalar products is given in [19].Then, we may write(25)LjtLt=ψi00ψi11⋮ψiM-1M-1=KLiLtwhereKLi is a M×M square matrix.Then it comes(26)Lt⊗Lt=KL0KL1⋮KLM-1=KLLtwhereKL∈RM2×M is the Kronecker product operational matrix of Legendre polynomials.Based on relations (26) and (23), we define(27)ht⊗ht=KL0⊗E1,1N×NKL1⊗E1,1N×N⋮KLM-1⊗E1,1N×N⋯KL0⊗E2,2N×NKL1⊗E2,2N×N⋮KLM-1⊗E2,2N×N⋯⋮⋯KL0⊗EN,NN×NKL1⊗EN,NN×N⋮KLM-1⊗EN,NN×N=KhhtwhereKh∈Rw2×w is the Kronecker product operational matrix of hybrid functions.


## 2.1. Operational Matrix of Integration

The integration ofh(t) can be approximated by [26](4)∫0thνdν≅PhtwhereP is the integration operational matrix of order w×w(5)Pw×w=THH⋯H0TH⋯H00T⋯H⋮⋮⋮⋱⋮000⋯Twhere(6)HM×M=tfN100⋯0000⋯0000⋯0⋮⋮⋮⋮⋮000⋯0and(7)TM×M=tf2N.1100⋯000-130130⋯0000-15015⋯000⋮⋮⋮⋮⋮⋮⋮⋮0000⋯-12N-3012N-30000⋯0-12N-10


## 2.2. Delay Modeling with Hybrid Functions

A vector functiong(t) of r dimensional components which are square integrable in [0,tf] can be approximated by a block pulse series as(8)gt≅∑i=1Ngiϕit=GTΦtwhereG=[g1,g2,…,gN]T.For anr component delay vector variable g(t-τ) with(9)gt=ζt∀t∈-τ,0the block pulse series approximation ofg(t-τ) may be defined as [31](10)gt-τ≅∑i=1Ng∗τϕit=G∗τΦtwhere(11)gi∗τ=NT∫itf/Ni+1tf/Ngt-τdt=ζiτfori<μgi-μfori≥μwith(12)ζiτ=NT∫itf/Ni+1tf/Nζt-τdtfori<μandμ is the number of block pulse functions considered over 0≤t≤τ, and G∗(τ)=[g1∗(τ),g2∗(τ),…,gN∗(τ)]T.Let(13)ζμτ=ζ1∗τ,ζ2∗τ,…,ζμ-1∗τTThen, it comes [31](14)vecG∗τ=Er,μvecζμτ+Dr,μvecGE ( r , μ ) and D(r,μ) are called the shift operational matrices, given by(15)Er,μ=Irμ×rμ⋯0rN-μ×rμand(16)Dr,μ=0rμ×rN-μ⋮0rμ×rμ⋯⋮⋯IrN-μ×rN-μ⋮0rN-μ×rμIt is worth noticing that the Shift operational matrices of hybrid functions could be derived forwardly from those of block pulse functions by(17)Eh=IM⊗Er,μDh=IM⊗Dr,μwhere⊗ stands for the Kronecker product.However, it should be noticed that block pulse functions are fundamental for delay modeling. The choice ofN depends on τ and tf, which issue had been addressed in [26]. In this framework, we propose to choose N as follows:(18)N=a.ninttfτ,a∈N∗wherea is a nonnegative integer, to be chosen bigger than one if possible in order to improve approximation and nint(.) denotes the nearest integer function [22] (implemented by round routine in MATLAB).


## 2.3. The Integration of the Cross Product

The integration of the cross product of two hybrid functions vectorsh(t) can be obtained as [26](19)C=∫0tfhthTtdt=L0M⋯0M0ML⋯0M⋮⋮⋱⋮0M0M⋯LwhereC is an w×w matrix. 0M stands for the zero M×M matrix and L is an M×M diagonal matrix that is given by(20)L=tfN10⋯0013⋯0⋮⋮⋱⋮00⋯12M-1


## 2.4. Kronecker Product Operational Matrix

It would be interesting, for bilinear systems, as it will be proven later, to investigate the Kronecker product operational matrix for hybrid functions. This particular matrix operator derivation, as it is the case for the integration, cross product, and delay operators, is highly inspired of the Kronecker operational matrix of both Legendre polynomials and block pulse functions.For the block pulse functions, we can state(21)Φt⊗Φt=ϕ1tΦtϕ2tΦt⋮ϕNtΦt=KΦΦtwith(22)KΦ=E1,1N×NE2,2N×N⋮EN,NN×NwhereKΦ∈RN2×N is the Kronecker product operational matrix of block pulse functions and the matrix Ei,jN×N is defined in Appendix.On the other hand, the product of two Shifted Legendre PolynomialsLi(t) and Lj(t) can be expressed by(23)LitLjt≅∑k=0M-1ψijkLktwith(24)ψijk=2k+1tf∫0tfLitLjtLktdtA practical implementation of the latter scalar products is given in [19].Then, we may write(25)LjtLt=ψi00ψi11⋮ψiM-1M-1=KLiLtwhereKLi is a M×M square matrix.Then it comes(26)Lt⊗Lt=KL0KL1⋮KLM-1=KLLtwhereKL∈RM2×M is the Kronecker product operational matrix of Legendre polynomials.Based on relations (26) and (23), we define(27)ht⊗ht=KL0⊗E1,1N×NKL1⊗E1,1N×N⋮KLM-1⊗E1,1N×N⋯KL0⊗E2,2N×NKL1⊗E2,2N×N⋮KLM-1⊗E2,2N×N⋯⋮⋯KL0⊗EN,NN×NKL1⊗EN,NN×N⋮KLM-1⊗EN,NN×N=KhhtwhereKh∈Rw2×w is the Kronecker product operational matrix of hybrid functions.


## 3. Numerical Solution of the Nonlinear Time Delay Optimal Control Problem

### 3.1. Description of the Studied System

We consider the nonlinear continuous system which can be represented by the following state space representation:(28)x.t=fxt,τ+gxt,τutxt-τ=x0,∀t∈0,τwherex(t)∈Rn is the state vector, x(t-τ)∈Rn is the delayed state where τ denotes the time delay, u(t)∈Rm is the control vector, and f(x(t,τ)) from Rn into Rn and g(x(t,τ)) from Rn into Rn×m are nonlinear analytic functions of x(t) and x(t-τ).Note that, any functionsf(x(t,τ)) and g(x(t,τ)) could be approached using truncated series of the Kronecker power of x(t) and x(t-τ) as follows:(29)fxt,τ≈∑i=1pFixit+∑j=1qF¯jxjt-τ+∑i=1p∑j=1qΓijxit⊗xjt-τwhereFi∈Rn×nj, F¯j∈Rn×nj and Γij∈Rn×ni+j are constant matrices. x[i](t) denotes the ith Kronecker power of the state vector (see Appendix).And(30)gxt,τ=g1xt,τ⋮g2xt,τ⋮⋯⋮gmxt,τwheregk(x(t,τ))∈Rn, for k=1,…,m, are defined by(31)gkxt,τ≈∑i=0rGkixit+∑j=1sG¯kjxjt-τ+∑i=1r∑j=1sΥkijxit⊗xjt-τwhereGki∈Rn×nj, G¯kj∈Rn×nj, and Υkij∈Rn×ni+j are constant matrices.Notice thatg(x(t,τ)) is composed of three terms, the first is a function of x[i](t), the second depends on x[j](t-τ), and the third is a function of x[i](t)⊗x[j](t-τ). Consider the first term of g(x(t,τ)); we note it g(x(t)). It could be written as follows:(32)gxt=∑i=0rG1ixit⋮G2ixit⋮⋯⋮Gmixit=∑i=0rG1i⋮G2i⋮⋯⋮Gmi︸Gixit0xit⋱0xit︸Im⊗xitTheng(x(t,τ)) could be generalized to the following expression:(33)gxt,τ=∑i=0rGiIm⊗xit+∑j=1sG¯jIm⊗xjt-τ+∑i=1r∑j=1sΥijIm⊗xit⊗xjt-τ


### 3.2. Statement of the Problem

Consider the system defined by (28), (29), and (33) with an initial condition x(0)=x0. Our objective is firstly to find the optimal open loop control u∗(t), which minimizes the performance index:(34)J=12xTtfSxtf+12∫0tfxTQx+uTRudtwhereQ and S are positive semidefinite matrices and R is symmetric positive definite with appropriate dimensions.The direct approach presented in this paper is based on expanding system equations (28), (29), and (33) as well as objective function (34) to be minimized over an hybrid functions basis. Hence, the main purpose is to transform the optimal control problem under dynamic constraints to a nonlinear programming problem. To this end, each of the state and control variables is approximated by a finite length of unknown parameters as follows:(35)xt≈XThtut≈UThtwhereX and U are unknown state and control parameters, respectively. Applying the vec operator (see Appendix) and related Kronecker product property [32] yields(36)vecxt≈hTt⊗InvecXTvecut≈hTt⊗ImvecUTwhereIn and Im are n×n and m×m identity matrices.Moreover, at the initial time,t0=0, the initial state could be written(37)x0≅X0Thtwhere(38)X0T=x00⋯0︸Mx00⋯0︸M⋯x00⋯0︸Mis anw constant vector.For clarity purpose, let us denotez as the whole unknown parameters vector. zx=vec(XT) and zu=vec(UT) are, respectively, the state parameters and the control ones, such that(39)z=zxzuandzx0=vec(X0T).According to (14) and (17), the delayed state coefficients are given by(40)vecX∗Tτ=Ehzx0+Dhzx


### 3.3. Optimal Control Problem Reformulation Using Hybrid Functions

The cost function (33) is composed of two parts. The first is the terminal penalty of the state, while the second is known to be the running cost.


#### 3.3.1. Cost of the Final State Approximation

At the final time,tf, the state approximation could be written(41)xtf≈XThtfIt is important to mention here that hybrid functions inherit an important property from Legendre polynomials(Li(tf)=1,∀i=0,…,M-1). In fact, the subset hNj(t) verifies(42)hNjtf=1,∀j=0,1,⋯M-1The rest of hybrid functions are null att=tf.The cross product of two hybrid functions at the final time is given by(43)htfhTtf=Qf=0M0M⋯0M0M⋱⋱⋮⋮⋱0M0M0M⋯0M1MwhereQf is NM×NM matrix. 1M stands for the all-ones M×M matrix.Hence, the terminal penalty of the state could be approximated as follows:(44)xTtfSxtf≈zxThtf⊗InShTtf⊗Inzx≈zxTQf⊗Szx


#### 3.3.2. Cost of the State Trajectory Approximation

The integral term∫0tf‍xT(t)Qx(t)+uT(t)Ru(t)dt in the criterion is approached now as(45)∫0tfzxTht⊗InQhTt⊗Inzx+zuTht⊗ImRhTt⊗Imzudtwhich is equivalent to(46)∫0tfzxThthTt⊗Qzx+zuThthTt⊗RzudtUsing the integral of the cross operational matrixC, it reduces to(47)zxTC⊗Qzx+zuTC⊗Rzu


#### 3.3.3. System Path Approximation

The expansion of the system state over a hybrid basis requires the development of functionsf(x(t,τ)) and g(x(t,τ)) over that basis. To this end, several preliminary lemmas need to be introduced.Lemma 1.
The development of theith Kronecker power of the state vector over a hybrid basis h(t) gives (48)xit=XiT.ht
where(49)XiT=XTi.κi
with(50)κi=κi-1⊗Iw.Kh,fori=3,4,…and(51)κ2=Kh,κ1=κ0=1
We recall thatXT[i] denotes the ith Kronecker power of the state coefficients XT, with Kh being the operational matrix of the Kronecker product.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma1 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., x[j](t-τ)=X[j]∗T.h(t)) and express it in terms of decision variable zx by the mean of relation (32):(52)Xj∗T=X∗Tj.κjExpansion of f ( x ( t , τ ) ) over the Hybrid Basis. The third term of (29) could be developed over the hybrid basis as follows:(53)xit⊗xjt-τ≃XTi⊗X∗Tjκi,jKh.ht,∀i=1,2,…where we notateκ[i,j]=κ[i]⊗κ[j].Now,f(x(t,τ)) could be approached as follows:(54)fxt,τ≃FT.htwhere(55)FT=∑i=1pFiXTiκi+∑j=1qF¯jX∗Tjκj+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhExpansion of g ( x ( t , τ ) ) u ( t ) over the Hybrid Basis. Notice that the first term of g(x(t,τ))u(t) under the sum could be expanded over the hybrid basis as follows:(56)Im⊗xitut≃UT⊗XTiκiKhht≃UT⊗XTiIw⊗κiKhht,∀i=1,2,…while the second one could be derived similarly.The third term could be approached as(57)Im⊗xit⊗xjt-τut≃UT⊗XTi⊗X∗Tjκi,jKh.Khht≃UT⊗XTi⊗X∗TjIw⊗κi,jKhKhhtwherei and j belong to N∗. Then it comes(58)gxt,τut≃GT.htwith(59)GT=G0UT+∑i=1rGiUT⊗XTiIw⊗κiKh+∑j=1sG¯jUT⊗X∗TjIw⊗κjKh+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhExpansion of System Equation over the Hybrid Basis. The integration of the system equation by introducing the operational matrix of integration P with respect to notations (58) and (54) gives(60)XT-X0T=FTP+GTPOur objective is to express the constraint (60) in terms of decision variables zx and zu; to this end we apply the vec operator to (60). That allows us to state(61)zx-zx0=vec∑i=1pFiXTiκiP+∑j=1qF¯jX∗TjκjP+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rGiUT⊗XTiκi.KhP+∑j=1sG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPUsing the linearity property of thevec operator it comes(62)zx-zx0=∑i=1pvecFiXTiκiP+∑j=1qvecF¯jX∗TjκjP+∑i=1p∑j=1qvecΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rvecGiUT⊗XTiIw⊗κiKhP+∑j=1svecG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1svecΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPwhich is equivalent to(63)zx-zx0=∑i=1pPTκiT⊗FivecXTi+∑j=1qPTκjT⊗Fj¯vecX∗Tj+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijvecXTi⊗X∗Tj+G0UTP+∑i=1rPTKhTIw⊗κiT⊗GivecUT⊗XTi+∑j=1sPTKhTIw⊗κiT⊗G¯jvecUT⊗X∗Tj+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijvecUT⊗XTi⊗X∗TjLemma 2.
(64) v e c X T i = Π n i - 1 , w i - 1 m a t z x . Π n i - 2 , w i - 2 m a t z x … Π n , w m a t z x . z x = Δ i z x ,  f o r   i = 2,3 , …
where the matrixΠ(.,.) is defined in Appendix.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma2 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., X∗T[j]) and express it in terms of decision variable zx by the mean of relation (42). We note(65)vecX∗Ti=Δi∗zxApplying thevec operator to (XT[i]⊗X∗T[j]) yields(66)vecXTi⊗X∗Tj=Πnj,wjmatvecXTi.vecX∗Tj=Πnj,wjmatΔizxΔj∗zxSimilarly(67)vecUT⊗XTi=Πni,wimatzu.vecXTi=Πni,wimatzuΔizxFinally, the system path constraint could be implemented using the following equation:(68)zx-zx0=∑i=1pPTκiT⊗FiΔizx+∑j=1qPTκjT⊗F¯jΔi∗zx+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijΠni,wjmatΔizxΔj∗zx+PT⊗G0zu+∑i=1rPTKhTIw⊗κiT⊗GiΠni,wimatzuΔizx+∑j=1sPTKhTIw⊗κiT⊗G¯jΠnj,wjmatzuΔj∗zx+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijΠni+j,wi+jmatzuΠni,wjmatΔizxΔj∗zxNow, it could be noticed that the system path constraint is expressed properly in terms of unknown parameterszx and zu.


### 3.4. The Nonlinear Programming Problem

The optimal control problem has been approximated by a nonlinear programming problem and is given by the following: find the optimal vectorz of the unknown parameters zx and zu that minimizes(69)12zTΩzsubject to (68).One has(70)Ω=Qf+Q0nNM×mNM0mNM×nNMRThe mathematical programming problem can be solved by using available nonlinear programming solvers like IPOPT or the routinefmincon of the MATLAB Toolbox.After solving the latter nonlinear programming problem and determining the optimal value of the unknown parameters vectorz, these parameters are substituted back into (28) to determine the optimal state vector and the optimal control.


## 3.1. Description of the Studied System

We consider the nonlinear continuous system which can be represented by the following state space representation:(28)x.t=fxt,τ+gxt,τutxt-τ=x0,∀t∈0,τwherex(t)∈Rn is the state vector, x(t-τ)∈Rn is the delayed state where τ denotes the time delay, u(t)∈Rm is the control vector, and f(x(t,τ)) from Rn into Rn and g(x(t,τ)) from Rn into Rn×m are nonlinear analytic functions of x(t) and x(t-τ).Note that, any functionsf(x(t,τ)) and g(x(t,τ)) could be approached using truncated series of the Kronecker power of x(t) and x(t-τ) as follows:(29)fxt,τ≈∑i=1pFixit+∑j=1qF¯jxjt-τ+∑i=1p∑j=1qΓijxit⊗xjt-τwhereFi∈Rn×nj, F¯j∈Rn×nj and Γij∈Rn×ni+j are constant matrices. x[i](t) denotes the ith Kronecker power of the state vector (see Appendix).And(30)gxt,τ=g1xt,τ⋮g2xt,τ⋮⋯⋮gmxt,τwheregk(x(t,τ))∈Rn, for k=1,…,m, are defined by(31)gkxt,τ≈∑i=0rGkixit+∑j=1sG¯kjxjt-τ+∑i=1r∑j=1sΥkijxit⊗xjt-τwhereGki∈Rn×nj, G¯kj∈Rn×nj, and Υkij∈Rn×ni+j are constant matrices.Notice thatg(x(t,τ)) is composed of three terms, the first is a function of x[i](t), the second depends on x[j](t-τ), and the third is a function of x[i](t)⊗x[j](t-τ). Consider the first term of g(x(t,τ)); we note it g(x(t)). It could be written as follows:(32)gxt=∑i=0rG1ixit⋮G2ixit⋮⋯⋮Gmixit=∑i=0rG1i⋮G2i⋮⋯⋮Gmi︸Gixit0xit⋱0xit︸Im⊗xitTheng(x(t,τ)) could be generalized to the following expression:(33)gxt,τ=∑i=0rGiIm⊗xit+∑j=1sG¯jIm⊗xjt-τ+∑i=1r∑j=1sΥijIm⊗xit⊗xjt-τ


## 3.2. Statement of the Problem

Consider the system defined by (28), (29), and (33) with an initial condition x(0)=x0. Our objective is firstly to find the optimal open loop control u∗(t), which minimizes the performance index:(34)J=12xTtfSxtf+12∫0tfxTQx+uTRudtwhereQ and S are positive semidefinite matrices and R is symmetric positive definite with appropriate dimensions.The direct approach presented in this paper is based on expanding system equations (28), (29), and (33) as well as objective function (34) to be minimized over an hybrid functions basis. Hence, the main purpose is to transform the optimal control problem under dynamic constraints to a nonlinear programming problem. To this end, each of the state and control variables is approximated by a finite length of unknown parameters as follows:(35)xt≈XThtut≈UThtwhereX and U are unknown state and control parameters, respectively. Applying the vec operator (see Appendix) and related Kronecker product property [32] yields(36)vecxt≈hTt⊗InvecXTvecut≈hTt⊗ImvecUTwhereIn and Im are n×n and m×m identity matrices.Moreover, at the initial time,t0=0, the initial state could be written(37)x0≅X0Thtwhere(38)X0T=x00⋯0︸Mx00⋯0︸M⋯x00⋯0︸Mis anw constant vector.For clarity purpose, let us denotez as the whole unknown parameters vector. zx=vec(XT) and zu=vec(UT) are, respectively, the state parameters and the control ones, such that(39)z=zxzuandzx0=vec(X0T).According to (14) and (17), the delayed state coefficients are given by(40)vecX∗Tτ=Ehzx0+Dhzx


## 3.3. Optimal Control Problem Reformulation Using Hybrid Functions

The cost function (33) is composed of two parts. The first is the terminal penalty of the state, while the second is known to be the running cost.


### 3.3.1. Cost of the Final State Approximation

At the final time,tf, the state approximation could be written(41)xtf≈XThtfIt is important to mention here that hybrid functions inherit an important property from Legendre polynomials(Li(tf)=1,∀i=0,…,M-1). In fact, the subset hNj(t) verifies(42)hNjtf=1,∀j=0,1,⋯M-1The rest of hybrid functions are null att=tf.The cross product of two hybrid functions at the final time is given by(43)htfhTtf=Qf=0M0M⋯0M0M⋱⋱⋮⋮⋱0M0M0M⋯0M1MwhereQf is NM×NM matrix. 1M stands for the all-ones M×M matrix.Hence, the terminal penalty of the state could be approximated as follows:(44)xTtfSxtf≈zxThtf⊗InShTtf⊗Inzx≈zxTQf⊗Szx


### 3.3.2. Cost of the State Trajectory Approximation

The integral term∫0tf‍xT(t)Qx(t)+uT(t)Ru(t)dt in the criterion is approached now as(45)∫0tfzxTht⊗InQhTt⊗Inzx+zuTht⊗ImRhTt⊗Imzudtwhich is equivalent to(46)∫0tfzxThthTt⊗Qzx+zuThthTt⊗RzudtUsing the integral of the cross operational matrixC, it reduces to(47)zxTC⊗Qzx+zuTC⊗Rzu


### 3.3.3. System Path Approximation

The expansion of the system state over a hybrid basis requires the development of functionsf(x(t,τ)) and g(x(t,τ)) over that basis. To this end, several preliminary lemmas need to be introduced.Lemma 1.
The development of theith Kronecker power of the state vector over a hybrid basis h(t) gives (48)xit=XiT.ht
where(49)XiT=XTi.κi
with(50)κi=κi-1⊗Iw.Kh,fori=3,4,…and(51)κ2=Kh,κ1=κ0=1
We recall thatXT[i] denotes the ith Kronecker power of the state coefficients XT, with Kh being the operational matrix of the Kronecker product.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma1 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., x[j](t-τ)=X[j]∗T.h(t)) and express it in terms of decision variable zx by the mean of relation (32):(52)Xj∗T=X∗Tj.κjExpansion of f ( x ( t , τ ) ) over the Hybrid Basis. The third term of (29) could be developed over the hybrid basis as follows:(53)xit⊗xjt-τ≃XTi⊗X∗Tjκi,jKh.ht,∀i=1,2,…where we notateκ[i,j]=κ[i]⊗κ[j].Now,f(x(t,τ)) could be approached as follows:(54)fxt,τ≃FT.htwhere(55)FT=∑i=1pFiXTiκi+∑j=1qF¯jX∗Tjκj+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhExpansion of g ( x ( t , τ ) ) u ( t ) over the Hybrid Basis. Notice that the first term of g(x(t,τ))u(t) under the sum could be expanded over the hybrid basis as follows:(56)Im⊗xitut≃UT⊗XTiκiKhht≃UT⊗XTiIw⊗κiKhht,∀i=1,2,…while the second one could be derived similarly.The third term could be approached as(57)Im⊗xit⊗xjt-τut≃UT⊗XTi⊗X∗Tjκi,jKh.Khht≃UT⊗XTi⊗X∗TjIw⊗κi,jKhKhhtwherei and j belong to N∗. Then it comes(58)gxt,τut≃GT.htwith(59)GT=G0UT+∑i=1rGiUT⊗XTiIw⊗κiKh+∑j=1sG¯jUT⊗X∗TjIw⊗κjKh+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhExpansion of System Equation over the Hybrid Basis. The integration of the system equation by introducing the operational matrix of integration P with respect to notations (58) and (54) gives(60)XT-X0T=FTP+GTPOur objective is to express the constraint (60) in terms of decision variables zx and zu; to this end we apply the vec operator to (60). That allows us to state(61)zx-zx0=vec∑i=1pFiXTiκiP+∑j=1qF¯jX∗TjκjP+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rGiUT⊗XTiκi.KhP+∑j=1sG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPUsing the linearity property of thevec operator it comes(62)zx-zx0=∑i=1pvecFiXTiκiP+∑j=1qvecF¯jX∗TjκjP+∑i=1p∑j=1qvecΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rvecGiUT⊗XTiIw⊗κiKhP+∑j=1svecG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1svecΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPwhich is equivalent to(63)zx-zx0=∑i=1pPTκiT⊗FivecXTi+∑j=1qPTκjT⊗Fj¯vecX∗Tj+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijvecXTi⊗X∗Tj+G0UTP+∑i=1rPTKhTIw⊗κiT⊗GivecUT⊗XTi+∑j=1sPTKhTIw⊗κiT⊗G¯jvecUT⊗X∗Tj+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijvecUT⊗XTi⊗X∗TjLemma 2.
(64) v e c X T i = Π n i - 1 , w i - 1 m a t z x . Π n i - 2 , w i - 2 m a t z x … Π n , w m a t z x . z x = Δ i z x ,  f o r   i = 2,3 , …
where the matrixΠ(.,.) is defined in Appendix.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma2 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., X∗T[j]) and express it in terms of decision variable zx by the mean of relation (42). We note(65)vecX∗Ti=Δi∗zxApplying thevec operator to (XT[i]⊗X∗T[j]) yields(66)vecXTi⊗X∗Tj=Πnj,wjmatvecXTi.vecX∗Tj=Πnj,wjmatΔizxΔj∗zxSimilarly(67)vecUT⊗XTi=Πni,wimatzu.vecXTi=Πni,wimatzuΔizxFinally, the system path constraint could be implemented using the following equation:(68)zx-zx0=∑i=1pPTκiT⊗FiΔizx+∑j=1qPTκjT⊗F¯jΔi∗zx+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijΠni,wjmatΔizxΔj∗zx+PT⊗G0zu+∑i=1rPTKhTIw⊗κiT⊗GiΠni,wimatzuΔizx+∑j=1sPTKhTIw⊗κiT⊗G¯jΠnj,wjmatzuΔj∗zx+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijΠni+j,wi+jmatzuΠni,wjmatΔizxΔj∗zxNow, it could be noticed that the system path constraint is expressed properly in terms of unknown parameterszx and zu.


## 3.3.1. Cost of the Final State Approximation

At the final time,tf, the state approximation could be written(41)xtf≈XThtfIt is important to mention here that hybrid functions inherit an important property from Legendre polynomials(Li(tf)=1,∀i=0,…,M-1). In fact, the subset hNj(t) verifies(42)hNjtf=1,∀j=0,1,⋯M-1The rest of hybrid functions are null att=tf.The cross product of two hybrid functions at the final time is given by(43)htfhTtf=Qf=0M0M⋯0M0M⋱⋱⋮⋮⋱0M0M0M⋯0M1MwhereQf is NM×NM matrix. 1M stands for the all-ones M×M matrix.Hence, the terminal penalty of the state could be approximated as follows:(44)xTtfSxtf≈zxThtf⊗InShTtf⊗Inzx≈zxTQf⊗Szx


## 3.3.2. Cost of the State Trajectory Approximation

The integral term∫0tf‍xT(t)Qx(t)+uT(t)Ru(t)dt in the criterion is approached now as(45)∫0tfzxTht⊗InQhTt⊗Inzx+zuTht⊗ImRhTt⊗Imzudtwhich is equivalent to(46)∫0tfzxThthTt⊗Qzx+zuThthTt⊗RzudtUsing the integral of the cross operational matrixC, it reduces to(47)zxTC⊗Qzx+zuTC⊗Rzu


## 3.3.3. System Path Approximation

The expansion of the system state over a hybrid basis requires the development of functionsf(x(t,τ)) and g(x(t,τ)) over that basis. To this end, several preliminary lemmas need to be introduced.Lemma 1.
The development of theith Kronecker power of the state vector over a hybrid basis h(t) gives (48)xit=XiT.ht
where(49)XiT=XTi.κi
with(50)κi=κi-1⊗Iw.Kh,fori=3,4,…and(51)κ2=Kh,κ1=κ0=1
We recall thatXT[i] denotes the ith Kronecker power of the state coefficients XT, with Kh being the operational matrix of the Kronecker product.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma1 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., x[j](t-τ)=X[j]∗T.h(t)) and express it in terms of decision variable zx by the mean of relation (32):(52)Xj∗T=X∗Tj.κjExpansion of f ( x ( t , τ ) ) over the Hybrid Basis. The third term of (29) could be developed over the hybrid basis as follows:(53)xit⊗xjt-τ≃XTi⊗X∗Tjκi,jKh.ht,∀i=1,2,…where we notateκ[i,j]=κ[i]⊗κ[j].Now,f(x(t,τ)) could be approached as follows:(54)fxt,τ≃FT.htwhere(55)FT=∑i=1pFiXTiκi+∑j=1qF¯jX∗Tjκj+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhExpansion of g ( x ( t , τ ) ) u ( t ) over the Hybrid Basis. Notice that the first term of g(x(t,τ))u(t) under the sum could be expanded over the hybrid basis as follows:(56)Im⊗xitut≃UT⊗XTiκiKhht≃UT⊗XTiIw⊗κiKhht,∀i=1,2,…while the second one could be derived similarly.The third term could be approached as(57)Im⊗xit⊗xjt-τut≃UT⊗XTi⊗X∗Tjκi,jKh.Khht≃UT⊗XTi⊗X∗TjIw⊗κi,jKhKhhtwherei and j belong to N∗. Then it comes(58)gxt,τut≃GT.htwith(59)GT=G0UT+∑i=1rGiUT⊗XTiIw⊗κiKh+∑j=1sG¯jUT⊗X∗TjIw⊗κjKh+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhExpansion of System Equation over the Hybrid Basis. The integration of the system equation by introducing the operational matrix of integration P with respect to notations (58) and (54) gives(60)XT-X0T=FTP+GTPOur objective is to express the constraint (60) in terms of decision variables zx and zu; to this end we apply the vec operator to (60). That allows us to state(61)zx-zx0=vec∑i=1pFiXTiκiP+∑j=1qF¯jX∗TjκjP+∑i=1p∑j=1qΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rGiUT⊗XTiκi.KhP+∑j=1sG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1sΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPUsing the linearity property of thevec operator it comes(62)zx-zx0=∑i=1pvecFiXTiκiP+∑j=1qvecF¯jX∗TjκjP+∑i=1p∑j=1qvecΓijXTi⊗X∗Tjκi,jKhP+G0UTP+∑i=1rvecGiUT⊗XTiIw⊗κiKhP+∑j=1svecG¯jUT⊗X∗TjIw⊗κjKhP+∑i=1r∑j=1svecΥijUT⊗XTi⊗X∗TjIw⊗κi,jKhKhPwhich is equivalent to(63)zx-zx0=∑i=1pPTκiT⊗FivecXTi+∑j=1qPTκjT⊗Fj¯vecX∗Tj+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijvecXTi⊗X∗Tj+G0UTP+∑i=1rPTKhTIw⊗κiT⊗GivecUT⊗XTi+∑j=1sPTKhTIw⊗κiT⊗G¯jvecUT⊗X∗Tj+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijvecUT⊗XTi⊗X∗TjLemma 2.
(64) v e c X T i = Π n i - 1 , w i - 1 m a t z x . Π n i - 2 , w i - 2 m a t z x … Π n , w m a t z x . z x = Δ i z x ,  f o r   i = 2,3 , …
where the matrixΠ(.,.) is defined in Appendix.Proof.
The proof of this lemma needs only a few manipulations.Notice that results of Lemma2 could be applied to the jth Kronecker power of the delayed state coefficients (i.e., X∗T[j]) and express it in terms of decision variable zx by the mean of relation (42). We note(65)vecX∗Ti=Δi∗zxApplying thevec operator to (XT[i]⊗X∗T[j]) yields(66)vecXTi⊗X∗Tj=Πnj,wjmatvecXTi.vecX∗Tj=Πnj,wjmatΔizxΔj∗zxSimilarly(67)vecUT⊗XTi=Πni,wimatzu.vecXTi=Πni,wimatzuΔizxFinally, the system path constraint could be implemented using the following equation:(68)zx-zx0=∑i=1pPTκiT⊗FiΔizx+∑j=1qPTκjT⊗F¯jΔi∗zx+∑i=1p∑j=1qPTKhTκi,jT⊗ΓijΠni,wjmatΔizxΔj∗zx+PT⊗G0zu+∑i=1rPTKhTIw⊗κiT⊗GiΠni,wimatzuΔizx+∑j=1sPTKhTIw⊗κiT⊗G¯jΠnj,wjmatzuΔj∗zx+∑i=1r∑j=1sPTKhTIw⊗KhTκi,jT⊗ΥijΠni+j,wi+jmatzuΠni,wjmatΔizxΔj∗zxNow, it could be noticed that the system path constraint is expressed properly in terms of unknown parameterszx and zu.


## 3.4. The Nonlinear Programming Problem

The optimal control problem has been approximated by a nonlinear programming problem and is given by the following: find the optimal vectorz of the unknown parameters zx and zu that minimizes(69)12zTΩzsubject to (68).One has(70)Ω=Qf+Q0nNM×mNM0mNM×nNMRThe mathematical programming problem can be solved by using available nonlinear programming solvers like IPOPT or the routinefmincon of the MATLAB Toolbox.After solving the latter nonlinear programming problem and determining the optimal value of the unknown parameters vectorz, these parameters are substituted back into (28) to determine the optimal state vector and the optimal control.


## 4. Suboptimal Feedback Control

Once the optimal open loop results are obtained by solving the nonlinear programming problem given by (69)-(68), let us note(71)z∗=zx∗zu∗,the optimal state and control coefficients.We are interested now, based on previous results, to synthesize the following nonlinear state feedback control law:(72)ut=-∑i=1lKixit-∑j=1vK¯jxjt-τThe idea is to find control matricesKi and K¯j such that the optimal vector (71) verifies the control equation (72).Expanding (72) over the hybrid basis yields(73)UT=-∑i=1lKiXTiκi-∑j=1vK¯jX∗TjκjSubstituting the control and state coefficients with their optimal values and applying thevec operator give(74)zu∗=-∑i=1lαivecKi-∑j=1vβjvecK¯jwithαi=(κ[i]T(mat(zx∗))[i]T⊗Im) and βj=(κ[j]T(mat(Ehzx0+Dhzx∗))[j]T⊗Im).Finding control parameters could be then reduced to solving, in the least square sense, the following problem:(75)A.θ=Bwhere(76)A=α1⋮α2⋮⋯⋮αl⋮β1⋮β2⋮⋯⋮βv,B=-zu∗,θ=vecK1vecK2⋮vecKlvecK¯1vecK¯2⋮vecK¯v


## 5. Computational Results

### 5.1. Example 1

Consider the system [33](77)x˙1t=x1tx2t+2x1t-τx2t-x1t-τ+utx˙2t=-x1t+x2t-τWhenu=0, the above system has two equilibria, one which is at the origin and the other one at (1/3,1/3).In this example, we aim to minimize the following criterion∫010‍(x12(t)+x22(t)+u2(t))dt in order to find optimal states and open loop control, x∗(t) and u∗(t), then a suboptimal control u(t)=-Kx(t)-K¯x(t-τ) is characterized.Nonlinear system (77) could be written under a polynomial form (28), (29), and (32) with(78)F1=00-10,F¯1=-100-1,G0=10,F2=01000000,Γ11=00200000The time delay is considered asτ=2.


#### 5.1.1. Open Loop Study

The development presented above is implemented in this subsection by using both hybrid functions (HFs) and block pulse functions (BPFs).Table1 summarizes considered parameters for simulations below and obtained performances indexes with the different bases. It is then clear that hybrid basis is superior over the piecewise constant one, both with the same number of elementary functions.Table 1
Hybrid functions and block pulse functions direct approach performance analysis for Example 1.
Method Parameters J Hybrid Functions N=5, M=5 1.9955 Block Pulse Functions N=25 2.1679Simulation results for the above open loop controlled system initialized with(0.5,0,5) are given in Figures 1 and 2 based on hybrid and block pulse functions.Figure 1
Optimal states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure 2
Optimal control signals with hybrid functions (HFs) and block pulse functions (BPFs).


#### 5.1.2. Closed Loop Framework

The state feedback control gains designed based on hybrid and block pulse open loop frameworks are as follows:(79)KHFs=3.0813-1.3873,K¯HFs=-1.1017-1.7575KBPFs=2.9425-0.6949K¯BPFs=-0.9131-2.1364Controlled state trajectories, obtained with determined gains, are depicted in Figure3. It is shown that the hybrid functions technique is also better in closed loop.Figure 3
Closed loop states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure4 shows the optimal states trajectories obtained by minimizing the formulated NLP problem, by using the hybrid of block pulse and Legendre polynomials basis, over a finite horizon tf=10. Controlled states with obtained suboptimal feedback are drawn on the same figure over a simulation time 20s. It could be seen that system states converge to the origin equilibrium with respect to imposed criterion.Figure 4
Optimal and suboptimal states trajectories.Figure5 exposes optimal control and suboptimal state feedback control signals using HFs.Figure 5
Optimal and suboptimal control signals.


### 5.2. Example 2: A Two-Stage Chemical Reactor

In this section, we consider a cascade chemical system with two reactors [34](80)x˙1t=-k1x1t-1θ1x1t-1θ1x1t-d+1-R2V1x2t+δ1x1t-dx˙2t=-k2x2t-1θ2x22t+R1V2x1t-d-1θ2x2t+R2V2x2t-d+FV2ut+δ2x2t-dwherexi,i=1,2, are the compositions, d is a known time delay, Ri are the recycle flow rates, θi are the reactor residence times, ki are the reaction constants, F is the feed rate, Vi are reactor volumes, and δi are nonlinear functions for describing the system uncertainties and external disturbances. Note that (80) is the transformed reactor model given by [34]. The original one may have a nonzero equilibrium point and the compositions xi actually denote the deviations from the equilibrium point.LetRi=0.5, Vi=0.5, d=0.5, ki=0.5, and F=0.5. The uncertainties δi are the functions δ1(x1(t-d))=θ3x1(t-d) and δ2(x2(t-d))=0.5θ4x22(t-d), respectively, with θi=1.The simulation is done by taking(x1(t0),x2(t0))=(1,-2) for t0∈[-0.5,0] and results are given in Figures 3 and 4.We consider in this example the criterion∫03‍(100x12(t)+100x22(t)+u2(t))dt in order to find optimal states and open loop control, x∗(t) and u∗(t), then suboptimal control gains such that u(t)=-K1x(t)-K2x[2](t)-K¯1x(t-τ)-K¯2x[2](t-τ) are investigated.Nonlinear system (80) could be written under a polynomial form (28), (29), and (32) with(81)
F
1
=
-
k
1
-
1
θ
1
1
-
R
2
V
1
0
-
k
2
-
1
θ
2
,
F
¯
1
=
-
1
θ
1
+
θ
3
0
R
1
V
2
R
2
V
2
,
G
0
=
0
F
V
2
,
F
2
=
0
0
0
0
0
0
0
-
1
θ
2
,
F
¯
2
=
0
0
2
0
0
0
0
0.5
1
θ
2Table2 gives simulations parameters and obtained Performances indexes with the two bases utilized in Example 1. Also, state feedback gains are included.Table 2
Hybrid functions and block pulse functions open loop and closed loop results for Example 2.
Method Parameters J State feedback gains Hybrid Functions N=6, M=5 57.9349 K 1 = [ 47.8366     10.7607 ], K 2 = [ - 65.4651     - 26.375     0     - 6.0281 ], K ¯ 1 =   - 17.7816     0.5511 ], K ¯ 2 = [ 9.6281     0     - 2.1138     - 0.7305 ] Block Pulse Functions N=30 69.2999 K 1 = [ 90.9235     - 1.9396 ], K 2 = [ - 131.5504     0     - 30.3858     - 8.1110 ], K ¯ 1 =   - 36.3910     - 9.3258 ], K ¯ 2 = [ 22.8197     14.7181     0     2.3220 ]It is worth noting that for this particular system, control gains determined by the mean of block pulse functions are not stabilizing. While hybrid functions results are illustrated on Figures6 and 7, controlled states with obtained suboptimal feedback designed over a finite horizon tf=3s are drawn on the same figure over a simulation time 5s. It could be seen that suboptimal system states coincide perfectly with optimal solution.Figure 6
Optimal and suboptimal states trajectories of the chemical reactor.Figure 7
Optimal and suboptimal control signals for the chemical reactor.Figure7 exposes optimal control and suggested suboptimal nonlinear state feedback control signals. The proposed nonlinear feedback reproduces sharply the optimal open loop control.


## 5.1. Example 1

Consider the system [33](77)x˙1t=x1tx2t+2x1t-τx2t-x1t-τ+utx˙2t=-x1t+x2t-τWhenu=0, the above system has two equilibria, one which is at the origin and the other one at (1/3,1/3).In this example, we aim to minimize the following criterion∫010‍(x12(t)+x22(t)+u2(t))dt in order to find optimal states and open loop control, x∗(t) and u∗(t), then a suboptimal control u(t)=-Kx(t)-K¯x(t-τ) is characterized.Nonlinear system (77) could be written under a polynomial form (28), (29), and (32) with(78)F1=00-10,F¯1=-100-1,G0=10,F2=01000000,Γ11=00200000The time delay is considered asτ=2.


### 5.1.1. Open Loop Study

The development presented above is implemented in this subsection by using both hybrid functions (HFs) and block pulse functions (BPFs).Table1 summarizes considered parameters for simulations below and obtained performances indexes with the different bases. It is then clear that hybrid basis is superior over the piecewise constant one, both with the same number of elementary functions.Table 1
Hybrid functions and block pulse functions direct approach performance analysis for Example 1.
Method Parameters J Hybrid Functions N=5, M=5 1.9955 Block Pulse Functions N=25 2.1679Simulation results for the above open loop controlled system initialized with(0.5,0,5) are given in Figures 1 and 2 based on hybrid and block pulse functions.Figure 1
Optimal states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure 2
Optimal control signals with hybrid functions (HFs) and block pulse functions (BPFs).


### 5.1.2. Closed Loop Framework

The state feedback control gains designed based on hybrid and block pulse open loop frameworks are as follows:(79)KHFs=3.0813-1.3873,K¯HFs=-1.1017-1.7575KBPFs=2.9425-0.6949K¯BPFs=-0.9131-2.1364Controlled state trajectories, obtained with determined gains, are depicted in Figure3. It is shown that the hybrid functions technique is also better in closed loop.Figure 3
Closed loop states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure4 shows the optimal states trajectories obtained by minimizing the formulated NLP problem, by using the hybrid of block pulse and Legendre polynomials basis, over a finite horizon tf=10. Controlled states with obtained suboptimal feedback are drawn on the same figure over a simulation time 20s. It could be seen that system states converge to the origin equilibrium with respect to imposed criterion.Figure 4
Optimal and suboptimal states trajectories.Figure5 exposes optimal control and suboptimal state feedback control signals using HFs.Figure 5
Optimal and suboptimal control signals.


## 5.1.1. Open Loop Study

The development presented above is implemented in this subsection by using both hybrid functions (HFs) and block pulse functions (BPFs).Table1 summarizes considered parameters for simulations below and obtained performances indexes with the different bases. It is then clear that hybrid basis is superior over the piecewise constant one, both with the same number of elementary functions.Table 1
Hybrid functions and block pulse functions direct approach performance analysis for Example 1.
Method Parameters J Hybrid Functions N=5, M=5 1.9955 Block Pulse Functions N=25 2.1679Simulation results for the above open loop controlled system initialized with(0.5,0,5) are given in Figures 1 and 2 based on hybrid and block pulse functions.Figure 1
Optimal states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure 2
Optimal control signals with hybrid functions (HFs) and block pulse functions (BPFs).


## 5.1.2. Closed Loop Framework

The state feedback control gains designed based on hybrid and block pulse open loop frameworks are as follows:(79)KHFs=3.0813-1.3873,K¯HFs=-1.1017-1.7575KBPFs=2.9425-0.6949K¯BPFs=-0.9131-2.1364Controlled state trajectories, obtained with determined gains, are depicted in Figure3. It is shown that the hybrid functions technique is also better in closed loop.Figure 3
Closed loop states trajectories with hybrid functions (HFs) and block pulse functions (BPFs).Figure4 shows the optimal states trajectories obtained by minimizing the formulated NLP problem, by using the hybrid of block pulse and Legendre polynomials basis, over a finite horizon tf=10. Controlled states with obtained suboptimal feedback are drawn on the same figure over a simulation time 20s. It could be seen that system states converge to the origin equilibrium with respect to imposed criterion.Figure 4
Optimal and suboptimal states trajectories.Figure5 exposes optimal control and suboptimal state feedback control signals using HFs.Figure 5
Optimal and suboptimal control signals.


## 5.2. Example 2: A Two-Stage Chemical Reactor

In this section, we consider a cascade chemical system with two reactors [34](80)x˙1t=-k1x1t-1θ1x1t-1θ1x1t-d+1-R2V1x2t+δ1x1t-dx˙2t=-k2x2t-1θ2x22t+R1V2x1t-d-1θ2x2t+R2V2x2t-d+FV2ut+δ2x2t-dwherexi,i=1,2, are the compositions, d is a known time delay, Ri are the recycle flow rates, θi are the reactor residence times, ki are the reaction constants, F is the feed rate, Vi are reactor volumes, and δi are nonlinear functions for describing the system uncertainties and external disturbances. Note that (80) is the transformed reactor model given by [34]. The original one may have a nonzero equilibrium point and the compositions xi actually denote the deviations from the equilibrium point.LetRi=0.5, Vi=0.5, d=0.5, ki=0.5, and F=0.5. The uncertainties δi are the functions δ1(x1(t-d))=θ3x1(t-d) and δ2(x2(t-d))=0.5θ4x22(t-d), respectively, with θi=1.The simulation is done by taking(x1(t0),x2(t0))=(1,-2) for t0∈[-0.5,0] and results are given in Figures 3 and 4.We consider in this example the criterion∫03‍(100x12(t)+100x22(t)+u2(t))dt in order to find optimal states and open loop control, x∗(t) and u∗(t), then suboptimal control gains such that u(t)=-K1x(t)-K2x[2](t)-K¯1x(t-τ)-K¯2x[2](t-τ) are investigated.Nonlinear system (80) could be written under a polynomial form (28), (29), and (32) with(81)
F
1
=
-
k
1
-
1
θ
1
1
-
R
2
V
1
0
-
k
2
-
1
θ
2
,
F
¯
1
=
-
1
θ
1
+
θ
3
0
R
1
V
2
R
2
V
2
,
G
0
=
0
F
V
2
,
F
2
=
0
0
0
0
0
0
0
-
1
θ
2
,
F
¯
2
=
0
0
2
0
0
0
0
0.5
1
θ
2Table2 gives simulations parameters and obtained Performances indexes with the two bases utilized in Example 1. Also, state feedback gains are included.Table 2
Hybrid functions and block pulse functions open loop and closed loop results for Example 2.
Method Parameters J State feedback gains Hybrid Functions N=6, M=5 57.9349 K 1 = [ 47.8366     10.7607 ], K 2 = [ - 65.4651     - 26.375     0     - 6.0281 ], K ¯ 1 =   - 17.7816     0.5511 ], K ¯ 2 = [ 9.6281     0     - 2.1138     - 0.7305 ] Block Pulse Functions N=30 69.2999 K 1 = [ 90.9235     - 1.9396 ], K 2 = [ - 131.5504     0     - 30.3858     - 8.1110 ], K ¯ 1 =   - 36.3910     - 9.3258 ], K ¯ 2 = [ 22.8197     14.7181     0     2.3220 ]It is worth noting that for this particular system, control gains determined by the mean of block pulse functions are not stabilizing. While hybrid functions results are illustrated on Figures6 and 7, controlled states with obtained suboptimal feedback designed over a finite horizon tf=3s are drawn on the same figure over a simulation time 5s. It could be seen that suboptimal system states coincide perfectly with optimal solution.Figure 6
Optimal and suboptimal states trajectories of the chemical reactor.Figure 7
Optimal and suboptimal control signals for the chemical reactor.Figure7 exposes optimal control and suggested suboptimal nonlinear state feedback control signals. The proposed nonlinear feedback reproduces sharply the optimal open loop control.


## 6. Conclusion

In this paper a practical approach is developed to solve the problem of finite time quadratic optimal control for polynomial time delay systems. The proposed method is based on the expansion of the system model on a complete set of orthogonal hybrid of block pulse and Legendre polynomials. Two types of optimal control laws have been investigated. In the first step, the method focuses on the determination of the open loop optimal control law. Thus, by defining a general NLP problem for the considered system class, in the second step, a nonlinear delay-depending state feedback control law has been derived in order to meet the optimal states trajectories. The developed results have been illustrated on different examples of nonlinear time delay systems; namely, a two-stage chemical reactor and the obtained results are significant. Note that the proposed method may be enhanced, by decreasing the polynomial system matrices order when introducing the nonredundant form. In addition, if the state information is not completely available, the output feedback is a good choice. The method is feasible for the output feedback control and could be even extended to handle dynamic nonlinear state/output feedback synthesis. Moreover, the presented development is limited to some class of analytical nonlinear systems. In future work, we intend to apply the presented method to the class of switched systems, where the used hybrid functions seems to be convenient to treat that models which are a mixture of smooth functions and piecewise constant signals.



---

*Source: 9596253-2019-04-02.xml*