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.

Keywords

Optimal control nonlinearity numerical analysis variational finite difference

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}} (x^{T} (t) Qx (t) + u^{T} (t) Ru (t)) dt,

(1)

subjected to the state equation

\overset{\cdot}{x} (t) = Ax (t) + Bu (t) + R_{nl} (x, u),

(2)

and initial/terminate conditions

x (t_{0}) = x_{t_{0}}, x (t_{f}) = x_{t_{f}},

(3)

where the state variables are

x = {[\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \end{matrix}]}_{m_{1} \times 1}, u = {[\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \end{matrix}]}_{m_{2} \times 1},

(4)

and ${[R_{nl}]}_{m_{1} \times 1}$ denotes the nonlinear vector in state equation. Also, ${[Q]}_{m_{1} \times m_{1}}$ , ${[R]}_{m_{2} \times m_{2}}$ signal the weighted matrices and ${[A]}_{m_{1} \times m_{1}}$ , ${[B]}_{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 $x$ and $u$ vectors.

Solution strategy

In order to solve the optimal control problem of finding $x$ and $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

t = {[\begin{matrix} t_{0} & t_{0} + Δ t & t_{0} + 2 Δ t & \dots & t_{f} \end{matrix}]}^{T}, t_{f} = t_{0} + (n - 1) Δ t .

(5)

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

x = {[\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \end{matrix}]}_{m_{1} n \times 1}, u = {[\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \end{matrix}]}_{m_{2} n \times 1},

(6)

\begin{matrix} x_{1} = {[\begin{matrix} x_{1} (t_{0}) & x_{1} (t_{0} + Δ t) & x_{1} (t_{0} + 2 Δ t) & \dots & x_{1} (t_{f}) \end{matrix}]}^{T}, \\ x_{2} = {[\begin{matrix} x_{2} (t_{0}) & x_{2} (t_{0} + Δ t) & x_{2} (t_{0} + 2 Δ t) & \dots & x_{2} (t_{f}) \end{matrix}]}^{T}, \dots, \end{matrix}

(7)

\begin{matrix} u_{1} = {[\begin{matrix} u_{1} (t_{0}) & u_{1} (t_{0} + Δ t) & u_{1} (t_{0} + 2 Δ t) & \dots & u_{1} (t_{f}) \end{matrix}]}^{T}, \\ u_{2} = {[\begin{matrix} u_{2} (t_{0}) & u_{2} (t_{0} + Δ t) & u_{2} (t_{0} + 2 Δ t) & \dots & u_{2} (t_{f}) \end{matrix}]}^{T}, \dots, \end{matrix}

(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

q = [\begin{matrix} x \\ u \end{matrix}] .

(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 (t)$ as vector ${[f]}_{n \times 1}$ , the r-th derivative of $f (t)$ in whole time domain is given by

\frac{d^{r} f}{d t^{r}} = D^{(r)} f,

(10)

where, ${[D^{(r)}]}_{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 ${[S]}_{1 \times n}$ (see Appendix B). It leads to

\int_{t_{0}}^{t_{f}} f (t) dt = S 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} S (Q_{ij} x_{i} ° x_{j} + R_{ij} u_{i} ° u_{j}),

(12)

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

{[X]}_{p_{1} \times p_{2}} ° {[Y]}_{p_{1} \times p_{2}} = {[\begin{matrix} X_{11} Y_{11} & X_{12} Y_{12} & \dots \\ X_{21} Y_{21} & X_{22} Y_{22} & \dots \\ ⋮ & ⋮ & ⋱ \end{matrix}]}_{p_{1} \times p_{2}} .

(13)

from the lemma $V_{1}^{T} (V_{2} ° V_{3}) = V_{2}^{T} V_{1} V_{3}$ for three arbitrary column vectors and considering that $〈 V 〉$ is the diagonal form of vector $V$ , equation (12) is rewritten as

J = \frac{1}{2} (x_{i}^{T} Q_{ij} 〈 S 〉 x_{j} + u_{i}^{T} R_{ij} 〈 S 〉 u_{j}),

(14)

which can also be given in the following total matrix representation

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

(15)

where $⊛$ is the symbol of Kronecker product defined as

{[X]}_{p_{1} \times p_{2}} ⊛ {[Y]}_{q_{1} \times q_{2}} = {[\begin{matrix} X_{11} Y & X_{12} Y & \dots \\ X_{21} Y & X_{22} Y & \dots \\ ⋮ & ⋮ & ⋱ \end{matrix}]}_{p_{1} q_{1} \times p_{2} q_{2}}

(16)

Also, discretization of equation (2) leads to

(I_{m_{1}} ⊛ D^{(1)}) x = (A ⊛ I_{n}) x + (B ⊛ I_{n}) u + R_{nl},

(17)

where $I_{p}$ denotes the p-by-p identity matrix and ${[R_{nl}]}_{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} q^{T} Sq,

(18)

C q - R_{nl} = 0,

(19)

E q - x_{0} = 0,

(20)

where

S = [\begin{matrix} Q ⊛ 〈 S 〉 & 0 \\ 0 & R ⊛ 〈 S 〉 \end{matrix}],

(21)

ℂ = [\begin{matrix} D - A & - B \end{matrix}], D = I_{m_{1}} ⊛ D^{(1)}, A = A ⊛ I_{n}, B = B ⊛ I_{n} .

(22)

and $E$ shows a matrix containing zero and one for the initial and terminate values denoted by the vector $x_{0}$ . It is worth mentioning that the discretized formulation of the vector $R_{nl}$ depends on the nature of nonlinearity in the problem. It should be also noted that $A$ and $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, $A$ and $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{matrix} minimize J (equation (18) : J = \frac{1}{2} q^{T} S q), \\ \begin{matrix} subject to equation (19) (C q - R_{nl} = 0) and \\ equation (2) 0 (E q - x_{0} = 0), \end{matrix} \end{matrix}

(23)

and the Lagrange function is given by

L (q, L_{1}, L_{2}) = \frac{1}{2} q^{T} S q + L_{1} (C q - R_{nl}) + L_{2} (E q - x_{0}) .

(24)

where $L_{1}$ is the first vector of Lagrange multipliers with size of $1 \times m_{1} n$ and $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 L (q, L_{1}, L_{2}) = 0$ ). Equivalently, the first variation of Lagrange function is calculated as

\begin{matrix} δ L & = (\frac{\partial L}{\partial q}) \cdot δ q + (\frac{\partial L}{\partial L_{1}}) \cdot δ L_{1}^{T} + (\frac{\partial L}{\partial L_{2}}) \cdot δ L_{2}^{T} \\ = δ q^{T} (S q + {(L_{1} (C - \partial R_{nl} / \partial q))}^{T} + {(L_{2} E)}^{T}) \\ + δ L_{1} (C q - R_{nl}) + δ L_{2} (E q - x_{0}) . \end{matrix}

(25)

The necessary condition of vanishing the first variation of the Lagrange function, that is, $δ L = 0$ , is achieved by letting the coefficients of independent vectors $δ_{L}$ , $δ L_{1}$ , and $δ L_{2}$ set to zeros. It, therefore, results that

S q + {(L_{1} (C - \partial R_{nl} / \partial q))}^{T} + {(L_{2} E)}^{T} = 0,

(26)

C q - R_{nl} = 0,

(27)

E q - x_{0} = 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

(Q_{1} q) ° (Q_{2} q) = \frac{1}{2} (〈 Q_{1} q 〉 Q_{2} + 〈 Q_{2} q 〉 Q_{1}) q = \frac{1}{2} R_{1} q,

(29)

(ℚ_{3} q) \circ (ℚ_{4} q) \circ (ℚ_{5} q) = \frac{1}{3} (ℚ_{3} q ℚ_{4} q ℚ_{5} + ℚ_{3} q ℚ_{5} q ℚ_{4} + ℚ_{4} q ℚ_{5} q ℚ_{3}) q = \frac{1}{3} ℝ_{2} q,

(30)

where $R_{1}$ , $R_{2}$ , $Q_{1}$ , …, $Q_{5}$ are matrices of order $m_{1} n \times (m_{1} + m_{2}) n$ and one can generally write

R_{nl} = (\frac{1}{2} R_{1} + \frac{1}{3} R_{2}) q,

(31)

which results that (see Appendix C)

\frac{\partial R_{nl}}{\partial q} = R_{1} + R_{2} .

(32)

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

K_{s} d = F,

(33)

where

\begin{matrix} d = [\begin{matrix} q \\ L_{1}^{T} \\ L_{2}^{T} \end{matrix}], \end{matrix}

(34)

F = [\begin{matrix} 0 \\ 0 \\ x_{0} \end{matrix}],

(35)

and the stiffness matrix $K_{s}$ is given by

K_{s} = K_{l} + \frac{1}{2} K_{n 1} + \frac{1}{3} K_{n 2},

(36)

\begin{matrix} K_{l} = [\begin{matrix} S & C^{T} & E^{T} \\ C & 0 & 0 \\ E & 0 & 0 \end{matrix}], K_{n 1} = - [\begin{matrix} 0 & 2 R_{1}^{T} & 0 \\ R_{1} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], \\ K_{n 2} = - [\begin{matrix} 0 & 3 R_{2}^{T} & 0 \\ R_{2} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}] . \end{matrix}

(37)

$K_{l}$ is the linear matrix, $K_{n 1}$ and $K_{n 2}$ 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 $Δ δ L = δ d^{T} K_{t} Δ d$ in which $K_{t}$ is the tangent stiffness matrix. Since $Δ δ q = Δ δ L_{1} = Δ δ L_{2} = 0$ , the increment of variation of the Lagrange function becomes

\begin{matrix} Δ δ L & = δ q^{T} (S Δ q + C^{T} Δ L_{1}^{T} - Δ {{(R_{1} + R_{2})}^{T} L_{1}^{T}} + E^{T} Δ L_{2}^{T}) \\ + δ L_{1} (C Δ q - Δ {(\frac{1}{2} R_{1} + \frac{1}{3} R_{2}) q}) + δ L_{2} E Δ q, \end{matrix}

(38)

considering

Δ R_{1} = 〈 Q_{1} Δ q 〉 Q_{2} + 〈 Q_{2} Δ q 〉 Q_{1},

(39)

\begin{matrix} Δ R_{2} & = 〈 Q_{3} Δ q 〉 〈 Q_{4} q 〉 Q_{5} + 〈 Q_{3} q 〉 〈 Q_{4} Δ q 〉 Q_{5} \\ + 〈 Q_{3} Δ q 〉 〈 Q_{5} q 〉 Q_{4} + 〈 Q_{3} q 〉 〈 Q_{5} Δ q 〉 Q_{4} \\ + 〈 Q_{4} Δ q 〉 〈 Q_{5} q 〉 Q_{3} + 〈 Q_{4} q 〉 〈 Q_{5} Δ q 〉 Q_{3}, \end{matrix}

(40)

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

\begin{matrix} Δ {{(R_{1} + R_{2})}^{T} L_{1}^{T}} & = {(Δ R_{1} + Δ R_{2})}^{T} L_{1}^{T} + {(R_{1} + R_{2})}^{T} Δ L_{1}^{T} \\ = ({\tilde{R}}_{1} + {\tilde{R}}_{2}) Δ q + {(R_{1} + R_{2})}^{T} Δ L_{1}^{T}, \end{matrix}

(41)

\begin{matrix} Δ {(\frac{1}{2} R_{1} + \frac{1}{3} R_{2}) q} & = (\frac{1}{2} Δ R_{1} + \frac{1}{3} Δ R_{2}) q \\ + (\frac{1}{2} R_{1} + \frac{1}{3} R_{2}) Δ q = (R_{1} + R_{2}) Δ q, \end{matrix}

(42)

where

{\tilde{R}}_{1} = Q_{1}^{T} 〈 L 〉_{1}^{T} Q_{2} + Q_{2}^{T} 〈 L 〉_{1}^{T} Q_{1},

(43)

\begin{matrix} {\tilde{R}}_{2} & = Q_{3}^{T} 〈 L_{1}^{T} 〉 〈 Q_{4} q 〉 Q_{5} + Q_{3}^{T} 〈 L_{1}^{T} 〉 〈 Q_{5} q 〉 Q_{4} \\ + Q_{4}^{T} 〈 L_{1}^{T} 〉 〈 Q_{3} q 〉 Q_{5} + Q_{4}^{T} 〈 L_{1}^{T} 〉 〈 Q_{5} q 〉 Q_{3} \\ + Q_{5}^{T} 〈 L_{1}^{T} 〉 〈 Q_{3} q 〉 Q_{4} + Q_{5}^{T} 〈 L_{1}^{T} 〉 〈 Q_{4} q 〉 Q_{3} . \end{matrix}

(44)

So, the tangent stiffness matrix is achieved as follows

K_{t} = K_{l} + K_{t 1} + K_{t 2},

(45)

K_{t 1} = - [\begin{matrix} {\tilde{R}}_{1} & R_{1}^{T} & 0 \\ R_{1} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], K_{t 2} = - [\begin{matrix} {\tilde{R}}_{2} & R_{2}^{T} & 0 \\ R_{2} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}] .

(46)

and as a result, the solution of equation (33) is obtained iteratively as $d_{i + 1} = d_{i} - K_{t}^{- 1} (d_{i}) V (d_{i})$ where the residual is given by $V = K_{s} d - 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{matrix} J = \frac{1}{2} \int_{0}^{\infty} (4 x_{1}^{2} (t) + R u^{2} (t)) dt, \\ {\overset{\cdot}{x}}_{1} (t) = x_{2} (t), \\ {\overset{\cdot}{x}}_{2} (t) = u (t), \\ x_{1} (0) = 10, x_{2} (0) = 0 . \end{matrix}

Comparing with the represented matrix formulation, one has

\begin{matrix} t_{0} = 0, t_{f} = \infty, m_{1} = 2, m_{2} = 1, \\ Q = [\begin{matrix} 4 & 0 \\ 0 & 0 \end{matrix}], R = R, A = [\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}], B = [\begin{matrix} 0 \\ 1 \end{matrix}], R_{nl} = [\begin{matrix} 0 \\ 0 \end{matrix}] . \end{matrix}

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 $E$ is a 2-by- $(m_{1} + m_{2}) n$ matrix to satisfy $q (1) = 10$ , $q (n + 1) = 0$ . As a result, the second Lagrange multiplier vector $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.²⁴	$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.²⁴	$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.²⁴	$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.²⁴	$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{matrix} J = \frac{1}{2} \int_{0}^{1} (x^{2} (t) + u^{2} (t)) dt, \\ \overset{\cdot}{x} (t) = tx (t) + u (t), \\ x (0) = 1 . \end{matrix}

So, the parameters are

\begin{matrix} t_{0} = 0, t_{f} = 1, m_{1} = 1, m_{2} = 1, \\ Q = 1, R = 1, A = t, B = 1, R_{nl} = 0 . \end{matrix}

According to what mentioned earlier, the discretized form of matrix $A$ is defined as $A = 〈 t 〉$ where the discretized time vector $t$ is given in equation (5). Also, $E$ is a 1-by- $(m_{1} + m_{2}) n$ matrix in which $q (1) = 1$ and $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.²⁴	$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{matrix} J = \frac{1}{2} \int_{0}^{10} (x_{1}^{2} (t) + x_{2}^{2} (t) + u^{2} (t)) dt, \\ {\overset{\cdot}{x}}_{1} (t) = (1 - x_{2}^{2} (t)) x_{1} (t) - x_{2} (t) + u (t), \end{matrix}

\begin{matrix} {\overset{\cdot}{x}}_{2} (t) = x_{1} (t), \\ x_{1} (0) = 0, x_{2} (0) = 3, \\ x_{1} (10) = 0, x_{2} (10) = 0 . \end{matrix}

So, one has

\begin{matrix} t_{0} = 0, t_{f} = 10, m_{1} = 2, m_{2} = 1, \\ Q = [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}], R = 1, A = [\begin{matrix} 1 & - 1 \\ 1 & 0 \end{matrix}], B = [\begin{matrix} 1 \\ 0 \end{matrix}], R_{nl} = [\begin{matrix} - x_{1} x_{2}^{2} \\ 0 \end{matrix}] . \end{matrix}

The discretized form of $R_{nl}$ and its derivative are obtained by taking the following matrices into account

Q_{1} = Q_{2} = 0, Q_{3} = [\begin{matrix} - I_{n} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], Q_{4} = Q_{5} = [\begin{matrix} 0 & I_{n} & 0 \\ 0 & 0 & 0 \end{matrix}] .

From boundary conditions, $E$ is a 4-by- $(m_{1} + m_{2}) n$ matrix so that $q (1) = q (n) = q (2 n) = 0$ , $q (n + 1) = 3$ and thus, $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 ( $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.²⁴	$depth = 5$		$depth = 7$			$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{matrix} J = \frac{1}{2} \int_{0}^{100} (u_{1}^{2} (t) + u_{2}^{2} (t) + u_{3}^{2} (t)) dt, \\ {\overset{\cdot}{x}}_{1} (t) = - \frac{(I_{3} - I_{2})}{I_{1}} x_{2} x_{3} + \frac{u_{1}}{I_{1}}, \end{matrix}

\begin{matrix} {\overset{\cdot}{x}}_{2} (t) = - \frac{(I_{1} - I_{3})}{I_{2}} x_{1} x_{3} + \frac{u_{2}}{I_{2}}, \\ {\overset{\cdot}{x}}_{3} (t) = - \frac{(I_{2} - I_{1})}{I_{3}} x_{1} x_{2} + \frac{u_{3}}{I_{3}}, \\ x_{1} (0) = 0.01, x_{2} (0) = 0.005, x_{3} (0) = 0.001, \\ x_{1} (100) = x_{2} (100) = x_{3} (100) = 0 . \end{matrix}

where $I_{1} = 86.24$ , $I_{2} = 85.07$ and $I_{3} = 113.59$ . The parameters of matrix formulation are

\begin{matrix} t_{0} = 0, t_{f} = 100, m_{1} = 3, m_{2} = 3, \\ \begin{matrix} Q = [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], R = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}], A = [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], \\ B = [\begin{matrix} 1 / I_{1} & 0 & 0 \\ 0 & 1 / I_{2} & 0 \\ 0 & 0 & 1 / I_{3} \end{matrix}], \end{matrix} \\ R_{nl} = [\begin{matrix} - ((I_{3} - I_{2}) / I_{1}) x_{2} x_{3} \\ - ((I_{1} - I_{3}) / I_{2}) x_{1} x_{3} \\ - ((I_{2} - I_{1}) / I_{3}) x_{1} x_{2} \end{matrix}] . \end{matrix}

Also, $R_{nl}$ is characterized by

\begin{matrix} Q_{1} = [\begin{matrix} 0 & - ((I_{3} - I_{2}) / I_{1}) I_{n} & 0 & 0 & 0 & 0 \\ - ((I_{1} - I_{3}) / I_{2}) I_{n} & 0 & 0 & 0 & 0 & 0 \\ - ((I_{2} - I_{1}) / I_{3}) I_{n} & 0 & 0 & 0 & 0 & 0 \end{matrix}], \\ Q_{2} = [\begin{matrix} 0 & 0 & I_{n} & 0 & 0 & 0 \\ 0 & 0 & I_{n} & 0 & 0 & 0 \\ 0 & I_{n} & 0 & 0 & 0 & 0 \end{matrix}], Q_{3} = Q_{4} = Q_{5} = 0 . \end{matrix}

and $E$ is a 6-by- $(m_{1} + m_{2}) n$ matrix in which $q (n) = q (2 n) = q (3 n) = 0$ , $q (1) = 0.01$ , $q (n + 1) = 0.005$ , $q (2 n + 1) = 0.001$ and $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⁶	$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.

Footnotes

Appendix A

To evaluate derivatives of function $f (t)$ with respect to t, different approximation schemes of the finite difference method can be used. Considering the equally-space discretization with interval h and the first- to fourth-order derivations, the coefficients of the central as well as forward/backward approximations are given in the following tables.

For instance

(A-1)

f^{(1)} (t) = \frac{f (t + h) - f (t - h)}{2 h} + O (h^{2}) .

For instance

(A-2)

f^{(1)} (t) = \frac{f (t + h) - f (t)}{h} + O (h),

(A-3)

f^{(1)} (t) = \frac{f (t) - f (t - h)}{h} + O (h) .

So, one can write $d^{r} f / d t^{r} = D^{(r)} f,$ where $f$ denotes the discretized vector of $f (t)$ and $D^{(r)}$ is the r-th order derivative operator matrix produced from Tables A.1 and A.2.

Appendix B

According to Shojaei and Ansari (28), the Taylor series representation of a function is used to evaluate the integral operator at the same points on which the differentiation is calculated. In a neighborhood of a typical point $t_{i}$ , where $i = 1, \dots, n$ , function $f (t)$ becomes

(B-1)

f (t) = \sum_{r = 0}^{\infty} \frac{f^{(r)} (t_{i})}{r!} {(t - t_{i})}^{r} .

Since $f^{(r)} (t_{i})$ can be evaluated by the numerical differentiation introduced in Appendix A, it is convenient to integrate the function in arbitrary domain; namely

(B-2)

\begin{matrix} \int_{(t_{i - 1} + t_{i}) / 2}^{(t_{i} + t_{i + 1}) / 2} f (t) dt = \sum_{r = 0}^{\infty} \frac{f^{(r)} (t_{i})}{2^{r + 1} (r + 1)!} \\ ({(t_{i + 1} - t_{i})}^{r + 1} - {(t_{i - 1} - t_{i})}^{r + 1}) . \end{matrix}

Considering the whole intervals of discretization scheme, it results that $\int_{t_{0}}^{t_{f}} f (t) dt = S f$ in which $f$ is the discretized vector of $f (t)$ and the integral operator $S$ is given by

(B-3)

S = \sum_{r = 0}^{\infty} T^{r} D^{(r)},

where, $D^{(r)}$ is the introduced differential operator of finite difference method and

(B-4)

\begin{matrix} T^{r} = \\ [\begin{matrix} \frac{{(t_{2} - t_{1})}^{r + 1}}{2^{r + 1} (r + 1)!} & \dots & \frac{{(t_{i + 1} - t_{i})}^{r + 1} - {(t_{i - 1} - t_{i})}^{r + 1}}{2^{r + 1} (r + 1)!} & \dots & - \frac{{(t_{n - 1} - t_{n})}^{r + 1}}{2^{r + 1} (r + 1)!} \end{matrix}] . \end{matrix}

Appendix C

The variation of discretized nonlinear state vector (equation (31)) is calculated as

(C-1)

δ R_{nl} = (\frac{1}{2} δ R_{1} + \frac{1}{3} δ R_{2}) q + (\frac{1}{2} R_{1} + \frac{1}{3} R_{2}) δ q .

Considering equation (29), equation (30) and the lemma $〈 X 〉 Y = 〈 Y 〉 X$ for two column vectors, one has

(C-2)

\begin{matrix} δ R_{1} q & = 〈 Q_{1} δ q 〉 Q_{2} q + 〈 Q_{2} δ q 〉 Q_{1} q = 〈 Q_{1} q 〉 Q_{2} δ q \\ + Q_{2} q Q_{1} δ q = R_{1} δ q, \end{matrix}

\begin{matrix} δ R_{2} q & = 〈 Q_{3} δ q 〉 〈 Q_{4} q 〉 Q_{5} q + 〈 Q_{3} q 〉 〈 Q_{4} δ q 〉 Q_{5} q \\ + 〈 Q_{3} δ q 〉 〈 Q_{5} q 〉 Q_{4} q + 〈 Q_{3} q 〉 〈 Q_{5} δ q 〉 Q_{4} q \\ + 〈 Q_{4} δ q 〉 〈 Q_{5} q 〉 Q_{3} q + 〈 Q_{4} q 〉 〈 Q_{5} δ q 〉 Q_{3} q \\ = 2 〈 Q_{3} q 〉 〈 Q_{4} q 〉 Q_{5} δ q + 2 〈 Q_{4} q 〉 〈 Q_{5} q 〉 Q_{3} δ q \\ + 2 〈 Q_{3} q 〉 〈 Q_{5} q 〉 Q_{4} δ q = 2 R_{2} δ q, \end{matrix}

and thus, it is shown that $δ R_{nl} = (R_{1} + R_{2}) δ q$ .

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Amir Norouzzadeh

Data availability

The datasets generated and analyzed during the current study are not publicly available due the data also forms part of an ongoing study but are available from the corresponding authors on reasonable request.

References

Lewis

Vrabie

Syrmos

VL.

Optimal control. Hoboken, NJ: Wiley, 2012.

Athans

Falb

PL.

Optimal control: an introduction to the theory and its applications. Chelmsford, MA: Courier Corporation, 2007.

Zhu

Trélat

Cerf

Geometric optimal control and applications to aerospace. Pac J Math Ind 2017; 9: 1–41.

Kalman

RE.

Contributions to the theory of optimal control. Bol Soc Mat Mexicana 1960; 5(2): 102–119.

Tsai

J-H

Shieh

L-S.

Discretized quadratic optimal control for continuous-time two-dimensional systems. IEEE Trans Circuits Syst I Fundam Theor Appl 2002; 49(1): 116–125.

Jaddu

Direct solution of nonlinear optimal control problems using quasilinearization and Chebyshev polynomials. J Franklin Inst 2002; 339(4–5): 479–498.

Huang

C-Y.

A three-dimensional optimal control problem in determining the boundary control heat fluxes. Heat Mass Transf 2003; 39: 589–598.

Gerdts

Greif

Pesch

HJ.

Numerical optimal control of the wave equation: optimal boundary control of a string to rest in finite time. Math Comp Simul 2008; 79(4): 1020–1032.

Tricaud

Chen

An approximate method for numerically solving fractional order optimal control problems of general form. Comp Math Appl 2010; 59(5): 1644–1655.

10.

Casas

Tröltzsch

Numerical analysis of some optimal control problems governed by a class of quasilinear elliptic equations. ESAIM Control Optim Calc Var 2011; 17(3): 771–800.

11.

Yousefi

Lotfi

Dehghan

The use of a Legendre multiwavelet collocation method for solving the fractional optimal control problems. J Vib Control 2011; 17(13): 2059–2065.

12.

Neitzel

Tröltzsch

. Numerical analysis of state-constrained optimal control problems for PDEs. In: Leugering

(ed.) Constrained optimization and optimal control for partial differential equations. Cham: Springer, 2012, pp.467–482.

13.

Foroozandeh

Shamsi

Solution of nonlinear optimal control problems by the interpolating scaling functions. Acta Astronaut 2012; 72: 21–26.

14.

Kammann

Tröltzsch

Volkwein

A posteriori error estimation for semilinear parabolic optimal control problems with application to model reduction by POD. ESAIM Math Model Numeric Anal 2013; 47(2): 555–581.

15.

Bhrawy

Doha

Tenreiro Machado

, et al. An efficient numerical scheme for solving multi-dimensional fractional optimal control problems with a quadratic performance index. Asian J Control 2015; 17(6): 2389–2402.

16.

Mamehrashi

Yousefi

SA.

Numerical solution of a class of two-dimensional quadratic optimal control problems by using Ritz method. Optimal Control Appl Methods 2016; 37(4): 765–781.

17.

Keshavarz

Ordokhani

Razzaghi

A numerical solution for fractional optimal control problems via Bernoulli polynomials. J Vib Control 2016; 22(18): 3889–3903.

18.

Mamehrashi

Yousefi

SA.

A numerical method for solving a nonlinear 2-D optimal control problem with the classical diffusion equation. Int J Control 2017; 90(2): 298–306.

19.

Lotfi

A combination of variational and penalty methods for solving a class of fractional optimal control problems. J Optimiz Theor Appl 2017; 174: 65–82.

20.

Bourdin

Trélat

Linear–quadratic optimal sampled-data control problems: convergence result and Riccati theory. Automatica 2017; 79: 273–281.

21.

Chen

Shi

Zhang

Optimal control via neural networks: a convex approach. arXiv preprint arXiv: 180511835, 2018.

22.

Hajiloo

Nariman-Zadeh

Moeini

Pareto optimal robust design of fractional-order PID controllers for systems with probabilistic uncertainties. Mechatronics 2012; 22(6): 788–801.

23.

Yıldırım Aksoy

Hào

Yagub

. Finite difference method for an optimal control problem for a nonlinear time-dependent Schrödinger equation. Numeri Function Anal Optimiz 2017; 38(6): 799–817.

24.

Mohammadi

Nariman-Zadeh

Jamali

The archived-based genetic programming for optimal design of linear/non-linear controllers. Trans Inst Meas Control 2020; 42(8): 1475–1491.

25.

Diehl

Ferreau

Haverbeke

. Efficient numerical methods for nonlinear MPC and moving horizon estimation. In: Magni

Raimondo

Allgöwer

(eds.) Nonlinear model predictive control: towards new challenging applications. Heidelberg: Springer, 2009, pp.391–417.

26.

Zhang

Liao

Near-optimal control of nonlinear dynamical systems: a brief survey. Ann Rev Control 2019; 47: 71–80.

27.

Vukov

Gros

Horn

, et al. Real-time nonlinear MPC and MHE for a large-scale mechatronic application. Control Eng Pract 2015; 45: 64–78.

28.

Shojaei

Ansari

Variational differential quadrature: a technique to simplify numerical analysis of structures. Appl Math Modell 2017; 49: 705–738.

29.

Ansari

Hasrati

Shakouri

, et al. Nonlinear large deformation analysis of shells using the variational differential quadrature method based on the six-parameter shell theory. Int J Non-Linear Mech 2018; 106: 130–143.

30.

Hasrati

Ansari

Rouhi

Elastoplastic postbuckling analysis of moderately thick rectangular plates using the variational differential quadrature method. Aerosp Sci Technol 2019; 91: 479–493.

31.

Hassani

Ansari

Rouhi

A VDQ-based multifield approach to the 2D compressible nonlinear elasticity. Int J Numer Methods Eng 2019; 118(6): 345–370.

32.

Ansari

Torabi

Norouzzadeh

An integral nonlocal model for the free vibration analysis of Mindlin nanoplates using the VDQ method. Europ Phys J Plus 2020; 135(2): 1–16.

33.

Kirk

DE.

Optimal control theory: an introduction. Chelmsford, MA: Courier Corporation, 2004.

34.

Radhoush

Samavat

Vali

MA.

Optimal control of linear time-varying systems using the Chebyshev wavelets (a comparative approach). Syst Sci Control Eng 2014; 2(1): 691–698.

35.

Maher

Mohamed

MJ.

An enhanced genetic programming algorithm for optimal controller design. Intell Control Automat 2013; 4(1): 1–8.