# An Efficient Pseudospectral Method for Solving a Class of Nonlinear Optimal Control Problems

Authors: Emran Tohidi; Atena Pasban; A. Kilicman; S. Lotfi Noghabi
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper gives a robust pseudospectral scheme for solving a class of nonlinear optimal control problems (OCPs) governed by differential inclusions. The basic idea includes two major stages. At the first stage, we linearize the nonlinear dynamical system by an interesting technique which is called linear combination property of intervals. After this stage, the linearized dynamical system is transformed into a multi domain
dynamical system via computational interval partitioning. Moreover, the integral form of this multidomain dynamical system is considered. Collocating these constraints at the Legendre Gauss Lobatto (LGL) points together with using the Legendre Gauss Lobatto quadrature rule for approximating the involved integrals enables us to transform the basic OCPs into the associated nonlinear programming problems (NLPs). In all parts of this procedure, the associated control and state functions are approximated by piecewise constants and piecewise polynomials, respectively. An illustrative example is provided for confirming the accuracy and applicability of the proposed idea.


---

## Body

## 1. Introduction

Optimal control problems (OCPs) have received considerable attention during the last four decades because of their applications. Such problems arise in many areas of science and engineering and play an important role in the modeling of real-life phenomena in other fields of science. The principal difficulty in studying OCPs via traditional and classical methods lies in their special nature. Obviously, most of OCPs cannot be solved by the well-known indirect methods [1, 2]. Therefore, it is highly desirable to design accurate direct numerical approaches to approximate the solutions of OCPs [3].Among all of the numerical techniques for solving smooth OCPs, orthogonal functions and polynomials have been applied in a huge size of research works. High accuracy and ease of applying these polynomials and functions for OCPs are two important advantages which have encouraged many authors to use them for different types of problems. For solving smooth OCPs, there exist a broad class of methods based on orthogonal polynomials which were presented by famous applied mathematics scientists such as [4, 5]. The fundamental idea of these methods is based upon pseudospectral (or spectral collocation) operational matrices of differentiation. However, Legendre spectral operational matrix of differentiation was used in [6] (for other applications of spectral operational matrices of differentiation see [7]). The best property of the spectral operational matrices of differentiation is the sparsity, while the pseudospectral ones are relatively filled matrices. Another computational approach for solving OCPs which is based on high order Gauss quadrature rules was presented in [8]. However, high order of accuracy may be obtained by this method, but suitable preconditionings should be explored because of its ill-conditioning of the associated algebraic system.In many real mathematical models, the controller should be restricted. In other words, the control functions of OCPs are bounded in many cases. According to the classical theory of optimal control [9], if the control functions are bounded and appear linearly in the cost functionals and dynamical systems, the resulting problem is a Bang-Bang OCP. In this case, the control functions are discontinuous. Therefore, we deal with a nonsmooth OCP. For dealing with such nonsmooth OCPs, some new numerical methods have been proposed in the literature such as [10, 11]. These approaches are based on finite difference methods (FDMs). Simplicity of the discretization by FDMs is usually easy to handle, but lower order of accuracy may make them unsuccessful. Therefore, we should look at high order numerical methods such as spectral or pseudospectral techniques. But, as it is mentioned in the literature, spectral schemes are the best tools just for the problems with smooth solutions and data. In other words, if we apply these methods for approximating nonsmooth functions we usually observe the Gibbs phenomena. The following example illustrates this fact.Example 1 (see [9]).
We consider the following OCP:(1)MinJ=∫02‍(3u(τ)-2y(τ))dτs.t.11y˙(τ)=y(τ)+u(τ),0≤τ≤2,s.t.11y(0)=4,y(2)=39.392,s.t.11u(τ)∈[0,2],0≤τ≤2.
Since the computational interval is[0,2], we should change it into [-1,1] by a simple transformation as follows:
(2)MinJ=∫-11‍(3u(t)-2y(t))dts.t.1y˙(t)=y(t)+u(t),-1≤t≤1,s.t.1y(-1)=4,y(1)=39.392,s.t.1u(t)∈[0,2],-1≤t≤1.
The optimal control of the above-mentioned problem is given in the following form:(3)u*(t)={2-1≤t≤0.096,00.096≤t≤1.
For approximating the control function of this problem, we use the classical spectral method [6]. As it is depicted in Figures 1 and 2, the desired optimal control cannot be obtained in a good manner. From these Figures one can observe that not only the exact value of switching point (i.e., t=0.096) is not detected with a high accuracy, but also the obtained solutions have additional jumps in the boundary of domain. These are the disadvantages of the applying the classical spectral methods for solving nonsmooth problems.Figure 1
Approximate optimal control history of Example1 for N=21.Figure 2
Approximate optimal control history of Example1 for N=23.To delete these mentioned disadvantages, a robust spectral method is presented in [12] for solving a class of nonsmooth OCPs that has some fundamental differences with the classical spectral techniques. First, the computational interval is partitioned into subintervals where the size of each subinterval is considered as an unknown parameter, and this enables us to compute the switching times more efficiently. Second, in contrast with the classical spectral schemes, the integral form of the dynamical system is considered. This equivalent form is found by integrating the differential dynamics and adding the initial conditions.Our fundamental goal of this paper is to extend a new idea which was introduced in [12] to approximate the control and state functions of the following nonsmooth OCP:
(4)MinJ=∫0tf‍f(t,y(t))dts.t.1y˙(t)∈d(t),t∈[0,tf],s.t.1(y(0),y(tf))∈S,
where d(t) is a set of continuous functions on [0,tf] and S is a set which contains boundary points of state variable y(t). Also, y(t)∈ℜn, u(t)∈ℜn, and f(t,y(t))∈ℜ. According to discussions in [11], we can assume that
(5)d(t)={D(t,u(t)):u(t)∈U},t∈[0,tf],
where U⊂ℜn is a compact set and D(t,u(t))=(D1(t,u(t)),D2(t,u(t)),…,Dn(t,u(t)))T is a continuous function on [0,tf]×U. Thus, OCP (4) can be rewritten in the form
(6)MinJ=∫0tf‍f(t,y(t))dts.t.1y˙(t)=D(t,u(t)),u(t)∈U,s.t.1y(0)=y0,y(tf)=yf.It should be noted that the dynamical system of (6) is nonlinear in terms of control u(t). For handling OCP (6) in a proper manner, we first linearize the nonlinear dynamical system by an interesting technique which is called linear combination property of intervals (LCPI). After this stage, the linearized dynamical system is transformed into a multidomain dynamical system via computational interval partitioning. Collocating these constraints at the Legendre Gauss Lobatto (LGL) points together with using the Legendre Gauss Lobatto quadrature rule for approximating the involved integrals enables us to transform the basic OCPs into the associated nonlinear programming problems (NLPs).The paper is organized as follows. Section2 is devoted to linearize the nonlinear dynamical system by using LCPI. In Section 3, we design our basic idea which is based on approximation of the associated control and state functions by piecewise constant and piecewise polynomials, respectively. It should be noted that Legendre Gauss Lobatto points are used for collocating the linearized dynamical system. In Section 4, we present a numerical example, demonstrating the efficiency of the suggested numerical algorithm. Concluding remarks are given in Section 5.


## 2. Dynamical System Linearization

SinceD:[0,tf]×U→ℜn is continuous and [0,tf]×U is a compact and connected subset of ℜn+1, then {D(t,u(t)):u∈U} is a closed set in ℜn. Thus, {Di(t,u(t)):u∈U} for i=1,2,…,n is closed in ℜ. Now, suppose that the lower and upper bounds of the {Di(t,u(t)):u∈U} are li(t) and vi(t), respectively. Therefore,
(7)li(t)≤Di(t,u(t))≤vi(t),t∈[0,tf].In other words,(8)li(t)=minu{Di(t,u(t)):u∈U},t∈[0,tf],vi(t)=maxu{Di(t,u(t)):u∈U},t∈[0,tf].By using LCPI, which was first introduced in [10], Di(t,u(t)) can be approximated as a convex linear combination of its minimum li(t) and maximum vi(t) in the following form:
(9)Di(t,u(t))≈λi(t)vi(t)+(1-λi(t))li(t)=λi(t)αi(t)+li(t),
where αi(t)=vi(t)-li(t) and λi(t)∈[0,1]. It should be mentioned that all the λi(t) are the new associated control variables. Now, the main problem (6) is approximated by the following OCP:
(10)Min∫0tf‍f(t,y(t))dts.t.y˙(t)=A(t)Λ(t)+l(t),s.t.Λ(t)∈[0,1]×[0,1]×⋯×[0,1]︷ntimes,t∈[0,tf],s.t.y(0)=y0,y(tf)=yf,
where A(t)=diag(α1(t),α2(t),…,αn(t))n×n, Λ(t)=(λ1(t),λ2(t),…,λn(t))n×1, and l(t) = (l1(t),l2(t),…,ln(t)))n×1. Note that problem (10) is a Bang-Bang OCP, because in this problem the new control Λ(t) has lower 0 and upper 1 bounds and appears linearly in the dynamical equations. As soon as the controls are assumed to be bang-bang, the problem of finding the required controls becomes one of finding the switching times.


## 3. Discretization of the New OCP Containing Linearized Dynamical System

In the sequel and for simplicity in the discretization procedure, we assume thatn=1 (in other words, Λ(t)=λ(t)) and suppose that problem (10) has an optimal solution with m≥1 switching points denoted by t1,t2,…,tm. So if we set t0=0 and tm+1=tf, then the interval [0,tf] breaks into m+1 subintervals. That is,
(11)[0,tf]=[t0,t1]∪[t1,t2]∪⋯∪[tm,tm+1],
where on each subinterval, λ(t) is constant. We denote λ(t) in kth subinterval with constant bk. Since 0≤λ(t)≤1, we have
(12)0≤bk≤1,k=1,2,…,m+1.Moreover, we take the restriction ofy(t) to the kth subinterval by yk(t). By considering (11), the dynamical system of (10) is conveyed as
(13)22y˙k(t)=A(t)bk+l(t),tk-1≤t≤tk,k=1,2,…,m+1,(14)y1(0)=y0,(15)yk(tk-1)=yk-1(tk-1),k=2,3,…,m+1.It should be noted that (15) is assumed to guarantee the continuity of state functions. Integrating (13) gives rise to the following dynamic equations:
(16)yk(t)=ck-1+∫tk-1t‍(A(s)bk+l(s))ds,k=1,…,m+1,
where tk-1≤t≤tk and
(17)ck={y0k=0,yk(tk)k=1,2,…,m.The final conditiony(tf)=yf is imposed only on ym+1(t) as
(18)ym+1(tm+1)=yf.Therefore, problem (10) is transformed into the following optimization problem:
(19)MinJ=∑k=1s+1‍∫tk-1tk‍f(s,yk(s))dss.t.yk(t)=ck-1+∫tk-1t‍(A(s)bk+l(s))ds,s.t.tk-1≤t≤tkk=1,2,…,m+1,s.t.ys+1(ts+1)=yf,s.t.0≤bk≤1,k=1,2,…,m+1.
For discretizing (19), we assume s^ik, i=0,1,…,N, to be the shifted LGL nodes to subinterval [tk-1,tk]; that is, s^ik=(si)((tk-tk-1)/2)+((tk+tk-1)/2). By using Lagrange interpolation, we approximate yk(t) by
(20)yk(t)≈∑i=0N‍yikL^ik(t),
where yik=yk(s^ik) and
(21)L^ik(t)=Li(2tk-tk-1t-tk+tk-1tk-tk-1).It should be noted thatLi(t) is the ith Lagrange basis function. Since yk(t) is approximated, therefore f can be approximated in the kth subinterval as follows:
(22)f(s,yk(s))≈∑j=0N‍f(s^jk,yk(s^jk))L^jk(t)=∑j=0N‍f(s^jk,yjk)L^jk(t).Now by substituting approximations (20) and (22) into (19), we get
(23)∑j=0N‍yjkL^jk(t)=ck-1+∫tk-1t‍(A(s)bk+l(s))ds,tk-1≤t≤tk.From (17) and (20) for k=1,…,m, we have ck=yNk. Now if we set yN0=y0, then we obtain
(24)ck=yNk,k=0,…,m+1.Collocating (23) at s^ik,  i=0,…,N,   k=1,…,m+1, yields
(25)yik=yNk-1+bk∫tk-1s^ik‍A(s)ds+∫tk-1s^ik‍l(s)ds,i=0,1,…,N.Now, by applying a simple linear transformation, we transform the interval[tk-1,s^ik] into [-1,1] as follows:
(26)s=s^ik-tk-12η+s^ik+tk-12,i=0,1,…,N.Therefore, the Legendre Gauss Lobatto quadrature rule can be applied in the following form:(27)yik=yNk-1+s^ik-tk-12(bk∫-11‍A^(η)dη+∫-11‍l^(η)dη)≈yNk-1+s^ik-tk-12{∑q=0N‍wq(bkA^(sq)+l^(sq))},0≤i≤N,
where A^(η)=A(((s^ik-tk-1)/2)η+((s^ik+tk-1)/2)),  l^(η)=l(((s^ik-tk-1)/2)η+(s^ik+tk-1)/2),  wq = (2/N(N+1))(1/PN2(sq)) for q=0,1,…,N are the LGL weights and PN(x) is the Nth degree Legendre Polynomial.So by consideringym+1(tf)=yNm+1, problem (19) is discretized to the following NLP:
(28)MinJN,m=∑k=1m+1‍∑j=0N‍∑q=0N‍tk-tk-12f(s^jk,yjk)wqLj(sq)s.t.yik-yNk-1-s^ik-tk-12{∑q=0N‍wq(bkA^(sq)+l^(sq))}s.t.=0s.t.i=0,1,…,N,k=1,2,…,m+1,s.t.yNs+1-yf=0,s.t.0≤bk≤1,k=1,2,⋯,m+1.Here,bk, yik, i=0,…,N, k=1,…,m, and parameters t1,…,ts,tf are unknown variables in the NLP. Note that yN0 is known and yN0=y0.In the above discretization procedure, the number of switching points,s, is considered as a known parameter. So at first we should guess the number of switching points. This is the disadvantage of the proposed method. It should be noted that if the number of switching points, s, is chosen correctly, then the resulting value of bk is equal to its lower or upper bounds; furthermore, bk changes in each switching point.


## 4. Numerical Example

We now apply the proposed idea for solving a nonlinear OCP governed by differential inclusion. This example was first introduced in [11]. In the mentioned work, the authors used the simplest form of FDMs. One of the advantages of [11] is that we finally solve a Linear Programming (LP) problem. However, this method has other disadvantages such as needing higher values of approximations (i.e., N), and this leads to ill-conditioning of the associated discrete problem. Our presented ideas do not contains the difficulties of the classical spectral methods and FDMs for solving nonsmooth OCPs and also achieve superior results with respect to at least 3 other methods. These advantages confirm the efficiency of this modern spectral approximation. The following example is modeled using the mathematical software package MAPLE, and the corresponding nonlinear programming problem is solved using the command NLPSolve. It should be noted that if the NLP is univariate and unconstrained except for finite bounds, quadratic interpolation method may be used. If the problem is unconstrained and the gradient of the objective function is available, the preconditioned conjugate gradient (PCG) method may be used. Otherwise, the sequential quadratic programming (SQP) method can be used. According to the structure of our NLP, the SQP method is used.Example 2.
We consider the following nonlinear OCP governed by differential inclusion:(29)MinJ=∫01‍sin(3πt)y(t)dts.t.y˙(t)∈{-tan(π8u3(t)+t):u(t)∈[0,1]},s.t.y(0)=1,y(1)=0.
According to discussions in [11], the above OCP can be rewritten in the following form:
(30)MinJ=∫01‍sin(3πt)y(t)dts.t.y˙(t)=-tan(π8u3(t)+t),s.t.u(t)∈[0,1],y(0)=1,y(1)=0.In this problem, the control function appears nonlinearly, and we should linearize the initial dynamical system. According to the idea of LCPI, the above OCP can be reduced to a linear OCP which is Bang-Bang. Here,D(t,u(t))=-tan((π/8)u3(t)+t). Thus,
(31)l(t)=Minu{-tan(π8u3(t)+t):u(t)∈[0,1]}=-tan(π8+t),v(t)=Maxu{-tan(π8u3(t)+t):u(t)∈[0,1]}=-tan(t),
and hence, α(t)=v(t)-l(t)=tan((π/8)+t)-tan(t). Therefore, D(t,u(t)) can be approximated as follows:
(32)D(t,u(t))≈α(t)λ(t)+l(t)≈(tan(π8+t)-tan(t))λ(t)-tan(π8+t).It should be noted thatλ(t)∈[0,1] is the new control function, which is called associated control. By considering this approximation for D(t,u(t)), the basic OCP is approximated by the following Bang-Bang OCP:
(33)MinJ=∫01‍sin(3πt)y(t)dts.t.y˙(t)=(tan(π8+t)-tan(t))λ(t)-tan(π8+t),s.t.λ(t)∈[0,1],y(0)=1,y(1)=0.According to our experiences in [11], we assume that the number of switching points is s=3. Since by applying this assumption we reach to the exact results in which the new associated control λ(t) is switched between its lower and upper bounds, the numerical results related to the values of switching points and objective function for different values of N are provided in Table 1. Moreover, the associated control λ(t), control u(t), and optimal state y(t) are depicted in Figures 3, 4, and 5, respectively. Moreover, in Table 2 comparisons of the numerical results of the proposed method with respect to the methods of [6, 10, 13] are given. From this table one can see the efficiency and applicability of the suggested method for solving nonlinear OCPs governed by differential inclusions.Table 1
Numerical results of Example2.
N
t
1
t
2
t
3
J
N
6
0.2218
0.4538
0.8874
0.0970
8
0.2214
0.4653
0.8683
0.0960
10
0.2214
0.4317
0.8341
0.0968
12
0.2214
0.4591
0.8890
0.0963
14
0.2214
0.4473
0.8762
0.0969
16
0.2214
0.4481
0.8971
0.0962Table 2
Comparisons of the methods in evaluating the objective functionJ*.
N
Proposed  method
Method  of [6]
Method  of [13]
Method  of [10]
6
0.0970
0.1196
0.1371
—
8
0.0960
0.1058
0.1289
—
10
0.0968
0.0964
0.1152
—
11
0.0960
0.1009
0.1094
—
100
—
—
—
0.0985Figure 3
Approximate optimal associated controlλ(t) history of Example 2 for N=16.Figure 4
Approximate optimal controlu(t) history of Example 2 for N=16.Figure 5
Approximate optimal statey(t) history of Example 2 for N=16.


## 5. Concluding Remarks

In this study, a robust numerical technique has been used for solving a class of optimal control problems (OCPs) governed by differential inclusions. The proposed idea includes linearizing the dynamical system in which the resulting problem is a Bang-Bang OCP. After obtaining this nonsmooth OCP, we use the general idea of [12] for dealing with such problems in the best manner. As observed in the numerical example, the proposed scheme has superior results with regard to at least 3 methods which confirm the applicability of the method. One of the disadvantages of our method is more sensitivity to initial guess in comparison with the classical spectral schemes. However, our idea is terminated successfully by considering an initial guess from the solution of the traditional spectral techniques, even for small values of N.



---

*Source: 357931-2013-09-17.xml*