A novel unified variational finite difference (VFD) solution method for optimal control problems

Abstract

This article is aimed to propose a simple yet efficient unified numerical strategy for solving both linear and non-linear optimal control problems. To do so, the general form of quadratic performance index function and nonlinear state equations are considered first. Then, the idea of variational differential quadrature method is used to convert the integral/differential equations to the equivalent algebraic form. Since time is the only independent variable is this research, the finite difference method with an equally-spaced discretization scheme would be a more appropriate technique rather than the differential quadrature approach. So, the implemented numerical solution is called now as the variational finite difference method. The method of Lagrange multipliers is then utilized for minimization purpose and, as a result, the final set of nonlinear algebraic equations are obtained. Finally, the quadratic and triadic forms of non-linearity are considered and an explicit formulation is represented for the residual and Jacobian of the Newton-based iterative solution procedure. To demonstrate the accuracy and efficiency of the proposed approach in the quadratic optimal control area, several benchmark problems involving linear time-invariant, linear time-variant and nonlinear examples are successfully solved and the results are confirmed with those existed in the literature.

Introduction

Optimal control is an active research area in the control theory with several multi-disciplinary applications. It is used in various fields including bioengineering, economics, robotics, management, process control, finance, aerospace and communication networks.^1–3 The well-known theories of linear-quadratic and Gaussian optimal control are dated back to 1960s represented by Kalman. ⁴ The optimal control theory is mainly utilized to transfer the system to the final desired situation by minimizing an objective function and satisfying certain conditions.

Many researchers have addressed various optimal control problems, and several solution approaches have been proposed in their works. Generally speaking, the reason of enormous attempts is the difficulty of developing an efficient solution technique to find the optimum value of a cost function (in integral form) subjected to the nonlinear state equations (in differential form) and certain boundary conditions (in algebraic form). Tsai et al. ⁵ developed an approach to design an optimal digital regulator for continuous-time two-dimensional systems described by linear partial differential equations (PDEs). A new state vector was introduced to directly convert the original continuous-time two-dimensional quadratic cost function into a decoupled discretized form. A numerical method was developed by Jaddu ⁶ to solve nonlinear optimal control problems with terminal state constraints, control inequality constraints and simple bounds on the state variables. Huang and Li ⁷ used the conjugate gradient method with adjoint equation to solve a general three-dimensional optimal boundary control problem with the aim of determining the optimal boundary control functions. Gerdts et al. ⁸ studied optimal control problems for the wave equation. To accomplish this aim, they employed a full discretization method based upon finite differences through which the PDE constrained optimal control problem is transformed into a nonlinear programming problem. On the other hand, Tricaud and Chen ⁹ proposed a numerical approximate method to solve fractional-order optimal control problems in general form. Casas and Tröltzsch ¹⁰ analyzed a distributed optimal control problem governed by a quasilinear elliptic equation of non-monotone type using a numerical approach. Yousefi et al. ¹¹ suggested an approximate solution for the fractional optimal control problems by using the Legendre multiwavelet in conjunction with the collocation method. Also, state-constrained optimal control problems for PDEs were numerically analyzed by Neitzel and Tröltzsch. ¹² Based on interpolating scaling functions, Foroozandeh and Shamsi ¹³ presented a numerical approach to solve nonlinear optimal control problems including state and control inequality constraints. Kammann et al. ¹⁴ presented a posteriori error estimation for semilinear parabolic optimal control problems with application to model reduction by proper orthogonal decomposition (POD). Bhrawy et al. ¹⁵ utilized the shifted Legendre orthogonal polynomials to numerically solve the multi-dimensional fractional optimal control problem with a quadratic performance index.

More recently, by employing the Ritz method, Mamehrashi and Yousefi ¹⁶ obtained numerical solution for a class of two-dimensional quadratic optimal control problems. Keshavarz et al. ¹⁷ derived a general formulation for the Bernoulli operational matrix of fractional order integration and multiplication. Mamehrashi and Yousefi ¹⁸ also presented a numerical solution approach based upon shifted Legendre polynomials for solving a nonlinear two-dimensional optimal control problem. Furthermore, based on the combination of epsilon penalty and variational methods, Lotfi ¹⁹ developed an approximate method to address a class of multidimensional fractional optimal control problems. Bourdin and Trélat ²⁰ perused the linear-quadratic optimal sampled-data control problems. The system under consideration involves the continuous state variables and piecewise constant control in time. Also, Chen et al. ²¹ presented an investigation on the system identification and control of convex problems by respectively employing the neural networks and optimal controller design. Hajiloo et al. ²² presented a multi-objective optimization approach for Pareto optimum design of robust integer-order and fractional-order PID controllers for both integer-order and fractional-order plants with parametric uncertainties. Yıldırım Aksoy et al. ²³ utilized the discretization scheme of the finite difference method for the analysis of nonlinear optimal control problems subjected to Schrödinger equations with complex coefficients. They provided an analytical study to prove the stability and convergence of the finite difference scheme. Recently, Mohammadi et al. ²⁴ investigated the linear and nonlinear optimal control problems through the archived-based genetic programming (AGP). In addition of addressing the quadratic case, they also studied a new class of cost functional named as the absolute performance index (API). On the other hand, readers are invited to read the comprehensive review paper of Diehl et al. ²⁵ to find detailed data on the available algorithms for numerically solving the optimal control problems based on the nonlinear model predictive control (MPC) and moving horizon estimation (MHE). Also, Zhang et al. ²⁶ presented a review to study the near-optimal control problems in the area of nonlinear dynamical systems through different nonlinear programming approaches including the methods via control parameterization, model predictive control and the method via Taylor expansion. The computational demand of real-time nonlinear MPC and MHE approaches in mechatronic applications is examined by Vukov et al. ²⁷

According to the foregoing literature review on the optimal control solution strategies, one may find that solving the nonlinear problems requires an extensive computational cost based on the available approaches. For instance, refer to Jaddu, ⁶ in which the nonlinear optimal control problem is converted to a sequence of constrained linear-quadratic problems; Mamehrashi and Yousefi, ¹⁶ in which the Ritz method based on the certain approximations of shifted Legendre polynomials was utilized; Chen et al. ²¹ and Mohammadi et al. ²⁴ in which time-consuming solution processes of deep neural networks and genetic algorithm were employed, respectively. So, developing an accurate and efficient solution methodology could be promising and majorly contributes to this research area. Accordingly, a novel numerical method is proposed in present paper to solve both linear and nonlinear optimal control problems.

The method developed in this work is based on the variational differential quadrature (VDQ) approach which has been originally established in the area of computational mechanics. ²⁸ In the short time of its emergence, VDQ has been employed in different engineering problems which are classically known as the realm of powerful numerical techniques such as the finite element analysis. Readers can refer to the references^29–32 to find the examples. Using the differential operators of the generalized differential quadrature (GDQ) method, VDQ represents an appropriate discretization scheme for the problems with space-domain parameters. In the case of present study, with time-dependent state variables, it would be an efficient choice to consider the differential operators of the finite difference method. And, the integral operators are correspondingly derived in the similar way of VDQ, that is, through the Taylor series expansion. Hence, the present approach is called as the variational finite difference (VFD) method.

Here, the basic idea is converting the performance index, state equations and boundary conditions to a functional via the Lagrange method which must be minimized. In this regard, section “Problem statement” is dedicated to defining the corresponding minimization problem. Presented in section “Solution strategy,” accurate integral and differential operational matrices of the VFD method are introduced in order to discretize the obtained weak form equation. Since the majority of nonlinear case studies in optimal control theory can be modeled via second- and third-order polynomial terms, an explicit formulation is developed for quadratic and triadic nonlinearities. Finally, a system of nonlinear algebraic equations is obtained which can be solved by means of the iterative methods. Several illustrative test examples are presented in section “Selected case studies” to demonstrate the validity and general applicability of the proposed method for the linear time-invariant, linear time-variant and nonlinear optimal control problems. Also, it is shown that the computational cost of analyzing different problems is significantly reduced by employing the proposed VFD solution strategy in comparison with the AGP method in recent relevant paper. ²⁴

Problem statement

The optimal control theory is aimed to minimize the performance index function constrained by a specific state path for dynamic system. In order to take positive and negative deviations from the optimal control signal into account, one may consider the absolute or quadratic forms.^24,33 However, it would be mathematically straightforward to deal with the quadratic formulation. So, it is aimed to minimize the following quadratic form of the cost function

$$J=\frac{1}{2}{\int }_{t_0}^{t_f}({\mathbf{x}}^{\mathrm{T}}\left(t\right)\mathbf{Qx}\left(t\right)+{\mathbf{u}}^{\mathrm{T}}\left(t\right)\mathbf{Ru}\left(t\right))\mathrm{dt},$$

(1)

subjected to the state equation

$$\stackrel{·}{\mathbf{x}}\left(t\right)=\mathbf{Ax}\left(t\right)+\mathbf{Bu}\left(t\right)+{\mathbf{R}}_{\mathrm{nl}}(\mathbf{x},\mathbf{u}),$$

(2)

and initial/terminate conditions

$$\mathbf{x}\left(t_0\right)={\mathbf{x}}_{t_0},\mathbf{x}\left(t_f\right)={\mathbf{x}}_{t_f},$$

(3)

where the state variables are

$$\mathbf{x}={\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\\ \multicolumn{1}{c}{⋮}\end{array}\right]}_{m_1\times 1},\mathbf{u}={\left[\begin{array}{c}{u}_1\hfill \\ \multicolumn{1}{c}{u_2}\\ \multicolumn{1}{c}{⋮}\end{array}\right]}_{m_2\times 1},$$

(4)

and ${\left[{\mathbf{R}}_{\mathrm{nl}}\right]}_{m_1\times 1}$ denotes the nonlinear vector in state equation. Also, ${\left[\mathbf{Q}\right]}_{m_1\times m_1}$ , ${\left[\mathbf{R}\right]}_{m_2\times m_2}$ signal the weighted matrices and ${\left[\mathbf{A}\right]}_{m_1\times m_1}$ , ${\left[\mathbf{B}\right]}_{m_1\times m_2}$ are the system matrices, in which $m_1$ and $m_2$ respectively stand for the number of state variables in $\mathbf{x}$ and $\mathbf{u}$ vectors.

Solution strategy

In order to solve the optimal control problem of finding $\mathbf{x}$ and $\mathbf{u}$ , equations (1) and (2), which respectively are in integral and differential forms, must be converted to the equivalent algebraic form. In this regard, a numerical approach consists of appropriate integral and differential operators is introduced in the following.

First, the time domain is discretized by n equally-space points

$$\mathbb{t}={\left[\begin{array}{ccccc}{t}_0\hfill & t_0+\mathrm{\Delta }t\hfill & t_0+2\mathrm{\Delta }t\hfill & \dots \hfill & t_f\hfill \end{array}\right]}^{\mathrm{T}},t_f=t_0+(n-1)\mathrm{\Delta }t.$$

(5)

So, the state and control variables will have the following discretized forms

$$\mathbb{x}={\left[\begin{array}{c}{\mathbb{x}}_1\hfill \\ \multicolumn{1}{c}{{\mathbb{x}}_2}\\ \multicolumn{1}{c}{⋮}\end{array}\right]}_{m_{1}n\times 1},\mathbb{u}={\left[\begin{array}{c}{\mathbb{u}}_1\hfill \\ \multicolumn{1}{c}{{\mathbb{u}}_2}\\ \multicolumn{1}{c}{⋮}\end{array}\right]}_{m_{2}n\times 1},$$

(6)

$$\begin{array}{c}{\mathbb{x}}_1={\left[\begin{array}{ccccc}{x}_1\left(t_0\right)\hfill & x_1(t_0+\mathrm{\Delta }t)\hfill & x_1(t_0+2\mathrm{\Delta }t)\hfill & \dots \hfill & x_1\left(t_f\right)\hfill \end{array}\right]}^{\mathrm{T}},\\ {\mathbb{x}}_2={\left[\begin{array}{ccccc}{x}_2\left(t_0\right)\hfill & x_2(t_0+\mathrm{\Delta }t)\hfill & x_2(t_0+2\mathrm{\Delta }t)\hfill & \dots \hfill & x_2\left(t_f\right)\hfill \end{array}\right]}^{\mathrm{T}},\dots ,\end{array}$$

(7)

$$\begin{array}{c}{\mathbb{u}}_1={\left[\begin{array}{ccccc}{u}_1\left(t_0\right)\hfill & u_1(t_0+\mathrm{\Delta }t)\hfill & u_1(t_0+2\mathrm{\Delta }t)\hfill & \dots \hfill & u_1\left(t_f\right)\hfill \end{array}\right]}^{\mathrm{T}},\\ {\mathbb{u}}_2={\left[\begin{array}{ccccc}{u}_2\left(t_0\right)\hfill & u_2(t_0+\mathrm{\Delta }t)\hfill & u_2(t_0+2\mathrm{\Delta }t)\hfill & \dots \hfill & u_2\left(t_f\right)\hfill \end{array}\right]}^{\mathrm{T}},\dots ,\end{array}$$

(8)

and one has the total unknown parameters in the form of the augmented vector of both states and control efforts of order $(m_1+m_2)n$ as follows

$$\mathbb{q}=\left[\begin{array}{c}\mathbb{x}\hfill \\ \multicolumn{1}{c}{\mathbb{u}}\end{array}\right].$$

(9)

Now, the numerical differentiation and integration operators can be introduced. According to the time discretization, the finite difference method could be an appropriate candidate to evaluate derivatives in the state equation. If, we consider the discretized form of a typical function $f\left(t\right)$ as vector ${\left[\mathbb{f}\right]}_{n\times 1}$ , the r-th derivative of $f\left(t\right)$ in whole time domain is given by

$$\frac{d^r\mathbb{f}}{d{t}^r}={\mathbf{D}}^{\left(r\right)}\mathbb{f},$$

(10)

where, ${\left[{\mathbf{D}}^{\left(r\right)}\right]}_{n\times n}$ is the r-th order derivative operator matrix of the finite difference method (see Appendix A).

In order to calculate the integral value at the same points of time discretization, the integral operator must be proportional to the differential one. ²⁸ The idea of variational differential quadrature approach is used herein to develop the variational finite difference solution methodology and tackle the optimal control problems. Thus, the corresponding approach can be adopted here to obtain the integral operator vector ${\left[\mathbf{S}\right]}_{1\times n}$ (see Appendix B). It leads to

$${\int }_{t_0}^{t_f}f\left(t\right)\mathrm{dt}=\mathbf{S}\mathbb{f}.$$

(11)

Therefore, using the discretization of the quadratic form of the integrand, the algebraic form of equation (1) is accordingly obtained as

$$J=\frac{1}{2}\mathbf{S}({\mathbf{Q}}_{\mathrm{ij}}{\mathbb{x}}_i°{\mathbb{x}}_j+{\mathbf{R}}_{\mathrm{ij}}{\mathbb{u}}_i°{\mathbb{u}}_j),$$

(12)

where i and j stand for the vector indexes and ° shows the Hadamard product with the following definition

$${\left[\mathbf{X}\right]}_{p_1\times p_2}°{\left[\mathbf{Y}\right]}_{p_1\times p_2}={\left[\begin{array}{ccc}{X}_{11}{Y}_{11}\hfill & X_{12}{Y}_{12}\hfill & \dots \hfill \\ \multicolumn{1}{c}{X_{21}{Y}_{21}}& X_{22}{Y}_{22}\hfill & \dots \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill \end{array}\right]}_{p_1\times p_2}.$$

(13)

from the lemma ${\mathbf{V}}_1^{\mathrm{T}}({\mathbf{V}}_2°{\mathbf{V}}_3)={\mathbf{V}}_2^{\mathrm{T}}{\mathbf{V}}_1{\mathbf{V}}_3$ for three arbitrary column vectors and considering that $〈\mathbf{V}〉$ is the diagonal form of vector $\mathbf{V}$ , equation (12) is rewritten as

$$J=\frac{1}{2}({\mathbb{x}}_i^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{ij}}〈\mathbf{S}〉{\mathbb{x}}_j+{\mathbb{u}}_i^{\mathrm{T}}{\mathbf{R}}_{\mathrm{ij}}〈\mathbf{S}〉{\mathbb{u}}_j),$$

(14)

which can also be given in the following total matrix representation

$$J=\frac{1}{2}({\mathbb{x}}^{\mathrm{T}}(\mathbf{Q}⊛〈\mathbf{S}〉)\mathbb{x}+{\mathbb{u}}^{\mathrm{T}}(\mathbf{R}⊛〈\mathbf{S}〉)\mathbb{u}),$$

(15)

where $\mathrm{⊛}$ is the symbol of Kronecker product defined as

$${\left[\mathbf{X}\right]}_{p_1\times p_2}⊛{\left[\mathbf{Y}\right]}_{q_1\times q_2}={\left[\begin{array}{ccc}{X}_{11}\mathbf{Y}\hfill & X_{12}\mathbf{Y}\hfill & \dots \hfill \\ \multicolumn{1}{c}{X_{21}\mathbf{Y}}& X_{22}\mathbf{Y}\hfill & \dots \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill \end{array}\right]}_{p_1{q}_1\times p_2{q}_2}$$

(16)

Also, discretization of equation (2) leads to

$$\left({\mathbf{I}}_{m_1}⊛{\mathbf{D}}^{\left(1\right)}\right)\mathbb{x}=(\mathbf{A}⊛{\mathbf{I}}_n)\mathbb{x}+(\mathbf{B}⊛{\mathbf{I}}_n)\mathbb{u}+{\mathbb{R}}_{\mathrm{nl}},$$

(17)

where ${\mathbf{I}}_p$ denotes the p-by-p identity matrix and ${\left[{\mathbb{R}}_{\mathrm{nl}}\right]}_{m_{1}n\times 1}$ is the discretized form of the nonlinear state vector.

Therefore, the cost function, state equation, and the corresponding boundary conditions in equation (3) are, respectively, represented in the following forms

$$J=\frac{1}{2}{\mathbb{q}}^{\mathrm{T}}\mathbb{Sq},$$

(18)

$$\mathbb{C}\mathbb{q}-{\mathbb{R}}_{\mathrm{nl}}=\mathbf{0},$$

(19)

$$\mathbb{E}\mathbb{q}-{\mathbb{x}}_0=\mathbf{0},$$

(20)

where

$$\mathbb{S}=\left[\begin{array}{cc}\mathbf{Q}⊛〈\mathbf{S}〉\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{R}⊛〈\mathbf{S}〉\hfill \end{array}\right],$$

(21)

$$ℂ=\left[\begin{array}{cc}\mathbb{D}-\mathbb{A}& -\mathbb{B}\end{array}\right],\mathbb{D}={\mathbf{I}}_{m_1}⊛{\mathbf{D}}^{\left(1\right)},\mathbb{A}=\mathbf{A}⊛{\mathbf{I}}_n,\mathbb{B}=\mathbf{B}⊛{\mathbf{I}}_n.$$

(22)

and $\mathbb{E}$ shows a matrix containing zero and one for the initial and terminate values denoted by the vector ${\mathbb{x}}_0$ . It is worth mentioning that the discretized formulation of the vector ${\mathbb{R}}_{\mathrm{nl}}$ depends on the nature of nonlinearity in the problem. It should be also noted that $\mathbb{A}$ and $\mathbb{B}$ are written for the time-invariant problems. However, such idea can be readily extended to the case of time-variant state equations, that is, $\mathbf{A}$ and $\mathbf{B}$ are the function of time, in which each component should be discretized in its diagonal area (see section “Selected case studies”).

As a result of the developed VFD formulation and instead of the primary integro-differential-algebraic equations in equations (1)–(3), the design of the optimal control problem is now observed as that a system of $(m_1+m_2)n$ algebraic equations in terms of the augmented vector in equation (9) should be minimized with the $m_{1}n$ constraints given in equation (19) and the ones in equation (20). At this stage, by using the method of Lagrange multipliers, the problem is restated as

$$\begin{array}{c}\mathrm{minimize}J(\mathrm{equation}(18):J=\frac{1}{2}{\mathbb{q}}^{\mathrm{T}}\mathbb{S}\mathbb{q}),\\ \begin{array}{cc}\multicolumn{1}{c}{}& \mathrm{subject}\mathrm{to}\mathrm{equation}(19)(\mathbb{C}\mathbb{q}-{\mathbb{R}}_{\mathrm{nl}}=\mathbf{0})\mathrm{and}\hfill \\ \multicolumn{1}{c}{}& \mathrm{equation}(2)0(\mathbb{E}\mathbb{q}-{\mathbb{x}}_0=\mathbf{0}),\hfill \end{array}\end{array}$$

(23)

and the Lagrange function is given by

$$\mathcal{L}(\mathbb{q},{\mathbb{L}}_1,{\mathbb{L}}_2)=\frac{1}{2}{\mathbb{q}}^{\mathrm{T}}\mathbb{S}\mathbb{q}+{\mathbb{L}}_1(\mathbb{C}\mathbb{q}-{\mathbb{R}}_{\mathrm{nl}})+{\mathbb{L}}_2(\mathbb{E}\mathbb{q}-{\mathbb{x}}_0).$$

(24)

where ${\mathbb{L}}_1$ is the first vector of Lagrange multipliers with size of $1\times m_{1}n$ and ${\mathbb{L}}_2$ is the second one.

Setting the gradient of the Lagrange function equal to zero, the governing equations would be determined (i.e. $\nabla \mathcal{L}(\mathbb{q},{\mathbb{L}}_1,{\mathbb{L}}_2)=\mathbf{0}$ ). Equivalently, the first variation of Lagrange function is calculated as

$$\begin{array}{cc}\multicolumn{1}{c}{\delta \mathcal{L}}& =\left(\frac{\partial \mathcal{L}}{\partial \mathbb{q}}\right)·\delta \mathbb{q}+\left(\frac{\partial \mathcal{L}}{\partial {\mathbb{L}}_1}\right)·\delta {\mathbb{L}}_1^{\mathrm{T}}+\left(\frac{\partial \mathcal{L}}{\partial {\mathbb{L}}_2}\right)·\delta {\mathbb{L}}_2^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& =\delta {\mathbb{q}}^{\mathrm{T}}(\mathbb{S}\mathbb{q}+{\left({\mathbb{L}}_1(\mathbb{C}-\partial {\mathbb{R}}_{\mathrm{nl}}/\partial \mathbb{q})\right)}^{\mathrm{T}}+{\left({\mathbb{L}}_2\mathbb{E}\right)}^{\mathrm{T}})\hfill \\ \multicolumn{1}{c}{}& +\delta {\mathbb{L}}_1(\mathbb{C}\mathbb{q}-{\mathbb{R}}_{\mathrm{nl}})+\delta {\mathbb{L}}_2(\mathbb{E}\mathbb{q}-{\mathbb{x}}_0).\hfill \end{array}$$

(25)

The necessary condition of vanishing the first variation of the Lagrange function, that is, $\delta \mathcal{L}=0$ , is achieved by letting the coefficients of independent vectors ${\delta }_{\mathbb{L}}$ , $\delta {\mathbb{L}}_1$ , and $\delta {\mathbb{L}}_2$ set to zeros. It, therefore, results that

$$\mathbb{S}\mathbb{q}+{\left({\mathbb{L}}_1(\mathbb{C}-\partial {\mathbb{R}}_{\mathrm{nl}}/\partial \mathbb{q})\right)}^{\mathrm{T}}+{\left({\mathbb{L}}_2\mathbb{E}\right)}^{\mathrm{T}}=\mathbf{0},$$

(26)

$$\mathbb{C}\mathbb{q}-{\mathbb{R}}_{\mathrm{nl}}=\mathbf{0},$$

(27)

$$\mathbb{E}\mathbb{q}-{\mathbb{x}}_0=\mathbf{0}.$$

(28)

Such obtained nonlinear set of equations can now be solved in a straightforward manner by utilizing the iterative solution procedure of Newton-based methods. Here, an explicit form of nonlinearity, if any, is considered to complete the formulation. In this way, it is assumed that the nonlinear vector is consisted of second- and third-order terms so that the quadratic and triadic nonlinearities are, respectively, given by

$$\left({\mathbb{Q}}_1\mathbb{q}\right)°\left({\mathbb{Q}}_2\mathbb{q}\right)=\frac{1}{2}(〈{\mathbb{Q}}_1\mathbb{q}〉{\mathbb{Q}}_2+〈{\mathbb{Q}}_2\mathbb{q}〉{\mathbb{Q}}_1)\mathbb{q}=\frac{1}{2}{\mathbb{R}}_1\mathbb{q},$$

(29)

$$\left({\mathbb{Q}}_3\mathbb{q}\right)\circ \left({\mathbb{Q}}_4\mathbb{q}\right)\circ \left({\mathbb{Q}}_5\mathbb{q}\right)=\frac{1}{3}\left({\mathbb{Q}}_3\mathbb{q}{\mathbb{Q}}_4\mathbb{q}{\mathbb{Q}}_5+{\mathbb{Q}}_3\mathbb{q}{\mathbb{Q}}_5\mathbb{q}{\mathbb{Q}}_4+{\mathbb{Q}}_4\mathbb{q}{\mathbb{Q}}_5\mathbb{q}{\mathbb{Q}}_3\right)\mathbb{q}=\frac{1}{3}{ℝ}_2\mathbb{q},$$

(30)

where ${\mathbb{R}}_1$ , ${\mathbb{R}}_2$ , ${\mathbb{Q}}_1$ , …, ${\mathbb{Q}}_5$ are matrices of order $m_{1}n\times (m_1+m_2)n$ and one can generally write

$${\mathbb{R}}_{\mathrm{nl}}=(\frac{1}{2}{\mathbb{R}}_1+\frac{1}{3}{\mathbb{R}}_2)\mathbb{q},$$

(31)

which results that (see Appendix C)

$$\frac{\partial {\mathbb{R}}_{\mathrm{nl}}}{\partial \mathbb{q}}={\mathbb{R}}_1+{\mathbb{R}}_2.$$

(32)

Thus, from equations (26)–(28) to equations (31)–(32), the final equation reads

$${\mathbb{K}}_s\mathbb{d}=\mathbb{F},$$

(33)

where

$$\begin{array}{c}\multicolumn{1}{c}{\mathbb{d}=\left[\begin{array}{c}\mathbb{q}\hfill \\ \multicolumn{1}{c}{{\mathbb{L}}_1^{\mathrm{T}}}\\ \multicolumn{1}{c}{{\mathbb{L}}_2^{\mathrm{T}}}\end{array}\right],}\end{array}$$

(34)

$$\mathbb{F}=\left[\begin{array}{c}\mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}\\ \multicolumn{1}{c}{{\mathbb{x}}_0}\end{array}\right],$$

(35)

and the stiffness matrix ${\mathbb{K}}_s$ is given by

$${\mathbb{K}}_s={\mathbb{K}}_l+\frac{1}{2}{\mathbb{K}}_{n1}+\frac{1}{3}{\mathbb{K}}_{n2},$$

(36)

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mathbb{K}}_l=\left[\begin{array}{ccc}\mathbb{S}\hfill & {\mathbb{C}}^{\mathrm{T}}\hfill & {\mathbb{E}}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{\mathbb{C}}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbb{E}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],{\mathbb{K}}_{n1}=-\left[\begin{array}{ccc}\mathbf{0}\hfill & 2{\mathbb{R}}_1^{\mathrm{T}}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbb{R}}_1}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],\hfill \\ \multicolumn{1}{c}{}& {\mathbb{K}}_{n2}=-\left[\begin{array}{ccc}\mathbf{0}\hfill & 3{\mathbb{R}}_2^{\mathrm{T}}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbb{R}}_2}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right].\hfill \end{array}$$

(37)

${\mathbb{K}}_l$ is the linear matrix, ${\mathbb{K}}_{n1}$ and ${\mathbb{K}}_{n2}$ respectively denote the second- and third-order nonlinear parts.

In addition of the above-mentioned formulation, the Jacobi of nonlinear algebraic equations is required in Newton approach. Deriving the incremental form of equation (25), one has $\mathrm{\Delta }\delta \mathcal{L}=\delta {\mathbb{d}}^{\mathrm{T}}{\mathbb{K}}_t\mathrm{\Delta }\mathbb{d}$ in which ${\mathbb{K}}_t$ is the tangent stiffness matrix. Since $\mathrm{\Delta }\delta \mathbb{q}=\mathrm{\Delta }\delta {\mathbb{L}}_1=\mathrm{\Delta }\delta {\mathbb{L}}_2=0$ , the increment of variation of the Lagrange function becomes

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{\Delta }\delta \mathcal{L}}& =\delta {\mathbb{q}}^{\mathrm{T}}(\mathbb{S}\mathrm{\Delta }\mathbb{q}+{\mathbb{C}}^{\mathrm{T}}\mathrm{\Delta }{\mathbb{L}}_1^{\mathrm{T}}-\mathrm{\Delta }\left\{{({\mathbb{R}}_1+{\mathbb{R}}_2)}^{\mathrm{T}}{\mathbb{L}}_1^{\mathrm{T}}\right\}+{\mathbb{E}}^{\mathrm{T}}\mathrm{\Delta }{\mathbb{L}}_2^{\mathrm{T}})\hfill \\ \multicolumn{1}{c}{}& +\delta {\mathbb{L}}_1(\mathbb{C}\mathrm{\Delta }\mathbb{q}-\mathrm{\Delta }\left\{(\frac{1}{2}{\mathbb{R}}_1+\frac{1}{3}{\mathbb{R}}_2)\mathbb{q}\right\})+\delta {\mathbb{L}}_2\mathbb{E}\mathrm{\Delta }\mathbb{q},\hfill \end{array}$$

(38)

considering

$$\mathrm{\Delta }{\mathbb{R}}_1=〈{\mathbb{Q}}_1\mathrm{\Delta }\mathbb{q}〉{\mathbb{Q}}_2+〈{\mathbb{Q}}_2\mathrm{\Delta }\mathbb{q}〉{\mathbb{Q}}_1,$$

(39)

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{\Delta }{\mathbb{R}}_2}& =〈{\mathbb{Q}}_3\mathrm{\Delta }\mathbb{q}〉〈{\mathbb{Q}}_4\mathbb{q}〉{\mathbb{Q}}_5+〈{\mathbb{Q}}_3\mathbb{q}〉〈{\mathbb{Q}}_4\mathrm{\Delta }\mathbb{q}〉{\mathbb{Q}}_5\hfill \\ \multicolumn{1}{c}{}& +〈{\mathbb{Q}}_3\mathrm{\Delta }\mathbb{q}〉〈{\mathbb{Q}}_5\mathbb{q}〉{\mathbb{Q}}_4+〈{\mathbb{Q}}_3\mathbb{q}〉〈{\mathbb{Q}}_5\mathrm{\Delta }\mathbb{q}〉{\mathbb{Q}}_4\hfill \\ \multicolumn{1}{c}{}& +〈{\mathbb{Q}}_4\mathrm{\Delta }\mathbb{q}〉〈{\mathbb{Q}}_5\mathbb{q}〉{\mathbb{Q}}_3+〈{\mathbb{Q}}_4\mathbb{q}〉〈{\mathbb{Q}}_5\mathrm{\Delta }\mathbb{q}〉{\mathbb{Q}}_3,\hfill \end{array}$$

(40)

and the lemma $〈\mathbf{X}\mathbf{〉}\mathbf{Y}=\mathbf{〈}\mathbf{Y}\mathbf{〉}\mathbf{X}$ for two arbitrary vectors, one has

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{\Delta }\left\{{({\mathbb{R}}_1+{\mathbb{R}}_2)}^{\mathrm{T}}{\mathbb{L}}_1^{\mathrm{T}}\right\}}& ={(\mathrm{\Delta }{\mathbb{R}}_1+\mathrm{\Delta }{\mathbb{R}}_2)}^{\mathrm{T}}{\mathbb{L}}_1^{\mathrm{T}}+{({\mathbb{R}}_1+{\mathbb{R}}_2)}^{\mathrm{T}}\mathrm{\Delta }{\mathbb{L}}_1^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& =({\tilde{\mathbb{R}}}_1+{\tilde{\mathbb{R}}}_2)\mathrm{\Delta }\mathbb{q}+{({\mathbb{R}}_1+{\mathbb{R}}_2)}^{\mathrm{T}}\mathrm{\Delta }{\mathbb{L}}_1^{\mathrm{T}},\hfill \end{array}$$

(41)

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{\Delta }\left\{(\frac{1}{2}{\mathbb{R}}_1+\frac{1}{3}{\mathbb{R}}_2)\mathbb{q}\right\}}& =(\frac{1}{2}\mathrm{\Delta }{\mathbb{R}}_1+\frac{1}{3}\mathrm{\Delta }{\mathbb{R}}_2)\mathbb{q}\hfill \\ \multicolumn{1}{c}{}& +(\frac{1}{2}{\mathbb{R}}_1+\frac{1}{3}{\mathbb{R}}_2)\mathrm{\Delta }\mathbb{q}=({\mathbb{R}}_1+{\mathbb{R}}_2)\mathrm{\Delta }\mathbb{q},\hfill \end{array}$$

(42)

where

$${\tilde{\mathbb{R}}}_1={\mathbb{Q}}_1^{\mathrm{T}}〈\mathbb{L}{〉}_1^{\mathrm{T}}{\mathbb{Q}}_2+{\mathbb{Q}}_2^{\mathrm{T}}〈\mathbb{L}{〉}_1^{\mathrm{T}}{\mathbb{Q}}_1,$$

(43)

$$\begin{array}{cc}\multicolumn{1}{c}{{\tilde{\mathbb{R}}}_2}& ={\mathbb{Q}}_3^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_4\mathbb{q}〉{\mathbb{Q}}_5+{\mathbb{Q}}_3^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_5\mathbb{q}〉{\mathbb{Q}}_4\hfill \\ \multicolumn{1}{c}{}& +{\mathbb{Q}}_4^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_3\mathbb{q}〉{\mathbb{Q}}_5+{\mathbb{Q}}_4^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_5\mathbb{q}〉{\mathbb{Q}}_3\hfill \\ \multicolumn{1}{c}{}& +{\mathbb{Q}}_5^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_3\mathbb{q}〉{\mathbb{Q}}_4+{\mathbb{Q}}_5^{\mathrm{T}}〈{\mathbb{L}}_1^{\mathrm{T}}〉〈{\mathbb{Q}}_4\mathbb{q}〉{\mathbb{Q}}_3.\hfill \end{array}$$

(44)

So, the tangent stiffness matrix is achieved as follows

$${\mathbb{K}}_t={\mathbb{K}}_l+{\mathbb{K}}_{t1}+{\mathbb{K}}_{t2},$$

(45)

$${\mathbb{K}}_{t1}=-\left[\begin{array}{ccc}{\tilde{\mathbb{R}}}_1\hfill & {\mathbb{R}}_1^{\mathrm{T}}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbb{R}}_1}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],{\mathbb{K}}_{t2}=-\left[\begin{array}{ccc}{\tilde{\mathbb{R}}}_2\hfill & {\mathbb{R}}_2^{\mathrm{T}}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbb{R}}_2}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right].$$

(46)

and as a result, the solution of equation (33) is obtained iteratively as ${\mathbb{d}}_{i+1}={\mathbb{d}}_i-{\mathbb{K}}_t^{-1}\left({\mathbb{d}}_i\right)\mathbb{V}\left({\mathbb{d}}_i\right)$ where the residual is given by $\mathbb{V}={\mathbb{K}}_s\mathbb{d}-\mathbb{F}$ and i is the iteration number. The flowchart of the proposed VFD solution for the nonlinear quadratic optimal control problems can be seen in Figure 1.

Figure 1.

Flowchart of the solution strategy.

Selected case studies

In this section, the accuracy and efficiency of the proposed VFD method for the analyzing of both linear and nonlinear optimal control problems are examined through several case studies. In different examples, the convergence of results based on the VFD method is investigated. Also, to examine the validity of present study, the obtained results are compared with those of available in open literature, that is, with the paper of Jaddu ⁶ and Mohammadi et al. ²⁴ It should be noted that for $n=140$ the convergence is achieved in all cases, as the low executed time for the analysis demonstrated the significant efficiency of this approach. Such a claim is supported by comparing the run-times of the presented VFD method with the recently developed archived-based genetic algorithm. ²⁴ The reported values are obtained from the parallel computation in MATLAB software with system properties as Core i7-4712MQ CPU @ 2.3 GHz, 8 GB RAM, 64 – bit. It should be mentioned that the linear solution is considered as the initial guess of Newton’s method and the relative error of convergence is ${10}^{-6}$ .

Linear time-invariant system (manned spacecraft)

As the first case study, consider the following cost function, state equations and conditions^24,33

$$\begin{array}{c}J=\frac{1}{2}{\int }_0^{\infty }(4x_1^2\left(t\right)+R{u}^2\left(t\right))\mathrm{dt},\\ {\stackrel{·}{x}}_1\left(t\right)=x_2\left(t\right),\\ {\stackrel{·}{x}}_2\left(t\right)=u\left(t\right),\\ x_1\left(0\right)=10,x_2\left(0\right)=0.\end{array}$$

Comparing with the represented matrix formulation, one has

$$\begin{array}{c}{t}_0=0,t_f=\infty ,m_1=2,m_2=1,\\ \mathbf{Q}=\left[\begin{array}{cc}4\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right],\mathbf{R}=R,\mathbf{A}=\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right],\mathbf{B}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\right],{\mathbf{R}}_{\mathrm{nl}}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right].\end{array}$$

It is assumed that $R=0.1,1,10,50$ . Practically, one can choose the terminate time as $t_f>10$ and thus, $t_f=15$ is considered in present numerical simulations. ²⁴ Also, from the boundary conditions, it reveals that $\mathbb{E}$ is a 2-by- $(m_1+m_2)n$ matrix to satisfy $\mathbb{q}\left(1\right)=10$ , $\mathbb{q}(n+1)=0$ . As a result, the second Lagrange multiplier vector ${\mathbb{L}}_2$ is a $1\times 2$ vector.

The produced numerical results are given in Table 1 and Figure 2. Considering different number of discretization points for the finite difference method, the values of performance index are reported in Table 1. Accordingly, the appropriate convergence behavior and validity of the proposed solution approach can be clearly observed. Corresponding to different values of R, time responses of the state variables are displayed in Figure 2. It is revealed that the same results of references are obtained via the VFD method. The average run time of this case is about 0.09 s.

Table 1.

Convergence and comparison study of J in the first numerical example (linear time-invariant system).

$R=0.1$	Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
		$J=215.853$	$J=218.525$	$J=220.992$	$J=221.281$	$J=221.287$
	Mohammadi et al. 24	$J=224.938$
$R=1$	Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
		$J=393.615$	$J=395.649$	$J=397.319$	$J=397.553$	$J=397.555$
	Mohammadi et al. 24	$J=399.999$
$R=10$	Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
		$J=707.038$	$J=708.377$	$J=709.137$	$J=709.614$	$J=709.623$
	Mohammadi et al. 24	$J=711.353$
$R=50$	Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
		$J=1060.28$	$J=1061.27$	$J=1061.94$	$J=1062.20$	$J=1062.27$
	Mohammadi et al. 24	$J=1064.54$

Figure 2.

Transient responses of the first case study (linear time-invariant system).

Linear time-variant system

The second numerical example is dedicated to the analysis of a time-variant system. The problem is identified with^24,34

$$\begin{array}{c}J=\frac{1}{2}{\int }_0^1(x^2\left(t\right)+u^2\left(t\right))\mathrm{dt},\\ \stackrel{·}{x}\left(t\right)=\mathrm{tx}\left(t\right)+u\left(t\right),\\ x\left(0\right)=1.\end{array}$$

So, the parameters are

$$\begin{array}{c}{t}_0=0,t_f=1,m_1=1,m_2=1,\\ \mathbf{Q}=1,\mathbf{R}=1,\mathbf{A}=t,\mathbf{B}=1,{\mathbf{R}}_{\mathrm{nl}}=0.\end{array}$$

According to what mentioned earlier, the discretized form of matrix $\mathbf{A}$ is defined as $\mathbb{A}=〈\mathbb{t}〉$ where the discretized time vector $\mathbb{t}$ is given in equation (5). Also, $\mathbb{E}$ is a 1-by- $(m_1+m_2)n$ matrix in which $\mathbb{q}\left(1\right)=1$ and ${\mathbb{L}}_2$ is a $1\times 1$ vector (scalar).

Table 2 shows the rich convergence of the proposed solution method to solve the time-variant problem. In Figure 3, the variation of the parameters x and u are plotted versus time. Also, a verification study is made with the reported results of literature. The time of analysis is 0.04 s which is far less than the previous attempts such as the ones based on genetic algorithm^24,35 or MPC.^25–27

Table 2.

Convergence study of J in the second numerical example (linear time-variant system).

Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
	$J=0.9683$	$J=0.9684$	$J=0.9684$	$J=0.9684$	$J=0.9684$
Mohammadi et al. 24	$J=2\times 0.4843=0.9686$

Figure 3.

Transient responses of the second case study (linear time-variant system).

Nonlinear system (Van der Pol problem)

In this case, the following optimal control problem with nonlinear state equations is considered^24,35

$$\begin{array}{c}J=\frac{1}{2}{\int }_0^{10}(x_1^2\left(t\right)+x_2^2\left(t\right)+u^2\left(t\right))\mathrm{dt},\\ {\stackrel{·}{x}}_1\left(t\right)=(1-x_2^2\left(t\right))x_1\left(t\right)-x_2\left(t\right)+u\left(t\right),\end{array}$$

$$\begin{array}{c}{\stackrel{·}{x}}_2\left(t\right)=x_1\left(t\right),\\ x_1\left(0\right)=0,x_2\left(0\right)=3,\\ x_1\left(10\right)=0,x_2\left(10\right)=0.\end{array}$$

So, one has

$$\begin{array}{c}{t}_0=0,t_f=10,m_1=2,m_2=1,\\ \mathbf{Q}=\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right],\mathbf{R}=1,\mathbf{A}=\left[\begin{array}{cc}1\hfill & -1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right],\mathbf{B}=\left[\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right],{\mathbf{R}}_{\mathrm{nl}}=\left[\begin{array}{c}-x_1{x}_2^2\hfill \\ \multicolumn{1}{c}{0}\end{array}\right].\end{array}$$

The discretized form of ${\mathbf{R}}_{\mathrm{nl}}$ and its derivative are obtained by taking the following matrices into account

$${\mathbb{Q}}_1={\mathbb{Q}}_2=\mathbf{0},{\mathbb{Q}}_3=\left[\begin{array}{ccc}-{\mathbf{I}}_n\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],{\mathbb{Q}}_4={\mathbb{Q}}_5=\left[\begin{array}{ccc}\mathbf{0}\hfill & {\mathbf{I}}_n\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right].$$

From boundary conditions, $\mathbb{E}$ is a 4-by- $(m_1+m_2)n$ matrix so that $\mathbb{q}\left(1\right)=\mathbb{q}\left(n\right)=\mathbb{q}\left(2n\right)=0$ , $\mathbb{q}(n+1)=3$ and thus, ${\mathbb{L}}_2$ is a $1\times 4$ vector.

Now, like the previous case studies, the convergence behavior and transient responses of the system are studied. From Table 3, the converged performance index can be compared with different values reported in Mohammadi et al. ²⁴ It is observed that there is a very good agreement with the more precise case ( $\mathrm{depth}=9$ ). Illustrated in Figure 4 is the variation of three state variables with time which also verified the present model once again. Furthermore, the run time of nonlinear analysis is 0.52 s corresponding to seven iterations in the Newton’s method. It should be mentioned that the developed analytical Jacobian (tangent stiffness matrix) plays a key role in increasing such an efficiency.

Table 3.

Convergence and comparison study of J in the third numerical example (first nonlinear system).

Present	$n=80$	$n=100$		$n=120$	$n=140$		$n=200$
	$J=21.4299$	$J=21.4191$		$J=21.4160$	$J=21.4152$		$J=21.4152$
Mohammadi et al. 24	$\mathrm{depth}=5$		$\mathrm{depth}=7$			$\mathrm{depth}=9$
	$J=21.4305$		$J=21.4247$			$J=21.4181$

Figure 4.

Transient responses of the third case study (first nonlinear system).

Nonlinear system (Rigid asymmetric spacecraft)

As the last example, the VFD method is used to solve the following nonlinear system ⁶

$$\begin{array}{c}J=\frac{1}{2}{\int }_0^{100}(u_1^2\left(t\right)+u_2^2\left(t\right)+u_3^2\left(t\right))\mathrm{dt},\\ {\stackrel{·}{x}}_1\left(t\right)=-\frac{(I_3-I_2)}{I_1}{x}_2{x}_3+\frac{u_1}{I_1},\end{array}$$

$$\begin{array}{c}{\stackrel{·}{x}}_2\left(t\right)=-\frac{(I_1-I_3)}{I_2}{x}_1{x}_3+\frac{u_2}{I_2},\\ {\stackrel{·}{x}}_3\left(t\right)=-\frac{(I_2-I_1)}{I_3}{x}_1{x}_2+\frac{u_3}{I_3},\\ x_1\left(0\right)=0.01,x_2\left(0\right)=0.005,x_3\left(0\right)=0.001,\\ x_1\left(100\right)=x_2\left(100\right)=x_3\left(100\right)=0.\end{array}$$

where $I_1=86.24$ , $I_2=85.07$ and $I_3=113.59$ . The parameters of matrix formulation are

$$\begin{array}{c}{t}_0=0,t_f=100,m_1=3,m_2=3,\\ \begin{array}{cc}\multicolumn{1}{c}{}& \mathbf{Q}=\left[\begin{array}{ccc}\mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],\mathbf{R}=\left[\begin{array}{ccc}1\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& 1\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & 1\hfill \end{array}\right],\mathbf{A}=\left[\begin{array}{ccc}\mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],\hfill \\ \multicolumn{1}{c}{}& \mathbf{B}=\left[\begin{array}{ccc}1/I_1\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& 1/I_2\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & 1/I_3\hfill \end{array}\right],\hfill \end{array}\\ {\mathbf{R}}_{\mathrm{nl}}=\left[\begin{array}{c}-((I_3-I_2)/I_1)x_2{x}_3\hfill \\ \multicolumn{1}{c}{-((I_1-I_3)/I_2)x_1{x}_3}\\ \multicolumn{1}{c}{-((I_2-I_1)/I_3)x_1{x}_2}\end{array}\right].\end{array}$$

Also, ${\mathbb{R}}_{\mathrm{nl}}$ is characterized by

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mathbb{Q}}_1=\left[\begin{array}{cccccc}\mathbf{0}\hfill & -((I_3-I_2)/I_1){\mathbf{I}}_n\hfill & 0\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{-((I_1-I_3)/I_2){\mathbf{I}}_n}& \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{-((I_2-I_1)/I_3){\mathbf{I}}_n}& \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],\hfill \\ \multicolumn{1}{c}{}& {\mathbb{Q}}_2=\left[\begin{array}{cccccc}\mathbf{0}\hfill & \mathbf{0}\hfill & {\mathbf{I}}_n\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& \mathbf{0}\hfill & {\mathbf{I}}_n\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{\mathbf{0}}& {\mathbf{I}}_n\hfill & 0\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right],{\mathbb{Q}}_3={\mathbb{Q}}_4={\mathbb{Q}}_5=\mathbf{0}.\hfill \end{array}$$

and $\mathbb{E}$ is a 6-by- $(m_1+m_2)n$ matrix in which $\mathbb{q}\left(n\right)=\mathbb{q}\left(2n\right)=\mathbb{q}\left(3n\right)=0$ , $\mathbb{q}\left(1\right)=0.01$ , $\mathbb{q}(n+1)=0.005$ , $\mathbb{q}(2n+1)=0.001$ and ${\mathbb{L}}_2$ is a $1\times 6$ vector.

The numerical results of this example involving convergence and comparative studies as well as time responses can be found in Table 4 and Figure 5. Comparing the performance index and transient solutions for the six state parameters, it is observed that the obtained results are in excellent agreement with Jaddu. ⁶ Also, the time of nonlinear analysis of the optimal control problem is 0.67 s for two iterations. Like the previous nonlinear case study, it is seen that the present solution methodology is reliable and efficient. Comparing with the computational demand of the available nonlinear programming algorithms, the executed run time of nonlinear case studies could be promising for future attempts in the area of optimal control theory.

Table 4.

Convergence and comparison study of $J\times 100$ in the fourth numerical example (second nonlinear system).

Present	$n=80$	$n=100$	$n=120$	$n=140$	$n=200$
	$J=0.46858$	$J=0.46863$	$J=0.46867$	$J=0.46869$	$J=0.46870$
Jaddu 6	$J=0.46878$

Figure 5.

Transient responses of the fourth case study (second nonlinear system).

Computational efficiency

The straightforward mathematical formulation as well as the accuracy of VFD method has been shown so far. Here, some comparisons are made to show the computational efficiency of the variational finite difference approach in context of the optimal control theory. In this regard, the numerical case studies are solved via the archived-based genetic programming as the recent attempt in this research area. ²⁴ To be ensure of the accuracy of the run-times reported for the proposed VFD method, several runs were performed and the average values are considered. In addition, tabulated in Table 5 are the executed run-times obtained from AGP solution of the single-objective optimal control problems with 200 generations. Results are quite promising and clearly reveal the efficiency of the suggested solution strategy.

Table 5.

Comparison of run-times between VFD and AGP methods.

Case study	Run-time (seconds)
	VFD	AGP
Linear time-invariant system (manned spacecraft)	0.09	5262
Linear time-variant system	0.04	5115
Nonlinear system (Van der Pol problem)	0.52	5338
Nonlinear system (Rigid asymmetric spacecraft)	0.67	5398

Concluding remarks

A novel solution was proposed in the context of nonlinear quadratic optimal control problems using variational time discretized finite difference approach inspired from the variational differential quadrature method. Also, with the aim of straightforward implementation of solution strategy, a matrix-vector formulation was conveniently developed for the performance index, state equations, and boundary conditions usually represented in optimal control problems. The resulted algebraic equations of the minimization problem were then solved with the help of the Lagrange multiplier method. Due to utilizing consistent integral/differential operators and the analytical tangent stiffness matrix rather than the numerical Jacobian, the developed methodology is robust and more efficient in comparison with the available solutions in literature. Such a claim was supported with successfully analyzing different case studies including the linear time-invariant, linear time-variant and nonlinear systems. The significant efficiency of the developed VFD method could be find from the fact that the run time of analysis was less than 1 s even in case of nonlinear optimal control problems. In future communications, the proposed approach will be extended for the case of multi-dimensional nonlinear optimal control problems.