A Laguerre spectral method for quadratic optimal control of nonlinear systems in a semi-infinite interval

This paper presents a Laguerre homotopy method for quadratic optimal control problems in semi-infinite intervals (LaHOC), with particular interests given to nonlinear interconnected large-scale dynamic systems. In LaHOC, the spectral homotopy analysis method is used to derive an iterative solver for the nonlinear two-point boundary value problem derived from Pontryagin's maximum principle. A proof of local convergence of the LaHOC is provided. Numerical comparisons are made between the LaHOC, Matlab BVP5C generated results and results from the literature for two nonlinear optimal control problems. The results show that LaHOC is superior in both accuracy and efficiency.


Introduction
Large-scale systems are found in many practical applications, such as power systems and physical plants. During the past several years, the problem of analysis and synthesis for dynamic large-scale systems has received considerable attention. Based on the characteristics of large-scale systems, many results have been proposed, such as modelling, stability, robust control, decentralized, and so on [1][2][3][4][5][6].
The optimal control of nonlinear large-scale systems has been widely investigated in recent decades. For instance, a new successive approximation approach (SAA) was proposed in [7]. In this approach, instead of directly solving the nonlinear large-scale two-point boundary value problem (TPBVP), derived from the maximum principle, a sequence of non-homogeneous linear time-varying TPBVPs is solved iteratively. Also, in [8] a new technique, called the modal series method, has been extended to solve a class of infinite horizon OCPs of nonlinear interconnected large-scale dynamic systems, where the cost function is assumed to be quadratic and decoupled. This method provides the solution of autonomous nonlinear systems in terms of fundamental and interacting modes. Conventional methods of optimal control are generally impractical for many nonlinear large-scale systems because of the dimensionality problem and high complexity in calculations. One example is the state-dependent Riccati equation (SDRE) method [9]. Although this scheme has been widely used in many applications, its major limitation is that it needs to solve a sequence of matrix Riccati algebraic equations at each sample state along the trajectory. This property may take a long computing time and large memory space. Therefore, developing new methods is necessary for solving nonlinear large-scale optimal control problems (OCPs) [10].
The use of spectral methods for optimal control problems usually leads to a more efficient method than finite element or finite difference approaches. Chebyshev's and Legendre's methods are commonly used for problems in finite intervals [11,12]. For infinite or semi-infinite intervals, there are several choices for the approximation bases: Hermite polynomials/functions [13], Laguerre polynomials/functions [14], mapped Jacobi bases [15][16][17]. Furthermore, one class of very important applications of OCP in unbounded intervals is the so-called Minimum Action Method (MAM) [18] used in finding the most probable transition path in phase transition phenomena. Using MAM to study spatial extended transitions, such as fluid instability transition, is usually equivalent to solve a large-scale nonlinear optimal control problem [18].
The homotopy analysis method is an analytical technique for solving nonlinear differential equations. The HAM [19,20] was first proposed by Liao in 1992 to solve lots of nonlinear problems. This method has been successfully applied to many nonlinear problems, such as physical models with an infinite number of singularities [21], nonlinear eigenvalue problems [22], fractional Sturm-Liouville problems [23], optimal control problems [24,25], Cahn-Hilliard initial value problem [26], semi-linear elliptic boundary value problems [27] and so on [28]. The HAM contains a certain auxiliary parameter which provides us with a simple way to adjust and control the convergence region and rate of convergence of the series solution. Moreover, by means of the so-called -curve, it is easy to determine the valid regions of to gain a convergent series solution. The HAM, however, suffers from a number of restrictive measures, such as the requirement that the solution sought ought to conform to the so-called rule of solution expression and the rule of coefficient ergodicity. These HAM requirements are meant to ensure that the implementation of the method results in a series of differential equations can be solved analytically.
Recently, Motsa et al. [29][30][31][32] proposed a spectral modification of the homotopy analysis method, the spectral-homotopy analysis method (SHAM). The SHAM approach imports some of the ideas of the HAM such as the use of the convergence controlling auxiliary parameter. In the implementation of the SHAM, the sequence of the so-called deformation differential equations is converted into a matrix system by applying the Chebyshev or Legendre pseudospectral method [31]. But so far, to our knowledge, there is no work concerning the combination of Laguerre polynomials [33] with the HAM. This paper presents a spectral homotopy analysis method based on modified Laguerre-Radau interpolation to solve nonlinear large-scale optimal control problems (OCPs). This process has several advantages. First, it possesses spectral accuracy [34,35]. Next, it is easier to be implemented, especially for nonlinear systems. Furthermore, it is applicable to long-time calculations.
The paper is organized as follows. The nonlinear interconnected OCP and optimality conditions are described in Section 2. In Section 3, we propose the new algorithm by using the modified Laguerre polynomials. The convergence of the proposed method is proved in Section 4. We present the numerical results in Section 5, which demonstrate the spectral accuracy of the proposed methods. The final section is for concluding remarks.

The nonlinear interconnected OCP
Consider a nonlinear interconnected large-scale dynamic system which can be decomposed into N interconnected subsystems. The ith subsystem for i = 1, 2, . . . , N is described bẏ with x i ∈ R n i denoting the state vector, u i ∈ R m i the control vector of the ith subsystem, respectively, x = is a nonlinear analytic vector function where F i (0) = 0, and x i 0 ∈ R n i is the initial state vector. Also, A i and B i are constant matrices of appropriate dimensions such that the pair (A i , B i ) is completely controllable [8]. Furthermore, the infinite horizon quadratic cost function to be minimized is given by where Q i ∈ R n i ×n i and R i ∈ R m i ×m i are positive semidefinite and positive definite matrices, respectively. Note that quadratic cost function (2) is assumed to be decoupled as a superposition of the cost functions of the subsystems. According to Pontryagin's maximum principle, the optimality conditions are obtained as the following nonlinear TPBVP: where λ i (t) ∈ R n i is the co-state vector, λ = (λ T 1 , λ T 2 , . . . λ T N ) T , and i (x(t), λ(t)) = N j=1 (∂f j (x(t))/ ∂x i (t))λ j (t). Also the optimal control law of the ith subsystem is given by Unfortunately, problem (3) is a nonlinear large-scale TPBVP which is decomposed into N interconnected subproblems. In general, it is extremely difficult to solve this problem analytically or even numerically, except in a few simple cases. In order to overcome this difficulty, we will present the LaHOC method in the next section.

Laguerre polynomials and spectral homotopy analysis method
In this section, we give a brief description of the basic idea of the Laguerre homotopy method for solving nonlinear boundary value problems. At first, we take into account the following properties of the modified Laguerre polynomials.

Properties of the modified Laguerre polynomials
Let ω β (t) = e −βt , β > 0, and define the weighted space L 2 ω β (0, ∞) as usual, with the following inner product and norm [14]: The modified Laguerre polynomial of degree l is defined by They satisfy the recurrence relation The set of Laguerre polynomials is a complete L 2 ω β (0, ∞)-orthogonal system, namely, where δ l,m is the Kronecker symbol. Thus, for any v ∈ where the coefficientsv l are given bŷ Now, let N be any positive integer, and P N (0, ∞) the set of all algebraic polynomials of degree at most N. We denote by t N β,j , 0 ≤ j ≤ N, the nodes of modified Laguerre-Radau interpolation. Indeed, t N β,0 = 0 and t N β,j , 1 ≤ j ≤ N, are the distinct zeros of (d/dt)L β N+1 (t). By using (7), the corresponding Christoffel numbers are as follows: , For Next, we define the following discrete inner product and norm, For any , ψ ∈ P N (0, ∞),

Spectral homotopy analysis method
In this section, we give a description of the SHAM with the Laguerre polynomials basis. This will be followed by a description of the new version of the SHAM algorithm [29]. To this end, consider a general n-dimensional initial value problem described aṡ We make the usual assumption that f is sufficiently smooth for linearization techniques to be valid. If z = (z 1 , z 2 , . . . , z n ), we can apply the SHAM by rewriting Equation (15) aṡ subject to the initial conditions where z 0 r are the given initial conditions, σ r,k are known constant parameters and g r is the nonlinear component of the rth equation.
The SHAM approach imports the conventional ideas of the standard homotopy analysis method by defining the following zeroth-order deformation equations: (19) where q ∈ [0, 1] is an embedding parameter,z r (t; q) are unknown functions, and r is a convergence controlling parameter. The operators L r and N r are defined as Using the ideas of the standard HAM approach [20], we differentiate zeroth-order Equations (19) m times with respect to q and then set q = 0 and finally divide the resulting equations by m! to obtain the following equations, which are referred to as the mth order (or higher order) deformation equations: subject to where and After obtaining solutions for Equation (22), the approximate solution for each z r (t) is determined as the series solution A HAM solution is said to be of order M if the above series is truncated at m = M, that is, if A suitable initial guess to start off the SHAM algorithm is obtained by solving the linear part of (17) subject to the given initial conditions, that is, we solve If Equation (28) cannot be solved exactly, the spectral collocation method is used as a means of solution. The solution z r,0 (t) of Equation (28) is then fed to (22) which is iteratively solved for z r,m (t) (for m = 1, 2, 3, . . . , M ).
In this paper, we use the Laguerre pseudo-spectral method to solve Equations (22)- (24). The pseudospectral derivative D N (z) of a continuous function z is defined by that is, D N (z) is the derivative of the interpolating polynomial of z. Moreover, D N can be expressed in terms of a matrix, the pseudo-spectral derivation matrix D β : Indeed, given the nodes {x N of an unknown function and {(h β ) j }, the Lagrange interpolation polynomials associated with the points x j , differentiating m times the expression We now state two important results. The first ensures that it is sufficient to compute the first-order differentiation matrix, and the second gives the general expression of its entries.

Lemma 3.1 ([33]):
Next we have: j } N j=0 have the following form: Applying the Laguerre spectral collocation method in Equations (22)-(24) gives where R m−1 is an (N + 1)n × 1 vector corresponding to R r,m−1 when evaluated at the collocation points and The matrix A is an (N + 1)n × (N + 1)n matrix that is derived from transforming the linear operator L r using the derivative matrix D β (we omit subscript β for simplicity) and is defined as where I is an identity matrix of order N + 1. Thus, starting from the initial approximation, recurrence formula (33) can be used to obtain the solution z r (t).

Convergence analysis of LaHOC
To analyse the convergence of LaHOC, we first recall the mth order (or higher order) deformation equation subject to the initial condition where H(t) = 0 is an auxiliary function, where z r,m , L r and N r in (22) are the rth components of z m−1 and operators L and N , respectively. Let us define the nonlinear operator N and the sequence {Z m } ∞ m=0 as Therefore, we have from (39) we have subject to the initial condition Consequently, the collocation method is based on a solution Z N (t) ∈ P N+1 (0, ∞), for (41) such that subject to the initial condition From (43), we have Theorem 4.1: Assume that for any k = 0, 1, . . . , N, for some constant L f > 0. Then for any initial n-vector Z N 0 (t N β,k ), Z k converges to someẐ(t N β,k ) which is the exact solution of (17), at any GLR point, t N β,k , if Proof: 1. Using (5) and integrating by parts yield that then we have by (51) and from the Cauchy inequality, we obtain that from where 2. Taking the discrete weighted inner product of (47) withZ N m (t N β,k ), we have Therefore, a combination with Cauchy inequality and (51) leads to Then by using inverse inequality of Laguerre polynomial and (48), we get Hence, we have Then for any m ≥ m ≥ 1, Since γ ∈ [0, 1), Z N m − Z N m ω β → 0 as m, m → ∞. Thus Z k is a Cauchy sequence, and since R n is a Banach space, Z k has a limitẐ(t N β,k ). Taking limit m → ∞ in (43) yields Thus,Ẑ(t N β,k ) is the exact solution of (17) at any GLR point t N β,k . Also, by noticing the definition ofẐ N (t), it is easy to verifyẐ N (t N β,k ) =Ẑ(t N β,k ), and hence, the proof is completed.

Numerical experiments
To demonstrate the applicability of the LaHOC algorithm as an appropriate tool for solving infinite horizon optimal control for nonlinear large-scale dynamical systems, we apply the proposed algorithm to several test problems.
Test problem 3.1. Consider the two-order nonlinear composite system described by [7]: The quadratic cost functional to be minimized is given by In this example, we have Then, according to optimal control theory (3), the optimality conditions can be written aṡ Also the optimal control laws are u 1 (t) = −λ 1 , u 2 (t) = −λ 2 . In this example, the parameters used in the LaHOC algorithms are With these definitions, the LaHOC algorithm gives Because the right-hand side of Equation (73) is known, the solution can easily be obtained by using methods for solving a linear system of equations. Table 1 gives a comparison between the present LaHOC results for N = 100 and = −0.6 and the numerically generated BVP5C, at selected values of time t. It can be seen from the table that there is good agreement between the two results. Moreover, our calculations show the better accuracy of LaHOC. In comparison with the BVP5C, it is noteworthy that the LaHOC controls the error bounds while preserving the CPU time. The CPU time of LaHOC is 0.606532 s, and BVP5C is 1.109817 s. Figure 1 and 2 show the suboptimal states and control for m = 20 iterations of LaHOC, compared to MATLAB built-in function BVP5C. The convergence of LaHOC iteration is depicted in Figure 3. Also, Figure 4 presents that the minimum objective functional |J j − J jj |, j = 1, 2, . . . , 11 converges to J jj , where j = 20 : 10 : 120 and jj = 11. Table 1. Comparison between the LaHOC solution when N = 100 and = −0.6 and BVP4C solution.  The results obtained with the present method are in good agreement with the results of the successive approximation method used by Tang and Sun [7].
Because the right-hand side of Equation (92) is known, the solution can easily be obtained by using methods for solving linear systems of equations. Tables 2 and 3 give a comparison between the present SHAM results for N = 50 and = −1 and the numerically generated BVP5C at selected values of time t. It can be seen from the tables that there is good agreement between the two results. Moreover, our calculations show the better accuracy of LaHOC. In comparison with the BVP5C, it is noteworthy that the LaHOC controls the error bounds while preserving the CPU time. The CPU time of LaHOC is 1.009860 s, and BVP5C is 4.514071 s. Figures 5-9 show the suboptimal states and control for m = 20 iterations of LaHOC compared to MATLAB the built-in function BVP5C. The convergence of the LaHOC iteration is depicted in Figure 10.
The obtained optimal trajectories and optimal controls are identical to those obtained by Jajarmi et al. [8].

Conclusion
In this paper, an effective method based upon the spectral homotopy method with Laguerre basis (LaHOC) is proposed for finding the numerical solutions of the infinite horizon optimal control problem of nonlinear interconnected large-scale dynamic systems. A modified Laguerre method is used to discretize the equation of optimal condition, while a homotopy method is used to construct an iterative scheme. Two illustrative examples demonstrated that LaHOC has spectral accuracy and very good efficiency, which is comparable to wellestablished numerical methods such as the MATLAB BVP5C solver. The second example shows when the multi-components have different time and amplitude scales, one need to use adaptive rescaling technique in the Laguerre bases to improve accuracy, which deserves a further study.

Disclosure statement
No potential conflict of interest was reported by the author(s).