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

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.

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: Min𝐽=²₀(𝜏)(𝜏)𝑦̇(𝜏)(𝜏)(𝜏)𝑦(0)(2)𝑢(𝜏)∈[]3𝑢−2𝑦𝑑𝜏s.t.=𝑦+𝑢,0≤𝜏≤2,=4,𝑦=39.392,0,2,0≤𝜏≤2. (1)
Since the computational interval is [0,2], we should change it into [−1,1] by a simple transformation as follows: Min𝐽=¹₋₁(𝑡)(𝑡)𝑦̇(𝑡)(𝑡)(𝑡)𝑦()(1)𝑢(𝑡)∈[]3𝑢−2𝑦𝑑𝑡s.t.=𝑦+𝑢,−1≤𝑡≤1,−1=4,𝑦=39.392,0,2,−1≤𝑡≤1. (2)
The optimal control of the above-mentioned problem is given in the following form: 𝑢^∗(𝑡)=2−1≤𝑡≤0.096,00.096≤𝑡≤1. (3)
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., 𝑡=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 Example 1 for 𝑁=21.

Figure 2 

Approximate optimal control history of Example 1 for 𝑁=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: Min𝐽=^𝑡𝑓₀𝑓(𝑡)𝑦̇(𝑡)(𝑡)𝑡,𝑦𝑑𝑡s.t.∈𝑑,𝑡∈0,𝑡_𝑓,𝑦(0),𝑦𝑡_𝑓∈𝑆, (4) where 𝑑(𝑡) is a set of continuous functions on [0,𝑡_𝑓] and 𝑆 is a set which contains boundary points of state variable 𝑦(𝑡). Also, 𝑦(𝑡)∈ℜ^𝑛, 𝑢(𝑡)∈ℜ^𝑛, and 𝑓(𝑡,𝑦(𝑡))∈ℜ. According to discussions in [11], we can assume that 𝑑(𝑡)={𝐷(𝑡,𝑢(𝑡))∶𝑢(𝑡)∈𝑈},𝑡∈0,𝑡_𝑓, (5) where 𝑈⊂ℜ^𝑛 is a compact set and 𝐷(𝑡,𝑢(𝑡))=(𝐷₁(𝑡,𝑢(𝑡)),𝐷₂(𝑡,𝑢(𝑡)),…,𝐷_𝑛(𝑡,𝑢(𝑡)))^𝑇 is a continuous function on [0,𝑡_𝑓]×𝑈. Thus, OCP (4) can be rewritten in the form Min𝐽=^𝑡𝑓₀𝑓(𝑡)𝑦̇(𝑡)((𝑡(𝑡)𝑦(0)𝑡,𝑦𝑑𝑡s.t.=𝐷𝑡,𝑢)),𝑢∈𝑈,=𝑦₀,𝑦𝑡_𝑓=𝑦_𝑓. (6)

It should be noted that the dynamical system of (6) is nonlinear in terms of control 𝑢(𝑡). 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. Section 2 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

Since 𝐷∶[0,𝑡_𝑓]×𝑈→ℜ^𝑛 is continuous and [0,𝑡_𝑓]×𝑈 is a compact and connected subset of ℜ^𝑛+1, then {𝐷(𝑡,𝑢(𝑡))∶𝑢∈𝑈} is a closed set in ℜ^𝑛. Thus, {𝐷_𝑖(𝑡,𝑢(𝑡))∶𝑢∈𝑈} for 𝑖=1,2,…,𝑛 is closed in ℜ. Now, suppose that the lower and upper bounds of the {𝐷_𝑖(𝑡,𝑢(𝑡))∶𝑢∈𝑈} are 𝑙_𝑖(𝑡) and 𝑣_𝑖(𝑡), respectively. Therefore, 𝑙_𝑖(𝑡)≤𝐷_𝑖(𝑡,𝑢(𝑡))≤𝑣_𝑖(𝑡),𝑡∈0,𝑡_𝑓. (7)

In other words, 𝑙_𝑖(𝑡)=min_𝑢𝐷_𝑖(𝑡,𝑢(𝑡))∶𝑢∈𝑈,𝑡∈0,𝑡_𝑓,𝑣_𝑖(𝑡)=max_𝑢𝐷_𝑖(𝑡,𝑢(𝑡))∶𝑢∈𝑈,𝑡∈0,𝑡_𝑓. (8)

By using LCPI, which was first introduced in [10], 𝐷_𝑖(𝑡,𝑢(𝑡)) can be approximated as a convex linear combination of its minimum 𝑙_𝑖(𝑡) and maximum 𝑣_𝑖(𝑡) in the following form: 𝐷_𝑖(𝑡,𝑢(𝑡))≈𝜆_𝑖(𝑡)𝑣_𝑖(𝑡)+1−𝜆_𝑖(𝑡)𝑙_𝑖(𝑡)=𝜆_𝑖(𝑡)𝛼_𝑖(𝑡)+𝑙_𝑖(𝑡), (9) where 𝛼_𝑖(𝑡)=𝑣_𝑖(𝑡)−𝑙_𝑖(𝑡) and 𝜆_𝑖(𝑡)∈[0,1]. It should be mentioned that all the 𝜆_𝑖(𝑡) are the new associated control variables. Now, the main problem (6) is approximated by the following OCP: Min^𝑡𝑓₀𝑓(𝑡)𝑦̇(𝑡)(𝑡)Λ(𝑡)(𝑡),Λ(𝑡)∈𝑡,𝑦𝑑𝑡s.t.=𝐴+𝑙_𝑛times[0,1]×[0,1]×⋯×[0,1],𝑡∈0,𝑡_𝑓,𝑦(0)=𝑦₀,𝑦𝑡_𝑓=𝑦_𝑓, (10) where 𝐴(𝑡)=diag(𝛼₁(𝑡),𝛼₂(𝑡),…,𝛼_𝑛(𝑡))_𝑛×𝑛, Λ(𝑡)=(𝜆₁(𝑡),𝜆₂(𝑡),…,𝜆_𝑛(𝑡))_𝑛×1, and 𝑙(𝑡) = (𝑙₁(𝑡),𝑙₂(𝑡),…,𝑙_𝑛(𝑡)))_𝑛×1. Note that problem (10) is a Bang-Bang OCP, because in this problem the new control Λ(𝑡) 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 that 𝑛=1 (in other words, Λ(𝑡)=𝜆(𝑡)) and suppose that problem (10) has an optimal solution with 𝑚≥1 switching points denoted by 𝑡₁,𝑡₂,…,𝑡_𝑚. So if we set 𝑡₀=0 and 𝑡_𝑚+1=𝑡_𝑓, then the interval [0,𝑡_𝑓] breaks into 𝑚+1 subintervals. That is, 0,𝑡_𝑓=𝑡₀,𝑡₁∪𝑡₁,𝑡₂∪⋯∪𝑡_𝑚,𝑡_𝑚+1, (11) where on each subinterval, 𝜆(𝑡) is constant. We denote 𝜆(𝑡) in 𝑘th subinterval with constant 𝑏^𝑘. Since 0≤𝜆(𝑡)≤1, we have 0≤𝑏^𝑘≤1,𝑘=1,2,…,𝑚+1. (12)

Moreover, we take the restriction of 𝑦(𝑡) to the 𝑘th subinterval by 𝑦^𝑘(𝑡). By considering (11), the dynamical system of (10) is conveyed as 𝑦̇^𝑘(𝑡)=𝐴(𝑡)𝑏^𝑘(𝑡),𝑡+𝑙_𝑘−1≤𝑡≤𝑡_𝑘𝑦,𝑘=1,2,…,𝑚+1,₁(0)=𝑦₀,𝑦^𝑘𝑡_𝑘−1=𝑦^𝑘−1𝑡_𝑘−1,𝑘=2,3,…,𝑚+1. (13)(14)(15)

It should be noted that (15) is assumed to guarantee the continuity of state functions. Integrating (13) gives rise to the following dynamic equations: 𝑦^𝑘(𝑡)=𝑐^𝑘−1+^𝑡_𝑡𝑘−1𝐴(𝑠)𝑏^𝑘+𝑙(𝑠)𝑑𝑠,𝑘=1,…,𝑚+1, (16) where 𝑡_𝑘−1≤𝑡≤𝑡_𝑘 and 𝑐^𝑘=𝑦₀𝑦𝑘=0,_𝑘𝑡_𝑘𝑘=1,2,…,𝑚. (17)

The final condition 𝑦(𝑡_𝑓)=𝑦_𝑓 is imposed only on 𝑦^𝑚+1(𝑡) as 𝑦^𝑚+1𝑡_𝑚+1=𝑦_𝑓. (18)

Therefore, problem (10) is transformed into the following optimization problem: Min𝐽=^𝑠+1_𝑘=1^𝑡𝑘_𝑡𝑘−1𝑓𝑠,𝑦^𝑘(𝑠)𝑑𝑠s.t.𝑦^𝑘(𝑡)=𝑐^𝑘−1+^𝑡_𝑡𝑘−1𝐴(𝑠)𝑏^𝑘(𝑠)𝑡+𝑙𝑑𝑠,_𝑘−1≤𝑡≤𝑡_𝑘𝑦𝑘=1,2,…,𝑚+1,_𝑠+1𝑡_𝑠+1=𝑦_𝑓,0≤𝑏^𝑘≤1,𝑘=1,2,…,𝑚+1. (19) For discretizing (19), we assume 𝑠^𝑘_𝑖̂, 𝑖=0,1,…,𝑁, to be the shifted LGL nodes to subinterval [𝑡_𝑘−1,𝑡_𝑘]; that is, 𝑠^𝑘_𝑖̂=(𝑠_𝑖)((𝑡_𝑘−𝑡_𝑘−1)/2)+((𝑡_𝑘+𝑡_𝑘−1)/2). By using Lagrange interpolation, we approximate 𝑦^𝑘(𝑡) by 𝑦^𝑘(𝑡)≈^𝑁_𝑖=0𝑦^𝑘_𝑖𝐿^𝑘_𝑖(𝑡), (20) where 𝑦^𝑘_𝑖=𝑦^𝑘(𝑠^𝑘_𝑖̂) and 𝐿^𝑘_𝑖(𝑡)=𝐿_𝑖2𝑡_𝑘−𝑡_𝑘−1𝑡𝑡−^𝑘+𝑡_𝑘−1𝑡_𝑘−𝑡_𝑘−1. (21)

It should be noted that 𝐿_𝑖(𝑡) is the 𝑖th Lagrange basis function. Since 𝑦^𝑘(𝑡) is approximated, therefore 𝑓 can be approximated in the 𝑘th subinterval as follows: 𝑓𝑠,𝑦^𝑘(𝑠)≈^𝑁_𝑗=0𝑓𝑠^𝑘_𝑗̂,𝑦^𝑘𝑠^𝑘_𝑗̂𝐿^𝑘_𝑗(𝑡)=^𝑁_𝑗=0𝑓𝑠^𝑘_𝑗̂,𝑦^𝑘_𝑗𝐿^𝑘_𝑗(𝑡). (22)

Now by substituting approximations (20) and (22) into (19), we get ^𝑁_𝑗=0𝑦^𝑘_𝑗𝐿^𝑘_𝑗(𝑡)=𝑐^𝑘−1+^𝑡_𝑡𝑘−1𝐴(𝑠)𝑏^𝑘(𝑠)𝑡+𝑙𝑑𝑠,_𝑘−1≤𝑡≤𝑡_𝑘. (23)

From (17) and (20) for 𝑘=1,…,𝑚, we have 𝑐^𝑘=𝑦^𝑘_𝑁. Now if we set 𝑦⁰_𝑁=𝑦₀, then we obtain 𝑐^𝑘=𝑦^𝑘_𝑁,𝑘=0,…,𝑚+1. (24)

Collocating (23) at 𝑠^𝑘_𝑖̂,  𝑖=0,…,𝑁,   𝑘=1,…,𝑚+1, yields 𝑦^𝑘_𝑖=𝑦_𝑁^𝑘−1+𝑏^𝑘𝑠̂^𝑘_{𝑖𝑡𝑘−1}𝐴(𝑠)𝑠̂𝑑𝑠+^𝑘𝑖_𝑡𝑘−1𝑙(𝑠)𝑑𝑠,𝑖=0,1,…,𝑁. (25)

Now, by applying a simple linear transformation, we transform the interval [𝑡_𝑘−1,𝑠^𝑘_𝑖̂] into [−1,1] as follows: 𝑠̂𝑠=^𝑘𝑖−𝑡_𝑘−12𝑠̂𝜂+^𝑘𝑖+𝑡_𝑘−12,𝑖=0,1,…,𝑁. (26)

Therefore, the Legendre Gauss Lobatto quadrature rule can be applied in the following form: 𝑦^𝑘_𝑖=𝑦_𝑁^𝑘−1+𝑠^𝑘_𝑖̂−𝑡_𝑘−12𝑏^𝑘1₋₁𝐴𝜂𝑑𝜂+¹₋₁𝑙̂𝜂𝑑𝜂≈𝑦_𝑁^𝑘−1+𝑠^𝑘_𝑖̂−𝑡_𝑘−12^𝑁_𝑞=0𝑤_𝑞𝑏^𝑘𝐴𝑠_𝑞+𝑙̂𝑠_𝑞,0≤𝑖≤𝑁, (27) where 𝐴(𝜂)=𝐴(((𝑠^𝑘_𝑖̂−𝑡_𝑘−1)/2)𝜂+((𝑠^𝑘_𝑖̂+𝑡_𝑘−1)/2)),  𝑙(̂𝜂)=𝑙(((𝑠^𝑘_𝑖̂−𝑡_𝑘−1)/2)𝜂+(𝑠^𝑘_𝑖̂+𝑡_𝑘−1)/2),  𝑤_𝑞 = (2/𝑁(𝑁+1))(1/𝑃²_𝑁(𝑠_𝑞)) for 𝑞=0,1,…,𝑁 are the LGL weights and 𝑃_𝑁(𝑥) is the 𝑁th degree Legendre Polynomial.

So by considering 𝑦^𝑚+1(𝑡_𝑓)=𝑦_𝑁^𝑚+1, problem (19) is discretized to the following NLP: Min𝐽_𝑁,𝑚=^𝑚+1𝑁_𝑘=1^𝑁_𝑗=0𝑞=0𝑡_𝑘−𝑡_𝑘−12𝑓𝑠^𝑘_𝑗̂,𝑦^𝑘_𝑗𝑤_𝑞𝐿_𝑗𝑠_𝑞s.t.𝑦^𝑘_𝑖−𝑦_𝑁^𝑘−1−𝑠^𝑘_𝑖̂−𝑡_𝑘−12^𝑁_𝑞=0𝑤_𝑞𝑏^𝑘𝐴𝑠_𝑞+𝑙̂𝑠_𝑞𝑦=0𝑖=0,1,…,𝑁,𝑘=1,2,…,𝑚+1,_𝑁𝑠+1−𝑦_𝑓=0,0≤𝑏^𝑘≤1,𝑘=1,2,⋯,𝑚+1. (28)

Here, 𝑏^𝑘, 𝑦^𝑘_𝑖, 𝑖=0,…,𝑁, 𝑘=1,…,𝑚, and parameters 𝑡₁,…,𝑡_𝑠,𝑡_𝑓 are unknown variables in the NLP. Note that 𝑦⁰_𝑁 is known and 𝑦⁰_𝑁=𝑦₀.

In the above discretization procedure, the number of switching points, 𝑠, 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, 𝑠, is chosen correctly, then the resulting value of 𝑏^𝑘 is equal to its lower or upper bounds; furthermore, 𝑏^𝑘 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., 𝑁), 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: Min𝐽=¹₀()𝑦(𝑡)𝑦̇(𝑡)∈𝜋sin3𝜋𝑡𝑑𝑡s.t.−tan8𝑢³(𝑡)(𝑡)∈[,𝑦(0)(1)+𝑡∶𝑢0,1]=1,𝑦=0. (29)
According to discussions in [11], the above OCP can be rewritten in the following form: Min𝐽=¹₀()𝑦(𝑡)𝑦̇(𝑡)𝜋sin3𝜋𝑡𝑑𝑡s.t.=−tan8𝑢³(𝑡),𝑢(𝑡)∈[](0)(1)+𝑡0,1,𝑦=1,𝑦=0. (30)

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, 𝐷(𝑡,𝑢(𝑡))=−tan((𝜋/8)𝑢³(𝑡)+𝑡). Thus, 𝑙(𝑡)=Min_𝑢𝜋−tan8𝑢³(𝑡)(𝑡)∈[𝜋+𝑡∶𝑢0,1]=−tan8,𝑣(𝑡)+𝑡=Max_𝑢𝜋−tan8𝑢³(𝑡)(𝑡)∈[(𝑡),+𝑡∶𝑢0,1]=−tan (31) and hence, 𝛼(𝑡)=𝑣(𝑡)−𝑙(𝑡)=tan((𝜋/8)+𝑡)−tan(𝑡). Therefore, 𝐷(𝑡,𝑢(𝑡)) can be approximated as follows: 𝐷((𝑡(𝑡)𝜆(𝑡)(𝑡)≈𝜋𝑡,𝑢))≈𝛼+𝑙tan8(𝑡)𝜆(𝑡)𝜋+𝑡−tan−tan8+𝑡. (32)

It should be noted that 𝜆(𝑡)∈[0,1] is the new control function, which is called associated control. By considering this approximation for 𝐷(𝑡,𝑢(𝑡)), the basic OCP is approximated by the following Bang-Bang OCP: Min𝐽=¹₀()𝑦(𝑡)𝑦̇(𝑡)=𝜋sin3𝜋𝑡𝑑𝑡s.t.tan8(𝑡)𝜆(𝑡)𝜋+𝑡−tan−tan8,𝜆(𝑡)∈[](0)(1)+𝑡0,1,𝑦=1,𝑦=0. (33)

According to our experiences in [11], we assume that the number of switching points is 𝑠=3. Since by applying this assumption we reach to the exact results in which the new associated control 𝜆(𝑡) 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 𝑁 are provided in Table 1. Moreover, the associated control 𝜆(𝑡), control 𝑢(𝑡), and optimal state 𝑦(𝑡) 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 Example 2.

Table 2 

Comparisons of the methods in evaluating the objective function 𝐽^∗.

Figure 3 

Approximate optimal associated control 𝜆(𝑡) history of Example 2 for 𝑁=16.

Figure 4 

Approximate optimal control 𝑢(𝑡) history of Example 2 for 𝑁=16.

Figure 5 

Approximate optimal state 𝑦(𝑡) history of Example 2 for 𝑁=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 𝑁.

Conflict of Interests

The authors declare that they do not have any conflict of interests in their submitted paper.

Acknowledgment

The authors thank the referee for his or her valuable comments and helpful suggestions which led to the improved version of the paper.