# Robustness of Operational Matrices of Differentiation for Solving State-Space Analysis and Optimal Control Problems

Authors: Emran Tohidi; F. Soleymani; Adem Kilicman
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The idea of approximation by monomials together with the collocation technique over a uniform mesh for solvingstate-space analysis and optimal control problems (OCPs) has been proposed in this paper. After imposing the Pontryagins maximum principle to the main OCPs, the problems reduce to a linear or nonlinear boundary value problem. In the linear case we propose a monomial collocation matrix approach, while in the nonlinear case, the general collocation method has been applied. We also show the efficiency of the operational matrices of differentiation with respect to the operational matrices of integration in our numerical examples. These matrices of integration are related to the Bessel, Walsh, Triangular, Laguerre, and Hermite functions.


---

## Body

## 1. Introduction

In the last four decades, numerical methods which are based on the operational matrices of integration (especially for orthogonal polynomials and functions) have received considerable attention for dealing with a huge size of applied mathematics problems such asstate-space analysis and optimal control. The key idea of these methods is based on the integral expression
(1)∫0tΦ(τ)dτ≈Φ(t)P,
where Φ(t)=[Φ1(t),Φ2(t),…,ΦN(t)] is an arbitrary basis vector and P is a (N+1)×(N+1) constant matrix, called the operational matrix of integration. The matrix P has already been determined for many types of orthogonal (or nonorthogonal) bases such as Walsh functions [1–3], block-pulse functions [4], Laguerre polynomials [5], Chebyshev polynomials [6], Legendre polynomials [7], Hermite polynomials [8], Fourier series [9], Bernstein polynomials [10], and Bessel functions [11]. As a primary research work which was based on the operational matrices of integration, one can refer to the work of Corrington [1]. In [1], the author proposed a method of solving nonlinear differential and integral equations using a set of Walsh functions as the basis. His method is aimed at obtaining piecewise constant solutions of dynamic equations and requires previously prepared tables of coefficients for integrating Walsh functions. To alleviate the need for such tables, Chen and Hsiao [2, 3] introduced an operational matrix to perform integration of Walsh functions. This operational matrix approach has been applied to various problems such as time-domain analysis and synthesis of linear systems, and piecewise constant-feedback-gain determination for optimal control of linear systems and for inverting irrational Laplace transforms.On the other hand, since the beginning of 1994, the Bernoulli, Chebyshev, Laguerre, Bernstein, Legendre, Taylor, Hermite, and Bessel matrix methods have been used in the works [12–24] to solve high-order linear and nonlinear differential (including hyperbolic partial differential equations) Fredholm Volterra integrodifferential difference delay equations and their systems. The main characteristic of these approaches is based on the operational matrices of differentiation instead of integration. The best advantage of these techniques with respect to the integation methods is that, in the fundamental matrix relations, there is not any approximation symbol, meanwhile in the integration forms such as (1) the approximation symbol could be seen obviously. In other words
(2)Φ′(τ)=Φ(τ)B,
where B is the operational matrix of differentiation for any selected basis such as the previously mentioned polynomials, functions, and truncated series. The readers can see that there is no approximation symbol in (2), meanwhile this can be seen in (1) by using operational matrices of integration. For justifying this expression, one can refer to this subject that after differentiating an Nth degree polynomial we usually reach to a polynomial which has less than Nth degree. However, in the integration processes the degree of polynomials would be increased.In this paper, we generalize a new collocation matrix method that was applied for solving a huge size of applied mathematics models (see for instance [16] and the references therein), to several special classes of systems of ordinary differential equations (ODEs). Two important classes of such systems of ODEs are(i)
State space analysis,(ii)
Hamiltonian system,which are the necessary (and also are sufficient in several special cases) conditions for optimality of the solutions of OCPs, originate from the PMP, and have considerable importance in optimal control and calculus of variation.We again emphasized that the methods that are based on the operational matrices of differentiation are more accurate and effective with regard to the integration ones. We illustrate this fact through several examples for dealing with the previously mentioned systems in the section of numerical examples. It should be noted that one of the best tools for the integration approaches is using high accurate Gauss quadrature rules such as the method of [25, 26]. However, more CPU times are required for using such quadrature rules, and also the matrix coefficient associated to these methods is ill-conditioned usually and should be preconditioned.The remainder of this paper is organized as follows. In Section2, the considered problems such as state-space analysis and Hamiltonian system are introduced. In Section 3, the fundamental matrix relations together with the method of obtaining approximate solutions are described. In Section 4, several numerical examples are provided for confirming high accuracy of the proposed method. The last Section is devoted to the conclusions.


## 2. Problems Statement

In this section two types of problems are considered. In the first subsection, we show that how theHamiltonian systems can be obtained in both linear and nonlinear forms. In the second subsection, we introduce a general form of state-space analysis problems.


### 2.1. Hamiltonian Systems

#### 2.1.1. Linear Quadratic Optimal Control Problems

In this part, we consider the following linear optimal control problem (OCP):(3)minJ=12∫t0tf(xTPx+2xTQu+uTRu)dts.t.x˙=A(t)x(t)+B(t)u(t),x(t0)=x0,
where x∈ℝn, u∈ℝm, A(·)∈ℝn×n and B(·)∈Rm×n. The control u(t) is an admissible control if it is piecewise continuous in t for t∈[t0,tf]. Its values belong to a given closed subset U of ℝm. The input u(t) is derived by minimizing the quadratic performance index J, where P∈ℝn×n is positive semidefinite matrix and R∈ℝm×m is a positive definite matrix. We consider Hamiltonian for system (3) as
(4)H(x,u,λ,t)=12(xTPx+2xTQu+uTRu)+λT(A(t)x+B(t)u),
where λ∈ℝn is the costate vector.According to the Pontryagin’s maximum principle, we have [27]
(5)λ˙=-∂H∂x=-Px-Qu-A(t)Tλ,∂H∂u=QTx+Ru+B(t)Tλ=0.
The optimal control is computed by [27]
(6)u*=-R-1QTx-R-1B(t)Tλ,
where λ and x are the solution of the Hamiltonian system:
(7)x˙=[A(t)-B(t)R-1QT]x-B(t)R-1B(t)Tλ,λ˙=[-P+QR-1QT]x+[QR-1B(t)T-A(t)T]λ,x(t0)=x0,λ(tf)=0.


#### 2.1.2. Nonlinear Quadratic Optimal Control Problems

Consider the nonlinear dynamical system(8)x˙(t)=f(t,x(t))+g(t,x(t))u(t),t∈[t0,tf],x(t0)=x0,
with x(t)∈ℝn denoting the state variable, u(t)∈ℝm the control variable, and x0 is the given initial state at t0. Moreover, f(t,x(t))∈ℝn and g(t,x(t))∈ℝn×m are two continuously differentiable functions in all arguments. Our aim is to minimize the quadratic objective functional
(9)J[x,u]=12∫t0tf(xT(t)Qx(t)+uT(t)Ru(t))dt,
subject to the nonlinear system (8), for Q∈ℝn×n, R∈ℝm×m positive semidefinite and positive definite matrices, respectively. Since the performance index (9) is convex, the following extreme necessary conditions are also sufficient for optimality [28]:
(10)x˙=f(t,x)+g(t,x)u*,λ˙=-Hx(x,u*,λ),u*=argminuH(x,u,λ),x(t0)=x0,λ(tf)=0,
where H(x,u,λ)=(1/2)[xTQx+uTRu]+λT[f(t,x)+g(t,x)u] is referred to the Hamiltonian. Equivalently, (10) can be written in the form of
(11)x˙=f(t,x)+g(t,x)[-R-1gT(t,x)λ]λ˙=-(Qx+(∂f(t,x)∂x)Tλ=.+∑i=1nλi[-R-1gT(t,x)λ]T∂gi(t,x)∂x)x(t0)=x0,λ(tf)=0,
where λ(t)∈ℝn is the costate vector with the ith component λi(t), i=1,…,n and g(t,x)=[g1(t,x),…,gn(t,x)]T with gi(t,x)∈ℝm,  i=1,…,n.Also the optimal control law is obtained by(12)u*=-R-1gT(t,x)λ.For solving such a two-point boundary value problem (TPBVP) in (11), we apply a similar collocation method that was proposed in [29].


### 2.2. State Space Analysis Problems

In this part, we consider the followingstate space analysis problem:
(13)x˙(t)=A(t)x(t)+B(t)u(t),x(0)=x0,
where A(t)∈ℝn×n, B(t)∈ℝn, and u(t)∈ℝ are known, meanwhile x(t)∈ℝn is unknown. The goal is to obtain the approximation of x(t) in (13). The previously mentioned system (13) is similar to Hamiltonian system (7) and the scheme of their solutions is the same.Remark 1.
We recall that the main goal of this paper is to approximate the solution of the systems (7), (11), and (13) by applying a new matrix method which is based on the operational matrix of differentiation and also the uniform collocation scheme in the parts of Hamiltonian systems and state space analysis problems.


## 2.1. Hamiltonian Systems

### 2.1.1. Linear Quadratic Optimal Control Problems

In this part, we consider the following linear optimal control problem (OCP):(3)minJ=12∫t0tf(xTPx+2xTQu+uTRu)dts.t.x˙=A(t)x(t)+B(t)u(t),x(t0)=x0,
where x∈ℝn, u∈ℝm, A(·)∈ℝn×n and B(·)∈Rm×n. The control u(t) is an admissible control if it is piecewise continuous in t for t∈[t0,tf]. Its values belong to a given closed subset U of ℝm. The input u(t) is derived by minimizing the quadratic performance index J, where P∈ℝn×n is positive semidefinite matrix and R∈ℝm×m is a positive definite matrix. We consider Hamiltonian for system (3) as
(4)H(x,u,λ,t)=12(xTPx+2xTQu+uTRu)+λT(A(t)x+B(t)u),
where λ∈ℝn is the costate vector.According to the Pontryagin’s maximum principle, we have [27]
(5)λ˙=-∂H∂x=-Px-Qu-A(t)Tλ,∂H∂u=QTx+Ru+B(t)Tλ=0.
The optimal control is computed by [27]
(6)u*=-R-1QTx-R-1B(t)Tλ,
where λ and x are the solution of the Hamiltonian system:
(7)x˙=[A(t)-B(t)R-1QT]x-B(t)R-1B(t)Tλ,λ˙=[-P+QR-1QT]x+[QR-1B(t)T-A(t)T]λ,x(t0)=x0,λ(tf)=0.


### 2.1.2. Nonlinear Quadratic Optimal Control Problems

Consider the nonlinear dynamical system(8)x˙(t)=f(t,x(t))+g(t,x(t))u(t),t∈[t0,tf],x(t0)=x0,
with x(t)∈ℝn denoting the state variable, u(t)∈ℝm the control variable, and x0 is the given initial state at t0. Moreover, f(t,x(t))∈ℝn and g(t,x(t))∈ℝn×m are two continuously differentiable functions in all arguments. Our aim is to minimize the quadratic objective functional
(9)J[x,u]=12∫t0tf(xT(t)Qx(t)+uT(t)Ru(t))dt,
subject to the nonlinear system (8), for Q∈ℝn×n, R∈ℝm×m positive semidefinite and positive definite matrices, respectively. Since the performance index (9) is convex, the following extreme necessary conditions are also sufficient for optimality [28]:
(10)x˙=f(t,x)+g(t,x)u*,λ˙=-Hx(x,u*,λ),u*=argminuH(x,u,λ),x(t0)=x0,λ(tf)=0,
where H(x,u,λ)=(1/2)[xTQx+uTRu]+λT[f(t,x)+g(t,x)u] is referred to the Hamiltonian. Equivalently, (10) can be written in the form of
(11)x˙=f(t,x)+g(t,x)[-R-1gT(t,x)λ]λ˙=-(Qx+(∂f(t,x)∂x)Tλ=.+∑i=1nλi[-R-1gT(t,x)λ]T∂gi(t,x)∂x)x(t0)=x0,λ(tf)=0,
where λ(t)∈ℝn is the costate vector with the ith component λi(t), i=1,…,n and g(t,x)=[g1(t,x),…,gn(t,x)]T with gi(t,x)∈ℝm,  i=1,…,n.Also the optimal control law is obtained by(12)u*=-R-1gT(t,x)λ.For solving such a two-point boundary value problem (TPBVP) in (11), we apply a similar collocation method that was proposed in [29].


## 2.1.1. Linear Quadratic Optimal Control Problems

In this part, we consider the following linear optimal control problem (OCP):(3)minJ=12∫t0tf(xTPx+2xTQu+uTRu)dts.t.x˙=A(t)x(t)+B(t)u(t),x(t0)=x0,
where x∈ℝn, u∈ℝm, A(·)∈ℝn×n and B(·)∈Rm×n. The control u(t) is an admissible control if it is piecewise continuous in t for t∈[t0,tf]. Its values belong to a given closed subset U of ℝm. The input u(t) is derived by minimizing the quadratic performance index J, where P∈ℝn×n is positive semidefinite matrix and R∈ℝm×m is a positive definite matrix. We consider Hamiltonian for system (3) as
(4)H(x,u,λ,t)=12(xTPx+2xTQu+uTRu)+λT(A(t)x+B(t)u),
where λ∈ℝn is the costate vector.According to the Pontryagin’s maximum principle, we have [27]
(5)λ˙=-∂H∂x=-Px-Qu-A(t)Tλ,∂H∂u=QTx+Ru+B(t)Tλ=0.
The optimal control is computed by [27]
(6)u*=-R-1QTx-R-1B(t)Tλ,
where λ and x are the solution of the Hamiltonian system:
(7)x˙=[A(t)-B(t)R-1QT]x-B(t)R-1B(t)Tλ,λ˙=[-P+QR-1QT]x+[QR-1B(t)T-A(t)T]λ,x(t0)=x0,λ(tf)=0.


## 2.1.2. Nonlinear Quadratic Optimal Control Problems

Consider the nonlinear dynamical system(8)x˙(t)=f(t,x(t))+g(t,x(t))u(t),t∈[t0,tf],x(t0)=x0,
with x(t)∈ℝn denoting the state variable, u(t)∈ℝm the control variable, and x0 is the given initial state at t0. Moreover, f(t,x(t))∈ℝn and g(t,x(t))∈ℝn×m are two continuously differentiable functions in all arguments. Our aim is to minimize the quadratic objective functional
(9)J[x,u]=12∫t0tf(xT(t)Qx(t)+uT(t)Ru(t))dt,
subject to the nonlinear system (8), for Q∈ℝn×n, R∈ℝm×m positive semidefinite and positive definite matrices, respectively. Since the performance index (9) is convex, the following extreme necessary conditions are also sufficient for optimality [28]:
(10)x˙=f(t,x)+g(t,x)u*,λ˙=-Hx(x,u*,λ),u*=argminuH(x,u,λ),x(t0)=x0,λ(tf)=0,
where H(x,u,λ)=(1/2)[xTQx+uTRu]+λT[f(t,x)+g(t,x)u] is referred to the Hamiltonian. Equivalently, (10) can be written in the form of
(11)x˙=f(t,x)+g(t,x)[-R-1gT(t,x)λ]λ˙=-(Qx+(∂f(t,x)∂x)Tλ=.+∑i=1nλi[-R-1gT(t,x)λ]T∂gi(t,x)∂x)x(t0)=x0,λ(tf)=0,
where λ(t)∈ℝn is the costate vector with the ith component λi(t), i=1,…,n and g(t,x)=[g1(t,x),…,gn(t,x)]T with gi(t,x)∈ℝm,  i=1,…,n.Also the optimal control law is obtained by(12)u*=-R-1gT(t,x)λ.For solving such a two-point boundary value problem (TPBVP) in (11), we apply a similar collocation method that was proposed in [29].


## 2.2. State Space Analysis Problems

In this part, we consider the followingstate space analysis problem:
(13)x˙(t)=A(t)x(t)+B(t)u(t),x(0)=x0,
where A(t)∈ℝn×n, B(t)∈ℝn, and u(t)∈ℝ are known, meanwhile x(t)∈ℝn is unknown. The goal is to obtain the approximation of x(t) in (13). The previously mentioned system (13) is similar to Hamiltonian system (7) and the scheme of their solutions is the same.Remark 1.
We recall that the main goal of this paper is to approximate the solution of the systems (7), (11), and (13) by applying a new matrix method which is based on the operational matrix of differentiation and also the uniform collocation scheme in the parts of Hamiltonian systems and state space analysis problems.


## 3. Fundamental Matrix Relations and Method of the Solution

In this section, by using the collocation points and the matrix relations between the monomials{1,t,t2,…,tN} and their derivatives, we will find the approximate solution of the system (7) expressed in the truncated monomial series form (assuming that x(t) and λ(t)∈ℝ and also A(t) together with B(t) is independent of time t, that is, A=A(t) and B=B(t))
(14)x(t)≈xN(t)=∑n=0Na1,ntn,λ(t)≈λN(t)=∑n=0Na2,ntn,
so that a1,n and a2,n; n=0,1,2,…,N are the unknown coefficients.Let us consider the desired solutionsx(t) and λ(t), of (7) defined by the truncated monomial series (14). We can write the approximate solutions, which are given in relation (14) in the matrix form
(15)xN(t)=X(t)A1,λN(t)=X(t)A2,
where X(t)=[1t⋯tN] and Ai=[ai,0ai,1⋯ai,N]T, i=1,2.The matrix form of the relation between the matrixX(t) and its kth derivative X(k)(t) is
(16)X(k)(t)=X(t)(BT)k
so that
(17)BT=[010⋯0002⋯0⋮⋮⋮⋱⋮000⋯N000⋯0],(BT)0=[I](N+1)×(N+1)istheidentitymatrix.
By using the relations (15) and (16), we have the following relations:
(18)xN(i)(t)=X(t)(BT)iA1,λN(i)(t)=X(t)(BT)iA2,i=0,1.
Thus, we can express the matrices y(t) and y(1)(t) as follows:
(19)yN(i)(t)=X-(t)(B-)iA,i=0,1,
where
(20)yN(i)(t)=[xN(i)(t)λN(i)(t)],A=[A1A2],X-(t)=[X(t)00X(t)],B-=[BT00BT].Now, we can restate the system (7) in the matrix form
(21)yN(1)(t)-MyN(t)=O,
where
(22)yN(1)(t)=[xN(1)(t)λN(1)(t)],M=[A-BR-1QT-BR-1BT-P+QR-1QTQR-1BT-AT],yN(t)=[xN(t)λN(t)],O=[00].
Applying the following collocation points in (21):
(23)ts=t0+tf-t0Ns,s=0,1,…,N
yields to N+1 equations as follows:
(24)yN(1)(ts)-MyN(ts)=O,s=0,1,…,N.
All of the these equations can be written in the following matrix form:
(25)Y(1)-M-Y=O,
where
(26)Y(1)=[yN(1)(t0)yN(1)(t1)⋮yN(1)(tN)],M-=I⊗M=kron(I,M),Y=[yN(t0)yN(t1)⋮yN(tN)],
where ⊗ denotes the Kronecker product and I is the identity matrix of dimension N+1.With the aid of relation (19) and the collocation points (23), we gain
(27)yN(ts)=X-(ts)A,yN(1)(ts)=X-(ts)B-A,s=0,1,…,N,
which can be written as
(28)Y=XA,Y(1)=XB-A,
where
(29)X=[X-(t0)X-(t1)⋮X-(tN)],X-(ts)=[X(ts)00X(ts)],s=0,1,…,N.
If the relation (28) is substituted into (25), the fundamental matrix equation is obtained as
(30){XB--M-X}A=O.
Thus, the fundamental matrix equation (30) corresponding to (7) can be written in the form
(31)WA=Oor[W;O],
which corresponds to a linear system of 2(N+1) algebraic equations in 2(N+1) the unknown monomial coefficients so that
(32)W={XB--M-X}=[wp,q],p,q=1,2,…,2(N+1).
By the aid of the relation (19), the matrix form for the boundary conditions which are given in (7) can be written as
(33)UA=[R-],or[U;R-].Finally, by replacing the rows of the matrices[U;R-] by the last rows of the matrices [W;O], we obtain the new augmented matrix
(34)W~A=O~.
The unknown monomials coefficients which exist in A are determined by solving this linear system, and hence ai,0,ai,1,…,ai,N, (i=1,2) are substituted in (14). Therefore, we find the approximated solutions
(35)xN(t)=∑n=0Na1,ntn,λN(t)=∑n=0Na2,ntn.We can easily check the accuracy of the method. Since the truncated monomial series (14) are the approximate solutions of (7), when the functions xN(t), λN(t) and their derivatives are substituted in (7), the resulting equation must be satisfied approximately; that is, for t=tq∈[t0,tf],q=0,1,2,…(36)E1(tq)=|x(1)(tq)-[A(t)-B(t)R-1QT]x(tq)z+B(t)R-1B(t)Tλ(tq)|≅0,E2(tq)=|λ(1)(tq)-[-P+QR-1QT]x(tq)z-[QR-1B(t)T-A(t)T]λ(tq)|≅0,
and Ei,N(tq)≤10-kq, i=1,2 (kq positive integer).Ifmax10-kq=10-k (k positive integer) is prescribed, then the truncation limit N is increased until the difference Ei,N(tq) at each of the points becomes smaller than the prescribed 10-k, see [24].Remark 2.
We must recall that a similar approach can be applied for thestate space analysis problem (13). Moreover, as we say before, for solving a general nonlinear system of ODEs such as (11), we apply a generalization of the collocation method that was proposed in [29].


## 4. Numerical Examples

In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed method. All calculations are designed in MAPLE 13 and run on a Pentium 4 PC Laptop with 2 GHz of CPU and 2 GB of RAM. In this regard, in tables and figures, we report the absolute error functions associated to the trajectory and control variables and also the approximated values of performance index. In the first example, we provide an OCP that was recently considered by a new method (which is based on the operational matrix of integration of Triangular functions) [30] and reach to more accurate results. Also, in the second example, we consider another OCP (or Hamiltonian system) with time variant dynamical system, in which our results have more accuracy and credit with regard to methods [30, 31]. Moreover, we consider a nonlinear OCP as our third numerical illustration. In the fourth example, we provide a state space analysis problem together with a full comparison with the methods that are based on the operational matrices of integration such as Bessel [11] and Laguerre [32].Example 3 (see [30] linear Hamiltonian system).
Consider the problem of minimizing(37)J=12∫01(2x2(t)+u2(t))dt,
subject to
(38)x˙(t)=-12x(t)+u(t),x(0)=1.
The purpose is to find the optimal control u(t) which minimizes (37) subject to (38). The Optimal value of performance index for this problem is J*=0.463566653481105 and also exact solutions have been given in [30] as
(39)x*(t)=12+e-3(2e-(3t/2)+e-3+(3t/2)),u*(t)=12+e-3(-e-(3t/2)+e-3+(3t/2)).
Since the objective function of this OCP is convex, therefore the following necessary conditions (i.e., linear Hamiltonian system) for optimality are also sufficient:
(40)x˙(t)=-12x(t)-λ(t),λ˙(t)=-2x(t)+12λ(t),x(0)=1,λ(1)=0.
Hence, we need to solve the previous system of differential equations such that the obtained numerical solution is the optimal solution of problem (37)-(38). It should be noted that according to (6) the optimal control is computed by u*(t)=-λ(t), where λ(t) is the solution of the previous system.
We solve this problem by using our proposed method in the cases ofN=4,5,6,7, and 8. The approximated solutions corresponding to these values of N are provided below
(41)x4(t)=0.130043t4-0.4921311t3+1.1171671t2-1.4266579t+1.0,u4(t)=-0.1020785t4+0.5247564t3-1.0326910t2+1.5366710t-0.9266579,x5(t)=-0.0317935t5+0.1886526t4-0.5272247t3+1.123919t2-1.427164t+1.0,u5(t)=0.04038275t5-0.1772002t4+0.5697892t3-1.042226t2+1.536418t-0.9271638,x6(t)=0.009145352t6-0.05359579t5+0.2075571t4-0.5343086t3+1.124911t2-1.4271300t+1.0,u6(t)=-0.006951979t6+0.05688585t5-0.1914114t4+0.5750789t3-1.042907t2+1.536435t-0.9271298,x7(t)=-0.001606958t7+0.01380419t6-0.05883564t5+0.210416t4-0.5350703t3+1.124991t2-1.427133t+1.0,u7(t)=0.002084015t7-0.01300909t6+0.06371979t5-0.1951577t4+0.5760817t3-1.043019t2+1.536433t-0.9271335,x8(t)=0.0003561777t8-0.002824592t7+0.01548942t6-0.06004358t5+0.2108896t4-0.5351673t3+1.124999t2-1.427133t+1.0,u8(t)=-0.0002657152t8+0.00299044t7-0.01426032t6+0.06461366t5-0.1955065t4+0.5761526t3-1.043024t2+1.536433t-0.9271333.
The associated performance indexes for the selected values of N are J4=0.463550469, J5=0.4635575131, J6=0.4635665731, J7=0.463566618, and J8=0.4635666532. We provide the exN=max0≤t≤1|xN(t)-x*(t)|, eJN=|JN-J*| associated to our proposed method (PM) and a new method that is based on the operational matrix of integration of Triangular functions [30] for different values of N in Table 1. It can be seen from this table that our obtained results for such considered values of N (i.e., 4, 5, 6, 7, and 8) are the same and equal to the obtained results of [30] for higher values of N such as 4,8,16,32, and 64 in computation of exN. Moreover, our results corresponding to the eJN are more accurate with regard to the method of [30] even by choosing lower values of N.Table 1
Comparison results of Example3.
N
e
x
N of PM
e
J
N of PM
N
e
x
N of TFM [30]
e
J
N of TFM [30]
4
1.8441
e
-
003
1.6184
e
-
005
4
7.68816
e
-
03
4.27347
e
-
03
5
2.2886
e
-
004
9.1404
e
-
006
8
3.19260
e
-
03
1.06808
e
-
03
6
3.3538
e
-
005
8.0381
e
-
008
16
1.74696
e
-
04
2.67002
e
-
04
7
4.5138
e
-
005
3.5481
e
-
008
32
4.36429
e
-
05
6.67494
e
-
05
8
3.0718
e
-
005
2.8111
e
-
010
64
1.09092
e
-
05
1.66873
e
-
05Example 4 (see [30, 31] linear Hamiltonian system).
Consider the linear time-varying system(42)x˙(t)=tx(t)+u(t),
with the cost functional
(43)J=12∫01(x2(t)+u2(t))dt.
The problem is to obtain the optimal control u*(t) which minimizes (43) subject to (42). The optimal control is
(44)u*(t)=-K(t)x*(t),
where K(t) is the feedback controller gain matrix and the solution of the Riccati equation [30]
(45)K(t)=K2(t)-2tK(t)-1,K(1)=0.
According to the optimality conditions (5) and (6) we have
(46)x˙=tx(t)-λ(t),λ˙=-x(t)-tλ(t),x(0)=1,λ(1)=0,u(t)=-λ(t).
We first solve the previous system and obtain the numerical solutions xN(t) and λN(t) for N=4,5,6, and then solve (45) by ODE solver commands which exist in MAPLE 13. Since, K(t)=-(u(t)/x(t)) then our numerical results of K(t) are equal to -(uN(t)/xN(t)), that is, KN(t)=-(uN(t)/xN(t)). The numerical results of system (46), which are obtained by the proposed method could be deduced as
(47)x4(t)=0.20203304t4-0.31263421t3+1.0000762t2-0.96700804t+1.0,u4(t)=0.011211016t4-0.055485904t3+0.011282926t2+1.0t-0.96700804,x5(t)=0.021615359t5+0.16221315t4-0.28939672t3+0.99513749t2-0.96857084t+1.0,u5(t)=0.047811786t5-0.078776302t4-0.00037310025t3-0.000091547293t2+1.0t-0.96857084,x6(t)=0.043743856t6-0.083312898t5+0.25396668t4-0.32429532t3+1.0001826t2-0.96853474t+1.0,u6(t)=-0.002379348t6+0.052612074t5-0.081954336t4+0.00027888358t3-0.000022529694t2+1.0t-0.96853474.
Also, the exact solution K(t) of the Riccati equation (45) at the uniform mesh in the interval (0,1) are K(0)=9.6854e-001, K(0.1)=9.5147e-001, K(0.2)=9.1063e-001, K(0.3)=8.4416e-001, K(0.4)=7.5241e-001, K(0.5)=6.3856e-001, K(0.6)=5.0873e-001, K(0.7)=3.7127e-001, K(0.8)=2.3540e-001, K(0.9)=1.0955e-001, and K(1)=0. In Table 2, we provide the absolute values of errors at the selected points for the previously considered values of N together with the same errors associated with other methods [30, 31]. Again, we can see the accuracy of method with regard to the methods that are based on operational matrices of integration.Table 2
The error histories at the selected points for Example4.
t
PM forN=4
PM forN=5
PM forN=6
TFM [30] N=6
MGL [31] N=6
TFM [30] N=64
0.0
1.5274
e
-
003
3.5407
e
-
005
6.8533
e
-
007
7.4890
e
-
002
9.3400
e
-
002
7.6980
e
-
003
0.1
1.9116
e
-
003
6.7725
e
-
005
1.4752
e
-
006
2.2140
e
-
002
1.0020
e
-
001
1.6770
e
-
003
0.2
2.3352
e
-
003
9.5549
e
-
005
1.7583
e
-
006
5.2150
e
-
002
1.0660
e
-
001
5.6050
e
-
003
0.3
2.6500
e
-
003
9.9318
e
-
005
1.8482
e
-
006
7.1120
e
-
002
1.5530
e
-
001
5.8040
e
-
003
0.4
2.8209
e
-
003
1.0115
e
-
004
2.2289
e
-
006
1.5220
e
-
002
2.3530
e
-
001
2.0900
e
-
003
0.5
2.9118
e
-
003
1.1223
e
-
004
2.5894
e
-
006
9.9270
e
-
002
2.9390
e
-
001
1.0200
e
-
002
0.6
2.9909
e
-
003
1.1855
e
-
004
2.3877
e
-
006
2.8060
e
-
002
2.7180
e
-
001
1.9860
e
-
003
0.7
3.0063
e
-
003
1.0774
e
-
004
2.1119
e
-
006
5.7040
e
-
002
1.7540
e
-
001
5.9630
e
-
003
0.8
2.7289
e
-
003
9.2586
e
-
005
2.4315
e
-
006
6.5130
e
-
002
7.5600
e
-
002
5.5260
e
-
003
0.9
1.8221
e
-
003
8.1824
e
-
005
1.6928
e
-
006
1.5240
e
-
002
1.8200
e
-
002
1.7450
e
-
003Example 5 (nonlinear Hamiltonian system).
As our third illustration, consider the following nonlinear optimal control problem:(48)minJ=∫01u2(t)dts.t.x˙=12x2(t)sinx(t)+u(t),t∈[0,1]s...zx(0)=0,x(1)=0.5.
Trivially f(t,x(t))=(1/2)x2(t)sinx(t), g(t,x(t))=1, Q=0, R=1, t0=0, and tf=1. As mentioned in Section 2.2, we solve the following system of ordinary differential equations:
(49)x˙=12x2(t)sinx(t)-12λ(t),λ˙=-λ(t)x(t)sinx(t)-12λ(t)x2(t)cosx(t),-λ(t)x(t)==zzz0(t)zzz(t)…x2t∈[0,1],x(0)=0,x(1)=0.5,λ(1)=0.
Also the optimal control law is given by
(50)u*(t)=-12λ(t).
Similar to the linear cases, we suppose that the state and costate variables could be written in terms of linear combination of monic polynomials which are defined in Section3, with the unknown monomial coefficients. These coefficients will be determined after imposing the previous system of differential equations at the uniform mesh in the interval (0,1). In other words, applying these collocation points to the main system together with the considered boundary conditions on x(t) and λ(t) transforms the basic problem to the corresponding system of nonlinear algebraic equations. By assuming different values of N such as 5,  7, and 9, we solve the previously mentioned system. In Table 3, we provide the approximated performance index JN, which is obtained by our proposed method and also the difference between the approximated xN(1)-x(1)=xN(1)-(1/2) for the considered values of N.Table 3
Numerical results of Example5.
N
J
N of PM
x
N
(
1
)
-
x
(
1
) of PM
5
0.2426752
0
7
0.2426750
1.0000
e
-
007
9
0.2426740
2.0000
e
-
015Example 6 (see [11, 32] state-space analysis).
We consider a linear-time invariant state equation(51)[x˙1(t)x˙2(t)]=A[x1(t)x2(t)]+bu(t),t∈(0,1),
where
(52)A=[-312-4],b=[11].
We are given that the input u(t) is the unit step function in the interval (0,1) and the initial state is
(53)x0=[x1(0)x2(0)]=[11].
The exact solution for (51) is
(54)x1(t)=x2(t)=1+e-2t2.
We solve this problem by using our proposed method in the cases of N=4,5,6,7,8,9,10,12, and 20. The approximated solutions corresponding to the N=4,5,6,7,8,9, and 10 are provided later
(55)x1,4(t)=x2,4(t)=0.1658t4-0.58031t3+0.98446t2-1.0t+1.0,x1,5(t)=x2,5(t)=-0.06235t5+0.28058t4-0.64844t3+0.99761t2-1.0t+1.0,x1,6(t)=x2,6(t)=0.019951t6-0.10973t5+0.32143t4-0.66365t3+0.99969t2-1.0t+1.0,x1,7(t)=x2,7(t)=-0.0055374t7+0.035993t6-0.12776t5+0.33126t4-0.66626t3+0.99997t2-1.0t+1.0,x1,8(t)=x2,8(t)=0.0013547t8-0.010161t7+0.042378t6-0.13232t5+0.33304t4-0.66662t3+1.0t2-1.0t+1.0,x1,9(t)=x2,9(t)=-0.00029605t9+0.0025165t8-0.012061t7+0.044057t6-0.13318t5+0.3333t4-0.66666t3+1.0t2-1.0t+1.0,x1,10(t)=x2,10(t)=0.000058425t10-0.00055504t9+0.003006t8-0.012576t7+0.044386t6-0.13331t5+0.33333t4-0.66667t3+1.0t2-1.0t+1.0.
An interesting comparison between our presented method (PM) and the method of [11] in the absolute value of the errors at the uniform mesh for N=8,12, and 20 is considered in Table 4. Moreover, the error histories in the computational interval (0,1) associated with our method for N=4 and 5 together with the error history of the method [11] (Bessel Integration) for N=4 and also the error history of the method [32] (Laguerre Integration) for N=5 are depicted in Figure 1. From this figure one can see the applicability and high accuracy of the presented method with respect to the methods which are based on operational matrices of integration such as [11, 32].Table 4
The error histories at the selected points for Example6.
t
Bessel IntN=8
Bessel IntN=12
Bessel IntN=20
PM forN=8
PM forN=12
PM forN=20
0.0
9.6890
e
-
003
6.6498
e
-
003
4.1214
e
-
003
0
0
0
0.1
2.3974
e
-
002
1.0755
e
-
002
1.6289
e
-
003
6.5649
e
-
009
6.8945
e
-
014
6.1062
e
-
014
0.2
7.2706
e
-
003
1.6306
e
-
003
4.6111
e
-
003
4.5248
e
-
009
5.0071
e
-
014
4.8073
e
-
014
0.3
2.5780
e
-
003
4.6009
e
-
003
8.9332
e
-
003
3.6167
e
-
009
4.1522
e
-
014
3.7859
e
-
014
0.4
2.2427
e
-
003
7.2766
e
-
003
1.1642
e
-
002
3.2895
e
-
009
3.3973
e
-
014
3.0309
e
-
014
0.5
6.3782
e
-
003
1.0152
e
-
002
1.3160
e
-
002
2.2422
e
-
009
2.7311
e
-
014
2.4425
e
-
014
0.6
6.9801
e
-
003
1.1190
e
-
002
1.3796
e
-
002
2.4272
e
-
009
2.2649
e
-
014
1.9873
e
-
014
0.7
9.5374
e
-
003
1.1427
e
-
002
1.3807
e
-
002
1.1916
e
-
009
1.8097
e
-
014
1.6431
e
-
014
0.8
8.6910
e
-
003
1.1990
e
-
002
1.3398
e
-
002
1.8602
e
-
009
1.0991
e
-
014
1.3878
e
-
014
0.9
1.0169
e
-
002
1.1068
e
-
002
1.2680
e
-
002
4.0343
e
-
009
6.6391
e
-
014
1.8874
e
-
014
1.0
5.8131
e
-
003
8.4105
e
-
003
1.0482
e
-
002
2.2065
e
-
007
3.9901
e
-
012
8.3072
e
-
012Figure 1
The error histories of our method together with the methods [11, 32] in Example 6.


## 5. Conclusions

The aim of this paper is to present an indirect approach for solving optimal control problems using truncated monomial series together with the collocation method on a uniform mesh. Our method applies an operational matrix of differentiation which has more efficiency with respect to the integration ones. Operational matrices of differentiation have several specific properties that other integration ones do not have them. One of them is that through using differentiation ones, we solve our problem directly and do not need to integrate the dynamical system. Another property is that the fundamental relations, which are based on differentiation matrices, are the exact relations, meanwhile the methods which are based on integration matrices impose the approximation to the main problem. These properties are shown through several numerical examples such as state-space analysis and specially optimal control problems.



---

*Source: 535979-2013-04-04.xml*