Application of Optimal HAM for Finding Feedback Control of Optimal Control Problems

Abstract

An optimal homotopy-analysis approach is described for Hamilton-Jacobi-Bellman equation (HJB) arising in nonlinear optimal control problems. This optimal approach contains at most three convergence-control parameters and is computationally rather efficient. A kind of averaged residual error is defined. By minimizing the averaged residual error, the optimal convergence-control parameters can be obtained. This optimal approach has general meanings and can be used to get fast convergent series solutions of different types of equations with strong nonlinearity. The closed-loop optimal control is obtained using the Bellman dynamic programming. Numerical examples are considered aiming to demonstrate the validity and applicability of the proposed techniques and to compare with the existing results.

1. Introduction

Many real-world problems in engineering, economics, and biomedicine can be formulated as optimal control problem of the form minimum𝐽=ℎ𝑥𝑡_𝑓,𝑡_𝑓+^𝑡𝑓_𝑡0𝑔(𝑥(𝑠)(𝑠))𝑥̇(𝑡)(𝑥(𝑡)(𝑡)),𝑥𝑡,𝑢,𝑠𝑑𝑠,subjectto=𝑓,𝑢,𝑡₀=𝑥₀,𝑥𝑡_𝑓=𝑥_𝑓, (1) where 𝑢(⋅) is the control function, 𝑥(⋅) is the state function, constant parameters, 𝑥₀ and 𝑥_𝑓 are the initial and the final states, ℎ(⋅) is the salvage value or scrap value function, 𝑓(⋅) is the vector-valued transition function, 𝑔(⋅) is the integrand function, and 𝑠 is the dummy variable of integration rather than 𝑡.

There are various efficient methods such as those reported in [1–3] for the computation of open-loop optimal controls. However, feedback controls are much preferred in many engineering applications. In order to determine the optimal control law, there is an approach using dynamic programming [4]. This approach leads to the HJB equation that is hard to solve in most cases. An excellent literature review on the methods for solving the HJB equation is provided in [5], where a successive Galerkin approximation (SGA) method is also considered. In the SGA, a sequence of generalized HJB equations is solved iteratively to obtain a sequence of approximations approaching the solution of HJB equation. However, the proposed sequence may converge very slowly or even diverge. Another approach is to treat the problem with a measure theory approach [6]. This changes the nonlinear OCP to a linear programming and gives a piecewise constant control law. In the past two decades, the indirect methods have been extensively developed. It is well known that the nonlinear OCP leads to a nonlinear two-point boundary value problem (TPBVP) or a Hamilton-Jacobi-Bellman partial differential equation. Many recent researches have been devoted to solve these two problems.

In particular, we will derive the fundamental first-order partial differential equation obeyed by the optimal value function, known as the HJB equation. This shift in our attention, moreover, will lead us to a different form for the optimal value of the control vector, namely, the feedback or closed-loop form of the control. This form of the optimal control typically gives the optimal value of the control vector as a function of the current date, the current state, and the parameters of the control problem. In contrast, the form of the optimal control vector derived via the necessary conditions of optimal control theory is termed open loop and in general gives the optimal value of the control vector as a function of the independent variable time, the parameters, and the initial and/or terminal values of the planning horizon and the state vector. Essentially, the feedback form of the optimal control is a decision rule, for it gives the optimal value of the control for any current period and any admissible state in the current period that may arise. In contrast, the open-loop form of the optimal control is a curve, for it gives the optimal values of the control as the independent variable time varies over the planning horizon. Thus, it is required to solve an HJB equation by numerical-analytic-approximate method. In [7], the authors used the basic ADM for an example of HJB equation and then finding optimal control signal. In [8–10], standard HPM, piecewise HPM and differential transform method (DTM) were used for finding feedback control of optimal control problems.

One of the semiexact methods for solving linear and nonlinear equation which does not need small/large parameters is homotopy analysis method (HAM), first proposed by Liao in 1992 [10–13]. Since Liao’s book [14] for the homotopy analysis method was published in 2003, more and more researchers have been successfully applying this method to various nonlinear problems in science and engineering, such as the viscous flows of non-Newtonian fluids [15], the KdV-type equations [16], and finance problems [17].

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. Thus, through HAM, explicit analytic solutions of nonlinear problems are possible.

However, as illustrated in [18], the usual HAM has only one convergence-control parameter 𝑐₀, and it is a pity that curves for convergence-control parameter (i.e., 𝑐₀-curves) cannot tell us which value of 𝑐₀∈ℝ gives the fastest convergent series. Recently, to overcome this shortcoming, Liao [18] proposed an optimal HAM with more than one convergence-control parameter. In this optimal method, Liao introduced the so-called averaged residual error to get the optimal convergence-control parameters efficiently. In general, the optimal HAM can greatly modify the convergence of homotopy-series solution.

The results obtained in this paper show that the solutions given by the optimal HAM give much better approximations and convergence much faster than those given by the usual HAM. In this paper, we employed a new proposal of Liao [18, 19], namely, optimal homotopy-analysis approach to solve an HJB equation obtained from (2)–(4). We use here special deformation functions, which are determined completely by only one characteristic parameter |𝑐₁|<1 and |𝑐₂|<1, respectively. In this way, there exist at most only three convergence-control parameters 𝑐₀,𝑐₁, and 𝑐₂ at any order of approximations.

2. Nonlinear Time-Variant HJB Equation

Consider a process described by the following system of nonlinear differential equations, which is called the equation of motion, on a fixed interval [𝑡₀,𝑡_𝑓]𝑥̇(𝑡)=𝑓(𝑥(𝑡),𝑢(𝑡),𝑡), (2) where 𝑥∈ℝ^𝑛 is the state vector. Let 𝑈(⋅)∈ℝ^𝑛 be a closed set. A piecewise continuous function 𝑢∶[𝑡₀,𝑡_𝑓]→ℝ^𝑛 is said to be an admissible control if 𝑢(𝑡)∈𝑈(⋅). Let 𝒰 be the class of such admissible controls. The function 𝑓 is a real-valued continuously differentiable function and has continuous first partial derivative with respect to 𝑥. The initial condition for (2) is𝑥𝑡₀=𝑥₀. (3) Along with this controlled process, we have a cost functional of the form 𝐽𝑥₀,𝑢=Φ𝑥𝑡_𝑓,𝑡_𝑓+^𝑡𝑓_𝑡0𝐿(𝑥(𝜏),𝑢(𝜏),𝜏)𝑑𝜏. (4) Here, 𝐿(𝑥(𝜏),𝑢(𝜏),𝜏) is the running cost, and Φ(𝑥(𝑡_𝑓),𝑡_𝑓) is the terminal cost. This cost functional depends on the initial position (𝑥₀,𝑡₀) and the choice of control 𝑢(⋅). The optimization problem is therefore to minimize 𝐽(𝑥₀,𝑢), for each (𝑥₀,𝑢), over all controls 𝑢(𝑡)∈𝑈(⋅). The pair (𝑥₀,𝑢) which achieves this minimum is called an optimal control. In fact, the optimization problem with performance index as in (4) is called a Bolza problem.

In dynamic programming, the minimum value of the performance index is considered as a function of this initial point. This function is called the value function. Whenever the value function is differentiable, it satisfies a nonlinear first-order hyperbolic partial differential equation called the partial differential equation of dynamic programming. This equation is used for constructing a nonlinear optimal feedback control law. If we consider a family of optimization problems with differential initial condition (𝑥,𝑡), we consider the dependence of the value of these optimization problems on their initial conditions. Thus, a value function is defined by 𝑉(𝑥,𝑡)=𝐽𝑥,𝑢^∗=min_𝑢∈𝒰Φ𝑥𝑡_𝑓,𝑡_𝑓+^𝑡𝑓_𝑡𝐿(𝑥(𝜏),𝑢(𝜏),𝜏)𝑑𝜏. (5)

Theorem 1 .. Let (𝑥,𝑡) be any interior point of the set {(𝑥,𝑡)∈ℝ^𝑛+1|𝒰≠0} at which the function 𝑉(𝑥,𝑡) is differentiable. Function 𝑉(𝑥,𝑡) satisfies the partial differential inequality 𝜕𝑉+𝜕𝑡𝜕𝑉𝑓()()𝜕𝑥𝑥,𝑢,𝑡+𝐿𝑥,𝑢,𝑡≥0, (6) for all 𝑢∈𝒰. If there is an optimal control 𝑢^∗∈𝒰, then the partial differential equation 𝜕𝑉𝜕𝑡+min_𝑢∈𝑈𝜕𝑉𝑓()(𝜕𝑥𝑥,𝑢,𝑡+𝐿𝑥,𝑢,𝑡)=0 (7) is satisfied [20].

This nonlinear time-variant differential equation is called “HJB equation.” We have the following boundary condition: 𝐽^∗𝑥𝑡_𝑓,𝑡_𝑓=𝑉𝑥𝑡_𝑓,𝑡_𝑓=Φ𝑥𝑡_𝑓,𝑡_𝑓, (8) and by introducing the Hamiltonian function 𝐻𝑥,𝑢,𝑉_𝑥,𝑡=𝜕𝑉𝑓()(),𝜕𝑥𝑥,𝑢,𝑡+𝐿𝑥,𝑢,𝑡 (9) we have 𝐻𝑥,𝑢^∗,𝑉_𝑥,𝑡=min_𝑢∈𝑈𝐻𝑥,𝑢,𝑉_𝑥,𝑡. (10) Therefore, by substitution of Hamiltonian function (10) in (7), we have 𝜕𝑉𝜕𝑡+𝐻𝑥,𝑢^∗,𝑉_𝑥,𝑡=0. (11)

3. Basic Idea of Optimal Homotopy Analysis Method

In this section, we give a brief description of optimal homotopy-analysis approaches in general cases. For a general nonlinear differential/integral equation as follows: 𝑁(𝑈(𝑟,𝑡))=0, (12) where 𝑈(𝑟,𝑡) is an unknown function, 𝑟 and 𝑡 denote, respectively, spatial and temporal independent variables, we can always choose a proper initial guess 𝑈₀(𝑟,𝑡) and a proper auxiliary linear operator 𝐿 to construct the so-called zeroth-order deformation equation 1−𝐵_𝑛′𝑞;𝑏𝐿𝜙𝑟,𝑡;𝑞−𝑈₀(𝑟,𝑡)=𝑐₀𝐴_𝑚′𝑞,𝑎𝑁𝜙𝑟,𝑡;𝑞,𝑞∈[0,1], (13) where 𝐴_𝑚′(𝑞,𝑎) is a deformation function with 𝑚^′(𝑚^′≥0) unknown convergence-control parameters denoted by 𝑐₁,𝑐₂,…,𝑐^′_𝑚, and 𝐵_𝑛′(𝑞,𝑏) is a deformation function with 𝑛^′(𝑛^′≥0) unknown convergence-control parameters denoted by 𝑐_𝑚′+1,𝑐_𝑚′+2,…,𝑐_{𝑚′+𝑛′}, respectively. Here, 𝑎=𝑐₁,𝑐₂,…,𝑐_𝑚′,𝑏=𝑐_𝑚′+1,𝑐_𝑚′+2,…,𝑐_{𝑚′+𝑛′}. (14) Put 𝑙=𝑚^′+𝑛^′. Then, we have 𝑙+1 unknown convergence-control parameters 𝑐₀,𝑐₁,𝑐₂,…,𝑐_𝑙. Assuming that the initial guess 𝑈₀(𝑟,𝑡), the auxiliary linear operator 𝐿, and the 𝑙+1 convergence-control parameters 𝑐₁,𝑐₂,…,𝑐_𝑙 are so properly chosen that the Taylor series 𝜙𝑟,𝑡;𝑞=𝑈₀(𝑟,𝑡)+^+∞_𝑛=1𝑈_𝑛(𝑟,𝑡)𝑞^𝑛 (15) converges at 𝑞=1, we have the homotopy-series solution 𝑈(𝑟,𝑡)=𝑈₀(𝑟,𝑡)+^+∞_𝑛=1𝑈_𝑛(𝑟,𝑡). (16) Substituting the series (15) into the zeroth-order deformation equation (13) and then equating the coefficients of the like power of the embedding parameter 𝑞, we have the high-order deformation equation. Consider𝐿𝑈_𝑚(𝑟,𝑡)−^𝑚−1_𝑘=1𝜎_𝑚−𝑘(𝑏)𝑈_𝑘(𝑟,𝑡)=𝑐₀^𝑚−1_𝑘=0𝜇_𝑚−𝑘(𝑎)𝛿_𝑘(𝑟,𝑡), (17) where 𝛿_𝑘()=1𝑟,𝑡𝜕𝑘!𝜕𝑞^𝑘𝑘𝑁^+∞_𝑛=1𝑈_𝑛(𝑟,𝑡)𝑞^𝑛_𝑞=0, (18) and 𝜇_𝑘(𝑎)and𝜎_𝑘(𝑏) are coefficients of the Taylor series 𝐴_𝑚′𝑞;𝑎=^+∞_𝑘=1𝜇_𝑘(𝑎)𝑞^𝑘,𝐵_𝑛′𝑞;𝑏=^+∞_𝑘=1𝜎_𝑘(𝑏)𝑞^𝑘. (19) Let 𝐴_𝑚′(𝑞;𝑎),𝐵_𝑛′(𝑞;𝑏) be complex functions analytic in the region |𝑞|≤1, which satisfy 𝐴_𝑚′(0;𝑎)=𝐵_𝑛′(0;𝑎)=0,𝐴_𝑚′(1;𝑎)=𝐵_𝑛′(1;𝑎)=1. (20) The special solution 𝑈^∗_𝑚(𝑟,𝑡) of (17) is given by 𝑈^∗_𝑚(𝑟,𝑡)=^𝑚−1_𝑘=1𝜎_𝑚−𝑘(𝑏)𝑈_𝑘(𝑟,𝑡)+𝑐₀^𝑚−1_𝑘=1𝜇_𝑚−𝑘(𝑎)𝑆_𝑘(𝑟,𝑡), (21) where 𝑆_𝑘(𝑟,𝑡)=𝐿⁻¹𝛿_𝑘(𝑟,𝑡), (22) and 𝐿⁻¹ is the inverse operator of 𝐿.

Now, we define a kind of averaged residual error [18] 𝐸_𝑚=1𝑘^𝑘_𝑗=0𝑁^𝑚_𝑖=0𝑈_𝑖𝑗Δ𝑡². (23)The value of 𝑘 is depends to problem, for Examples 2 and 3, we have Δ𝑡=1/10,𝑘=10, and Δ𝑡=1/20,𝑘=20, respectively.

Note that 𝐸_𝑚 contains 𝑙+1 unknown convergence-control parameters 𝑐₀,𝑐₁,…,𝑐_𝑙, whose optimal values are given by the minimum of 𝐸_𝑚, corresponding to a set of 𝑙+1 nonlinear algebraic equations 𝜕𝐸_𝑚𝜕𝑐_𝑗=0,𝑗=0,1,…,𝑙. (24) So, the above approach is called the 𝑙+1-parameter optimal homotopy-analysis approach. In general, the above-mentioned optimal homotopy-analysis approaches can greatly modify the convergence of homotopy-series solution. And the optimal homotopy-analysis approaches with one or two unknown convergence-control parameters are strongly suggested; an optimal approach with too many unknown convergence-control parameters is not efficient computationally.

4. Application

To illustrate the effectiveness of the optimal homotopy analysis method, we will consider two examples. Comparisons are made between the exact solution and optimal homotopy-analysis approaches.

Example 2. Consider the following purely mathematical optimal control problem: 𝑥̇=𝑥(𝑡)+𝑢(𝑡). (25) Suppose that we consider the following cost function for this system: 𝐽=𝑥²𝑡_𝑓+^𝑡𝑓₀𝑢²(𝑡)𝑑𝑡. (26) The Hamiltonian function will be 𝐻𝑥,𝑢,𝑉_𝑥,𝑡=𝑢²(𝑡)+𝑉_𝑥(𝑥,𝑡)[𝑥+𝑢]. (27) Since the Hamiltonian, which as you know is the right-hand-side expression to be minimum, is a convex function of the state and control variables, a solution of the necessary conditions is a solution of the optimal control. Thus, the necessary and sufficient condition for the above minimization problem in the HJB equation is given by 𝜕𝐻(𝑡)𝜕𝑢=2𝑢+𝑉_𝑥(𝑥,𝑡)=0, (28) which is easily solved to get 𝑢^∗(𝑡)1=−2𝑉_𝑥(𝑥,𝑡). (29) Because 𝜕²𝐻/𝜕𝑢²=2>0, 𝑢^∗ is a minimum and acceptable. Now, by substitution of 𝑢^∗ in HJB equation, we have the following: −𝑉_𝑡1=−4𝑉²_𝑥+𝑥𝑉_𝑥,𝑉𝑥𝑡_𝑓,𝑡_𝑓=𝑥²𝑡_𝑓. (30) In Chapter 19 of [21], the authors obtained the solution of the above HJB equation in the form 𝑉(𝑥,𝑡)=2𝑥²1+𝑒^{2(𝑡−𝑡𝑓)},𝑢^∗(𝑡)=−2𝑥1+𝑒^{2(𝑡−𝑡𝑓)}, (31) for simplicity, let 𝑢^∗(𝑥,𝑡)=𝑘(𝑡)𝑥(𝑡), where 𝑘(𝑡)=−21+𝑒^{2(𝑡−𝑡𝑓)}. (32) To solve (30) by means of optimal homotopy analysis method, we choose the initial approximation, final time, and auxiliary linear operator 𝑉₀(𝑥,𝑡)=𝑥²,𝑡_𝑓𝐿𝜙==1,𝑥,𝑡;𝑞𝜕𝜙𝑥,𝑡;𝑞.𝜕𝑡 (33) There are an infinite number of deformation functions satisfying the properties (19). For the sake of computation efficiency, we use here the following one-parameter deformation functions: 𝐴₁𝑞;𝑐₁=^+∞_𝑚=1𝜇_𝑚𝑐₁𝑞^𝑚,𝐵₁𝑞;𝑐₂=^+∞_𝑚=1𝜎_𝑚𝑐₂𝑞^𝑚, (34) where |𝑐₁|<1 and |𝑐₂|<1 are constants, called the convergence-control parameter, and 𝜇_𝑚𝑐₁=1−𝑐₁𝑐₁^𝑚−1𝜎,𝑚≥1,_𝑚𝑐₂=1−𝑐₂𝑐₂^𝑚−1,𝑚≥1. (35) The different values of 𝑐₁ give different paths of 𝐴₁(𝑞;𝑐₁), as shown in Figure 1.

Figure 1 

Deformation function 𝐴1(𝑞;𝑐1) defined by (34) and (35).

Let 𝑞∈[0,1] denote the embedding parameter, 𝑐₀≠0 an auxiliary parameter, called the convergence-control parameter, and 𝜙(𝑥,𝑡;𝑞) a kind of continuous mapping of 𝑉(𝑥,𝑡), respectively. Obviously, when 𝑞=0 and 𝑞=1, it holds 𝜙(𝑥,𝑡;0)=𝑉₀(𝑥,𝑡),𝜙(𝑥,𝑡;1)=𝑉(𝑥,𝑡), (36) respectively. Thus, as 𝑞 increases from 0 to 1, the solutions 𝜙(𝑥,𝑡;𝑞)vary from the initial guesses 𝑉₀(𝑥,𝑡) to the solutions 𝑉(𝑥,𝑡). We construct the zeroth-order deformation equation: 1−𝐵₁𝑞;𝑐₂𝐿𝜙𝑥,𝑡;𝑞−𝑉₀(𝑥,𝑡)=𝑐₀𝐴₁𝑞,𝑐₁𝑁𝜙𝑥,𝑡;𝑞,𝑞∈[0,1], (37) subject to the boundary condition 𝜙(𝑥(1),1;0)=𝑥²(1). (38) As mentioned, we have the 𝑚th-order deformation equation 𝐿𝑉_𝑚(𝑥,𝑡)−^𝑚−1_𝑘=1𝜎_𝑚−𝑘𝑐₂V_𝑘(𝑥,𝑡)=𝑐₀^𝑚−1_𝑘=0𝜇_𝑚−𝑘𝑐₁𝛿_𝑘(𝑥,𝑡), (39) where 𝛿_𝑘(𝑥,𝑡)=𝜕𝑉_𝑘−1𝜕𝑡4^𝑘_𝑗=0𝜕𝑉_𝑗𝜕𝑥𝜕𝑉_𝑘−𝑗𝜕𝑥+𝑥𝜕𝑉_𝑘𝜕𝑥 (40) and the coefficients 𝜇_𝑘(𝑐₁)and𝜎_𝑘(𝑐₂) are defined by (35). Let 𝑉^∗(𝑥,𝑡) denote a special solution of (39) and 𝐿⁻¹ the inverse operator of 𝐿, respectively. We have 𝑉^∗_𝑚(𝑥,𝑡)=∑^𝑚−1_𝑘=1𝜎_𝑚−𝑘𝑐₂𝑉_𝑘(𝑥,𝑡)+𝑐∑₀^𝑚−1_𝑘=0𝜇_𝑚−𝑘𝑐₁𝐿⁻¹𝛿_𝑘(𝑥,𝑡). (41) We now successively obtain 𝑉₁(𝑥,𝑡)=−𝑐₀−1+𝑐₁(𝑡−1)𝑥²,𝑉₂(𝑥,𝑡)=𝑐₀𝑡𝑥²−𝑐₀𝑥²−𝑐₀𝑐₂𝑡𝑥²+𝑐₀𝑐₂𝑥²+𝑐₀𝑐₂𝑐₁𝑡𝑥²−𝑐₀𝑐₂𝑐₁𝑥²−𝑐₀𝑐²₁𝑡𝑥²+𝑐₀𝑐²₁𝑥²+𝑐²₀𝑡𝑥²−𝑐²₀𝑥²−2𝑐²₀𝑐₁𝑡𝑥₂+2𝑐²₀𝑐₁𝑥²+𝑐²₀𝑐²₁𝑡𝑥²−𝑐²₀𝑐²₁𝑥²,𝑉₃(𝑥,𝑡)=𝑐₀𝑥²𝑡−𝑐₀𝑥²−𝑐₀𝑐₂𝑥²𝑡+𝑐₀𝑥²𝑐₂+𝑐₀𝑥²𝑐₂𝑐₁𝑡−𝑐₀𝑥²𝑐₂𝑐₁−𝑐₀𝑥²𝑐²₁𝑡+𝑐₀𝑥²𝑐²₁+𝑐²₀𝑥²𝑡−𝑐²₀𝑥²−2𝑐²₀𝑥²𝑐₁𝑡+2𝑐²₀𝑥²𝑐₁+𝑐²₀𝑥²𝑐²₁𝑡−𝑐²₀𝑥²𝑐²₁,⋮ (42)

It should be emphasized that 𝑉_𝑚(𝑥,𝑡) contains at most three unknown convergence-control parameters 𝑐₀,𝑐₁, and 𝑐₂, which determine the convergence region and rate of the homotopy-series solution. Obviously, if the convergence-control parameters 𝑐₀,𝑐₁, and 𝑐₂ are properly chosen, the homotopy-series solution may converge fast. So, we should find out the good enough values of 𝑐₀,𝑐₁, and 𝑐₂, so that the homotopy-series solution converges fast enough.

4.1. Comparisons of Different Approaches for Example 2

In this section, we will give optimal homotopy-analysis approaches with different numbers of unknown convergence-control parameters and compare them in detail.

4.1.1. Optimal 𝑐₀ in Case of 𝑐₁=𝑐₂=0

In this case, only one convergence-control parameter 𝑐₀ is unknown. For given order of approximation 𝑚, the optimal value of 𝑐₀ is given by the minimum of 𝐸_𝑚, corresponding to a nonlinear algebraic equation 𝑑𝐸_𝑚/𝑑𝑐₀=0.

The curves of 𝐸_𝑚 versus 𝑐₀ at a different order of approximation 𝑚=5,7, and 9 are shown in Figure 2. In case of 𝑐₁=𝑐₂=0,  there is only one unknown convergence-control parameter 𝑐₀; thus the optimal value of 𝑐₀ is determined by the minimum of 𝐸₉, corresponding to the nonlinear algebraic equation 𝑑𝐸₉/𝑑𝑐₀=0. According to Table 1, the value minimum of 𝐸₉ is obtained at −0.8296. It should be noted that, when we fixed 𝑐₀=−1 and 𝑐₁=𝑐₂=0, automatically the problem turns into the so-called homotopy perturbation method (HPM). So, the HPM is only a special case of the HAM when 𝑐₀=−1.

Table 1 

Comparison of averaged residual error given by different 𝑐₀ in case of 𝑐₁=𝑐₂=0 and 𝑥=1.

Figure 2 

Averaged residual error 𝐸𝑚 versus 𝑐0 in case of 𝑐1=𝑐2=0.

In case of 𝑐₁=𝑐₂=0, one has the plain deformation functions 𝐴₁(𝑞,𝑐₁)=𝐵₁(𝑞,𝑐₂)=𝑞, which is exactly the same as that used in traditional HAM.

4.1.2. Optimal 𝑐₁=𝑐₂ in Case of 𝑐₀=−1

In this case, 𝐸₉ has the value minimum of 2.8203×10⁻⁶ at the optimal point 𝑐₁=𝑐₂=0.2458. This illustrates that the second one-parameter optimal homotopy-analysis approach is as good as the first one mentioned in Section 4.1.1. (See Table 2).

Table 2 

Comparison of averaged residual error given by different 𝑐₁=𝑐₂ in case of 𝑐₀=−1 and 𝑥=1.

4.1.3. Optimal 𝑐₀ in Case of 𝑐₁=𝑐₂

In this case, we obtain the unknown convergence-control parameters 𝑐₀ and 𝑐₁ in case of 𝑐₂=𝑐₁. The corresponding residual error 𝐸₉ is now a function of both 𝑐₀ and 𝑐₁, which has the value minimum of 1.4099×10⁻⁶ at the optimal point 𝑐₀=−1.0012, 𝑐₁=𝑐₂=0.0398. (See Table 3).

Table 3 

Comparison of averaged residual error given by different 𝑐₀ in case of 𝑐₁=𝑐₂.

4.1.4. Optimal 𝑐₀,𝑐₁, and 𝑐₂ in Case of 𝑐₁≠𝑐₂

Let us consider the optimal approach with the three unknown convergence-control parameters. The corresponding 𝐸₉ is now a function of 𝑐₀,𝑐₁, and 𝑐₂, which has the minimum 1.7716×10⁻⁵ at the optimal values 𝑐₀=−0.6932, 𝑐₁=−0.2440, and 𝑐₂=−0.1039. The corresponding homotopy-series solution converges much faster than homotopy-series solution in case of 𝑐₀=−1 and 𝑐₁=𝑐₂=0, as shown in Table 4.

Table 4 

Comparison of averaged residual error given by different approaches.

Based on the above calculations for this equation, we found that all optimal homotopy-analysis approaches can give much better approximations, which converge much faster than those without optimal convergence-control parameters in case of 𝑐₀=−1 and 𝑐₁=𝑐₂=0, as shown in Tables 1–4.

The absolute error of the 20th-order approximate solution and analytic solution (31) with 𝑐₀=−1, 𝑐₁=𝑐₂=0, and 𝑐₀=−0.8296, 𝑐₁=𝑐₂=0is shown in Figures 3 and 4, respectively. Comparison of the analytic solution (31) with the optimal homotopy analysis method solution by 𝑐₀=−0.8296 and 𝑐₁=𝑐₂=0 is given in Figure 5. The absolute errors |𝑘_analytic−𝑘_approximate| have been calculated for different cases of 𝑐₀,𝑐₁, and 𝑐₂ in Table 5.

Table 5 

The absolute error of approximate solution and exact solution of 𝑘(𝑡) by 20th-order approximate.

Figure 3 

|𝑘analytic−𝑘approximate| (𝑚=20), with 𝑐0=−1 and 𝑐1=𝑐2=0.

Figure 4 

|𝑘analytic−𝑘approximate| (𝑚=20), with 𝑐0=−0.8296 and 𝑐1=𝑐2=0.

Figure 5 

Comparison of the analytic solution (31) with the approximate solution of 𝑘(𝑡) (𝑚=20), by 𝑐0=−0.8296 and 𝑐1=𝑐2=0.

Example 3. Consider a single-input scalar system as follows: ̇(𝑡)(𝑡),1𝑥=−2𝑥+𝑢𝐽=2𝑥²(1)+12¹₀𝑥²(𝑡)+𝑢²(𝑡)𝑑𝑡. (43) The corresponding Hamiltonian function will be 𝐻𝑥,𝑢,𝑉_𝑥=1,𝑡2𝑥²(𝑡)+12𝑢²(𝑡)+𝜕𝑉(𝑥,𝑡)[(𝑡)(𝑡.𝜕𝑥−2𝑥+𝑢)] (44) For finding 𝑢^∗, we have 𝜕𝐻(𝑡)+𝜕𝑢=𝑢𝜕𝑉𝜕𝑥=0. (45) Therefore, we obtain 𝑢^∗(𝑡)=−𝜕𝑉.𝜕𝑥 (46) Because 𝜕²𝐻/𝜕𝑢²=1>0, 𝑢^∗ is a minimum and acceptable. Now, by substitution of 𝑢^∗ in HJB equation, we have the following: −𝜕𝑉=1𝜕𝑡2𝑥²(𝑡)+12−𝜕𝑉(𝑥,𝑡)𝜕𝑥²+𝜕𝑉(𝑥,𝑡)(𝑡)−()𝜕𝑥−2𝑥𝜕𝑉𝑥,𝑡,𝑉(𝑥(1))=1𝜕𝑥,12𝑥²(1) (47) or 𝜕𝑉1𝜕𝑡=−2𝑥²(𝑡)+12𝜕𝑉(𝑥,𝑡)𝜕𝑥²()+2𝑥𝜕𝑉𝑥,𝑡,𝑉(𝑥(1))=1𝜕𝑥,12𝑥²(1). (48) Note that, by (29), the optimal control law is 𝑢^∗(𝑡)=−𝑘(𝑡)𝑥(𝑡)=𝑓(𝑡)𝑥(𝑡), where the exact solution of 𝑘(𝑡) is [2] 𝑘(𝑡)=√5cosh√5(1−𝑡)−sinh√5(1−𝑡)√5cosh√5(1−𝑡)+3sinh√5().1−𝑡 (49) To solve (48) by means of optimal homotopy analysis method, we choose the initial approximation and auxiliary linear operator 𝑉₀()=1𝑥,𝑡2𝑥²,𝐿𝜙𝑥,𝑡;𝑞=𝜕𝜙𝑥,𝑡;𝑞.𝜕𝑡 (50) We construct the zeroth-order deformation equation 1−𝐵₁𝑞;𝑐₂𝐿𝜙𝑥,𝑡;𝑞−𝑉₀(𝑥,𝑡)=𝑐₀𝐴₁𝑞,𝑐₁𝑁𝜙𝑥,𝑡;𝑞,𝑞∈[0,1], (51) subject to the boundary condition 𝜙(𝑥(1))=1,1;02𝑥²(1). (52) As mentioned, we have the 𝑚th-order deformation equation 𝐿𝑉_𝑚(𝑥,𝑡)−^𝑚−1_𝑘=1𝜎_𝑚−𝑘𝑐₂𝑉_𝑘(𝑥,𝑡)=𝑐₀^𝑚−1_𝑘=0𝜇_𝑚−𝑘𝑐₁𝛿_𝑘(𝑥,𝑡), (53) where 𝛿_𝑘(𝑥,𝑡)=𝜕𝑉_𝑘−1𝜕𝑡2^𝑘_𝑗=0𝜕𝑉_𝑗𝜕𝑥𝜕𝑉_𝑘−𝑗𝜕𝑥−2𝑥𝜕𝑉_𝑘+1𝜕𝑥2𝑥² (54) and the coefficients 𝜇_𝑘(𝑐₁)and𝜎_𝑘(𝑐₂) are defined by (35). Let 𝑉^∗(𝑥,𝑡) denote a special solution of (53) and 𝐿⁻¹ the inverse operator of 𝐿, respectively. We have 𝑉^∗_𝑚(𝑥,𝑡)=^𝑚−1_𝑘=1𝜎_𝑚−𝑘𝑐₂𝑉_𝑘(𝑥,𝑡)+𝑐₀^𝑚−1_𝑘=0𝜇_𝑚−𝑘𝑐₁𝐿⁻¹𝛿_𝑘(𝑥,𝑡). (55) We now successively obtain 𝑉₁(𝑥,𝑡)=2𝑐₀−1+𝑐₁(𝑡−1)𝑥²,𝑉₂(𝑥,𝑡)=−2𝑐₀𝑡𝑥²+2𝑐₀𝑥²+2𝑐₀𝑐₂𝑡𝑥²−2𝑐₀𝑐₂𝑥²−2𝑐₀𝑐₂𝑐₁𝑡𝑥²+2𝑐₀𝑐₂𝑐₁𝑥²+6𝑐₀𝑡²𝑥²+8𝑐²₀𝑥²−12𝑐²₀𝑐₁𝑡²𝑥²−16𝑐²₀𝑐₁𝑥²+6𝑐²₀𝑐²₁𝑡²𝑥²+8𝑐²₀𝑐²₁𝑥²+2𝑐₀𝑐²₁𝑡𝑥²−2𝑐₀𝑐²₁𝑥²−14𝑐²₀𝑡𝑥²+28𝑐²₀𝑐₁𝑡𝑥²−14𝑐²₀𝑐²₁𝑡𝑥²,⋮ (56)

4.2. Comparisons of Different Approaches for Example 3

In this section, we will give optimal homotopy-analysis approaches with different numbers of unknown convergence-control parameters and compare them in detail.

4.2.1. Optimal 𝑐₀ in Case of 𝑐₁=𝑐₂=0

In this case, only one convergence-control parameter 𝑐₀ is unknown. We solved this example by standard HAM, and we investigated the influence of 𝑐₀ by means of the so-called 𝑐₀-curves. As pointed out by Liao [14], the valid region of 𝑐₀ is a horizontal line segment. So, we can just determine the possible valid region of 𝑐₀ as shown in Figure 6. However, usually the 𝑐₀-curves cannot tell us which value of 𝑐₀ gives the fastest convergent series.

Figure 6 

The 𝑐0-curve of the 9th-order approximation in case of 𝑐1=𝑐2=0.

To overcome this shortcoming, we determined the possible optimal value of 𝑐₀ by the minimum of averaged residual error 𝐸₉, corresponding to the nonlinear algebraic equation 𝐸^′₉=0. According to Table 6, the value minimum of 𝐸₉ is obtained at 𝑐₀=−0.44032. The curves of 𝐸_𝑚 versus 𝑐₀ at different order of approximation 𝑚=5,7, and 9 are shown in Figure 7.

Table 6 

Comparison of averaged residual error given by different 𝑐₀ in case of 𝑐₁=𝑐₂=0.

Figure 7 

Averaged residual error 𝐸𝑚 versus 𝑐0 in case of 𝑐1=𝑐2=0.

4.2.2. Optimal 𝑐₁=𝑐₂ in Case of 𝑐₀=−1

In this case, 𝐸₉ has the value minimum of 8.2039×10⁻² at the optimal point 𝑐₁=𝑐₂=0.0228. Our calculations showed that the corresponding 𝐸₉ in case of 𝑐₀=−1, 𝑐₁=𝑐₂=0 has the minimum 21.865. This illustrates that the second one-parameter optimal homotopy-analysis approach is not obviously better than the first one mentioned in Section 4.2.1. But this approach can give much better approximation than those without optimal convergence-control parameters in case of 𝑐₀=−1 and 𝑐₁=𝑐₂=0.

4.2.3. Optimal 𝑐₀ in Case of 𝑐₁=𝑐₂

In this case, we obtain the unknown convergence-control parameters 𝑐₀ and 𝑐₁ in case of 𝑐₂=𝑐₁. The corresponding residual error 𝐸₉ is now a function of both 𝑐₀ and 𝑐₁, which has the value minimum of 4.6895×10⁻⁴ at the optimal point 𝑐₀=−0.6084, 𝑐₁=𝑐₂=0.1611.

4.2.4. Optimal 𝑐₀,𝑐₁, and 𝑐₂in Case of 𝑐₁≠𝑐₂

Let us consider the optimal approach with the three unknown convergence-control parameters. The corresponding 𝐸₉ is now a function of 𝑐₀,𝑐₁, and 𝑐₂, which has the minimum 4.94×10⁻³ at the optimal values 𝑐₀=−0.1836, 𝑐₁=−0.2070, and 𝑐₂=−0.4554.

Based on the above calculations for this equation, we found that all optimal homotopy-analysis approaches can give much better approximations, which converge much faster than those without optimal convergence-control parameters in case of 𝑐₀=−1 and 𝑐₁=𝑐₂=0. But in this example, the first one-parameter optimal homotopy-analysis approach is better than other cases.

The absolute error and comparison of the 15th-order approximate solution and analytic solution (49) with 𝑐₀=−0.44032, 𝑐₁=𝑐₂ = 0 are shown in Figures 8 and 9, respectively. The absolute errors |𝑓_analytic−𝑓_approximate| have been calculated for different cases of 𝑐₀,𝑐₁, and 𝑐₂ in Table 7.

Table 7 

The absolute error of approximate solution and exact solution of 𝑓(𝑡) by 15th-order approximate.

Figure 8 

|𝑓analytic−𝑓approximate| with 𝑐0=−0.44032 and 𝑐1=𝑐2=0.

Figure 9 

Comparison of the exact solution with the approximate solution of 𝑓(𝑡) by 𝑐0=−0.44032 and 𝑐1=𝑐2=0.

5. Conclusion

In this paper, we have successfully developed an optimal homotopy-analysis approach for solving HJB equation arising in nonlinear optimal control problems. With the deformation functions defined, the method used here contains three auxiliary convergence-control parameters, which provide us with a simple way to adjust and control the convergence region of the solution. By minimizing the averaged square residual error, the optimal convergence-control parameters 𝑐₀,𝑐₁, and 𝑐₂ can be obtained.