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

and Applied Analysis 3 subject to the nonlinear system (8), for Q ∈ R, R ∈ Rm×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]: ?̇? = f (t, x) + g (t, x) u ∗ , ̇ λ = −H x (x, u ∗ , λ) , u ∗ = arg min u H(x, u, λ) , x (t 0 ) = x 0 , λ (t f ) = 0, (10) where H(x, u, λ) = (1/2)[xTQx + uTRu] + λT[f(t, x) + g(t, x)u] is referred to theHamiltonian. Equivalently, (10) can be written in the form of ?̇? = f (t, x) + g (t, x) [−R −1 g T (t, x) λ] ̇ λ = − (Qx + ( ∂f (t, x)


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 -space analysis and optimal control.The key idea of these methods is based on the integral expression where Φ() = [Φ 1 (), Φ 2 (), . . ., Φ  ()] 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][2][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 constantfeedback-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][13][14][15][16][17][18][19][20][21][22][23][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 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.

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.
where  ∈ R  is the costate vector.

State Space Analysis Problems.
In this part, we consider the following state space analysis problem: where () ∈ R × , () ∈ R  , and () ∈ R are known, meanwhile () ∈ R  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.

Fundamental Matrix Relations and Method of the Solution
In this section, by using the collocation points and the matrix relations between the monomials {1, ,  2 , . . .,   } and their derivatives, we will find the approximate solution of the system (7) expressed in the truncated monomial series form (assuming that () and () ∈ R and also () together with () is independent of time , that is,  = () and  = ()) 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 where The matrix form of the relation between the matrix X() and its th derivative X () () is By using the relations ( 15) and ( 16), we have the following relations: Thus, we can express the matrices y() and y (1) () as follows: where Now, we can restate the system (7) in the matrix form where y (1)   () = [ [  (1)   () Applying the following collocation points in (21): yields to  + 1 equations as follows: All of the these equations can be written in the following matrix form: where 1)   ( 0 ) y (1)   ( 1 ) . . .
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 y  (  ) = X (  ) A, y (1)   (  ) = X (  ) BA,  = 0, 1, . . ., , which can be written as where If the relation ( 28) is substituted into (25), the fundamental matrix equation is obtained as Thus, the fundamental matrix equation (30) corresponding to (7) can be written in the form which corresponds to a linear system of 2( + 1) algebraic equations in 2( + 1) the unknown monomial coefficients so that By the aid of the relation (19), the matrix form for the boundary conditions which are given in ( 7) can be written as 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 The unknown monomials coefficients which exist in A are determined by solving this linear system, and hence  ,0 ,  ,1 , . . .,  , , ( = 1, 2) are substituted in (14).Therefore, we find the approximated solutions 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 ,   ],  = 0, 1, 2, . . . 1 (  ) =       (1) and  , (  ) ≤ 10 −  ,  = 1, 2 (  positive integer).If max 10 −  = 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].

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 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 Since the objective function of this OCP is convex, therefore the following necessary conditions (i.e., linear Hamiltonian system) for optimality are also sufficient: 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 The associated performance indexes for the selected values of  are  4 = 0.463550469,  5 = 0.4635575131,  6 = 0.4635665731,  7 = 0.463566618, and  8 = 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 .
Example 5 (nonlinear Hamiltonian system).As our third illustration, consider the following nonlinear optimal control problem: = 0.5.
Trivially (, ()) = (1/2) 2 () sin (), (, ()) = 1,  = 0,  = 1,  0 = 0, and   = 1.As mentioned in Section 2.2, we solve the following system of ordinary differential equations: Also the optimal control law is given by 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 .
Example 6 (see [11,32] state-space analysis).We consider a linear-time invariant state equation where We are given that the input () is the unit step function in the interval (0, 1) and the initial state is The exact solution for (51) is 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

𝐸(𝑡)
Integration Laguerre ( = 5) Presented method ( = 4) Presented method ( = 5) Figure 1: The error histories of our method together with the methods [11,32] in Example 6.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].

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.

Table 1 :
Comparison results of Example 3.

Table 2 :
The error histories at the selected points for Example 4.

Table 3 :
Numerical results of Example 5.

Table 4 :
The error histories at the selected points for Example 6.