A pseudo-spectral approach for optimal control problems of variable-order fractional integro-differential equations

1.   Introduction

In recent decades, fractional calculus has been employed as a powerful tool to describe many engineering and physical phenomena ^[1,3,4,6,16]. A new concept that has been introduced recently is a variable-order fractional operator, which is an extension of classical fractional calculus that can vary in terms of time, place, or any other variable. These types of operators are of considerable importance due to their memory property, allowing many scientific problems to be modeled using differential equations based on these operators ^[5]. Therefore, it is important to find the approximate solution to variable-order fractional differential equations (VOFDEs). In the variable-order fractional optimal control problems (VOFOCPs), the VOFDEs are considered as the dynamic system of problems, which can be defined according to different definitions of fractional derivatives, with the the Riemann-Liouville and Caputo derivatives being the most important ones. In ^[7], the fractional-order Bessel wavelets were used by Dehestani et al. to solve optimal control problems under variable-order fractional dynamical systems. They proposed a collocation method and used pseudo-operational matrices of variable-order fractional derivatives, and dual operational matrices to solve the problem. In ^[10], Heydari and Avazzadeh used the Legendre wavelets to solve VOFOCPs. Then, using the operational matrix of the Riemann-Liouville fractional integration and the Legendre wavelet properties, they turned the performance index into a nonlinear algebraic equation and the dynamic system into a system of algebraic equations. Heydari ^[8] solved a class of VOFOCPs by introducing cardinal Chebyshev polynomials and obtaining their operational matrix corresponding to the Atangana-Baleanu-Caputo derivative.

In recent years, numerical techniques have been developed to propose approximate solutions for variable-order fractional integro-differential equations (VOFIDEs). As shown in ^[9], Heydari presented a method for computing nonlinear fractional quadratic integral equation solutions based on Chebyshev cardinal wavelets and a new operational matrix of variable-order fractional integration derived for the mentioned basis functions. Substituting the mentioned expansion into the intended problem results in a system of nonlinear algebraic equations. In ^[12], by applying piecewise integral quadratic spline interpolation to the estimation of fractional integral operators with variable order, new discretization techniques were proposed. In ^[13], using the second kind of Chebyshev polynomials, differential operational matrices, and integral operational matrices were derived. After adopting the collocation points, the original equation can be transformed into an algebraic system by combining two types of operational matrices. However, no studies have been published in the area of solving the optimal control (OC) problems governed by VOFIDEs. In this article, we attempt to find an approximate solution to such problems for the first time. We first solve these problems by using the collocation spectral method and a different definition of Lagrange polynomial, which is based on Chebyshev polynomials. By substituting the power series of the Chebyshev-modified polynomials, we then transform the dynamic system into a system of algebraic equations. By accurately calculating the fractional integral and derivative of the power series at collocation points, we can calculate the derivative matrix with precision. Also, we approximate the integral objective function using the Gauss quadrature rules.

The structure of the paper is as follows. In Section 2, some necessary preliminaries are given. A statement of the problem is presented in Section 3. In Section 4, a Chebyshev pseudo-spectral approach is given to solve the OC problem of VOFIDEs. Also, we present a convergence analysis for the method. In Section 5, some numerical examples are given to show the efficiency of the method. Finally, the conclusions and suggestions are presented in Section 6.

2.   Some preliminaries

2.1.   Fractional calculus

In this section, we present some basic definitions and mathematical preliminaries related to the fixed-order and variable-order fractional derivatives (see ^[2]).

Definition 2.1. Let $\chi(.)$ be defined on the interval $[0, T_{f}]$. The left and right Riemann-Liouville fractional integrals of fixed-order $\mu > 0$ are denoted by ${}_0I_\tau^\mu \chi(\tau)$ and ${}_\tau I_{T_f}^\mu \chi(\tau)$, respectively, and defined by

Definition 2.2. Consider $\chi(.)$ a function defined on interval $[0, T_{f}]$. The left and right Riemann-Liouville fractional derivatives of fixed-order $0 < \mu < 1$ are denoted by ${}_0D_\tau^\mu \chi(\tau)$ and ${}_\tau D_{T_f}^\mu \chi(\tau)$, respectively, and defined by

Definition 2.3. Let us suppose that function $\chi(.)$ is defined on the finite interval $[0, T_{f}]$. The left and right Caputo fractional derivatives of $\chi(.)$ of fixed-order $0 < \mu < 1$ are denoted by ${}_0^CD_\tau^\mu \chi(\tau)$ and ${}_\tau^CD_{T_f}^\mu \chi(\tau)$, respectively, and defined by

We now present the basic concepts of variable-order fractional calculus and take the fractional order in the derivative and integral as a continuous function on $(0, T_f)$. First, we introduce the generalization of a fixed-order fractional integral called the variable-order Riemann-Liouville integral.

Definition 2.4. Assuming that the continuously differentiable function $\chi$ is defined on $(0, T_f)$. The left and right Riemann-Liouville fractional integrals of order ${\mu (\tau)}$ are defined as follows:

and

Definition 2.5. Consider $\chi:[0, T_{f}] \to \mathbb{R}$ as a continuously differentiable function and suppose $\mu :[0, T_{f}] \to [0, 1]$ is a given function.

(1) The type I left and right Caputo variable-order fractional derivatives (VOFDs) of $\chi(\tau)$ of order $\mu (.)$, respectively, are defined by

(2) The type II left and right Caputo VOFDs of $\chi(\tau)$ of order $\mu (.)$, respectively, are given by

(3) The type III left and right Caputo VOFDs of $\chi(\tau)$ of order $\mu (\tau)$, respectively, are defined by

Lemma 2.1. Let $\chi(\tau) = (\tau-a)^{\beta }$ for $\tau\in[a, b]$ where $\beta > 0$. Then,

In this paper, we will focus on the type Ⅲ Caputo VOFDs.

2.2.   The shifted Chebyshev polynomials

We consider a special case of the Jacobi polynomials-Chebyshev polynomials (of the first kind), ${\tilde T_n}(t)$, which are proportional to Jacobi polynomials $J_n^{-\frac{1}{2}, -\frac{1}{2}}$ and are orthogonal with respect to the weight function $\omega(t) = \left(1-t^2\right)^{-\frac{1}{2}}$.

The three-term recurrence relation for the Chebyshev polynomial reads

For practical use of Chebyshev polynomials on the interval of interest $[0, T_{f}]$, it is necessary to change the defining domain by means of the following substitution:

So, the shifted Chebyshev polynomials $ST_{n}(\tau)$ on $[0, T_{f}]$ are obtained as follows:

The orthogonality condition for these shifted polynomials is

Shifted Chebyshev polynomials can be analytically written as follows:

3.   The statement of problem

In this paper, we consider the following OC problem of VOFIDE:

where $\chi_{0} \in \mathbb{R}$ is given, $f:\mathbb{R} \times {\mathbb{R}^n} \times {\mathbb{R}^m} \to \mathbb{R}$, and $g:\mathbb{R} \times {\mathbb{R}^n} \times {\mathbb{R}^m} \to \mathbb{R}$ are continuous functions, $\chi(\tau)$ and $\upsilon(\tau)$ are the state and control variables, respectively, $\mu_{1}, \mu_{2}:[0, T_{f}] \to[0, 1]$ are two given continuous functions, ${}_0^C \mathbb{D}_\tau^{\mu_{1}(\tau)}$ is the type Ⅲ Caputo VOFD operator, and ${}_0{I}_\tau^{\mu_{2}(\tau)}$ is the left Riemann-Liouville fractional integral.

4.   Numerical treatment of the OC problem of VOFIDE

In this section, we try to get an approximate solution to the optimal control problem (3.1)-(3.2). In subsection 4.1, we discuss the implementation of the method and present a new method for calculating the derivative and integral matrices. Also, we present a convergence analysis for the suggested method.

4.1.   The Chebyshev pseudo-spectral (CPS) method

For interpolating in the CPS method, the following Lagrange polynomials are utilized:

where

$\tau_{j} = {\frac{T_f}{2}}(t_{j} + 1)$ and ${t_j} = \cos ({ \dfrac{{\pi j}}{N}}), \; \; j = 0, 1, 2, ..., N,$ are roots of $(1 - {t }^2){\tilde{T}_N}(t)$.

The Lagrange polynomials have the useful property of the delta Kronecker, i.e.,

Now, we approximate the variables of problem (3.1)-(3.2) in terms of the Lagrange functions as follows:

Also, we have

Moreover, we approximate the control variable as

where $\bar{\chi} = ({{{\bar{\chi}}}_{1}}, {{{\bar{\chi}}}_{2}}, ..., {{{\bar{\chi}}}_{N}})$ and $\bar{\upsilon} = ({{{\bar{\upsilon}}}_{1}}, {{{\bar{\upsilon}}}_{2}}, ..., {{{\bar{\upsilon}}}_{N}})$ are unknown coefficients. By applying the interpolation property of the Lagrange polynomial, we get

Now, by relations (4.1) and (4.3), we have

Also, with the help of the Chebyshev polynomial analytical form shown in relation (2.18) in the previous section, we will have

By applying Lemma 2.1 on the above relation, we reach the following relation, which is the result of the fractional derivative and integral effect on the power function:

Lemma 4.1. Suppose $q:[-1, 1]\to \mathbb{R}$ is a continuous function. The following integral approximation is referred to as the Chebyshev-Gauss-Lobatto (CGL) quadrature rule:

where ${{\tau}_{j} = \cos (\tfrac{N-j}{N}\pi)}, \; j = 0, 1, ..., {N}$ are roots of $\left(1-{{\tau}^{2}} \right)\frac{d}{d\tau}{\tilde {T}_{N}}(\tau)$ and ${\tilde{T}_{N}}(\tau) = \cos (N{{\cos }^{-1}}(\tau))$ is the Chebyshev polynomial of order $N$ and ${{\bar \omega }_{j}} = \sqrt{1-\tau_{j}^2}\dfrac{\pi }{{{{\tilde{c}}}_{j}}N}, \; j = 0, 1, ..., N$ are the quadrature weights of the numerical approximation (4.10), where ${{{\tilde{c}}}_{0}} = {{{\tilde{c}}}_{N}} = 2, {{{\tilde{c}}}_{j}} = 1$ for $j = 1, ..., N-1$.

Now, using relations (4.6), (4.8) and (4.9) and the above lemma, we reach the following discrete system, which can be solved by optimization methods:

where

Here, $\bar \chi = (\bar \chi_{0}, \bar \chi_{1}, ..., \bar \chi_{N})$ and $\bar \upsilon = (\bar \upsilon_{0}, \bar \upsilon_{1}, ..., \bar \upsilon_{N})$ are unknown coefficients and $\mu_{j} = 2$ for $j = 0, N$ and $\mu_{j} = 1$ for $\; j = 1, 2, ..., N-1.$ Also, $\omega_j = \frac{T_f}{2}\tilde{\omega_j}$. Note that there is a new technique to calculate operation matrices of variable-order fractional derivatives and integrals. First, we rewrite the Lagrange polynomial as follows:

where ${{\eta }_{pj}}$ are unknown coefficients that are determined as below. Using the delta Kronecker property (4.2), we get

In matrix form, the above relation can be expressed as follows:

where

According to (4.2) and the invertibility of matrix $A$

finally,

Therefore, the components of the operational matrix of derivative can be obtained as follows:

and

According to (4.16) and (4.17), we get

and

4.2.   Convergence of the method

Here, we show the convergence of the method by applying an assumption.

We assume that the OC problem (3.1)-(3.2) has a Lagrange interpolating polynomial based on the shifted Chebyshev-Gauss-Lobatto (SCGL) points, which uniformly converges to it.

Theorem 4.1. Assume that $\{\big(\bar \chi _j^*, \bar v_j^*)\big\}_{j = 0}^{N}$ is an optimal solution of (4.11)-(4.12) and define $\bar{X}_{{N}}({\tau}) = \sum\limits_{j = 0}^{N} {{\bar \chi _j^*}{h_j}({\tau})}$ and $\bar U_{{N}}({\tau}) = \sum\limits_{j = 0}^{N} {{\bar \upsilon _j^*}{h_j}({\tau})}$ on $[0, T_{f}]$. Also, assume $\{\big(\bar X_{{N}}(.), \bar U_{{N}}(.)\big)\}_{{N} = {{N}_0}}^\infty$ uniformly converges to $\big(\bar \chi(.), \bar \upsilon(.)\big)$ such that $\bar \chi(.)$ and $\bar \upsilon(.)$ are continuously differentiable and $_{ 0}^C \mathbb{D}_{\tau} ^ {\mu_{1} ({\tau})}\bar \chi(.)$ and $_{ 0} \mathbb{I}_{\tau} ^ {\mu_{2} ({\tau})}\bar \chi(.)$ are in $C((0, T_{f}])$. Then, $\big(\bar X(.), \bar U(.) \big)$ is an optimal solution for the main OC problem of VOFIDE (3.1)-(3.2).

Proof. We first show that $\big(\bar X(.), \bar U(.)\big)$ is a feasible solution for the problem (3.1)-(3.2). Suppose that ${\tau}\in (0, T_f]$ is given. Since shifted CGL points $\left\{ {\tau}_k \right\}_{k = 0}^{N}$ with ${N}\to \infty$ is dense on $[0, { T_f}]$, there exists a subsequence $\left\{ {\tau}_{k_j} \right\}_{j = 0}^\infty$ such that $\lim_{j\to\infty}{k_j} = \infty$ and $\lim_{j\to\infty} { {\tau}_{k_j}} = {\tau}$. By continuity of functions $g(., ., .)$, $_{ 0}^C \mathbb{D}_{\tau} ^ {\mu_{1} ({\tau})}\bar \chi(.)$ and $_{ 0} \mathbb{I}_{\tau} ^ {\mu_{2} ({\tau})}\bar \chi(.)$, we get

Also, for the initial condition,

Now, we want to show that $\big(\bar X(.), \bar {U}(.) \big)$ is an optimal solution for the problem (3.1)-(3.2). By objective function (4.11), we define

Also, by objective functional (3.1) and replacing continuous function $q(.)$ in relation (4.10) with $f(., \bar X(.), \bar {U}(.))$, we gain

Moreover, since $\sum\limits_{k = 0}^{N} {\omega _k} = {T_f}$ and $\big(\bar X_{{N}}(.), \bar {U}_{{N}}(.) \big)$ is uniformly convergent to $\big(\bar X(.), \bar {U}(.) \big)$, we get

where $L_{1} > 0$ is the Lipschitz constant of continuously differentiable function $f(., ., .)$. Thus, by (4.20)–(4.22), we gain

Hence,

On the other hand, for any optimal solution $\big(X^{*}(.), {U}^{*}(.)\big)$ of problem (3.1)-(3.2), there exists a corresponding sequence $\{\big(\chi_j^{*}, \upsilon_j^{*}\big)\}_{i = 0}^\infty$ such that

Since $\big(X^{*}(.), {U}^{*}(.)\big)$ satisfies the constraint (3.2), sequence $\{\big(\chi_j^{*}, \upsilon_j^{*}\big)\}_{j = 0}^{N}$ with ${N} \to \infty$ satisfies constraint (4.12). Similar to the relation (4.23) and the process of achieving it, we can conclude

where $\chi^{*} = (\chi_0^{*}, \chi_1^{*}, ..., \chi_{N}^{*})$ and $\upsilon^{*} = (\upsilon_0^{*}, \upsilon_1^{*}, ..., \upsilon_{N}^{*})$. By relations (4.23) and (4.24), and optimality of pairs $\big(\bar \chi^*, \bar \nu^*\big)$ and $\big(X^{*}(.), {U}^{*}(.)\big)$, we achieve

which tends to $L(X^{*}(.), {U}^{*}(.)) = L(\bar X(.), \bar {U}(.))$. Thus, $(\bar X(.), \bar {U}(.))$ is an optimal solution for the OC problem (3.1)-(3.2). □

5.   Numerical examples

In this section, some examples are shown to depict the efficiency and practicability of the devised approximation method. MATLAB has been used for all examples. The absolute errors are computed as

and

where pairs $(\chi, \upsilon)$ and $(\chi_{N}, {\upsilon}_{N})$ are the exact and approximate solutions, respectively. We should also point out that the CPU time for program running in solving the discussed problems in all examples, requires less than 3 seconds for N = 5.

Example 5.1. Consider the OC problem

under the variable-order fractional dynamical system

and the initial condition is given by $\chi(0) = 0$, so that the optimal solutions are

Now, we solve the above problem for the following orders:

In Figure 1, we present the results obtained using the proposed technique as well as the exact solution for $\upsilon(\tau)$ and $\chi(\tau)$ when $N = 5$. Using the presented technique, Figure 2 shows the absolute errors for $\upsilon(\tau)$ and $\chi(\tau)$ at various $\mu_{1}^{i}(\tau)$ and $\mu_{2}^{i}(\tau), \; i = 1, 2, 3, 4$, and $N = 5$. The approximate values of the performance index $L$ with different $\mu_{1}^{i}(\tau)$ and $\mu_{2}^{i}(\tau), \; i = 1, 2, 3, 4$, and $N = 5$ are reported in Table 1. As the figures show, numerical solutions agree well with exact solutions when compared with the numerical results.

Example 5.2. Consider the OC problem

subject to the dynamical system

The optimal value for the performance index is ${L}^{*} = 0$, which is obtained by

This problem is solved by the proposed method with $N = 5$ for the following variable orders:

Figure 3 illustrates the behavior of the numerical solutions with $N = 5$ for state and control variables $\chi(\tau), \upsilon(\tau)$ for the above-mentioned $\mu _{1}^{1}(\tau), \mu _{2}^{1}(\tau)$. The absolute errors obtained by the suggested method in the state variable and control variable are reported in Figure 4. Table 2 contains the performance index with $N = 5$ for the variant $(\mu _{1}^{i}(\tau), \mu _{2}^{i}(\tau)), \; i = 1, 2, 3, 4$. In all these, it is evident that the proposed method yields numerical solutions that are highly accurate for all cases of orders of derivatives and integrals.

Example 5.3. Another OC problem of VOFIDE is as follows:

under the constraint

The solutions $\chi(\tau) = \tau^2$ and $\upsilon(\tau) = 2\tau+sin(\tau)$ minimize the performance index $L$ and $L^{*} = 0$. Using the method described in this paper, we are able to solve this problem for different variable orders $(\mu _{1}^{i}(\tau), \mu _{2}^{i}(\tau)), \; i = 1, 2, 3, 4$. The fractional orders used are as follows:

Considering Table 3 and Figures 5 and 6, it is evident that the proposed method achieves numerical solution with high accuracy.

Example 5.4. In this example, we consider an optimal control problem of the energy for a fractional RLC series electrical circuit. We solve this problem by using our method. This problem can be formulated as the following fractional form (where $R = 1(\Omega), L = 1(H)$ and $C = 1(F)$):

where state $\chi$ and control $\upsilon$ are the current and voltage in the RLC circuit, respectively. Note that if ${\mu _1}(t) = {\mu _2}(t) = 1$, then Eq (5.6) can be written as the following equivalent form:

which is the Kirchhoff's voltage law. Note that the integro-differential equation (5.8) has been discussed and analyzed by many researchers, for example, see relation (8.3) in ^[15]. Also, in relation (29) in ^[15], the fractional form of (5.8), i.e., integro-differential equation (5.6), was introduced and studied. Some other fractional form of a RLC electrical circuit can be seen in relation (1.3) in ^[14] and relation (9) in ^[11]. In fact, fractional derivatives and integrals play an important role in the modeling of electrical circuits that contain super resistants, super capacitors, and super inductors ^[17]. Moreover, fractional models provide a more efficient description and representation of real electrical systems. However, the goal of solving the minimum energy problem (5.5)-(5.7), is to move an electrical initial current $\chi(0) = {\chi_0}$ by a voltage in minimum energy to a desired final state $\chi(T) = {\chi_T}$. We assume ${\chi_0} = 1$ (ampere), ${\chi_T} = 0.25$ (ampere) and $T = 1$ (second) and solve problem (5.5)-(5.7) for different ${\mu _1}(t)$ and ${\mu _2}(t)$. The obtained optimal solutions are shown in Figures 7 and 8 for $N = 5$. Also, the obtained minimum energy L for different cases is given in Tables 4 and 5.

It can be seen that when ${\mu _1}(t)$ and ${\mu _2}(t)$ tend to 1, the obtained optimal state (current) and optimal control (voltage) go to the corresponding optimal solution with ${\mu _1}(t) = {\mu _2}(t) = 1$. This issue, regarding the obtained optimal value for the performance index L, can also be seen in the tables, which confirms the correctness of proposed method to solve this practical problem. Further, we can conclude that in fractional RLC model $(0 < \mu _1, \mu _2 < 1)$, more total energy (that is, L) is needed to bring the current in the circuit from an initial state to the desired state in a certain time compared with the real RLC model $(\mu _1 = \mu _2 = 1)$.

6.   Conclusions

In this paper, we investigated an approach in order to solve a nonlinear optimal control problem involving variable-order fractional integro-differential equations as the dynamic system. Pseudo-spectral collocation is the basis of this method. At first, by using the expansion of Lagrange polynomials in terms of Chebyshev polynomials and the power series of them, the problem was converted into an nonlinear programming problem, which was easier to solve. Then, variable-order fractional derivatives in the Caputo sense were represented by a new operational matrix, and fractional integrals were represented by an operational matrix. With the suggested method, the optimal control problem of the variable-order fractional integro-differential equation could easily be solved. Using the numerical results, we could see that the approximate and exact solutions are in good agreement and the method is efficient and accurate as well.

Author contributions

Zahra Pirouzeh, Mohammad Hadi Noori Skandari, Kamele Nassiri Pirbazari and Stanford Shateyi: Conceptualization, Methodology, Writing-review & editing, Software, Validation. All authors of this article have been contributed equally. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflicts of interest.