Convergence for the optimal control problems using collocation at Legendre-Gauss points

Abstract

The convergence of the novel Legendre-Gauss method is established for solving a continuous optimal control problem using collocation at Legendre-Gauss points. The method allows for changes in the number of Legendre-Gauss points to meet the error tolerance. The continuous optimal control problem is first discretized into a nonlinear programming problem at Gauss collocations by the Legendre-Gauss method. Subsequently, we prove the convergence of the Legendre-Gauss algorithm under the assumption that the continuous optimal control problem has a smooth solution. Compared with those of the shooting method, the single step method, and the general pseudospectral method, the numerical example shows that the Legendre-Gauss method has higher computational efficiency and accuracy in solving the optimal control problem.

Keywords

Optimal control Legendre-Gauss method nonlinear programming convergence

Introduction

The numerical methods for optimal control problems include indirect methods and direct methods. A pseudospectral method is a classic example of the direct method. The pseudospectral methods have been used extensively to achieve numerical solutions to optimal control problems by Garg et al. (2011), Chen et al. (2021), Im and Choi (2018), and Liu et al. (2014). Elnagar and Kazemi (1995) published a paper on the Legendre pseudospectral methods in IEEE Transactions on Automatic Control and explained the basic steps of discretizing the optimal control problem into an algebraic nonlinear programming problem, which marked the formal rise of the pseudospectral method. Ross and Fahroo (2003) proved the costate mapping theorem of the Legendre pseudospectral methods that provides an order-preserving transformation of the Lagrange multipliers associated with the discretized problem to the discrete covectors associated with the optimal control problem. The basis of the pseudospectral methods is the approximation of the state using a set of basis functions. Some differential-algebraic constraints are enforced at a specified set of collocation points, which are typically chosen based on orthogonal quadrature rules. By parameterizing the state variable and control variable and collocating the differential-algebraic equations using nodes obtained from a Gaussian quadrature, the pseudospectral methods are transformed into a nonlinear programming problem. The pseudospectral methods mainly include the Legendre method in Elnagar and Kazemi (1995) and Ross and Fahroo (2003), the Legendre-Gauss method in Todd (2007), and the Legendre-Radau method in Hou et al. (2015). The Legendre pseudospectral method is different from the Legendre-Gauss pseudospectral method in the location points interval, where the Legendre-Gauss pseudospectral and the Legendre pseudospectral location points lie on the open interval (-1,1), the closed interval [-1,1], respectively. Since endpoints 1 and +1 are collocation points, the Legendre pseudospectral method has no need for additional nodes at the endpoints. Hence, Garg et al. (2010) proved the state at the endpoints naturally appears in the discrete problem, which results in the oscillation of the curve of the costate variable obtained. The Legendre pseudospectral method no longer holds the costate mapping theorem at the endpoint, while the Legendre-Gauss pseudospectral method rule does not discretize the state equation at the endpoint, which is helpful to improve the accuracy of the solution. The Legendre-Gauss pseudospectral method was used by Zhang et al. (2019) to design controllers for the trajectory optimization of unmanned helicopter formations. The convergence rate of the Legendre-Radau pseudospectral method was established by Hager et al. (2015) and applied to the multiobjective optimization of the energy management strategy by Liu et al. (2019). Liu et al. (2015) proposed the $hp$ pseudospectral method for solving optimal control problems. The continuous-time nonlinear optimal control problem depends on the solution of the Hamilton Jacobi Bellman (HJB) equation, which is a nonlinear partial differential equation and difficult to obtain a global analytical solution. To calculate the solution of HJB, various numerical methods have been proposed. The successive Galerkin approximation algorithm for the HJB equations was first proposed by Beard (1995). Smears and Sueli (2014) have solved the completely nonlinear second-order elliptic HJB equation by using the $hp$ -version discontinuous Galerkin finite element.

Generally speaking, the existing results focus on three fundamental issues of pseudospectral methods, namely the state and costate approximation in Fahroo and Ross (2001), the existence of approximate solutions in Gong et al. (2006) and Kang et al. (2008), convergence in Dontchev et al. (1995), and the convergence rate. In recent years, convergence and convergence rate of pseudospectral methods for optimal control problems have been studied in the literature include the Legendre-Radau method by Hager et al. (2017), the Legendre-Gauss method by Francolin et al. (2015), Runge-Kutta method by Hager (2000) and Dontchev et al. (2001), and the segment widths and the polynomial degree (hp) collocation method by Elschner (1993).

The convergence of the discrete optimal solutions requires strong conditions, for instance, the normal form of feedback linearizable in Gong et al. (2006) or Lipschitz continuity of the inverse Karush-Kuhn-Tucker mapping in Gong et al. (2004). These conditions require information of the continuous optimal solution which is extremely difficult to get for nonlinear constrained systems. The convergence rate of the Legendre-Gauss method is determined by the number of Legendre-Gauss points and local minimizer for the nonlinear constrained in Hager et al. (2016) and Hager et al. (2015). However, as the number of Legende-Gauss points increases, the convergence rate of the Legende-Gauss method slows down, which increases the calculation time. We seek to carefully select the number of Legendre-Gauss points to reduce the calculation time under the condition of error tolerance. Therefore, based on the general Legendre-Gauss pseudospectral method, we propose a novel Legendre-Gauss algorithm for nonlinear constrained optimal control. Convergence of the novel Legendre-Gauss algorithm for optimal control problem is proved under the assumption that the continuous optimal control problem has a smooth solution in the paper. Different from previous works, the main innovations of this paper are summarized into the following three-fold:

We propose a novel Legendre-Gauss algorithm based on Gauss pseudospectral method to find the approximate numerical solution to the optimal feedback control problems for nonlinear systems. The novel Legendre-Gauss algorithm effectively reduces the computation time by selecting the fewer Legendre-Gauss points under satisfying the error tolerance.

The novel algorithm is designed to get the numerical computational solution for the optimal control problem. We first convert the continuous optimal control problem to discretized problem with collocation at Legendre-Gauss points, which is essentially transformed into a kind of nonlinear programming problem. Then the procedure of this method to compute the numerical solution to the optimal control problems is proposed.

Under the assumption that the continuous optimal control problem has a smooth solution, the convergence of the novel Legendre-Gauss algorithm is proved.

The rest of this paper is organized as follows. The problem description is presented in Section 2. In Section 3, we convert the optimal control problem to nonlinear programming and present the steps of the proposed algorithm. In Section 4, the algorithm of the optimal control problem is given, and the mathematical proof and demonstration of the proposed method are carried out in detail. And numerical example illustrates that convergence results of the Legendre-Gauss algorithm are more excellent than other methods in Section 5. Finally, a brief conclusion is present in Section 6.

Problem formula

The optimal control problem is considered as follows:

Minimize

C (x (t_{f})),

(2.1a)

subject to

\overset{\cdot}{x} (t) = f (x (t), u (t)), t \in [t_{0}, t_{f}],

(2.1b)

x (t_{0}) = x_{0}, (x (t), u (t)) \in C^{1} (R^{n}) \times C^{0} (R^{m}) .

(2.1c)

This problem is recorded as a continuous problem $A$ , where the state variable $x (t) \in R^{n}$ , the control variable $u (t) \in R^{m}$ , $f : R^{n} \times R^{m} \to R^{n}$ , $C : R^{n} \to R$ , $t_{0}$ indicates the initial time, $t_{f}$ indicates the final time and $x_{0}$ is the initial state variable, $C^{k}$ denotes the space of $k$ times continuously differentiable functions.

By introducing a new variable $τ$ and a linear transformation on time is

t = \frac{t_{f} - t_{0}}{2} τ + \frac{t_{f} + t_{0}}{2} .

(2.2)

By the linear transformation (2.2), the interval $[t_{f}, t_{0}]$ is mapped to the interval $[- 1, 1]$ and the continuous problem $A$ can be equivalently discretized into the following form:

Minimize

C (x (t_{f})),

(2.3a)

subject to

\overset{\cdot}{x} (τ_{i}) = \frac{t_{f} - t_{0}}{2} f (x (τ_{i}), u (τ_{i})), 1 \leq i \leq N,

(2.3b)

x (- 1) = x_{0}, x \in P_{N}^{n} .

(2.3c)

The problem is recorded as a discrete problem $B$ in this paper, where the collocation points $τ_{i}$ , 1 $\leq i \leq N$ , are Legendre-Gauss points, and $u (τ_{i})$ is the control approximation at $τ_{i}$ , 1 $\leq i \leq N$ . The sup-norm $∥ . ∥_{\infty}$ is given by

∥ x ∥_{\infty} = sup {| x (t) | : t \in [- 1, + 1]},

where |.| is the Euclidean norm.

Notation. The $∥ x ∥_{\infty}$ denotes the maximum of $| x (t) |$ over $t \in [- 1, + 1]$ . The notation $E_{i : j}$ denotes the Gauss pseudospectral differentiation matrix columns $i$ through $j$ and the notation $E_{i : j}^{†}$ denotes the adjoint differential approximation matrix columns $i$ through $j$ . In addition, $P_{N}$ denotes the space of polynomials of degrees at most defined on the interval $[- 1, + 1]$ and $P_{N}^{n}$ denotes the n-fold Cartesian product $P_{N} \times P_{N} \times . . . \times P_{N}$ . Finally, $0$ denotes zeros matrix and 0 denotes zero constant.

Legendre-Gauss method for optimal control

In numerical analysis and scientific computing, the Legendre-Gauss methods are a family of numerical methods for ordinary differential equations. More specifically, they are collocation methods based on the points of Legendre-Gauss quadrature.

Colloction at Legendre-Gauss points

Legendre-Gauss points are collocated and satisfy

- 1 < τ_{1} < τ_{2} < . . . < τ_{N} < + 1 .

In addition, two non-collocated points are used with $τ_{0} = - 1$ and $τ_{N + 1} = + 1$ . The state variables and control variables at the Legendre-Gauss points are discrete in the Legendre-Gauss method, where the Legendre-Gauss points are the zero point of the Legendre polynomial and the number of $N$ zero points are distributed over the interval $(- 1, 1)$ . The state variables are approximated by the Lagrange interpolation as follows

x (t) \approx X (τ) = \sum_{i = 0}^{N} X (τ_{i}) L_{i} (τ),

where $X (τ)$ is discrete at Legendre-Gauss points, and

L_{i} (τ) = Π_{j = 0, j \neq i}^{N} \frac{τ - τ_{j}}{τ_{i} - τ_{j}} .

Similarly, the control variables $u^{k} (τ)$ in $S_{k}$ is approximated by the Lagrange interpolation as follows

u^{k} (τ) \approx U^{k} (τ) = \sum_{j = 0}^{N} L_{j} (τ) U^{k} (τ_{j}) .

Differentiation of $X (τ)$ , then

\overset{\cdot}{X} (τ) = \sum_{i = 0}^{N} E_{j : i} X_{i},

where

E_{j : i} = {\overset{\cdot}{L}}_{i} (τ_{j}), 1 \leq i, j \leq N,

and

{\overset{\cdot}{L}}_{i} (τ) = \sum_{i = 0}^{N} Π_{j = 0, j \neq i}^{N} \frac{τ - τ_{j}}{τ_{i} - τ_{j}} .

$N \times (N + 1)$ matrix $E$ is called the Gauss pseudospectral differentiation matrix. From the Gauss integration formula, the terminal variable $X_{N + 1}$ can be expressed as

X_{N + 1} \approx X_{0} + \sum_{i = 1}^{N} ω_{i} f (X_{i}, U_{i}),

where $ω_{i}$ is Gauss weight. So the discrete problem $B$ can be reformulated as the nonlinear programming problem:

Minimize

C (X_{N + 1}),

(3.1a)

subject to

\sum_{i = 0}^{N} E_{j : i} X_{i} = \frac{t_{f} - t_{0}}{2} f (X_{i}, U_{i}), 1 \leq i, j \leq N,

(3.1b)

X_{N + 1} = X_{0} + \sum_{i = 1}^{N} ω_{i} f (X_{i}, U_{i}) .

(3.1c)

The problem is recorded as a discrete problem $C$ . And the nonlinear programming problem was solved by the solver sparse nonlinear optimize by Gill et al. (2005), which is a software package for solving large-scale optimization problems and employs a sequential quadratic programming method. For the Legendre-Gauss method, the costate variable $Λ (τ)$ should be discretized at the Legendre-Gauss points, so we have

λ (τ) \approx Λ (τ) = \sum_{i = 1}^{N + 1} L_{i}^{*} (τ) Λ (τ_{i}) = \sum_{i = 1}^{N + 1} L_{i}^{*} (τ) Λ_{i},

L_{i}^{*} (τ) = Π_{j = 1, j \neq i}^{N + 1} \frac{τ - τ_{j}}{τ_{i} - τ_{j}} .

The derivative of the costate variable can be expressed as

\overset{\cdot}{Λ} (τ) = \sum_{i = 1}^{N + 1} {\overset{\cdot}{L}}_{i}^{*} (τ) Λ (τ_{i}) = \sum_{i = 1}^{N + 1} E_{j : i}^{†} Λ_{i} .

Also, from Todd (2007), the following relationship exists between the adjoint differential approximation matrix and the Gauss pseudospectral differentiation matrix

E_{i : j} = - \frac{ω_{j}}{ω_{i}} E_{j : i}^{†} .

The novel Legendre-Gauss algorithm

Pseudospectral optimal control is a joint theoretical-computational method for solving optimal control problems. In a pseudospectral method, the continuous functions are approximated at a set of carefully selected quadrature nodes. The quadrature nodes are determined by the corresponding orthogonal polynomial basis used for the approximation. Mathematically, quadrature nodes can achieve high accuracy with a small number of points. The novel Legendre-Gauss algorithm is proposed for an optimal control problem as follows (Table 1).

Table 1.

The novel Legendre-Gauss algorithm.

INPUT

: the number of Legendre-Gauss points

N

, the errortolerance

ϵ

, and the initial state variable

x_{0}

OUTPUT

: the state variable

x^{I} (τ_{i})

and the control variable

u^{I} (τ_{i})

.
1. Get Legendre-Gauss points

τ_{i}, 1 \leq i \leq N

.
2. Set Gauss weights

ω_{i}

.
3. Obtain Lagrange polynomials

L_{i} (τ)

.
4. Compute the Gauss pseudospectral differentiation matrix

E_{j : i}

.
5. Obtain

{\overset{\cdot}{x}}^{I} (τ_{i})

and

{\overset{\cdot}{x}}^{*} (τ_{i})

.
6.

While

∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} > ϵ

do

N = N + 1

.
8. Update

τ_{i}, ω_{i}, L_{i} (τ), E_{j : i}

.
9. Update

{\overset{\cdot}{x}}^{I} (τ_{i})

and

{\overset{\cdot}{x}}^{*} (τ_{i})

.
10.

End

while

11. Convert the problem to NLP.
12. Verify the KKT condition of NLP.
13. Solve the NLP problem.
14.

Return

x^{I} (τ_{i})

and

u^{I} (τ_{i})

Where $x^{*} (τ_{i})$ is a state variable local minimizer for the continuous problem A, $x^{I} (τ_{i})$ and $u^{I} (τ_{i})$ are variable minimizers for the discrete problem C, NLP means the nonlinear programming problem of the optimal control problem. And the convergence of the Legendre-Gauss algorithm under the infinite norm $∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty}$ will be discussed in the next section.

Convergence

As a matter of the fact, the control variable and state variable in problem A have the following relationship

u (t) = g (x (t)),

which shows that the control variable $u (t)$ is the function of the state variable $x (t)$ . So the nonlinear programming problem of the optimal control problem can be transformed into the following problem:

Minimize

C (X_{N + 1}),

subject to

\sum_{i = 0}^{N} E_{j : i} X_{i} = \frac{t_{f} - t_{0}}{2} f (X_{i}), 1 \leq i, j \leq N,

X_{N + 1} = X_{0} + \sum_{i = 1}^{N} ω_{i} f (X_{i}) .

The Lagrangian of the nonlinear programming problem is

\begin{matrix} ℓ = C (X_{N + 1}) + < Λ_{i}, \sum_{i = 0}^{N} E_{j : i} X_{i} - \frac{t_{f} - t_{0}}{2} f (X_{i}) > \\ + < Λ_{N + 1}, X_{0} + \sum_{i = 1}^{N} ω_{i} f (X_{i}) - X_{N + 1} >, \end{matrix}

where $Λ_{i}, 1 \leq i \leq N + 1$ , is costate variable. The Karush-Kuhn-Tucker optimality conditions for the nonlinear programming problem are given as

\begin{matrix} \sum_{i = 1}^{N + 1} E_{j : i}^{†} Λ_{i} = - \nabla_{x} H (X_{i}, Λ_{i}), 1 \leq i, j \leq N, \\ Λ_{N + 1} = \nabla_{x} C (X_{N + 1}) . \end{matrix}

The discrete costate $Λ_{i}$ and the KKT multipliers $λ_{i}$ have the following relationship

ω_{i} Λ_{i} = λ_{i} + ω_{i} λ_{N + 1}, 1 \leq i \leq N,

and

Λ_{N + 1} = λ_{N + 1} .

The first-order necessary condition for the nonlinear programming problem can be written as $F (X, Λ) = 0$ , where $F = (F_{1}, F_{2}, F_{3}, F_{4})$ , which means

(F_{1}, F_{2}, F_{3}, F_{4}) (X, Λ) \in R^{nN} \times R^{n} \times R^{nN} \times R^{n} .

Their expressions are as follows

\begin{matrix} F_{1 i} (X, Λ) = \sum_{i = 0}^{N} E_{j : i} X_{i} - \frac{t_{f} - t_{0}}{2} f (X_{i}), 1 \leq i, j \leq N, \end{matrix}

\begin{matrix} F_{2 i} (X, Λ) = X_{N + 1} - X_{0} - \sum_{i = 1}^{N} ω_{i} f (X_{i}), 1 \leq i \leq N, \end{matrix}

\begin{matrix} F_{3 i} (X, Λ) = \sum_{i = 1}^{N + 1} E_{j : i}^{†} Λ_{i} + \nabla_{x} H (X_{i}, Λ_{i}), 1 \leq i, j \leq N, \end{matrix}

\begin{matrix} F_{4} (X, Λ) = Λ_{N + 1} - \nabla_{x} C (X_{N + 1}) . \end{matrix}

In Hager et al. (2016), studied the convergence rate of the Legendre-Gauss method but did not study the convergence. In the paper, the Legendre-Gauss algorithm is proposed and the convergence of the Legendre-Gauss algorithm is studied.

Lemma 1. Hager et al. (2016) Let $(x^{*}, u^{*})$ is a local minimizer for the continuous problem A with $(x^{*}, λ^{*})$ . If $(x^{*}, λ^{*}) \in C^{η}$ , $m = \min {η, N + 1}$ , and $N$ is the number of Legendre-Gauss points, then there exists a constant c independent of $N$ and $η$ such that

∥ F (X^{*}, Λ^{*}) ∥_{\infty} \leq {(\frac{c}{N})}^{m - 3} (∥ x^{* (m)} ∥_{\infty} + ∥ λ^{* (m)} ∥_{\infty}) .

(4.1)

Now, we analyze the convergence conditions for the Legendre-Gauss method applied to the optimal control. From Lemma 1, the convergence is equivalent to the fact that the right side of equation (4.1) tends to zero when $N$ approaches positive infinity. The mapping $F_{1 i}$ is considered, firstly refer to Garg et al. (2010), then

\sum_{i = 0}^{N} E_{j : i} X_{i}^{*} = {\overset{\cdot}{x}}^{I} (τ_{i}), 1 \leq i, j \leq N,

where $x^{I} \in P_{N}^{n}$ is the polynomial that passes through $x^{*} (τ_{i})$ for $0 \leq i \leq N$ . Since $x^{*}$ satisfies the dynamics of the continuous problem A, we have

f (X_{i}) = {\overset{\cdot}{x}}^{*} (τ_{i}) .

After the above analysis, $F_{1 i}$ can be converted into

F_{1 i} (X^{*}, Λ^{*}) = {\overset{\cdot}{x}}^{I} (τ_{i}) - \frac{t_{f} - t_{0}}{2} {\overset{\cdot}{x}}^{*} (τ_{i}) .

Next, we consider

F_{2 i} (X^{*}, Λ^{*}) = x^{*} (1) - x^{*} (- 1) - \sum_{i = 1}^{N} ω_{i} f (x^{*} (τ_{i})) .

(4.2)

By the fundamental theorem of calculus, we have

x^{I} (1) - x^{I} (- 1) - \int_{- 1}^{1} {\overset{\cdot}{x}}^{I} (t) dt = 0

Since $x^{I} (- 1) = x^{*} (- 1)$ , and the fact that $N$ -point Gauss quadrature is accurate for polynomials of degree up to $2 N - 1$ , we have

x^{I} (1) - x^{I} (- 1) - \int_{- 1}^{1} {\overset{\cdot}{x}}^{I} (t) dt = x^{I} (1) - x^{*} (- 1) - \sum_{i = 1}^{N} ω_{i} {\overset{\cdot}{x}}^{I} (τ_{i}),

then

F_{2 i} (X^{*}, Λ^{*}) = x^{*} (1) - x^{I} (1) + \sum_{i = 1}^{N} ω_{i} ({\overset{\cdot}{x}}^{I} (τ_{i}) - {\overset{\cdot}{x}}^{*} (τ_{i})) .

(4.3)

Now, we consider $F_{3 i}$ . It is derived by Lemma 1

\sum_{j = 1}^{N + 1} E_{i : j}^{†} Λ_{j}^{*} = {\overset{\cdot}{λ}}^{I} (τ_{i}), 1 \leq i, j \leq N,

where $λ^{I} \in P_{N}^{n}$ is the polynomial that passes through ${Λ^{*}}_{j} = λ (τ_{j})$ for $1 \leq j \leq N + 1$ . Since $λ^{*}$ satisfies

{\overset{\cdot}{λ}}^{*} (τ_{i}) = - \nabla_{x} H (X_{i}^{*}, {U^{*}}_{i}, Λ_{i}^{*}),

we have

F_{3 i} (X^{*}, Λ^{*}) = {\overset{\cdot}{λ}}^{I} (τ_{i}) - {\overset{\cdot}{λ}}^{*} (τ_{i}) .

Finally, we consider $F_{4}$ , it has come from Hager et al. (2016) and Garg et al. (2010)

λ_{N + 1}^{*} = \nabla_{x} C (X_{N + 1}),

and

λ_{N + 1}^{*} = Λ_{N + 1}^{*} .

Thus, we get

F_{4} (X, Λ) = Λ_{N + 1} - \nabla_{x} C (X_{N + 1}) = 0 .

Combining the analysis of $F_{1 i}$ to $F_{4 i}$ above, the norm $∥ F (X^{*}, Λ^{*}) ∥_{\infty}$ can be equivalently converted to norm $∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty}$ and $∥ {\overset{\cdot}{λ}}^{I} - {\overset{\cdot}{λ}}^{*} ∥_{\infty}$ . Then the forms of the norms $∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty}$ and $∥ {\overset{\cdot}{λ}}^{I} - {\overset{\cdot}{λ}}^{*} ∥_{\infty}$ are similar, so we only need to analyze the convergence of the norm $∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty}$ , and the convergence of the norm $∥ {\overset{\cdot}{λ}}^{I} - {\overset{\cdot}{λ}}^{*} ∥_{\infty}$ can be similarly obtained.

Theorem 1. If $(x^{*}, λ^{*})$ is a local minimizer for the continuous problem A, $(x^{I}, λ^{I})$ is a local minimizer for the discrete problem C, and $N$ is the number of Legendre-Gauss points, when $m = \min {η, N + 1} \geq 4$ , then

lim_{N \to + \infty} ∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} = 0,

and

lim_{N \to + \infty} ∥ {\overset{\cdot}{λ}}^{I} - {\overset{\cdot}{λ}}^{*} ∥_{\infty} = 0 .

More precisely, the Legendre-Gauss method applied to the optimal control converges when $m = \min {η, N + 1} \geq 4$ .

Proof. By Hager et al. (2017), then

\begin{matrix} ∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} \leq (1 + 2 N^{2}) inf_{q \in P_{N}^{n}} ∥ \overset{\cdot}{x} - \overset{\cdot}{q} ∥_{\infty} \\ + N^{2} (1 + L_{N}) inf_{q \in P_{N}^{n}} ∥ x - q ∥_{\infty} \end{matrix}

(4.4)

where $L_{N}$ is the Lebesgue constant at the Legendre-Gauss points set $τ_{i}$ , $0 \leq i \leq N$ . By Hager et al. (2017), we have

L_{N} = O (\sqrt{N}) .

(4.5)

When $N + 1 \geq η$ , by Elschner (1993), the new inequalities for some constant $c$ , independent of $N$ , we obtain

inf_{q \in P_{N}^{n}} ∥ x - q ∥_{\infty} \leq {(\frac{c}{N + 1})}^{η} ∥ x^{* (η)} ∥_{\infty},

(4.6)

and

inf_{q \in P_{N}^{n}} ∥ \overset{\cdot}{x} - \overset{\cdot}{q} ∥_{\infty} \leq {(\frac{c}{N})}^{η - 1} ∥ x^{* (η)} ∥_{\infty} .

(4.7)

Substituting (4.6) and (4.7) into (4.4), we have

\begin{matrix} ∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} \leq [(1 + 2 N^{2}) {(\frac{c}{N})}^{m - 1} + N^{2} (1 + L_{N}) {(\frac{c}{N})}^{η - 1}] \\ ∥ x^{* (η)} ∥_{\infty}, \end{matrix}

(4.8)

When $m = \min {η, N + 1} \geq 4$ , we obtain

lim_{N \to + \infty} (1 + 2 N^{2}) {(\frac{c}{N})}^{m - 1} ∥ x^{* (m)} ∥_{\infty} = 0 .

and

lim_{N \to + \infty} N^{2} (1 + L_{N}) {(\frac{c}{N})}^{m - 1} ∥ x^{* (m)} ∥_{\infty} = 0 .

Therefore, we have

lim_{N \to + \infty} ∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} = 0 .

Hence, the norm $∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty}$ converges when $m = \min {η, N + 1} \geq 4$ . Similarly, if $m = \min {η, N + 1} \geq 4$ , we obtain

lim_{N \to + \infty} ∥ {\overset{\cdot}{λ}}^{I} - {\overset{\cdot}{λ}}^{*} ∥_{\infty} = 0 .

Therefore, the Legendre-Gauss method applied to the optimal control converges when $m = \min {η, N + 1} \geq 4$ .

Theorem 2. If $m = \min {η, N + 1}$ , and $N$ is the number of Legendre-Gauss points, then

lim_{N \to + \infty} [∥ {\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i - 1}^{I} ∥_{\infty} - ∥ {\overset{\cdot}{x}}_{i}^{*} - {\overset{\cdot}{x}}_{i - 1}^{*} ∥_{\infty}] = 0 .

Proof: In fact, (4.8) is equal to the following inequality

∥ {\overset{\cdot}{x}}^{I} - {\overset{\cdot}{x}}^{*} ∥_{\infty} \leq {(\frac{c}{N})}^{m - 3} ∥ x^{* (m)} ∥_{\infty} .

(4.9)

By the triangle inequality, we have

\begin{matrix} ∥ {\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i - 1}^{I} ∥_{\infty} - ∥ {\overset{\cdot}{x}}_{i}^{*} - {\overset{\cdot}{x}}_{i - 1}^{*} ∥_{\infty} \leq ∥ ({\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i - 1}^{I}) - ({\overset{\cdot}{x}}_{i}^{*} - {\overset{\cdot}{x}}_{i - 1}^{*}) ∥_{\infty} \\ \leq ∥ {\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i}^{*} ∥_{\infty} + ∥ {\overset{\cdot}{x}}_{i - 1}^{I} - {\overset{\cdot}{x}}_{i - 1}^{*} ∥_{\infty} \end{matrix}

Combining with (4.9), we obtain

∥ {\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i - 1}^{I} ∥_{\infty} - ∥ {\overset{\cdot}{x}}_{i}^{*} - {\overset{\cdot}{x}}_{i - 1}^{*} ∥_{\infty} \leq {(\frac{c}{N})}^{m - 3} (∥ x_{i}^{* (m)} ∥_{\infty} + ∥ x_{i - 1}^{* (m)} ∥_{\infty}) .

For $∥ x_{i}^{* (m)} ∥_{\infty}$ and $∥ x_{i - 1}^{* (m)} ∥_{\infty}$ are constants, then

lim_{N \to + \infty} [∥ {\overset{\cdot}{x}}_{i}^{I} - {\overset{\cdot}{x}}_{i - 1}^{I} ∥_{\infty} - ∥ {\overset{\cdot}{x}}_{i}^{*} - {\overset{\cdot}{x}}_{i - 1}^{*} ∥_{\infty}] = 0 .

Numerical example

Numerical examples were performed in the software MATLAB R2018b, in which the computer was configured as CPU Intel(R) Core(TM) i5-7200U and the system was Windows 10 with a memory of 4.00GB. In the numerical example, the proposed Legendre-Gauss algorithm is compared to the shooting method by Bock and Plitt (1984), the single step method by Cui and Wang (2006), and the general pseudospectral method by Ross and Fahroo (2003), where the general pseudospectral method is operated in the software general pseudospectral optimal control software (GPOPS) by Rao et al. (2010). The single step method is a common numerical method for solving initial value problems of ordinary differential equations and has the advantages of high precision, unconditional stability, and no transcendence. The shooting method is an effective numerical method to solve the initial value problem of an ordinary differential equation. As a numerical solution of optimal control problems, the shooting method is presented in Bock and Plitt (1984). The general pseudospectral method takes the roots of orthogonal polynomials as collocation points and uses global interpolation polynomials as a finite basis to transform the optimal control problem into a nonlinear programming problem. In 5.1 the Bryson-Denham problem, we have shown that the proposed method achieves better performance than the shooting method and the single step method, and the calculation time is increased by adding the number of Legendre-Gauss points for the traditional pseudospectral method. Compared with GPOPS, the proposed method has higher computational efficiency and accuracy in solving the optimal control problem in 5.2 the moon soft landing problem.

The Bryson-Denham problem

Consider the Bryson-Denham optimal control problem taken from Bryon and Ho (1975):

Minimize the cost functional

\frac{1}{2} \int_{0}^{1} u^{2} dt,

subject to the dynamic constraints

\begin{matrix} \overset{\cdot}{x} (t) = v (t), \\ \overset{\cdot}{v} (t) = u (t), \end{matrix}

and the boundary conditions

{\begin{matrix} x (0) = x (1) = 0, \\ v (0) = - v (1) = 1, \end{matrix}

and the state inequality path constraint

{\begin{matrix} x (t) \leq ℓ, \\ ℓ \in (0, \frac{1}{6}) . \end{matrix}

Analytical solution of the Bryson-Denham optimal control problem as follows

u^{*} (t) = {\begin{matrix} - \frac{2}{3 ℓ} (1 - \frac{t}{3 ℓ}) & 0 \leq t < 3 ℓ, \\ 0 & 3 ℓ \leq t \leq 1 - 3 ℓ, \\ - \frac{2}{3 ℓ} (1 - \frac{1 - t}{3 ℓ}) & 1 - 3 ℓ \leq t \leq 1 . \end{matrix}

Note that the Bryson-Denham optimal control problem is the optimal control problem with the Lagrangian performance index. By introducing new state variables, the optimal control problem with the Lagrangian performance index can be equivalently transformed into the optimal Control problem with the endpoint cost. we introduce new state variables

\overset{\cdot}{z} (t) = \frac{1}{2} u^{2} (t),

(5.1)

and an initial value is

z (0) = 0 .

Integrating (5.1) from $t_{0}$ to $t$ , we have

z (t) = \int_{t_{0}}^{t} \frac{1}{2} u^{2} (t) dt .

Then, we have

z (t_{f}) = \int_{t_{0}}^{t_{f}} \frac{1}{2} u^{2} (t) dt,

where $t_{0} = 0$ and $t_{f} = 1$ . After introducing the new state variable, the Bryson-Denham problem is equivalent to the optimal control problem with the endpoint cost:

Minimize the cost functional

z (t_{f}) = \int_{t_{0}}^{t_{f}} \frac{1}{2} u^{2} (t) dt,

subject to the dynamic constraints

{\begin{matrix} \overset{\cdot}{x} (t) = v (t), \\ \overset{\cdot}{v} (t) = u (t), \\ \overset{\cdot}{z} (t) = \frac{1}{2} u^{2} (t), \end{matrix}

and the boundary conditions

{\begin{matrix} x (0) = x (1) = 0, \\ v (0) = - v (1) = 1, \\ z (0) = 0 . \end{matrix}

To demonstrate the effectiveness and usefulness of the proposed method for optimal control problems, the proposed method, single shooting method, and single step method are used to solve the problem. The solution results of each method are shown below Figure 1 and Figure 2 show the results obtained using the proposed method, Figure 3 and Figure 4 are solved by the single shooting method, and the results solved by the single step method are shown in Figure 5 and Figure 6. In these figures, the solid line refers to the true solution of the problem and the dotted line refers to the numerical solution of this problem. These three methods are discretized by the same transcription parameters where the Legendre-Gauss points are 30 in Figure 1 –Figure 6, Figure 7 and Figure 8 are discretized for the Legendre-Gauss points $N = 50$ .

To compare different computation effectiveness and performance in a same algorithm, we set number of node points to $N = 10, N = 30, N = 50$ , respectively. Table 2 lists the computation time for these three methods. It is shown from this table that the computation time of the single step method is the least where the time unit is second, but it is apparent that from Figure 3 and Figure 4 that the error of this method is maximal. To summarize, the proposed method achieves better performance than the other two methods.

Table 2.

Compare computation time with three methods.

Number TimeMethod	$N$ =10	$N$ =30	$N$ =50
Single shooting	9.1100	10.5435	12.0448
Single step	0.6939	2.4563	2.9946
Proposed method	0.5555	2.5634	8.9133

From Figure 1, Figure 2, Figure 7, and Figure 8, with the increase of Legendre-Gauss points, the accuracy is improved. Only adding the number of Legendre-Gauss points can improve the accuracy, but the calculation time is increased. Therefore, for the traditional pseudospectral method, it is important to select the least Legendre-Gauss points under the condition of error tolerance. The proposed Legendre-Gauss algorithm effectively reduces the calculation time by selecting the fewer Legendre Gauss points under satisfying the error tolerance. The more effective performance of the novel Legendre-Gauss algorithm compared with the traditional pseudospectral method is demonstrated in the next section.

The moon soft landing problem

In the moon soft landing problem numerical example, the proposed Legendre-Gauss algorithm is compared to the general pseudospectral method (GPOPS). Consider the moon soft landing optimal control problem taken from Sachan K, Padhi R. (2019):

Minimize the cost functional

J = - m (t_{f}),

subject to the dynamic constraints

{\begin{matrix} \overset{\cdot}{h} (t) = v (t), \\ \overset{\cdot}{v} (t) = - g + u (t), \\ \overset{\cdot}{m} (t) = - ku (t), \end{matrix}

and the boundary conditions

{\begin{matrix} h (0) = h_{0}, \\ v (0) = v_{0}, \\ m (0) = M + F, \\ h (t_{f}) = 0, \\ v (t_{f}) = 0, \end{matrix}

and the path constraint

0 \leq u (t) \leq a, \forall t \in [0, t_{f}],

where the mass of the spacecraft is $m (t)$ , the indicates the final time is $t_{f}$ , the altitude is $h (t)$ , the vertical velocity is $v (t)$ , the thrust of the engine is $u (t)$ , the max thrust of the engine is a constant $a$ , $k$ is a constant, the acceleration of the lunar gravity is a constant $g$ , the mass of the spacecraft without fuel is $M$ , and the total mass of the initial fuel is $F$ . Analytical solution of the moon soft landing optimal control problem as follows

u^{*} (t) = {\begin{matrix} a, λ_{v} (t) < k λ_{m} (t), \\ 0, λ_{v} (t) > k λ_{m} (t), \end{matrix}

where $λ_{v} (t)$ is a costate variable of the vertical velocity and $λ_{m} (t)$ is a costate variable of the mass. In the numerical example, the numerical parameters are given in Table 3.

Table 3.

The numerical parameters.

Parameters	Values
$g$	1.6
$k$	1
$a$	3
$h_{0}$	10
$m_{0}$	600
$v_{0}$	−2

We performed experiments on the error tolerance of $10^{- 4} ~ 10^{- 8}$ , respectively. Table 4 shows the number of Legendre-Gauss points of the proposed algorithm and GPOPS under different error tolerance. The calculation time of the proposed algorithm and GPOPS under different the error tolerance in Table 5. Figure 9 and Figure 10 show the change of control variables and state variables of the proposed algorithm when Legendre-Gauss points are $N = 22$ , respectively. Figure 11 and Figure 12 show the change of control variables and state variables of GPOPS when Legendre-Gauss points are $N = 20$ , respectively. The change images of the control variable and state variable of the proposed method when the error tolerance is $ϵ = 10^{- 7}$ in Figure 13 and Figure 14, respectively. Figure 15 and Figure 16 show the change of control variable and state variable of GPOPS at the error tolerance is $ϵ = 10^{- 7}$ , respectively. In Figures 9 –16, the yellow or red lines are the analytical solutions of the control variables or state variables. The discrete points circle or triangle is the numerical solution of the control variable or state variable. As shown in Figures 9 –16, the moon soft landing optimal control problem is a Bang-Bang control problem.

Table 4.

Compare the number of $N$ with GPOPS.

Error tolerance $N$ Method	$10^{- 4}$	$10^{- 5}$	$10^{- 6}$	$10^{- 7}$	$10^{- 8}$
GPOPS	20	33	46	58	68
Proposed method	15	22	37	39	48

Table 5.

Compare computation times with GPOPS.

Error tolerance $Times$ Method	$10^{- 4}$	$10^{- 5}$	$10^{- 6}$	$10^{- 7}$	$10^{- 8}$
GPOPS	1.7117	1.9015	2.0035	2.0296	2.0395
Proposed method	1.5324	1.5551	1.7556	1.7896	1.9312

As shown in Table. 4, when the error tolerance reaches $10^{- 4}$ , $10^{- 5}$ , $10^{- 6}$ , $10^{- 7}$ , and $10^{- 8}$ , the number of the Legendre-Gauss points required by GPOPS is 20, 33, 46, 58, and 68, respectively, while the Legendre-Gauss points required by the proposed algorithm is 15, 22, 37, 39, and 48 under the same error tolerance. The number of the Legendre-Gauss points of the proposed algorithm is less under the same error tolerance. As shown in Table. 5, compared with GPOPS, the calculation time of the proposed algorithm is less under the same error tolerance. For example, under the error tolerance $ϵ = 10^{- 5}$ , the calculation time of the proposed algorithm is 1.5551 s, while GPOPS is 1.9015 s.

In Figure 9 –Figure 12, the proposed algorithm with $N = 22$ reaches the error tolerance reaches $ϵ = 10^{- 5}$ . GPOPS with $N = 20$ reaches the error tolerance reaches $ϵ = 10^{- 4}$ . As shown in Figure 9 –Figure 12, the approximate of the proposed algorithm with $N = 22$ is better than that of GPOPS with $N = 20$ . From Table. 5, the calculation time of the proposed algorithm with $N = 22$ is 1.5551 s, while the calculation time of GPOPS with $N = 20$ is 1.7117 s.

In Figure 13 –Figure 16, under the error tolerance $ϵ = 10^{- 7}$ , the curve of the numerical solution of the proposed algorithm and GPOPS approximates the analytic solution completely. From Table 5, the calculation time of the proposed algorithm with $ϵ = 10^{- 7}$ is 1.7896 s, while the calculation time of GPOPS with $ϵ = 10^{- 7}$ is 2.0296 s. From Table 4, the number of the Legendre-Gauss points of the proposed algorithm with $ϵ = 10^{- 7}$ is 39, while GPOPS with $ϵ = 10^{- 7}$ is 58. Therefore, compared with GPOPS, the proposed Legendre-Gauss algorithm can effectively reduce the calculation time.

The numerical experiments show that the proposed Legendre-Gauss algorithm successfully has solved the moon soft landing problem and has obtained a very effective result. The proposed Legendre-Gauss algorithm effectively reduces the computation time by selecting fewer Legendre-Gauss points under satisfying the error tolerance. The numerical experiments show that the error tolerance decreases gradually with the increase of Gaussian nodes, which verifies the convergence theory.

Conclusion

The convergence for the optimal control problems using the Legendre-Gauss method is present in this paper. The precision of the Legendre-Gauss quadrature is higher than Gauss quadrature, and it is an advantage that the endpoint is one of the abscissas where the function is evaluated. The convergence theory of the Legendre-Gauss method applied to the optimal control problems is established. The proposed method can be used to find a numerical solution to some nonlinear optimal control problems, the numerical example shows the method performs effectively. The proposed method is not limited to application in the optimal control problems, but can also be applied to path planning problems in different fields, such as those of robots and spacecraft and the optimal path planning of autonomous parking.

Footnotes

Acknowledgements

The authors are very grateful to the anonymous referees and the editors for their valuable suggestions and comments.

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is partially supported by the National Natural Science Foundation of China (No.11671407).

ORCID iD

Heng Mai

References

Beard

(1995) Improving the closed-looped performance of nonlinear systems. Phd thesis, Rensselaer Polytechnic Institute, New York.

Bock

Plitt

(1984) A multiple shooting algorithm for direct solution of optimal control problems. In: IFAC World Congress, 1984, pp. 243–247.

Bryon

(1975) Applied Optimal Control. Taylor and Francis, Washington, D.C.

Chen

Mai

Zhang

(2021) A novel adaptive pseudospectral method for the optimal control problem of automatic car parking. Asian Journal Control. DOI:10.1002/asjc.2534 .

Cui

Wang

(2006) Research on a high order single step method for analytic expression of stiffness. Earthquake Engineering and Engineering Vibration 26(4): 63–67.

Dontchev

Hager

Poore

Yang

(1995) Optimality, stability, and convergence in nonlinear control. Applied Mathematics and Optimization 31(3): 297–326.

Dontchev

Hager

Veliov

(2001) Second-order runge-kutta approximations in control constrained optimal control. SIAM Journal on Numerical Analysis 38(1): 202–226.

Elnagar

Kazemi

(1995) The pseudospectral legendre method for discretizing optimal control problems. IEEE Transactions on Automatic Control 40(10): 1793–1796.

Elschner

(1993) The

h

p

-version of spline approximation methods for mellin convolution equations. Journal of Integral Equations and Applications 5(1): 47–73.

10.

Fahroo

Ross

(2001) Costate estimation by a Legendre pseudospectral method. Journal of Guidance Control and Dynamics 24(2): 270–277.

11.

Francolin

Benson

Hager

Rao

(2015) Costate approximation in optimal control using integral Gaussian quadrature orthogonal collocation methods. Optimal Control Applications and Methods 36(4): 381–397.

12.

Garg

Hager

Rao

(2011) Pseudospectral methods for solving infinite-horizon optimal control problems. Automatica 47(4): 829–837.

13.

Garg

Patterson

Hager

, et al. (2010) A unified framework for the numerical solution of optimal control problems using pseudospectral methods. Automatica 46(11): 1843–1851.

14.

Gill

Murray

Saunders

, et al. SNOPT: An SQP algorithm for large-scale constrained optimization[J]. SIAM Review, 2001, 47(1): 99–131.

15.

Gong

Kang

Ross

(2006) A pseudospectral method for the optimal control of constrained feedback linearizable systems. IEEE Transactions on Automatic Control 51(7): 1115–1129.

16.

Gong

Ross

Kang

Fahroo

(2004) Convergence of pseudospectral methods for constrained nonlinear optimal control problems. In: Proceedings of the IASTED International Conference on Intelligent Systems and Control, Honolulu, HI: Acta Press, Calgary, Canada; 209–214.

17.

Hager

(2000) Runge-Kutta methods in optimal control and the transformed adjoint system. Numerische Mathemati 87(2): 247–282.

18.

Hager

Hou

Rao

(2016) Convergence rate for a Gauss collocation method applied to unconstrained optimal control. Journal of Optimization Theory and Applications 169(3): 801–824.

19.

Hager

Hou

Rao

(2017) Lebesgue constants arising in a class of collocation methods. IMA Journal of Numerical Analysis 37(4): 1884–1901.

20.

Hager

Hou

Rao

(2015) Convergence rate for a radau collocation method applied to unconstrained optimal control. arXiv preprint arXiv:1508.03783.

21.

Hou

Hager

Rao

(2012) Convergence of a Gauss pseudospectral method for optimal control. In: AIAA Guidance, Navigation, and Control Conference and Exhibit. American Institute of Aeronautics and Astronautics, Minnesota, vol. 8, pp. 1–9.

22.

Choi

(2018) Helicopter rotor flow analysis using mapped Chebyshev pseudospectral method and overset mesh topology. Mathematical Problems in Engineering 2018 (5): 1–20.

23.

Kang

Ross

Gong

(2008) Pseudospectral Optimal Control and its Convergence Theorems, Analysis and design of nonlinear control systems. Springer, Berlin, Heidelberg: 109–124.

24.

Kapil

Radhakant

(2019) Waypoint constrained multi-phase optimal guidance of spacecraft for soft lunar landing. Unmanned Systems 7(2): 83–104.

25.

Liu

Hager

Rao

(2015) Adaptive mesh refinement method for optimal control using nonsmoothness detection and mesh size reduction. Journal of the Franklin Institute 352(10): 4081–4106.

26.

Liu

Lin

Zhao

Zhu

(2019) Multi-objective optimization of energy management strategy on hybrid energy storage system based on radau pseudospectral method. IEEE Access 7(5): 112483–112493.

27.

Liu

Zhu

Huang

Zheng

(2014) A Hermite-Lobatto pseudospectral method for optimal control. Asian Journal of Control 16(5): 1568–1575.

28.

Rao

Benson

Darby

, et al. (2010) GPOPS, A MATLAB software for solving multiple-phase optimal control problems using the Gauss pseudospectral method. ACM Transactions on Mathematical Software 37(2): 1–39.

29.

Ross

Fahroo

(2003) Legendre pseudospectral approximations of optimal control problems. In: Kang

Xiao

Borges

(eds) New Trends in Nonlinear Dynamics and Control, and Their Applications. Springer, Berlin, Heidelberg, 2003: 327–342.

30.

Smears

Sueli

(2014) Discontinuous Galerkin finite element approximation of Hamilton-Jacobi-Bellman equations with cordes coefficients. SIAM Journal On Numerical Analysis 52(2): 993–1016.

31.

Todd

(2007) Advancement and analysis of Gauss pseudospectral transcription for optimal control problems. Massachusetts Institute of Technology: Department of Aeronautics and Astronautics, 2007.

32.

Zhang

Zong

Dou

, et al. (2019) Trajectory optimization and finite-time control for unmanned helicopters formation. IEEE Access 7: 93023–93034.