Solution of Infinite Horizon Nonlinear Optimal Control Problems by Piecewise Adomian Decomposition Method

In this paper, a Piecewise Adomian Decomposition Method (PADM) is used to obtain the analytical approximate solution for a class of infinite horizon nonlinear optimal control problems (OCPs). The method is a new modification of the standard ADM, in which it is treated as an algorithm in a sequence of small intervals (i.e. with small time step) for finding accurate approximate solutions to the corresponding OCPs. Applying the PADM, the nonlinear two-point boundary value problem (TPBVP), derived from the application of Pontryagin’s maximum principle (PMP), is transformed into a sequence of linear time-invariant TPBVP’s. Through the finite iterations of algorithm, a suboptimal control law is obtained for the nonlinear optimal control problem. Comparing the methodology with some known techniques shows that the present approach is powerful and reliable. It is remarkable accuracy properties are finally demonstrated by two examples.


Introduction
The investigation of the optimal control is of great importance in modern control theory. One of the most active subjects in control theory is the optimal control which has a wide range of practical applications not only in all areas of physics but also in economy, aerospace, chemical engineering, robotic, etc. (see [6,17,20]). The theory and the application of optimal control for linear time-invariant systems have been developed perfectly. However, as for nonlinear systems, synthesis problems that are solved by classic control theory lead to difficult computations. It is well-known that the nonlinear optimal control problem can be reduced to a Two-Point Boundary Value Problem (TPBVP), implementing the Pontryagin's Maximum Principle (PMP). In general, this TPBVP cannot be solved exactly and most researches have been devoted to find an approximate solution, for nonlinear TPBVP's, see [4].
In order to find an analytical approximate solution, some methods have been proposed. One of these approaches is Successive Approximation Approach (SAA) which designs an suboptimal controller for a class of nonlinear systems with a quadratic performance index. In this approach, a sequence of nonhomogeneous linear time-varying TPBVP's is solved to produce a finitestep iteration of the nonlinear compensation sequence obtaining the suboptimal control law [20].
However, SAA needs to solve a linear time-varying TPBVP's which cannot be solved easily and thus, reduces the efficiency of this method. Recently, in [9], a novel method which implements Modal series to solve a class of nonlinear OCP's with quadratic performance index, has been proposed. This method which requires solving a sequence of linear time-invariant TPBVP's, has less efficiency for large scale problems. Other methods for solving the nonlinear TPBVP obtained from PMP, such as the Homotopy perturbation method and the Homotopy analisys methods are available in [5,7,8,12,13,15,16,22,23].
On the other hand, in the context of numerical analysis, the Adomian Decomposition Method (ADM) which was proposed originally by Adomian [1], has been proved by many authors to be a powerful mathematical tool for various kinds of linear and nonlinear ODE's or PDE's. Unlike the traditional numerical methods, ADM needs no discretization, linearization, transformation or perturbation. The method, has been widely applied to solve nonlinear problems, and different modifications are suggested to overcome the demerits arising in the solution procedure (see for example [1,2,10,18,19]).
In this work, we present a improved method known as PADM to solve a class of infinite horizon nonlinear OCPs. Using the PADM, we obtain the optimal control law and a suboptimal control law. Finally, simulation examples are employed to test the validity of the PADM.
The rest of the paper is organized as follows. In Section 2, the ADM are introduced. The nonlinear infinite horizon quadratic OCP and optimality conditions are presented in Section 3, also we give a brief description of PADM in this section. The numerical results are given in Section 4. Section 5 ends this paper with a brief conclusion.
In [1], G. Adomian developed a decomposition method for solving nonlinear (stochastic) differential equations using special polynomials A n , usually called Adomian polynomials. The A n 's are generated for each nonlinearity. Given a partial differential equation the application of the Adomian decomposition method separates the equation into linear and nonlinear parts. A solution is obtained in the form of a series whose terms are determined by a recursive relationship using the Adomian polynomials. The advantages of the Adomian's decomposition method are emphasized by many authors.
The usual numerical procedures such as finite differences, Galerkin method, Finite element method, etc., used to solve or find approximate solutions to nonlinear partial/ordinary differential equations, linearize the equation (or equations) or assume that the nonlinearities are relatively small, transforming the physical problem into a purely mathematical one with an available solution. Therefore, the use of these methods may change the real solution of the mathematical model which represents the physical reality. Generally, the numerical methods are based on discretization techniques, and allow only to calculate the approximate solutions for some values of time and space variables. Therefore, to find the value (or approximate value) in other points that do not belong to the grid, it is necessary to use interpolation. It is well known that the use of interpolation has several disadvantage since an error is related with the interpolation process and we could not overlooking some phenomenon related to the problem.
The Adomian's decomposition has many advantages: it does not require any kind of discretization, linearization or perturbation of the variables and of the equation, therefore it does not need any modification of the actual model that could change the solution; is efficient on providing an approximate or even exact solution in a closed form, to linear and nonlinear problems; provides a fast and accurate convergent series and therefore it is only necessary to calculate a few terms of the series in order to obtain a reliable approximate solution; the method depends only on the known function u 0 (x) and the algorithm is of simple implementation.
To illustrate the basic idea of this method, let us consider the equation with prescribed initial conditions, where u(t) is the unknown function, is a linear operator which is assumed to be invertible, L −1 exist, R is another linear differential operator, N (u) represents the nonlinear terms, and g is the source term. Applying the inverse operator L −1 to both sides of (2.1), and using the given conditions we obtain where the function f (t) represents the terms arising from integrating the source term g and from using the given initial or boundary conditions, all are assumed to be prescribed.
The standard Adomian method defines the solution u(t) of (2.1) as a series The decomposition method suggests that the zeroth component u 0 (t) is usually identified by the function f (t) defined above. The components u 0 , u 1 , u 2 , . . . are determined recursively, therefore the solution u(t) in a series form defined by (2.3) is readily determined. The closed form for the solution u(x), if exists, can be immediately obtained because of the rapid convergence presented by the method.
The last term in (2.3) can be computed by substituting A n (u 0 , u 1 , u 2 , . . . , u n ), (2.4) where A n are specially generated Adomian polynomials for the specific nonlinearity. They depend only on the u 0 , to u n components and form a rapidly convergent series. The A n , are given as , n = 1, 2, 3, . . . (2.5) or, Algorithms for formulating Adomian polynomials were investigated in [2]. Putting (2.3) and (2.4) Each term of series (2.3) is given by the recurrent relation For later numerical computation, we use the expression Φ n (t) = n−1 k=0 u k (t) to denote the n-term approximation to u(t). Now, we briefly describe how to apply the Adomian decomposition method to systems of ordinary differential equations. Let us consider a system of ordinary differential equations in the form . . .
Solution of Infinite Horizon Optimal Control Problems by Piecewise ADM 547 Applying the inverse operator L −1 to (2.6) we obtain the following canonical form Thus, by the Adomian decomposition method, each component of the solution of (2.6) can be expressed as a series of the form u j = ∞ i=0 f i,j and the integrands on the right side of (2.7), using (2.5), are expressed as where the A i,j are the Adomian polynomials corresponding to the nonlinear part f i . We should note that in order to solve system (2.6), we obtain a system of Volterra integral equations of the second kind, (2.7).

Nonlinear Infinite Horizon Quadratic OCP and Optimality Conditions
Consider an infinite horizon nonlinear OCP described bẏ with x(t) ∈ R n denoting the state variable, u(t) ∈ R m the control variable and x 0 the given initial state at t 0 , respectively. Moreover, F : R n → R n is a nonlinear analytic vector field, where F (0) = 0 (and hence x = 0 is an equilibrium point of the system), B is a constant matrix of appropriate dimension. Our aim is to minimize the quadratic objective functional subject to the nonlinear system (3.1), for Q ∈ R n×n , R ∈ R m×m positive semidefinite and positive definite matrices, respectively. Also, it is assumed that the pair (J F (0), B) is controllable and the pair (J F (0), Q 1 2 ) is observable, where J F is the Jacobian matrix of F , and Q 1 2 is the square root of matrix Q. These assumptions ensure the existence of a smooth optimal solution on a certain region of state space containing the equilibrium point, see [14].
According to Pontryagin's maximum principle, the optimality conditions are obtained as the following nonlinear TPBVP: . Also, the optimal control law is obtained by Remark 1 [Approximating the co-state initial guess]. It is necessary to transform the BVP in a IVP. Therefore, we must find an α ∈ R such that the condition λ(∞) = 0 can be replaced by the condition λ(t 0 ) = α. By considering this condition, and since the approximations of these methods are functions of both t and α, we have n k=0 λ k (∞, α) = 0. That is, α should be a real root of the equation Proof. Let us consider the solution of TPBVP (3.3) as (x(·), λ(·)). Because the initial value of λ is not known. Thus we rewrite the TPBVP in (3.3) as following:ẋ where α ∈ R is an unknown parameter. Since the nonlinear terms in Eq. (3.5) are analytic, (x(·), λ(·)), as the solution of IVP (3.5), is analytic with respect to x 0 , [3]. Thus, (x(·), λ(·)) as the solution of TPBVP (3.3) is analytic with respect to x 0 , and the proof is complete.

Piecewise Adomian decomposition method
In order to find the solution of (3.5), we will use piecewise Adomian decomposition method (PADM). The piecewise strategy is used to provided a steady state solution in whole time horizon rather than traditional ADM. The preference of the method lies in the fact that the piecewise strategy gives this relation in an arbitrary longtime interval, while the ADM gives the optimal solution just only in the neighborhood of the initial time. In order to overcome the difficulty, we modify the ADM in the following. We assume that equation (3.5) is defined on the time interval I = [t 0 , t f ]. We divide the time interval into M equal subintervals Now we employ the piecewise Adomian decomposition method for solving Equation (3.5) in each subinterval. In order to carry out the iterations in every subinterval of equal length ∆t, we consider the approximate solution of Equation (3.5) in the a sequence of intervals, which are subject to continuity conditions at the end points of each interval, and choose the corresponding initial approximations. But, in general, we do not have these information at our clearance except at the initial point t 0 . Let u 1 (t) be solution of Equation (3.5) in the subinterval I 1 . For 2 ≤ i ≤ M, u i (t) is solution of Equation (3.5) in the subinterval I i with initial conditions by obtaining the initial conditions from the interval I i−1 , The solution of Equation (3.5) for t ∈ [t 0 , t f ] is given by (3.8)

Suboptimal control design
Consider the OCP of nonlinear system (3.1) with quadratic cost function (3.2). Then, the N th order suboptimal trajectory-control pair is obtained as follows: (3.9) The integer N in (3.9) is generally determined according to a concrete control precision. As we will present in the next subsection, every time x i (t) and λ i (t) are obtained from the linear TPBVP sequence, we let N = i and calculate x (N ) (t) and λ (N ) (t) from (3.9). Then the following quadratic performance index (QPI) can be calculated:

Solution of Infinite Horizon Optimal Control Problems by Piecewise ADM 551
The N th order suboptimal trajectory-control pair in (3.9) has desirable accuracy if for a given positive constant > 0, the following condition holds: If the tolerance error bounds > 0 is small enough, the N th order suboptimal control law will be very close to the optimal control law u * (t), the value of QPI in (3.9) will be very close to its optimal value J * , and the boundary condition will be satisfied tightly.

Numerical Experiment
In this section, we consider two illustrative examples to illustrate the simplicity and efficiency of the proposed method. The codes are developed using symbolic computation software MATLAB 7.8 and the calculations are implemented on a machine with Intel Core 2 Due Processor 2.53 Ghz and 4 GB RAM to solve nonlinear quadratic OCP, that reduced to a TPBVP.
Test problem 4.1. As first numerical test problem we consider the optimal maneuvers of a rigid asymmetric spacecraft, [11]. Euler's equations for the angular velocities of the spacecraft are given bẏ where x 1 , x 2 , and x 3 are angular velocities of the spacecraft, u 1 , u 2 , and u 3 are control torques, I 1 = 86.24, I 2 = 85.07, and I 3 = 113.59 kg m 2 are the spacecraft principle inertia.
All the nonlinearities in (4.2) are of the form Φ(u, υ) = uυ. Thus, using (2.5) we have for Φ(u, υ) In order to apply the PADM, it is necessary to integrate (4.2) in order to the variable t in each I j = [t j−1 , t j [. Therefore, we obtain the following recurrence scheme: Solution of Infinite Horizon Optimal Control Problems by Piecewise ADM 553   We compared the results of PADM with the solutions of obtained using the standard ADM. In fact, Figs. 1, 2 confirm, the ADM is divergent in [0, 1000], even for large number of iterations. Therefore, it is reasonable to use the PADM instead. The following figures (Figs. [3][4][5] present the proposed solution obtained by the PADM, the modal series method, see [9], and the solution obtained by the MatLab package bvp4c.
In order to obtain an accurate enough suboptimal trajectory-control pair, we applied the strategy proposed in section 3, with the tolerance error bounds = 5 × 10 −5 . In this case, convergence is achieved after 2 iterations, i.e.
A minimum of J (2) = 0.005432169453226 is obtained. This confirms that the proposed method yields excellent results.
Test problem 4.2. Consider the two-order nonlinear composite system, see  [21], described by this is,ẋ

5)
x 2 (t) = −x 2 (t) + u 2 (t) + x 1 (t)x 2 (t) + x 3 2 (t), x 1 (0) = 0, x 2 (0) = 0.8. The quadratic cost functional to be minimized is given by since we considered Q = R = 1 0 0 1 . According to the optimal control theory (3.3), the optimality conditions can be written as  where α i ∈ R are unknown parameter. Also the optimal control law is given by  The Adomian polynomials for all the nonlinearities in (4.6) can be obtained from Φ(u, υ) = uυ, mentioned previously, or Γ (u, υ, z) = uυz whose Adomian polynomials are In order to obtain Γ 1 = u 3 it is only necessary to consider u = υ = z in (4.8). In order to obtain the nonlinearity for Γ 2 = u 2 z it is only necessary to consider u = υ in (4.8) and in order to obtain the Adomian polynomials for Ψ (u) = u 2 it is only necessary to consider u = υ in (4.3).
A minimum of J (3) = 0.5350 is obtained.
In this example, we considered a much small tolerance error. With one more iteration, the relative error to the J (n) is grater than the one in the previous example. This is due to the complexity of this example. This shows that the complexity of problem can influence the rate of convergence to the method.
The following figures (Figs. 6-9) present the proposed solution to problem obtained by the PADM, the SAA method and solution obtained by the MatLab package bvp4c.

Conclusions
In this paper, we have successfully developed ADM for a class of infinite horizon nonlinear optimal control problems (OCPs). Then we employed PADM to finding optimal control of infinite horizon nonlinear systems. This method can solve the TPBVP obtained from PMP. The results obtained by comparison with other methods justify the use of this procedure to obtain solutions to infinite horizon nonlinear OCP's problems.