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

Abstract

The idea of approximation by monomials together with the collocation technique over a uniform mesh for solving state-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.

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 as state-space analysis and optimal control. The key idea of these methods is based on the integral expression ^𝑡₀Φ(𝜏)𝑑𝜏≈Φ(𝑡)𝑃, (1) where Φ(𝑡)=[Φ₁(𝑡),Φ₂(𝑡),…,Φ_𝑁(𝑡)] is an arbitrary basis vector and 𝑃 is a (𝑁+1)×(𝑁+1) constant matrix, called the operational matrix of integration. The matrix 𝑃 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) where 𝐵 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 𝑁th degree polynomial we usually reach to a polynomial which has less than 𝑁th 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 Section 2, 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 the Hamiltonian 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): 1min𝐽=2^𝑡𝑓_𝑡0𝑥^𝑇𝑃𝑥+2𝑥^𝑇𝑄𝑢+𝑢^𝑇̇(𝑡)𝑥(𝑡)(𝑡)𝑢(𝑡)𝑡𝑅𝑢𝑑𝑡s.t.𝑥=𝐴+𝐵,𝑥₀=𝑥₀, (3) where 𝑥∈ℝ^𝑛, 𝑢∈ℝ^𝑚, 𝐴(⋅)∈ℝ^𝑛×𝑛 and 𝐵(⋅)∈𝑅^𝑚×𝑛. The control 𝑢(𝑡) is an admissible control if it is piecewise continuous in 𝑡 for 𝑡∈[𝑡₀,𝑡_𝑓]. Its values belong to a given closed subset 𝑈 of ℝ^𝑚. The input 𝑢(𝑡) is derived by minimizing the quadratic performance index 𝐽, where 𝑃∈ℝ^𝑛×𝑛 is positive semidefinite matrix and 𝑅∈ℝ^𝑚×𝑚 is a positive definite matrix. We consider Hamiltonian for system (3) as 𝐻()=1𝑥,𝑢,𝜆,𝑡2𝑥^𝑇𝑃𝑥+2𝑥^𝑇𝑄𝑢+𝑢^𝑇𝑅𝑢+𝜆^𝑇(𝐴(𝑡)𝑥+𝐵(𝑡)𝑢), (4) where 𝜆∈ℝ^𝑛 is the costate vector.

According to the Pontryagin’s maximum principle, we have [27] 𝜆̇=−𝜕𝐻(𝑡)𝜕𝑥=−𝑃𝑥−𝑄𝑢−𝐴^𝑇𝜆,𝜕𝐻𝜕𝑢=𝑄^𝑇𝑥+𝑅𝑢+𝐵(𝑡)^𝑇𝜆=0. (5) The optimal control is computed by [27] 𝑢^∗=−𝑅⁻¹𝑄^𝑇𝑥−𝑅⁻¹𝐵(𝑡)^𝑇𝜆, (6) where 𝜆 and 𝑥 are the solution of the Hamiltonian system: 𝑥̇=𝐴(𝑡)−𝐵(𝑡)𝑅⁻¹𝑄^𝑇𝑥−𝐵(𝑡)𝑅⁻¹𝐵(𝑡)^𝑇̇𝜆,𝜆=−𝑃+𝑄𝑅⁻¹𝑄^𝑇𝑥+𝑄𝑅⁻¹𝐵(𝑡)^𝑇−𝐴(𝑡)^𝑇𝑥𝑡𝜆,₀=𝑥₀,𝜆𝑡_𝑓=0. (7)

2.1.2. Nonlinear Quadratic Optimal Control Problems

Consider the nonlinear dynamical system 𝑥̇(𝑡)=𝑓(𝑡,𝑥(𝑡))+𝑔(𝑡,𝑥(𝑡))𝑢(𝑡),𝑡∈𝑡₀,𝑡_𝑓,𝑥𝑡₀=𝑥₀, (8) with 𝑥(𝑡)∈ℝ^𝑛 denoting the state variable, 𝑢(𝑡)∈ℝ^𝑚 the control variable, and 𝑥₀ is the given initial state at 𝑡₀. Moreover, 𝑓(𝑡,𝑥(𝑡))∈ℝ^𝑛 and 𝑔(𝑡,𝑥(𝑡))∈ℝ^𝑛×𝑚 are two continuously differentiable functions in all arguments. Our aim is to minimize the quadratic objective functional 𝐽[]=1𝑥,𝑢2^𝑡𝑓_𝑡0𝑥^𝑇(𝑡)𝑄𝑥(𝑡)+𝑢^𝑇(𝑡)𝑅𝑢(𝑡)𝑑𝑡, (9) subject to the nonlinear system (8), for 𝑄∈ℝ^𝑛×𝑛, 𝑅∈ℝ^𝑚×𝑚 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]: 𝑥̇=𝑓(𝑡,𝑥)+𝑔(𝑡,𝑥)𝑢^∗,𝜆̇=−𝐻_𝑥𝑥,𝑢^∗,𝑢,𝜆_∗=argmin_𝑢𝐻(𝑥,𝑢,𝜆),𝑥𝑡₀=𝑥₀,𝜆𝑡_𝑓=0, (10) where 𝐻(𝑥,𝑢,𝜆)=(1/2)[𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢]+𝜆^𝑇[𝑓(𝑡,𝑥)+𝑔(𝑡,𝑥)𝑢] is referred to the Hamiltonian. Equivalently, (10) can be written in the form of 𝑥̇=𝑓(𝑡,𝑥)+𝑔(𝑡,𝑥)−𝑅⁻¹𝑔^𝑇()𝜆̇()𝑡,𝑥𝜆=−𝑄𝑥+𝜕𝑓𝑡,𝑥𝜕𝑥^𝑇𝜆+^𝑛_𝑖=1𝜆_𝑖−R⁻¹𝑔^𝑇(𝑡,𝑥)𝜆^𝑇𝜕𝑔_𝑖(𝑡,𝑥)𝑥𝑡𝜕𝑥₀=𝑥₀,𝜆𝑡_𝑓=0, (11) where 𝜆(𝑡)∈ℝ^𝑛 is the costate vector with the 𝑖th component 𝜆_𝑖(𝑡), 𝑖=1,…,𝑛 and 𝑔(𝑡,𝑥)=[𝑔₁(𝑡,𝑥),…,𝑔_𝑛(𝑡,𝑥)]^𝑇 with 𝑔_𝑖(𝑡,𝑥)∈ℝ^𝑚,  𝑖=1,…,𝑛.

Also the optimal control law is obtained by 𝑢^∗=−𝑅⁻¹𝑔^𝑇(𝑡,𝑥)𝜆. (12)

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 following state space analysis problem: 𝑥̇(𝑡)=𝐴(𝑡)𝑥(𝑡)+𝐵(𝑡)𝑢(𝑡),𝑥(0)=𝑥₀, (13) where 𝐴(𝑡)∈ℝ^𝑛×𝑛, 𝐵(𝑡)∈ℝ^𝑛, and 𝑢(𝑡)∈ℝ are known, meanwhile 𝑥(𝑡)∈ℝ^𝑛 is unknown. The goal is to obtain the approximation of 𝑥(𝑡) 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,𝑡,𝑡²,…,𝑡^𝑁} and their derivatives, we will find the approximate solution of the system (7) expressed in the truncated monomial series form (assuming that 𝑥(𝑡) and 𝜆(𝑡)∈ℝ and also 𝐴(𝑡) together with 𝐵(𝑡) is independent of time 𝑡, that is, 𝐴=𝐴(𝑡) and 𝐵=𝐵(𝑡)) 𝑥(𝑡)≈𝑥_𝑁(𝑡)=^𝑁_𝑛=0𝑎_1,𝑛𝑡^𝑛,𝜆(𝑡)≈𝜆_𝑁(𝑡)=^𝑁_𝑛=0𝑎_2,𝑛𝑡^𝑛, (14) so that 𝑎_1,𝑛 and 𝑎_2,𝑛; 𝑛=0,1,2,…,𝑁 are the unknown coefficients.

Let us consider the desired solutions 𝑥(𝑡) and 𝜆(𝑡), 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) where 𝐗(𝑡)=[1𝑡⋯𝑡^𝑁] and 𝐀_𝑖=[𝑎_𝑖,0𝑎_𝑖,1⋯𝑎_𝑖,𝑁]^𝑇, 𝑖=1,2.

The matrix form of the relation between the matrix 𝐗(𝑡) and its 𝑘th derivative 𝐗^(𝑘)(𝑡) is 𝐗^(𝑘)(𝑡)=𝐗(𝑡)𝐁^𝑇𝑘 (16) so that 𝐁^𝑇=,𝐁010⋯0002⋯0⋮⋮⋮⋱⋮000⋯𝑁000⋯0_𝑇0=[𝐈]_{(𝑁+1)×(𝑁+1)}istheidentitymatrix. (17) By using the relations (15) and (16), we have the following relations: 𝑥_𝑁^(𝑖)(𝑡)=𝐗(𝑡)𝐁^𝑇𝑖𝐀₁,𝜆_𝑁^(𝑖)(𝑡)=𝐗(𝑡)𝐁^𝑇𝑖𝐀₂,𝑖=0,1. (18) Thus, we can express the matrices 𝐲(𝑡) and 𝐲⁽¹⁾(𝑡) as follows: 𝐲_𝑁^(𝑖)(𝑡)=𝐗(𝑡)𝐁^𝑖𝐀,𝑖=0,1, (19) where 𝐲_𝑁^(𝑖)(𝑡)=𝑥_𝑁^(𝑖)(𝑡)𝜆_𝑁^(𝑖)(𝑡)𝐀,𝐀=¹𝐀₂,𝐗(𝑡)=𝐗(𝑡)0(𝑡),0𝐗𝐁𝐁=^𝑇00𝐁^𝑇. (20)

Now, we can restate the system (7) in the matrix form 𝐲_𝑁⁽¹⁾(𝑡)−𝑀𝐲_𝑁(𝑡)=𝑂, (21) where 𝐲_𝑁⁽¹⁾(𝑡)=𝑥_𝑁⁽¹⁾(𝑡)𝜆_𝑁⁽¹⁾(𝑡),𝐌=𝐴−𝐵𝑅⁻¹𝑄^𝑇−𝐵𝑅⁻¹𝐵^𝑇−𝑃+𝑄𝑅⁻¹𝑄^𝑇𝑄𝑅⁻¹𝐵^𝑇−𝐴^𝑇,𝐲_𝑁(𝑡)=𝑥_𝑁(𝑡)𝜆_𝑁(𝑡)00,𝐎=. (22) Applying the following collocation points in (21): 𝑡_𝑠=𝑡₀+𝑡_𝑓−𝑡₀𝑁𝑠,𝑠=0,1,…,𝑁 (23) yields to 𝑁+1 equations as follows: 𝐲_𝑁⁽¹⁾𝑡_𝑠−𝑀𝐲_𝑁𝑡_𝑠=𝑂,𝑠=0,1,…,𝑁. (24) All of the these equations can be written in the following matrix form: 𝐘⁽¹⁾−𝐌𝐘=𝑂, (25) where 𝐘⁽¹⁾=𝐲_𝑁⁽¹⁾𝑡₀𝐲_𝑁⁽¹⁾𝑡₁⋮𝐲_𝑁⁽¹⁾𝑡_𝑁,(),𝐲𝐌=𝐼⊗𝑀=kron𝐼,𝑀𝐘=^𝑁𝑡₀𝐲_𝑁𝑡₁⋮𝐲_𝑁𝑡_𝑁, (26) where ⊗ denotes the Kronecker product and 𝐼 is the identity matrix of dimension 𝑁+1.

With the aid of relation (19) and the collocation points (23), we gain 𝐲_𝑁𝑡_𝑠=𝐗𝑡_𝑠𝐀,𝐲_𝑁⁽¹⁾𝑡_𝑠=𝐗𝑡_𝑠𝐁𝐀,𝑠=0,1,…,𝑁, (27) which can be written as 𝐘=𝐗𝐀,𝐘⁽¹⁾=𝐗𝐁𝐀, (28) where 𝐗=𝐗𝑡₀𝐗𝑡₁⋮𝐗𝑡_𝑁,𝐗𝑡_𝑠=𝐗𝑡_𝑠00𝐗𝑡_𝑠,𝑠=0,1,…,𝑁. (29) If the relation (28) is substituted into (25), the fundamental matrix equation is obtained as 𝐗𝐁−𝐌𝐗𝐀=𝐎. (30) Thus, the fundamental matrix equation (30) corresponding to (7) can be written in the form 𝐖𝐀=𝐎or[𝐖;𝐎], (31) which corresponds to a linear system of 2(𝑁+1) algebraic equations in 2(𝑁+1) the unknown monomial coefficients so that 𝐖=𝐗𝐁−𝐌𝐗=𝑤_𝑝,𝑞,𝑝,𝑞=1,2,…,2(𝑁+1). (32) By the aid of the relation (19), the matrix form for the boundary conditions which are given in (7) can be written as 𝐔𝐀=𝐑,or𝐔;𝐑. (33)

Finally, by replacing the rows of the matrices [𝐔;𝐑] by the last rows of the matrices [𝐖;𝐎], we obtain the new augmented matrix 𝐖^∼𝐀=𝐎^∼. (34) The unknown monomials coefficients which exist in 𝐀 are determined by solving this linear system, and hence 𝑎_𝑖,0,𝑎_𝑖,1,…,𝑎_𝑖,𝑁, (𝑖=1,2) are substituted in (14). Therefore, we find the approximated solutions 𝑥_𝑁(𝑡)=^𝑁_𝑛=0𝑎_1,𝑛𝑡^𝑛,𝜆_𝑁(𝑡)=^𝑁_𝑛=0𝑎_2,𝑛𝑡^𝑛. (35)

We can easily check the accuracy of the method. Since the truncated monomial series (14) are the approximate solutions of (7), when the functions 𝑥_𝑁(𝑡), 𝜆_𝑁(𝑡) and their derivatives are substituted in (7), the resulting equation must be satisfied approximately; that is, for 𝑡=𝑡_𝑞∈[𝑡₀,𝑡_𝑓],𝑞=0,1,2,…𝐸₁𝑡_𝑞=𝑥⁽¹⁾𝑡_𝑞−𝐴(𝑡)−𝐵(𝑡)𝑅⁻¹𝑄^𝑇𝑥𝑡_𝑞+𝐵(𝑡)𝑅⁻¹𝐵(𝑡)^𝑇𝜆𝑡_𝑞𝐸≅0,₂𝑡_𝑞=𝜆⁽¹⁾𝑡_𝑞−−𝑃+𝑄𝑅⁻¹𝑄^𝑇𝑥𝑡_𝑞−𝑄𝑅⁻¹𝐵(𝑡)^𝑇−𝐴(𝑡)^𝑇𝜆𝑡_𝑞≅0, (36) and 𝐸_𝑖,𝑁(𝑡_𝑞)≤10^−𝑘𝑞, 𝑖=1,2 (𝑘_𝑞 positive integer).

If max10^−𝑘𝑞=10^−𝑘 (𝑘 positive integer) is prescribed, then the truncation limit 𝑁 is increased until the difference 𝐸_𝑖,𝑁(𝑡_𝑞) at each of the points becomes smaller than the prescribed 10^−𝑘, see [24].

Remark 2. We must recall that a similar approach can be applied for the state 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 1𝐽=2¹₀2𝑥²(𝑡)+𝑢²(𝑡)𝑑𝑡, (37) subject to 𝑥̇(𝑡)1=−2𝑥(𝑡)+𝑢(𝑡),𝑥(0)=1. (38) The purpose is to find the optimal control 𝑢(𝑡) which minimizes (37) subject to (38). The Optimal value of performance index for this problem is 𝐽^∗=0.463566653481105 and also exact solutions have been given in [30] as 𝑥^∗(𝑡)=12+𝑒⁻³2𝑒^−(3𝑡/2)+𝑒^{−3+(3𝑡/2)},𝑢^∗(𝑡)=12+𝑒⁻³−𝑒^−(3𝑡/2)+𝑒^{−3+(3𝑡/2)}. (39) Since the objective function of this OCP is convex, therefore the following necessary conditions (i.e., linear Hamiltonian system) for optimality are also sufficient: 𝑥̇(𝑡)1=−2𝑥(𝑡)(𝑡),𝜆̇(𝑡)(𝑡)+1−𝜆=−2𝑥2𝜆(𝑡),𝑥(0)=1,𝜆(1)=0. (40) 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 𝑢^∗(𝑡)=−𝜆(𝑡), where 𝜆(𝑡) is the solution of the previous system.
We solve this problem by using our proposed method in the cases of 𝑁=4,5,6,7, and 8. The approximated solutions corresponding to these values of 𝑁 are provided below 𝑥₄(𝑡)=0.130043𝑡⁴−0.4921311𝑡³+1.1171671𝑡²𝑢−1.4266579𝑡+1.0,₄(𝑡)=−0.1020785𝑡⁴+0.5247564𝑡³−1.0326910𝑡²𝑥+1.5366710𝑡−0.9266579,₅(𝑡)=−0.0317935𝑡⁵+0.1886526𝑡⁴−0.5272247𝑡³+1.123919𝑡²𝑢−1.427164𝑡+1.0,₅(𝑡)=0.04038275𝑡⁵−0.1772002𝑡⁴+0.5697892𝑡³−1.042226𝑡²𝑥+1.536418𝑡−0.9271638,₆(𝑡)=0.009145352𝑡⁶−0.05359579𝑡⁵+0.2075571𝑡⁴−0.5343086𝑡³+1.124911𝑡²𝑢−1.4271300𝑡+1.0,₆(𝑡)=−0.006951979𝑡⁶+0.05688585𝑡⁵−0.1914114𝑡⁴+0.5750789𝑡³−1.042907𝑡²𝑥+1.536435𝑡−0.9271298,₇(𝑡)=−0.001606958𝑡⁷+0.01380419𝑡⁶−0.05883564𝑡⁵+0.210416𝑡⁴−0.5350703𝑡³+1.124991𝑡²𝑢−1.427133𝑡+1.0,₇(𝑡)=0.002084015𝑡⁷−0.01300909𝑡⁶+0.06371979𝑡⁵−0.1951577𝑡⁴+0.5760817𝑡³−1.043019𝑡²𝑥+1.536433𝑡−0.9271335,₈(𝑡)=0.0003561777𝑡⁸−0.002824592𝑡⁷+0.01548942𝑡⁶−0.06004358𝑡⁵+0.2108896𝑡⁴−0.5351673𝑡³+1.124999𝑡²𝑢−1.427133𝑡+1.0,₈(𝑡)=−0.0002657152𝑡⁸+0.00299044𝑡⁷−0.01426032𝑡⁶+0.06461366𝑡⁵−0.1955065𝑡⁴+0.5761526𝑡³−1.043024𝑡²+1.536433𝑡−0.9271333. (41) The associated performance indexes for the selected values of 𝑁 are 𝐽₄=0.463550469, 𝐽₅=0.4635575131, 𝐽₆=0.4635665731, 𝐽₇=0.463566618, and 𝐽₈=0.4635666532. We provide the 𝑒_𝑥𝑁=max_0≤𝑡≤1|𝑥_𝑁(𝑡)−𝑥^∗(𝑡)|, 𝑒_𝐽𝑁=|𝐽_𝑁−𝐽^∗| 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 𝑁 in Table 1. It can be seen from this table that our obtained results for such considered values of 𝑁 (i.e., 4, 5, 6, 7, and 8) are the same and equal to the obtained results of [30] for higher values of 𝑁 such as 4,8,16,32, and 64 in computation of 𝑒_𝑥𝑁. Moreover, our results corresponding to the 𝑒_𝐽𝑁 are more accurate with regard to the method of [30] even by choosing lower values of 𝑁.

Table 1 

Comparison results of Example 3.

Example 4 (see [30, 31] linear Hamiltonian system). Consider the linear time-varying system 𝑥̇(𝑡)=𝑡𝑥(𝑡)+𝑢(𝑡), (42) with the cost functional 1𝐽=2¹₀𝑥²(𝑡)+𝑢²(𝑡)𝑑𝑡. (43) The problem is to obtain the optimal control 𝑢^∗(𝑡) which minimizes (43) subject to (42). The optimal control is 𝑢^∗(𝑡)=−𝐾(𝑡)𝑥^∗(𝑡), (44) where 𝐾(𝑡) is the feedback controller gain matrix and the solution of the Riccati equation [30] 𝐾(𝑡)=𝐾²(𝑡)−2𝑡𝐾(𝑡)−1,𝐾(1)=0. (45) According to the optimality conditions (5) and (6) we have ̇(𝑡)(𝑡),̇(𝑡)(𝑡),𝑥(0)(1)(𝑡)(𝑡).𝑥=𝑡𝑥−𝜆𝜆=−𝑥−𝑡𝜆=1,𝜆=0,𝑢=−𝜆 (46) We first solve the previous system and obtain the numerical solutions 𝑥_𝑁(𝑡) and 𝜆_𝑁(𝑡) for 𝑁=4,5,6, and then solve (45) by ODE solver commands which exist in MAPLE 13. Since, 𝐾(𝑡)=−(𝑢(𝑡)/𝑥(𝑡)) then our numerical results of 𝐾(𝑡) are equal to −(𝑢_𝑁(𝑡)/𝑥_𝑁(𝑡)), that is, 𝐾_𝑁(𝑡)=−(𝑢_𝑁(𝑡)/𝑥_𝑁(𝑡)). The numerical results of system (46), which are obtained by the proposed method could be deduced as 𝑥₄(𝑡)=0.20203304𝑡⁴−0.31263421𝑡³+1.0000762𝑡²𝑢−0.96700804𝑡+1.0,₄(𝑡)=0.011211016𝑡⁴−0.055485904𝑡³+0.011282926𝑡²𝑥+1.0𝑡−0.96700804,₅(𝑡)=0.021615359𝑡⁵+0.16221315𝑡⁴−0.28939672𝑡³+0.99513749𝑡²𝑢−0.96857084𝑡+1.0,₅(𝑡)=0.047811786𝑡⁵−0.078776302𝑡⁴−0.00037310025𝑡³−0.000091547293𝑡²𝑥+1.0𝑡−0.96857084,₆(𝑡)=0.043743856𝑡⁶−0.083312898𝑡⁵+0.25396668𝑡⁴−0.32429532𝑡³+1.0001826𝑡²𝑢−0.96853474𝑡+1.0,₆(𝑡)=−0.002379348𝑡⁶+0.052612074𝑡⁵−0.081954336𝑡⁴+0.00027888358𝑡³−0.000022529694𝑡²+1.0𝑡−0.96853474. (47) Also, the exact solution 𝐾(𝑡) of the Riccati equation (45) at the uniform mesh in the interval (0,1) are 𝐾(0)=9.6854𝑒−001, 𝐾(0.1)=9.5147𝑒−001, 𝐾(0.2)=9.1063𝑒−001, 𝐾(0.3)=8.4416𝑒−001, 𝐾(0.4)=7.5241𝑒−001, 𝐾(0.5)=6.3856𝑒−001, 𝐾(0.6)=5.0873𝑒−001, 𝐾(0.7)=3.7127𝑒−001, 𝐾(0.8)=2.3540𝑒−001, 𝐾(0.9)=1.0955𝑒−001, and 𝐾(1)=0. In Table 2, we provide the absolute values of errors at the selected points for the previously considered values of 𝑁 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 Example 4.

Example 5 (nonlinear Hamiltonian system). As our third illustration, consider the following nonlinear optimal control problem: min𝐽=¹₀𝑢²(𝑡)̇1𝑑𝑡s.t.𝑥=2𝑥²(𝑡)(𝑡)(𝑡)[]𝑥(0)(1)sin𝑥+𝑢,𝑡∈0,1=0,𝑥=0.5. (48) Trivially 𝑓(𝑡,𝑥(𝑡))=(1/2)𝑥²(𝑡)sin𝑥(𝑡), 𝑔(𝑡,𝑥(𝑡))=1, 𝑄=0, 𝑅=1, 𝑡₀=0, and 𝑡_𝑓=1. As mentioned in Section 2.2, we solve the following system of ordinary differential equations: ̇1𝑥=2𝑥²(𝑡)(𝑡)−1sin𝑥2𝜆(𝑡),𝜆̇=−𝜆(𝑡)𝑥(𝑡)sin𝑥(𝑡)−12𝜆(𝑡)𝑥²(𝑡)(𝑡),[],𝑥(0)(1)(1)cos𝑥𝑡∈0,1=0,𝑥=0.5,𝜆=0. (49) Also the optimal control law is given by 𝑢^∗(𝑡)1=−2𝜆(𝑡). (50)
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 Section 3, 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 𝑥(𝑡) and 𝜆(𝑡) transforms the basic problem to the corresponding system of nonlinear algebraic equations. By assuming different values of 𝑁 such as 5,  7, and 9, we solve the previously mentioned system. In Table 3, we provide the approximated performance index 𝐽_𝑁, which is obtained by our proposed method and also the difference between the approximated 𝑥_𝑁(1)−𝑥(1)=𝑥_𝑁(1)−(1/2) for the considered values of 𝑁.

Table 3 

Numerical results of Example 5.

Example 6 (see [11, 32] state-space analysis). We consider a linear-time invariant state equation 𝑥̇₁(𝑡)𝑥̇₂(𝑡)𝑥=𝐴¹(𝑡)𝑥₂(𝑡)+𝑏𝑢(𝑡),𝑡∈(0,1), (51) where 11𝐴=−312−4,𝑏=. (52) We are given that the input 𝑢(𝑡) is the unit step function in the interval (0,1) and the initial state is 𝑥₀=𝑥₁(0)𝑥₂(0)=11. (53) The exact solution for (51) is 𝑥₁(𝑡)=𝑥₂(𝑡)=1+𝑒^−2𝑡2. (54) We solve this problem by using our proposed method in the cases of 𝑁=4,5,6,7,8,9,10,12, and 20. The approximated solutions corresponding to the 𝑁=4,5,6,7,8,9, and 10 are provided later 𝑥_1,4(𝑡)=𝑥_2,4(𝑡)=0.1658𝑡⁴−0.58031𝑡³+0.98446𝑡²𝑥−1.0𝑡+1.0,_1,5(𝑡)=𝑥_2,5(𝑡)=−0.06235𝑡⁵+0.28058𝑡⁴−0.64844𝑡³+0.99761𝑡²𝑥−1.0𝑡+1.0,_1,6(𝑡)=𝑥_2,6(𝑡)=0.019951𝑡⁶−0.10973𝑡⁵+0.32143𝑡⁴−0.66365𝑡³+0.99969𝑡²𝑥−1.0𝑡+1.0,_1,7(𝑡)=𝑥_2,7(𝑡)=−0.0055374𝑡⁷+0.035993𝑡⁶−0.12776𝑡⁵+0.33126𝑡⁴−0.66626𝑡³+0.99997𝑡²𝑥−1.0𝑡+1.0,_1,8(𝑡)=𝑥_2,8(𝑡)=0.0013547𝑡⁸−0.010161𝑡⁷+0.042378𝑡⁶−0.13232𝑡⁵+0.33304𝑡⁴−0.66662𝑡³+1.0𝑡²𝑥−1.0𝑡+1.0,_1,9(𝑡)=𝑥_2,9(𝑡)=−0.00029605𝑡⁹+0.0025165𝑡⁸−0.012061𝑡⁷+0.044057𝑡⁶−0.13318𝑡⁵+0.3333𝑡⁴−0.66666𝑡³+1.0𝑡²𝑥−1.0𝑡+1.0,_1,10(𝑡)=𝑥_2,10(𝑡)=0.000058425𝑡¹⁰−0.00055504𝑡⁹+0.003006𝑡⁸−0.012576𝑡⁷+0.044386𝑡⁶−0.13331𝑡⁵+0.33333𝑡⁴−0.66667𝑡³+1.0𝑡²−1.0𝑡+1.0. (55)
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 𝑁=8,12, and 20 is considered in Table 4. Moreover, the error histories in the computational interval (0,1) associated with our method for 𝑁=4 and 5 together with the error history of the method [11] (Bessel Integration) for 𝑁=4 and also the error history of the method [32] (Laguerre Integration) for 𝑁=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 Example 6.

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

Conflict of Interests

The authors declare that they do not have any conflict of interests in their submitted paper.

Acknowledgment

The authors are grateful to the comments of the referees, which have led to several improvements in the presentation of this paper.