Feedback Control Method Using Haar Wavelet Operational Matrices for Solving Optimal Control Problems

Abstract

Most of the direct methods solve optimal control problems with nonlinear programming solver. In this paper we propose a novel feedback control method for solving for solving affine control system, with quadratic cost functional, which makes use of only linear systems. This method is a numerical technique, which is based on the combination of Haar wavelet collocation method and successive Generalized Hamilton-Jacobi-Bellman equation. We formulate some new Haar wavelet operational matrices in order to manipulate Haar wavelet series. The proposed method has been applied to solve linear and nonlinear optimal control problems with infinite time horizon. The simulation results indicate that the accuracy of the control and cost can be improved by increasing the wavelet resolution.

1. Introduction

Optimal control is an important branch of mathematics and has been widely applied in a number of fields, including engineering, science, and economics. Although, the necessary and sufficient conditions for optimality have already been derived for 𝐻₂ and 𝐻_∞ optimal controls, they are only useful for finding analytical solutions for quite restricted cases. If we assume full-state knowledge, and if the optimal control problem is linear, then the optimal control is a linear feedback of the state, which is obtained by solving a matrix Riccati equation. However, if the system is nonlinear, then the optimal control is a state feedback function, which depends on the solution to a Hamilton-Jacobi-Bellman equation (HJB) or a Hamilton-Jacobi-Issac equation (HJI) for 𝐻₂ or 𝐻_∞ optimal control problem, respectively [1], and is usually difficult to solve analytically. Feng et al. [2] have solved an HJI equation iteratively by solving a sequence of HJB equation. In this paper, we are more concerned with approximate solution for HJB equation. Among numerous computational approach for solution of HJI equation, we refer in particular to [3–5]. Robustness of nonlinear state feedback is discussed in [6].

Broadly speaking, and in general, numerical methods for solving optimal control problem are divided into two categories: direct and indirect methods. The direct methods reduce optimal control problem to a nonlinear programming problem, by parameterizing or discretizing the infinite-dimensional optimal control problem, into finite-dimensional optimization problem. On the other hand, the indirect methods solve HJB equation or the first order necessary condition for optimality, which are obtained from Pontryagin minimum principle. Both these methods are important for solving optimal control problems; however, the difference between them is that the indirect methods are believed to yield more accurate result, whereas the direct methods tend to have better convergence properties. von Stryk and Bulirsch [7] have used both direct and indirect methods to solve optimal control problem for trajectory optimization in Apollo capsule. Beard et al. [8] have introduced Generalized Hamilton-Jacobi-Bellman equation to successively approximate solution of the HJB equation. Given an arbitrary stabilizing control law, their method can be used to improve the performance of the control. Moreover, Jaddu [9] has reported some numerical methods to solve unconstrained and constrained optimal control problems, by converting optimal control problems into quadratic programming problem. He has used a parameterization technique using the Chebyshev polynomials. Meanwhile, Beeler et al. [10] have performed a comparison study of five different methods for solving nonlinear control systems and studied the performance of the methods on several test problems. Park and Tsiotras [11] have proposed a successive wavelet collocation algorithm which used interpolating wavelets, to iteratively solve the Generalized Hamilton-Jacobi-Bellman equation and the corresponding optimal control law.

Wavelet basis that has compact support allows us to better represent functions with sharp spikes or edges than other bases. This property is advantageous in many applications in signal or image processing. In addition, the availability of fast transform makes it attractive as a computational tool. Numerical solutions of integral and differential equations have been discussed in many papers, which basically fall either in the class of spectral Galerkin and Collocation methods or finite element and finite difference methods.

Haar wavelet is the simplest orthogonal wavelet with a compact support. Chan and Hsiao [12] have used the Haar operational matrix method to solve lumped and distributed parameter systems. Hsiao and Wang [13] have solved optimal control of linear time-varying systems via Haar wavelets. Dai and Cochran Jr. [14] have considered a Haar wavelet technique to transform optimal control problems into nonlinear programming (NLP) parameters at collocation points. This NLP can be solved using nonlinear programming solver such as SNOPT.

In the present paper we have considered the method of Beard et al. [8] to successively approximate the solution of HJB equation. Instead of using the Galerkin method with polynomial basis, we have used collocation method with Haar wavelet basis to solve the Generalized Hamilton-Jacobi-Bellman equation. Galerkin method requires the computation of multidimensional integrals which makes the method impractical for higher order systems [15]. The main advantage of using collocation method in general is that computational burden of solving Generalized Hamilton-Jacobi-Bellman equation is reduced to matrix computation only. Our new successive Haar wavelet collocation method is used to solve linear and nonlinear optimal control problems. In the process of establishing the method we have to define new operational matrices of integration for a chosen stabilizing domain and new operational matrix for the product of two dimensions Haar wavelet functions.

2. Haar Wavelets

The orthogonal set of the Haar wavelets ℎ_𝑖(𝑥) is a group of square wave over the interval 𝑥∈[𝜏₁,𝜏₂) defined as follows: ℎ₀(𝑥)=1,𝜏₁≤𝑥<𝜏₂,ℎ0,elsewhere,₁(𝑥)=1,𝜏₁1≤𝑥<2𝜏₁+𝜏₂,1−1,2𝜏₁+𝜏₂≤𝑥<𝜏₂,0,elsewhere. (1) Other wavelets can be obtained by dilation and translation of the mother wavelet ℎ₁(𝑥). In general, ℎ_𝑖(𝑥)=ℎ₁(2^𝑗𝑥−𝑘), where 𝑖=2^𝑗+𝑘, 𝑗,𝑘∈𝑁∪{0}, and 0≤𝑘<2^𝑗.

Each 𝑓(𝑥)∈𝐿²([𝜏₁,𝜏₂)) can be expanded into Haar series of infinite terms: 𝑓(𝑥)=𝑐₀ℎ₀(𝑥)+𝑐₁ℎ₁(𝑥)+𝑐₂ℎ₂(𝑥)+⋯. (2) If 𝑓(𝑥) is approximated as piecewise constants then it can be decomposed as 𝑓(𝑥)=^𝑚−1_𝑖=0𝑐_𝑖ℎ_𝑖(𝑥), (3) where 𝑖=2^𝑗+𝑘, 𝑗=0,1,2,…,log₂𝑚, and 𝑘=0,1,2,…,2^𝑗−1.

The Haar coefficients that are 𝑐_𝑖=2^𝑗𝜏₂−𝜏₁^𝜏2_𝜏1𝑓(𝑥)ℎ_𝑖(𝑥)𝑑𝑥 (4) can be obtain by minimizing the integral square error ∫^𝜏2_𝜏1(𝑓(𝑥)−∑^𝑚−1_𝑖=0𝑐_𝑖ℎ_𝑖(𝑥))²𝑑𝑥.

The sum in (3) can be compactly written in the form 𝑓(𝑥)=𝑐^𝑇_𝑚𝐡_𝑚(𝑥), (5) where 𝑐^𝑇_𝑚=[𝑐₀𝑐₁⋯𝑐_𝑚−1] is called the coefficient vector and 𝐡_𝑚(𝑥)=[ℎ₀(𝑥)ℎ₁(𝑥)⋯ℎ_𝑚−1(𝑥)]^𝑇 is the Haar function vector.

At collocation points 𝑥_𝑗=(𝜏₁+((𝜏₂−𝜏₁)/2𝑚)(2𝑗−1)), 𝑗=1,2,3,…,𝑚, the Haar function vector can be expressed in matrix form as 𝐻_{𝑚𝑖,𝑗}=ℎ_𝑖𝑥_𝑗. (6) For instance, the fourth Haar wavelet matrix 𝐻₄ can be represented in matrix form as follows: 𝐻₄=111111−1−11−100001−1. (7)

3. Haar Wavelet Operational Matrices

The integration of ℎ_𝑖(𝑥) in the interval of [0,𝜏) can also be expanded into a Haar series, that is, ^𝑥₀𝐡_𝑚(𝑥)𝑑𝑥≅𝑃_𝑚𝐡_𝑚(𝑥), (8) where the 𝑚×𝑚 matrix 𝑃_𝑚 is called the operational matrix of integration obtain recursively as 𝑃_𝑚=12𝑚2𝑚𝑃_𝑚/2−𝜏𝐻_𝑚/2𝜏𝐻⁻¹_𝑚/20_𝑚/2,𝑃₁=𝜏2. (9) The formula in the interval of [0,1) was first given by Chen and Hsiao [12].

In order to solve nonlinear optimal control problem, it is essential to have the product of 𝐡(𝑥) and 𝐡^𝑇(𝑥). The product of two functions 𝑓(𝑥)=𝑐^𝑇𝐡(𝑥) and 𝑔(𝑥)=𝑑^𝑇𝐡(𝑥) can be expanded into a Haar series with a Haar coefficient matrix 𝑀_𝑚 as 𝑑^𝑇𝐡(𝑥)𝐡^𝑇(𝑥)𝑐=𝑑^𝑇𝑀_𝑚𝐡(𝑥), (10) where 𝑀_𝑚 is an 𝑚×𝑚 matrix referred to as the product operational matrix. It was first given by Hsiao and Wu [16] as 𝑀_𝑚=𝑀_𝑚/2𝐻_𝑚/2diag𝑐_𝑏diag𝑐_𝑏𝐻⁻¹_𝑚/2diag𝑐^𝑇_𝑎𝐻_𝑚/2, (11) where 𝑀₁=𝑐₀ and 𝑐_𝑎=[𝑐₀,…,𝑐_𝑚/2−1]^𝑇, 𝑐_𝑏=[𝑐_𝑚/2,…,𝑐_𝑚−1]^𝑇.

Two-dimensional Haar wavelets basis can be formed by taking a tensor product of 𝐡_𝑛(𝑥) and 𝐡_𝑚(𝑥). Let the basis be {ℎ_𝑖(𝑥₁)ℎ_𝑗(𝑥₂)}, 𝑖=1,2,…,𝑛, 𝑗=1,2,…,𝑚. Then the two dimensions Haar function vector can be expressed as 𝐇𝑥₁,𝑥₂=ℎ₁𝑥₁ℎ₁𝑥₂⋮ℎ₁𝑥₁ℎ_𝑚𝑥₂ℎ₂𝑥₁ℎ₁𝑥₂⋮ℎ_𝑛𝑥₁ℎ_𝑚𝑥₂. (12) Any function 𝑓∈𝐿²([−𝜏₁,𝜏₁)×[−𝜏₂,𝜏₂)) can be written as 𝑓𝑥₁,𝑥₂=𝐶^𝑇𝐇𝑥₁,𝑥₂, (13) where 𝐶^𝑇=[𝑐₁₁⋯𝑐_1𝑛𝑐₂₁⋯𝑐_2𝑛⋯𝑐_𝑚1⋯𝑐_𝑚𝑛]. Subsequently, we assume that 𝑛=𝑚 and 𝜏₁=𝜏₂=𝜏, so that the operation matrix will be a square matrix. Let 𝐶=vec(𝐶̌) where 𝐶̌ is a 𝑚×𝑚 matrix. By using the Haar wavelet matrix in (6), the coefficient 𝐶^𝑇 in (13) can be obtained from 𝐶̌ as follows: 𝐶̌=𝐻_𝑚^−1𝑇⋅𝑓⋅𝐻_𝑚⁻¹ (14) and 𝑓_𝑖,𝑗=[𝑓(𝑥_𝑖,𝑥_𝑗)], 𝑖,𝑗=1,2,…,𝑚.

The integration of two dimensions Haar function vectors in [−𝜏,𝜏)×[−𝜏,𝜏) is ^𝑥𝑖₀𝐇𝑥₁,𝑥₂𝑑𝑥_𝑖=𝑄_𝑖−𝜏𝐸_𝑖𝐇𝑥₁,𝑥₂, (15) where 𝑄_𝑖 and 𝐸_𝑖 for 𝑖=1,2 are the 𝑚²×𝑚² operational matrices given as follows: 𝑄₁=𝑃_𝑚⊗𝐼_𝑚,𝑄₂=𝐼_𝑚⊗𝑃_𝑚,𝐸₁=𝐴_𝑚⊗𝐼_𝑚,𝐸₂=𝐼_𝑚⊗𝐴_𝑚, (16) where ⊗ denotes the Kronecker product [17], 𝐼_𝑚 denotes 𝑚×𝑚 identity matrix, and 𝐴_{𝑚𝑖,𝑗}=1,𝑖=1,2,𝑗=1,0,otherwise. (17) As in (10), we also required the product of 𝐇(𝑥₁,𝑥₂) and 𝐇^𝑇(𝑥₁,𝑥₂). Let 𝐇𝑥₁,𝑥₂𝐇^𝑇𝑥₁,𝑥₂𝐶=𝑁_𝐶𝐇𝑥₁,𝑥₂. (18) The algorithm to obtain 𝑁_𝐶 is as follows.

Step 1. Let 𝐶̌ be a matrix of 𝐶, or equivalently 𝐶=vec(𝐶̌).

Step 2. Compute 𝑀_̌𝐶𝑖, 𝑖=1,2,…,𝑚 according to (11) using the column 𝐶̌_𝑖 as the coefficient vector.

Step 3. For 𝑖=1,2,…,𝑚, compute vec(𝑀_̌𝐶𝑖).

Step 4. Form a big matrix by concatenating all vectors from Step 3; that is, 𝑆=[vec(𝑀_̌𝐶1)vec(𝑀_̌𝐶2)⋯vec(𝑀_̌𝐶𝑚)].

Step 5. For each row 𝑘 of matrix 𝑆, compute 𝑁_𝑖,𝑗 according to (11) using the row 𝑆_𝑘 as the coefficient vector.

Step 6. Form the matrix 𝑁_̆𝐶 as follows: 𝑁_̌𝐶=𝑁₁₁𝑁₁₂…𝑁_1𝑚𝑁₂₁𝑁₂₂…𝑁_2𝑚𝑁⋮⋮…⋮_𝑚1𝑁_𝑚2…𝑁_𝑚𝑚. (19)

Step 7. End.

4. Problem Statement

The system to be controlled is given by the nonlinear differential equation of the form 𝑥̇(𝑡)=𝑓(𝑥)+𝑔(𝑥)𝑢(𝑥),𝑥(0)=𝑥₀, (20) where 𝑥(𝑡)∈Ω⊂ℝ^𝑛 is the state vector, 𝑢∶Ω→ℝ^𝑚 is the control, 𝑓∶Ω→ℝ^𝑛 and 𝑔∶Ω→ℝ^𝑛×𝑚 are continuously differentiable with respect to all its arguments, 𝑥₀ is the initial condition vector, and Ω is domain of attraction.

The problem is to find the optimal control 𝑢^∗(𝑥) that minimizes the following performance index: 𝐽𝑥₀,𝑢=^∞₀𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢𝑑𝑡, (21) where 𝑄∈ℝ^𝑛×𝑛 is a positive semidefinite matrix and 𝑅∈ℝ^𝑚×𝑚 is a positive definite matrix. Given an arbitrary control 𝑢, the performance of the control at 𝑥∈Ω⊂ℝ^𝑛 is given by a Lyapunov function for the system [8] 𝑉(𝑥,𝑢)=^∞₀𝑙(𝑥(𝑡))+‖𝑢(𝑥(𝑡))‖²_𝑅𝑑𝑡, (22) where, ‖𝑢‖²_𝑅=𝑢^𝑇𝑅𝑢 and 𝑙(𝑥)=𝑥^𝑇𝑄𝑥. The optimal controller in feedback form is presented as follows [8]: 𝑢^∗(𝑥)1=−2𝑅⁻¹𝑔^𝑇(𝑥)𝜕𝑉^∗(𝑥),𝜕𝑥 (23) where 𝑉^∗(𝑥) is the solution to the following Hamilton-Jacobi-Bellman (HJB) equation 𝜕𝑉^∗𝑇(𝑥)𝑓(𝑥)(𝑥)−1𝜕𝑥+𝑙4𝜕𝑉^∗𝑇(𝑥)𝑔(𝑥)𝑅𝜕𝑥⁻¹×𝑔(𝑥)^𝑇𝜕𝑉^∗(𝑥)𝜕𝑥=0 (24) with boundary condition 𝑉^∗(0)=0; that is 𝑉(𝑥^∗,𝑢^∗)≤𝑉(𝑥,𝑢) for all 𝑢, and 𝑥^∗(𝑡) is the solution of 𝑥̇=𝑓(𝑥)+𝑔(𝑥)𝑢^∗(𝑡). Basically, it is not so easy to solve the nonlinear partial differential equation in (24) for the purpose of obtaining 𝑉^∗(𝑥) and consequently 𝑢^∗(𝑥) from (23); rather the following two linear equations have been iterated by the algorithm proposed by [8] 𝜕𝑉^(𝑖)𝑇(𝑥)𝑓(𝑥)(𝑥)𝑢𝜕𝑥+𝑔^(𝑖)(𝑥)+𝑙(𝑥)+𝑢^(𝑖)(𝑥)²_𝑅=0 (25) with initial condition 𝑉^(𝑖)(0)=0 and 𝑢^(𝑖+1)(𝑥)1=−2𝑅⁻¹𝑔^𝑇(𝑥)𝜕𝑉^(𝑖)(𝑥).𝜕𝑥 (26) Equation (25) is called the Generalized Hamilton-Jacobi-Bellman (GHJB) equation in [8]. In case of moderate presumptions, it has been established in [8] that the iteration between the GHJB (25) and the control (26) coincide with original HJB equation solution (24). If we can find a stabilizing control 𝑢⁽⁰⁾(𝑥) to start off, it is possible to iteratively enhance the performance of this controller using (25), (26), and finally the optimal controller can be optimally approximated. Moreover, at each iteration step the controller 𝑢^(𝑖) is a stable control.

5. The Successive Haar Wavelet Collocation Method

The following section describes the successive Haar wavelet collocation method (SHWCM) used for obtaining the two dimensional numerical solution to the HJB equation. In every step of this algorithm, an approximate solution to the GHJB equation (25) has been identified, namely, 𝜕𝑉^(𝑖)/𝜕𝑥, 𝑉^(𝑖), and 𝑢^(𝑖); all can be approximately expressed in term of Haar wavelets. As 𝑖→∞, 𝑉^(𝑖) and 𝑢^(𝑖) will approach the optimal solution 𝑉^∗ and 𝑢^∗, respectively.

Let us consider the following two-dimensional optimal feedback control problem min𝑉𝑥₀,𝑢=^∞₀𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢𝑑𝑡 (27) subject to the dynamics 𝑥̇=𝑓(𝑥)+𝑔(𝑥)𝑢(𝑥),𝑥(0)=𝑥₀, (28) where 𝑥=^𝑥1_𝑥2, 𝑓(𝑥)=^{𝑓1(𝑥1,𝑥2)}_{𝑓2(𝑥1,𝑥2)}, 𝑔(𝑥)=^{𝑔1(𝑥1,𝑥2)}_{𝑔2(𝑥1,𝑥2)}, and 𝑢∶Ω→ℝ.

Without loss of generality, the domain of attraction has been selected as Ω=[−𝜏,𝜏]×[−𝜏,𝜏] for the sake of convenience. The following equations express the pair of GHJB equation and the control law: 𝜕𝑉^(𝑖)𝑇(𝑥)𝑓(𝑥)(𝑥)𝑢𝜕𝑥+𝑔^(𝑖)(𝑥)+𝑥^𝑇𝑄𝑥+𝑢^(𝑖)𝑇𝑅𝑢^(𝑖)=0, (29) with initial condition 𝑉^(𝑖)(0)=0 and 𝑢^(𝑖+1)(𝑥)1=−2𝑅⁻¹𝑔^𝑇(𝑥)𝜕𝑉^(𝑖)(𝑥).𝜕𝑥 (30) For (28), if initially 𝑢⁽⁰⁾ is a stabilizing control, then from (29) the solution to GHJB equation affiliated with 𝑢⁽⁰⁾ becomes a Lyapunov function for the system and equals to the cost associated with 𝑢⁽⁰⁾ as follows: 𝜕𝑉^(0)𝑇(𝑥)𝑓(𝑥)(𝑥)𝑢𝜕𝑥+𝑔⁽⁰⁾(𝑥)+𝑥^𝑇𝑄𝑥+𝑢^(0)𝑇𝑅𝑢⁽⁰⁾=0. (31) According to (13), function approximation for 𝑓₁(𝑥)+𝑔₁(𝑥)𝑢⁽⁰⁾(𝑥), 𝑓₂(𝑥)+𝑔₂(𝑥)𝑢⁰(𝑥) and 𝑥^𝑇𝑄𝑥+𝑢^(0)𝑇(𝑥)𝑅𝑢⁽⁰⁾(𝑥), can be written as 𝑓₁(𝑥)+𝑔₁(𝑥)𝑢⁽⁰⁾(𝑥)=𝜃^𝑇𝐇𝑥₁,𝑥₂,𝑓₂(𝑥)+𝑔₂(𝑥)𝑢⁽⁰⁾(𝑥)=𝜇^𝑇𝐇𝑥₁,𝑥₂,𝑥^𝑇𝑄𝑥+𝑢^(0)𝑇(𝑥)𝑅𝑢⁽⁰⁾(𝑥)=𝑘^𝑇𝐇𝑥₁,𝑥₂, (32) where the coefficient vectors, 𝜃^𝑇, 𝜇^𝑇, and 𝑘^𝑇, can be calculate from (14). Since it is not possible to differentiate Haar functions, and as (29) only involves first-order derivatives of 𝑉, we assume that second-order partial derivative of 𝑉 exists; that is, 𝜕²𝑉𝜕𝑥₁𝜕𝑥₂=𝜔^𝑇𝐇𝑥₁,𝑥₂ (33) for some coefficient vector 𝜔.

With the assumption 𝜕²𝑉𝜕𝑥₁𝜕𝑥₂=𝜕²𝑉𝜕𝑥₂𝜕𝑥₁, (34) the first-order partial derivative can be obtained by integrating (33), with respect to 𝑥₁ and 𝑥₂, respectively, 𝜕𝑉𝜕𝑥₁=𝜔^𝑇𝑄₂−𝜏𝐸₂𝐇𝑥₁,𝑥₂+𝛼^𝑇₁𝐇𝑥₁,𝑥₂𝜕𝑉𝜕𝑥₂=𝜔^𝑇𝑄₁−𝜏𝐸₁𝐇𝑥₁,𝑥₂+𝛼^𝑇₂𝐇𝑥₁,𝑥₂, (35) where 𝛼^𝑇₁=[𝛼₁₁,…,𝛼_1𝑚,0,…,0] and 𝛼^𝑇₂=[𝛼₂₁,0,…,𝛼₂₂,0,…,0,…,𝛼_2𝑚,0,…0].

It should be noted that 𝜔^𝑇 has 𝑚² unknown variables while 𝛼^𝑇₁ and 𝛼^𝑇₂ have only 𝑚 unknown variables each. Now substituting (32) and (35) into (29), we have 𝜔^𝑇𝑄₂−𝜏𝐸₂𝑁_𝜃+𝑄₁−𝜏𝐸₁𝑁_𝜇+𝛼^𝑇₁𝑁_𝜃+𝛼^𝑇₂𝑁_𝜇=−𝑘^𝑇. (36) Equation (36) is a system of underdetermined linear equations with 𝑚² equations and (𝑚²+2𝑚) unknown variables which can solve for the unknown vectors 𝜔^𝑇, 𝛼^𝑇₁, and 𝛼^𝑇₂ using Moore-Penrose pseudoinverse [18]. The underdetermined equation is expected because the Lyapunov function is not unique. The Moore-Penrose solution is the particular solution whose vector 2-norm is minimal.

By using the solution of GHJB equation (29), a feedback control law 𝑢⁽¹⁾ is constructed using (30), which improves the efficiency of 𝑢⁽⁰⁾. The solution of the Hamilton-Jacobi-Bellman equation is uniformly approximated by repeating the above process.

Knowing that 𝑉(𝑥)=^𝑥₀∇𝑉^𝑇𝑑𝑥 (37) depends only on the initial and final points, not on the path followed, we can calculate the Lyapunov function 𝑉(𝑥) by integrating parallel to the axes [19] as follows: 𝑉𝑥₁,𝑥₂=^𝑥1₀𝜕𝑉𝜕𝑥₁𝑥₁,0𝑑𝑥₁+^𝑥2₀𝜕𝑉𝜕𝑥₂𝑥₁,𝑥₂𝑑𝑥₂. (38) This gives 𝑉𝑥₁,𝑥₂=𝛽^𝑇𝑄₁−𝜏𝐸₁+𝜔^𝑇𝑄₁−𝜏𝐸₁𝑄₂−𝜏𝐸₂+𝛼^𝑇₂𝑄₂−𝜏𝐸₂𝐇𝑥₁,𝑥₂, (39) where 𝛽^𝑇=(𝜔^𝑇(𝑄₂−𝜏𝐸₂)+𝛼^𝑇₁)𝐇(𝑥₁,0).

6. Numerical Examples

To show the efficiency of the proposed method, we applied our method to a linear quadratic optimal control problem and two nonlinear quadratic optimal control problems.

Example 1. Consider the following linear quadratic regulator (LQR): 1𝐽=2^∞₀𝑥²₁(𝑡)+𝑢²(𝑡)𝑑𝑡 (40) subject to ̇01𝑥=0100𝑥+𝑢. (41)
To solve this problem we take the initial stabilizing control 𝑢⁽⁰⁾(𝑥)=−𝑥₁−𝑥₂. Tables 1 and 2 show sample iteration results for 𝑢^(𝑖) and 𝑉^(𝑖), respectively, when 𝑚=8, 𝑥₁=−1/8. The iteration is terminated when the difference between two successive controls is less than 𝜖=0.001. Subsequent, in order to display two-dimensional plots, we fix the value for 𝑥₁ at 𝑥₁[𝑚/2]=−𝜏/𝑚 and 𝑥₂∈[−1,1). Figure 1 shows that for the particular LQR problem, the usage of 𝑚=16 is enough to approximate the exact optimal feedback control 𝑢^∗(𝑥)=−𝑥₁−√2𝑥₂; however, to approximate the exact cost function we require higher value of 𝑚 as shown in Figure 2.

Table 1 

Iteration results 𝑢^(𝑖) for Example 1 when 𝑚=8 and 𝑥₁=−1/8.

Table 2 

Iteration results 𝑉^(𝑖) for Example 1 when 𝑚=8 and 𝑥₁=−1/8.

Figure 1 

Optimal feedback control for Example 1 via the SHWCM with 𝑚=8,16 and 𝑥1=−0.1250,−0.0625, respectively.

Figure 2 

Value cost function for Example 1 via the SHWCM with 𝑚=8,16,32 and 𝑥1=−0.1250,−0.0625,−0.0313, respectively.

Example 2. Consider the following nonlinear optimal control problem [15]: 𝐽=^∞₀𝑥²₂+𝑢²𝑑𝑡 (42) subject to ̇𝑥𝑥=²−𝑥₁𝜋2+tan⁻¹5𝑥₁−5𝑥²₁21+25𝑥²₁+4𝑥₂+03𝑢. (43) The optimum solution for this problem is 𝑢^∗(𝑥)=−3𝑥₂ and 𝑉^∗=𝑥²₁(𝜋/2+tan⁻¹(5𝑥₁))+𝑥²₂. To solve this nonlinear optimal control problem, we started with the initial stabilizing control 𝑢⁽⁰⁾(𝑥)=−1.8𝑥₂. Figure 3 shows approximate optimal feedback control law 𝑢^∗ for 𝑚=8,16, and 32. The graph for 𝑚=64 overlaps with the exact optimal feedback control, and Figure 4 shows that the approximate cost function converges to the exact cost function as we increase the resolution. Figure 5 compares the exact state trajectories with approximate trajectories.

Figure 3 

Optimal feedback control for Example 2 via the SHWCM with 𝑚=8,16,32 and 𝑥1=−0.1250,−0.0625,−0.0313, respectively.

Figure 4 

Value cost function for Example 2 via the SHWCM with 𝑚=8,16,32 and 𝑥1=−0.1250,−0.0625,−0.0313, respectively.

Figure 5 

State trajectory comparison for Example 2.

Example 3. Consider the following optimal control problem [8]: 𝐽=^∞₀𝑥²₁+𝑥²₂+𝑢²𝑑𝑡 (44) subject to 𝑥̇=−𝑥³₁−𝑥₂𝑥₁+𝑥₂+01𝑢. (45)
The initial stabilizing control 𝑢⁽⁰⁾(𝑥)=0.4142𝑥₁−1.3522𝑥₂ can be obtained using feedback linearization method as outlined in [20]. The optimal feedback control and cost function obtained using SHWCM for various resolution 𝑚 = 8, 16, and 32 are illustrated in Figures 6 and 7, respectively. We believe that, by increasing Haar wavelet resolution, the SHWCM will be capable of yielding more accurate results. Figure 8 shows simulation of the system trajectories.

Figure 6 

Optimal feedback control for Example 3 via the SHWCM with 𝑚=8,16,32 and 𝑥1=−0.1250,−0.0625,−0.0313, respectively.

Figure 7 

Value cost function for Example 3 via the SHWCM with 𝑚=8,16,32 and 𝑥1=−0.1250,−0.0625,−0.0313, respectively.

Figure 8 

Some state trajectories for Example 3.

7. Conclusion

In this paper we had proposed a new numerical method for solving the Hamilton-Jacobi-Bellman equation, which appears in the formulation of optimal control problems. Our approach uses a combination of successive Generalized Hamilton-Jacobi-Bellman equation and Haar wavelets operational matrix methods. The proposed approach is simple and stable and has been tested on linear and nonlinear optimal control problem in two-dimensional state space. Generally, by using our method, the approximate solutions for optimal feedback control require lower resolution, than the approximate solutions for the cost function. However, in both cases, it is clear that more accurate results can be obtained by increasing the resolution of Haar wavelet.

Acknowledgments

The authors are very grateful to the referees for their valuable comments and suggestions, which greatly improved the presentation of this paper. This research has been funded by University of Malaya, under Grant No. RG208-11AFR.