Solving distributed-order fractional optimal control problems via the Fibonacci wavelet method

Abstract

A new approach to finding the approximate solution of distributed-order fractional optimal control problems (D-O FOCPs) is proposed. This method is based on Fibonacci wavelets (FWs). We present a new Riemann–Liouville operational matrix for FWs using the hypergeometric function. Using this, an operational matrix of the distributed-order fractional derivative is presented. Implementing the mentioned operational matrix with the help of the Gauss–Legendre numerical integration, the problem converts to a system of algebraic equations. Error analysis is proposed. Finally, the validation of the present technique is checked by solving some numerical examples.

Introduction

There exist many real-life phenomena that can be modeled by different types of optimal control problems such as harmonic oscillator (Banks and Burns, 1978), chemotherapy for breast cancer (Newbury, 2007; Sabermahani and Ordokhani, 2019), electrical networks (Bahaa, 2017), and cancer treatment (Hassani et al., 2021). Also, some applications of OCPs appeared in Bohannan (2008); Sua’rez et al. (2008).

The existence of solutions to various classes of these problems is studied in Liu et al. (2021, 2020, 2022b, 2019). Moreover, in some cases of the mentioned problems, often in applications, it is not possible to find a solution for OCPs in analytical form. Then, various techniques are developed to solve them. The author in Liu et al. (2022a) investigated fractional optimal feedback control problems with history-dependent operators in separable reflexive Banach spaces. A sensitivity analysis for fractional optimal control problems is presented in Liu et al. (2018). Also, some other numerical techniques developed to find solutions for these problems are the fractional Boubaker wavelet method (Rabiei and Razzaghi, 2022), the Bernoulli wavelet method (Barikbin and Keshavarz, 2020), a method based on hybrid biorthogonal cubic Hermite spline multiwavelets (Mohammadzadeh and Lakestani, 2018), the radial basis functions and collocation method (Chen et al., 2020), Dickson polynomials (Chen et al., 2021), and so on.

The fractional derivative is a special case of the disturbance-order fractional derivative. In recent years, the D-O fractional derivative has been used to model some phenomena, which has led to a clearer model of the behavior of the system or the phenomenon in question and a more accurate analysis of the behavior of a phenomenon. The D-O FOCP is one of the important types of OCPs. Recently, the D-O FOCP has been considered by researchers, and studies have been conducted on these kinds of problems.

The authors in Gariboldi and Tarzia (2007) investigated the convergence of D-O optimal controls using the penalization method. Foderaro and Ferrari (2010) presented the necessary conditions for the optimality of a class of D-O FOCPs. A Legendre pseudo-spectral method is proposed to solve distributed optimal control of the viscous Burgers equations (Sabeh et al., 2016). Bezier curves are used to solve a class of D-O OCPs (Darehmiraki et al., 2016). Zaky proposed a Legendre collocation method to solve D-O FOCPs (Zaky, 2018).

In order to make a more comprehensive investigation of the studies done on the subject of this manuscript, we present a brief bibliometric analysis. For this approach, we consider the following main keyword for data collection in Scopus and Core collocation of Web of Science (WoS) (Sabermahani and Ordokhani, 2021b):

((TITLE(“optimal control problem*”) OR TITLE(“OCP*”)) AND (TITLE(“distributed”) OR TITLE(“distributed-order”) OR TITLE(“distributed order”)))

The title search for the search term on April 21, 2022, led to the retrieval of 142 documents from Scopus and 107 documents from WoS.

Figure 1 illustrates the publication trends in Scopus and WoS in the period 1964–2022.

Figure 1.

The number of publications in the research area from 1964 to 2022.

In WoS, H-index is “17” and the average per item is “13.05” which is an acceptable value for Mathematics. Due to the application of the mentioned problem, Figure 2 displays the subject area of documents in Scopus.

Figure 2.

The subject area of documents in Scopus.

Moreover, Figure 3 demonstrates a bibliometric map through VOSviewer. In this figure, the map is plotted as a network during the span of time.

Figure 3.

The relationships between achieved network data Scopus in considered span.

Due to the above discussion, we consider the D-O FOCP in the following form

$\min \mathfrak{J}(\mathsf{x},\mathsf{u})={\int }_0^1\mathcal{F}(t,\mathsf{x}(t),\mathsf{u}(t))\mathrm{d}t,$

(1)

with the dynamic system

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)={\mathcal{G}}_1\mathsf{x}(t)+{\mathcal{G}}_2\mathsf{u}(t),\mathrm{\varrho }\in (\mathrm{0,1}),$

(2)

and $\mathsf{x}(0)={\mathsf{x}}_0$ , where ${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)$ denotes the distributed-order (D-O) fractional derivative of order ξ(ϱ) and is defined by Pourbabaee and Saadatmandi (2019)

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)={\int }_0^1\xi (\mathrm{\varrho }){}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\mathsf{x}(t)\mathrm{d}\mathrm{\varrho },$

(3)

${\mathcal{G}}_1,{\mathcal{G}}_2$ are constant, and

${\int }_0^1\xi (\mathrm{\varrho })\mathrm{d}\mathrm{\varrho }<\infty .$

Remark 1

A generalization of the Caputo fractional derivative (see Sabermahani et al., 2020b) is the D-O fractional derivative. For ξ(ϱ) = δ(ϱ−ν) which is the Dirac delta function, ${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}$ becomes the Caputo fractional derivative of order ν ${(}_0^C{\mathfrak{D}}_t^{\nu })$ (Zaky, 2018).

In recent years, the important and efficient effects of wavelets have been proven in solving various problems with smooth and non-smooth solutions. In this study, the Fibonacci wavelet has been used as the main tool in providing an efficient computational method for solving D-O FOCPs. Fibonacci polynomials have some advantages which have been mentioned in different studies. For instance, the advantages Fibonacci polynomials have over shifted Legendre polynomials are listed in Sabermahani et al. (2020a):

• Fibonacci polynomials have fewer terms than the shifted Legendre polynomials, which reduces CPU time using Fibonacci polynomials.

• The coefficient of individual terms in Fibonacci polynomials is smaller than corresponding ones in shifted Legendre polynomials, which reduces CPU time using Fibonacci polynomials.

• The integration operational matrix of Fibonacci polynomials has less error than the same operational matrix for shifted Legendre polynomials.

These properties are also inherited from Fibonacci wavelets. Then, considering the feature of Fibonacci wavelets, we present a new computational method for solving an important class of fractional optimal control problems.

Here, a new fractional integration operational matrix based on the hypergeometric function is proposed for these wavelets. The properties of the operational matrix have reflected well in the process of the numerical scheme and created an approximate solution with high precision. Next, with the help of this operational matrix and the Gauss–Legendre integration formula, a new D-O fractional derivative operational matrix is presented for Fibonacci wavelets. Then, by applying these achievements, we solve the considered problem. On the other hand, most of the proposed numerical methods based on spectral methods take control variables from the existing dynamic system in the problem and place them in the performance index. Then, the numerical method is proposed to solve the resulting problem by approximating the state variable. However, in this paper, the control variable and the largest order derivative of the state variable are approximated using Fibonacci wavelets, and the proposed method is presented by using the mentioned approximations in the performance index and dynamic system, with the help of Lagrangian multipliers.

FW and approximation function

Fibonacci polynomials φ _ζ (t) are defined by the following formula (Falcon and Plaza, 2009)

${\phi }_{\zeta }(t)={\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\zeta }{2}\rfloor }(\begin{array}{c}\zeta -\iota \\ \iota \end{array})t^{\zeta -2\iota }.$

(4)

Theorem 1

If $\widehat{y}(t)$ is defined on [0, 1] where $\left|{\widehat{y}}^{(i)}(0)\right|\le {\mathfrak{M}}^i,i\ge 0$ and $\widehat{y}(t)={\displaystyle \sum }_{j=1}^{\infty }{c}_j{\phi }_j(t)$ , then we have (Abd-Elhameed and Youssri, 2016)

• $\left|c_j\right|{\mathfrak{M}}^{j-1}\cosh (\mathfrak{M})/2^{j-1}(j-1)!$ ;

• The series converges absolutely.

where $\mathfrak{M}$ is a positive constant.

Theorem 2

If $\widehat{y}(t)$ satisfies the hypothesis of Theorem 1, and ${\mathcal{E}}_i(t)={\displaystyle \sum }_{j=i+1}^{\infty }{c}_j{\phi }_j(t)$ , then we have the following error estimate (Abd-Elhameed and Youssri, 2016)

$\left|{\mathcal{E}}_i(t)\right|<\frac{\cosh (\mathfrak{M}){(\mathfrak{M}\sigma )}^{i+1}{e}^{\frac{\mathfrak{M}\sigma }{2}}}{2^{i+1}(i+1)!}.$

Also, the truncation error of the generalized Fibonacci polynomials is discussed in Abd-Elhameed and Youssri (2019).

For t ∈ [0, 1), $\mathsf{n}=\mathrm{1,2},\dots ,2^{\mathsf{k}-1}$ , and $\mathsf{m}=\mathrm{0,1},\dots ,M-1$ , FWs $({\psi }_{\mathsf{n},\mathsf{m}}(t))$ are defined as (Sabermahani et al., 2020a)

${\psi }_{\mathsf{n},\mathsf{m}}(t)=\{\begin{array}{cc}\frac{2^{\frac{\mathsf{k}-1}{2}}}{\sqrt{{\omega }_{\mathsf{m}}}}{\phi }_{\mathfrak{m}}(2^{\mathsf{k}-1}t-\mathsf{n}+1),& t\in [\frac{\mathsf{n}-1}{2^{\mathsf{k}-1}},\frac{\mathsf{n}}{2^{\mathsf{k}-1}}),\\ 0,& otherwise\end{array}$

(5)

subject to ${\omega }_{\mathrm{m}}={\int }_0^1{\phi }_{\mathfrak{m}}^2(t)\mathrm{d}t$ .

An arbitrary function $\mathsf{y}(t)$ , defined over the interval [0,1), can be expanded in FWs as

$\mathsf{y}(t)\simeq {\displaystyle \sum }_{\mathrm{n}=1}^{2^{\mathsf{k}-1}}{\displaystyle \sum }_{\mathrm{m}=0}^{M-1}{y}_{\mathrm{n},\mathrm{m}}{\psi }_{\mathrm{n},\mathrm{m}}(t)={\mathbf{y}}^T{\mathbf{\Psi }}_{\mathsf{k}}^M(t),$

(6)

where

$\mathbf{y}\triangleq {\left[y_{\mathrm{1,0}},\dots ,y_{1,M-1},\dots ,y_{2^{\mathsf{k}-1},0},\dots ,y_{2^{\mathsf{k}-1},M-1}\right]}^T,$

$\begin{array}{ccc}{\mathbf{\Psi }}_{\mathsf{k}}^M(t)& \triangleq & [{\psi }_{\mathrm{1,0}}(t),\dots ,{\psi }_{1,M-1}(t),\dots ,{\psi }_{2^{\mathsf{k}-1},0}(t),\\ & & {\cdots ,{\psi }_{2^{\mathsf{k}-1},M-1}(t)]}^T.\end{array}$

(7)

FW operational matrix of the Riemann–Liouville integral

Here, we consider the Riemann–Liouville integral operator of order ν > 0 for a continuous function, which is defined as follows (Sabermahani et al., 2020b)

${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }\mathrm{y}(t)=\{\begin{array}{cc}\frac{1}{\mathrm{\Gamma }(\nu )}{\int }_0^t{(t-s)}^{\nu -1}\mathrm{y}(s)\mathrm{d}s=\frac{t^{\nu -1}}{\mathrm{\Gamma }(\nu )}*\mathrm{y}(t),& \nu >0,\\ \mathrm{y}(t),& \nu =0,\end{array}$

(8)

and satisfies the following useful property

${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{t}^{\mathrm{m}}=\frac{\mathrm{\Gamma }(\mathsf{m}+1)}{\mathrm{\Gamma }(\mathsf{m}+\nu +1)}{t}^{\mathrm{m}+\nu },\hspace{1em}\mathrm{m}>-1.$

(9)

Then, a new operational matrix of the Riemann–Liouville fractional integration for FWs is proposed. Also, in this section, a new D-O fractional derivative operational matrix for FWs is presented. Due to the aforesaid equation and its property, and the hypergeometric function, we achieve the Riemann–Liouville operational matrix for FWs, where

${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{\mathbf{\Psi }}_{\mathsf{k}}^M(t)\simeq \mathfrak{P}(\nu ){\mathbf{\Psi }}_{\mathsf{k}}^M(t).$

For this aim, the following lemma is recalled.

Lemma 1

If ${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }$ is the Riemann–Liouville integral operator, the following relation is established (Rabiei and Razzaghi, 2021; Sabermahani et al., 2022)

${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }(t^{\mathrm{m}}{\chi }_{\gamma }(t))=\{\begin{array}{cc}\frac{\mathrm{\Gamma }(\mathrm{m}+1)t^{\mathrm{m}+\nu }}{\mathrm{\Gamma }(\mathrm{m}+\nu +1)}-\frac{t^{\nu -1}{\gamma }^{\mathrm{m}+1}}{\mathrm{\Gamma }(\nu )(\mathrm{m}+1)}& \\ \times {}_2{\mathfrak{F}}_1(1-\nu ,\mathrm{m}+1,\mathrm{m}+2;\frac{\gamma }{t}),& t\ge \gamma ,\\ 0,& t<\gamma ,\end{array}$

(10)

where χ _c is the unit step function, ${}_2{\mathfrak{F}}_1$ is the hypergeometric function, and

${}_2{\mathfrak{F}}_1(\alpha ,\beta ,\gamma ;\eta )={\displaystyle \sum }_{i=0}^{\infty }\frac{{(\alpha )}_i{(\beta )}_i}{{(\gamma )}_i}\cdot \frac{{\eta }^i}{i!},$

(11)

where $\left|\eta \right|<1$ and

${(p)}_i=\{\begin{array}{cc}1,& \text{if} i=0,\\ p(p+1)\cdots (p+i-1),& \text{if} i>0.\end{array}$

Theorem 3

For ν > 0, we have

${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{\psi }_{\mathrm{m},\mathrm{n}}(t)=\{\begin{array}{cc}0,& t<\frac{\mathrm{n}-1}{2^{\mathsf{k}-1}},\\ {\mathrm{\Lambda }}_1(t),& \frac{\mathrm{n}-1}{2^{\mathrm{k}-1}}\le t<\frac{\mathrm{n}}{2^{\mathrm{k}-1}},\\ {\mathrm{\Lambda }}_1(t)-{\mathrm{\Lambda }}_2(t),& t\ge \frac{\mathrm{n}}{2^{\mathrm{k}-1}},\end{array}$

(12)

subject to

${\mathrm{\Lambda }}_1(t)=\frac{2^{\frac{\mathsf{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathsf{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array}){(2^{\mathrm{k}-1})}^{\kappa }{(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }$

$\times [t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}{(\frac{n-1}{2^{k-1}})}^{(i+1)}$

${}_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n-1}{2^{k-1}}))],$

and

${\mathrm{\Lambda }}_2(t)=\frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathsf{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array}){(2^{\mathrm{k}-1})}^{\kappa }{(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }$

$\times [t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}{(\frac{n}{2^{k-1}})}^{(i+1)}$

${}_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n}{2^{k-1}}))].$

Proof

Due to equations (4) and (5), FWs can be rewritten as

$\begin{array}{ccc}{\psi }_{\mathsf{n},\mathsf{m}}(t)& =& \frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\phi }_{\mathrm{m}}(2^{\mathrm{k}-1}t-\mathrm{n}+1)({\xi }_{\frac{\mathrm{n}-1}{2^{\mathrm{k}-1}}}(t)-{\xi }_{\frac{\mathrm{n}}{2^{\mathrm{k}-1}}}(t))\end{array}$

$\begin{array}{ccc}& =& \frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array}){(2^{\mathrm{k}-1}t-\mathrm{n}+1)}^{\mathrm{m}-2\iota }\end{array}$

(13)

$\begin{array}{ccc}& \times & ({\xi }_{\frac{\mathrm{n}-1}{2^{\mathrm{k}-1}}}(t)-{\xi }_{\frac{\mathrm{n}}{2^{\mathrm{k}-1}}}(t))\end{array}$

$\begin{array}{ccc}& =& \frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array}){(2^{\mathrm{k}-1})}^{\kappa }{t}^{\kappa }\end{array}$

$\begin{array}{ccc}& \times & {(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }({\xi }_{\frac{\mathrm{n}-1}{2^{\mathrm{k}-1}}}(t)-{\xi }_{\frac{\mathrm{n}}{2^{\mathrm{k}-1}}}(t))\end{array}.$

Due to the linearity of the operator ${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }$ , we obtain

$\begin{array}{ccc}{}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{\psi }_{\mathrm{n},\mathrm{m}}(t)& =& \frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array})\end{array}$

(14)

$\begin{array}{ccc}& & \times {(2^{\mathrm{k}-1})}^{\kappa }{(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }\end{array}\begin{array}{ccc}& & \times ({}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }(t^{\kappa }{\xi }_{\frac{\mathsf{n}-1}{2^{\mathsf{k}-1}}}(t))\end{array}$

$\begin{array}{ccc}& & {-}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }(t^{\kappa }{\xi }_{\frac{\mathrm{n}}{2^{\mathrm{k}-1}}}(t)))\end{array}.$

Then, using Lemma 1 and equation (14), we get

$\begin{array}{ccc}{}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{\psi }_{\mathrm{n},\mathrm{m}}(t)& =& \frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array})\end{array}$

(15)

$\begin{array}{ccc}& & \times {(2^{\mathrm{k}-1})}^{\kappa }{(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }\end{array}$

$\begin{array}{ccc}& & ([t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}\end{array}$

$\begin{array}{ccc}& & \times {(\frac{n-1}{2^{k-1}})}^{(i+1)}\end{array}$

$\begin{array}{ccc}& & \times {}_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n-1}{2^{k-1}}))]\end{array}$

$\begin{array}{ccc}& -& [t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}\end{array}$

$\begin{array}{ccc}& & \times {(\frac{n}{2^{k-1}})}^{(i+1)}\end{array}$

$\begin{array}{ccc}& & {\times }_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n}{2^{k-1}}))])\end{array}$

$\begin{array}{ccc}& \triangleq & {\theta }_{\mathsf{n},\mathsf{m}}(t)\end{array},$

by setting

$\begin{array}{cc}{\mathrm{\Lambda }}_1(t)& =\frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array}){(2^{\mathrm{k}-1})}^{\kappa }\\ & \times {(1-\mathsf{n})}^{\mathrm{m}-2\iota -\kappa }\\ & \times [t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}{(\frac{n-1}{2^{k-1}})}^{(i+1)}\\ & \times {}_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n-1}{2^{k-1}}))],\end{array}$

and

$\begin{array}{cc}{\mathrm{\Lambda }}_2(t)& =\frac{2^{\frac{\mathrm{k}-1}{2}}}{\sqrt{{\omega }_{\mathrm{m}}}}{\displaystyle \sum }_{\iota =0}^{\lfloor \frac{\mathrm{m}}{2}\rfloor }{\displaystyle \sum }_{\kappa =0}^{\mathrm{m}-2\iota }(\begin{array}{c}\mathrm{m}-\iota \\ \iota \end{array})(\begin{array}{c}\mathrm{m}-2\iota \\ \kappa \end{array}){(2^{\mathrm{k}-1})}^{\kappa }\\ & \times {(1-\mathrm{n})}^{\mathrm{m}-2\iota -\kappa }\\ & \times [t^{i+\gamma }\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+\gamma +1)}-\frac{t^{\gamma -1}}{\mathrm{\Gamma }(\gamma )(i+1)}{(\frac{n}{2^{k-1}})}^{(i+1)}\\ & \times {}_2{\mathfrak{F}}_1(1-\gamma ,i+1,i+2;t(\frac{n}{2^{k-1}}))].\end{array}$

Then, the proof is complete.

Now, to present the FW operational matrix of the Riemann–Liouville integral, we approximate ${}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{\nu }{\psi }_{\mathsf{n},\mathsf{m}}(t)$ derived in equation (15) with terms of FWs in the following form

${\theta }_{\mathrm{n},\mathrm{m}}(t)\simeq {\displaystyle \sum }_{i=1}^{2^{\mathrm{k}-1}}{\displaystyle \sum }_{j=0}^{\mathrm{M}-1}{\stackrel{\sim }{\theta }}_{i,j}^{\mathrm{n},\mathrm{m}}{\psi }_{i,j}(t)={\mathbf{\Theta }}_{\mathrm{n},\mathrm{m}}^T{\mathbf{\Psi }}_{\mathsf{k}}^M(t),$

(16)

where $\mathsf{n}=\mathrm{1,2},\dots ,2^{\mathsf{k}-1},\mathsf{m}=\mathrm{0,1},\dots ,M-1$ , and ${\mathbf{\Theta }}_{\mathsf{n},\mathsf{m}}$ are the rows of the FW operational matrix of the Riemann–Liouville integral. For instant, the mentioned matrix for k = 2, M = 2 and ν = 0.2 is

$\mathfrak{P}(0.2)=(\begin{array}{cccc}0.5746& 0.2488& 0.2224& -0.0376\\ -0.0778& 0.8081& 0.1735& -0.0363\\ 0.& 0.& 0.7330& 0.1180\\ 0.& 0.& -0.0957& 0.8582\end{array})$

Description of the present method

According to Remark 1, we divide our numerical method into the following two cases and explain our method for each case.

* Problem (I): In the problem given in equations (1) and (2), we consider ξ(ϱ) = δ(ϱ−ν).

* Problem (II): In the problem given in equations (1) and (2), we consider ${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}$ is not fractional-order derivative.

Problem (I)

By considering problem given in equations (1) and (2), we suppose that

${}_0^C{\mathfrak{D}}_t^{\nu }\mathsf{x}(t)\simeq {\displaystyle \sum }_{\mathrm{n}=1}^{2^{\mathrm{k}-1}}{\displaystyle \sum }_{\mathrm{m}=0}^{\mathrm{M}-1}{x}_{\mathrm{n},\mathrm{m}}{\psi }_{\mathrm{n},\mathrm{m}}(t)={\mathbf{X}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}),$

(17)

and

$\mathsf{u}(t)\simeq {\displaystyle \sum }_{\mathrm{n}=1}^{2^{\mathrm{k}-1}}{\displaystyle \sum }_{\mathrm{m}=0}^{\mathrm{M}-1}{u}_{\mathrm{n},\mathrm{m}}{\psi }_{\mathrm{n},\mathrm{m}}(t)={\mathbf{U}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}).$

(18)

According to the Riemann–Liouville operational matrix presented in the previous section, we obtain

$\mathsf{x}(t)\simeq {\mathbf{X}}^T\mathfrak{P}(\nu ){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})+{\mathbf{E}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}),$

(19)

where

$\mathsf{x}(0)\simeq {\mathbf{E}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}).$

Substituting equations (17)–(19) into equation (2), we achieve

${\mathbf{X}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})={\mathcal{G}}_1({\mathbf{X}}^T\mathfrak{P}(\nu ){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})+{\mathbf{E}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}))+{\mathcal{G}}_2({\mathbf{U}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})),$

(20)

then, with the help of the above relation, we get

$\mathbf{Y}={\mathbf{X}}^T-{\mathcal{G}}_1({\mathbf{X}}^T\mathfrak{P}(\nu )+{\mathbf{E}}^T)-{\mathcal{G}}_2({\mathbf{U}}^T).$

(21)

Next, by substituting the approximation in equations (17)–(19) in the performance index $\mathfrak{J}$ and utilizing the above equation, we have

${\mathfrak{J}}^{*}=\mathfrak{J}+{\mathbf{S}}^T\mathbf{Y},$

(22)

where $\mathbf{S}={\left[{\mu }_{\mathrm{1,0}},\dots ,{\mu }_{2^{\mathsf{k}-1},M-1}\right]}^T$ is the unknown multiplier coefficient vector. To derive extremum of ${\mathfrak{J}}^{*}$ , the necessary conditions are that the following relations hold

$\frac{\partial {\mathfrak{J}}^{*}}{\partial \mathbf{X}}=0,\hspace{1em}\frac{\partial {\mathfrak{J}}^{*}}{\partial \mathbf{U}}=0,\hspace{1em}\frac{\partial {\mathfrak{J}}^{*}}{\partial \mathbf{S}}=0.$

(23)

Then, the relations in equation (23) can be solved by applying the “Find Root” package in “Mathematica software.” By computing the coefficient vectors $\mathbf{X}$ and $\mathbf{U}$ , the values of the performance index and state and control variables are obtained.

Problem (II)

In this case, by considering the problem given in equations (1) and (2), we suppose that

${\mathsf{x}}^{\prime }(t)\simeq {\mathbf{X}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}),$

(24)

and due to the Riemann–Liouville operational matrix for ν = 1, we have

$\mathsf{x}(t)\simeq {\mathbf{X}}^T\mathfrak{P}(1){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})+{\mathbf{E}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}).$

(25)

Theorem 4

Considering ${\mathbf{\Psi }}_{\mathsf{k}}^M(t)$ given in equation (7) and equations (24) and (25), we have

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)\simeq {\mathbf{X}}^T\mathfrak{D}{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}),$

(26)

where

$\mathfrak{D}=\left(\frac{1}{2}{\displaystyle \sum }_{j=1}^{\mathrm{s}}{\mathsf{w}}_j\left[\xi \left(\frac{{\mathsf{t}}_j+1}{2}\right)\mathfrak{P}(1-\frac{{\mathsf{t}}_j+1}{2})\right]\right).$

Proof

Using the Riemann–Liouville operational matrix, we have

${}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\mathsf{x}(t)\simeq {\mathbf{X}}^T\mathfrak{P}(1-\mathrm{\varrho }){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}),$

also, considering the definition of the D-O fractional derivative given in equation (3) and the above relation, we get

$\begin{array}{ccc}{}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)& =& {\int }_0^1\xi (\mathrm{\varrho }){}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\mathsf{x}(t)\mathrm{d}\mathrm{\varrho }\\ & \simeq & {\mathbf{X}}^T({\int }_0^1\xi (\mathrm{\varrho })\mathfrak{P}(1-\mathrm{\varrho })\mathrm{d}\mathrm{\varrho }){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}).\end{array}$

(27)

Utilizing the Gauss–Legendre quadrature scheme for evaluating the integral in equation (27), we achieve

${\int }_0^1\xi (\mathrm{\varrho })\mathfrak{P}(1-\mathrm{\varrho })\mathrm{d}\mathrm{\varrho }\simeq \frac{1}{2}{\displaystyle \sum }_{j=1}^{\mathrm{s}}{\mathrm{w}}_j\left[\xi \left(\frac{{\mathsf{t}}_j+1}{2}\right)\mathfrak{P}\left(1-\frac{{\mathsf{t}}_j+1}{2}\right)\right],$

(28)

in which ${\mathsf{w}}_j$ and ${\mathsf{t}}_j$ are weights and nodes of Gauss–Legendre. Substituting the above relation into equation (27), we get

$\begin{array}{ccc}{}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathsf{x}(t)& =& {\int }_0^1\xi (\mathrm{\varrho }){}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\mathsf{x}(t)\mathrm{d}\mathrm{\varrho }\\ & \simeq & {\mathbf{X}}^T({\int }_0^1\xi (\mathrm{\varrho })\mathfrak{P}(1-\mathrm{\varrho })\mathrm{d}\mathrm{\varrho }){\mathbf{\Psi }}_{\mathsf{k}}^M(t)\\ & \simeq & {\mathbf{X}}^T(\frac{1}{2}{\displaystyle \sum }_{j=1}^{\mathsf{s}}{\mathsf{w}}_j[\xi \left(\frac{{\mathsf{t}}_j+1}{2}\right)\\ & & \times \mathfrak{P}\left(1-\frac{{\mathsf{t}}_j+1}{2}\right)]){\mathbf{\Psi }}_{\mathsf{k}}^M(t)\\ & =& {\mathbf{X}}^T\mathfrak{D}{\mathbf{\Psi }}_{\mathsf{k}}^M(t),\end{array}$

(29)

where

$\mathfrak{D}=\left(\frac{1}{2}{\displaystyle \sum }_{j=1}^{\mathrm{s}}{\mathrm{w}}_j\left[\xi \left(\frac{{\mathsf{t}}_j+1}{2}\right)\mathfrak{P}\left(1-\frac{{\mathsf{t}}_j+1}{2}\right)\right]\right).$

Now, by substituting equations (18), (24)–(26) into equation (2), we derive

$y(t):={\mathbf{X}}^T\mathfrak{D}{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})-{\mathcal{G}}_1({\mathbf{X}}^T\mathfrak{P}(1){\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})\hspace{1em}$

(30)

$-{\mathbf{E}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t}))-{\mathcal{G}}_2({\mathbf{U}}^T{\mathbf{\Psi }}_{\mathrm{k}}^{\mathrm{M}}(\mathrm{t})),$

then, we can write

$\stackrel{\sim }{\mathbf{Y}}={\mathbf{X}}^T\mathfrak{D}-{\mathcal{G}}_1({\mathbf{X}}^T\mathfrak{P}(1)+{\mathbf{E}}^T)-{\mathcal{G}}_2({\mathbf{U}}^T).$

(31)

Next, by substituting the approximation in equations (18) and (25) in the performance index $\mathfrak{J}$ and utilizing the above equation, we have

${\mathfrak{J}}^{\star }=\mathfrak{J}+{\mathbf{S}}^T\stackrel{\sim }{\mathbf{Y}},$

(32)

where S is the unknown multiplier coefficient vector. Similar to Problem (I), the necessary conditions to extremum ${\mathfrak{J}}^{\star }$ are

$\frac{\partial {\mathfrak{J}}^{\star }}{\partial \mathbf{X}}=0,\hspace{1em}\frac{\partial {\mathfrak{J}}^{\star }}{\partial \mathbf{U}}=0,\hspace{1em}\frac{\partial {\mathfrak{J}}^{\star }}{\partial \mathbf{S}}=0.$

(33)

The “Find Root” package is used to solve the above system of algebraic equations.

Error analysis

In this section, we propose an error bound of the D-O fractional derivative defined in equation (3).

Lemma 2

Assume y ∈ H ^ϑ (0, 1), ϑ ≥ 0 and $\stackrel{\sim }{y}$ be the best approximation of y, then (Mashayekhi et al., 2016)

$\Vert y-\stackrel{\sim }{y}{\Vert }_{L^2(\mathrm{0,1})}\le \mathsf{c}{(M-1)}^{-\vartheta }{(2^{\mathsf{k}-1})}^{-\vartheta }\Vert y^{(\vartheta )}{\Vert }_{L^2(\mathrm{0,1})},$

(34)

and for v ≥ 1, we derive

$\Vert y-\stackrel{\sim }{y}{\Vert }_{H^v(\mathrm{0,1})}\le \mathsf{c}{(M-1)}^{2v-\frac{1}{2}-\vartheta }{(2^{\mathsf{k}-1})}^{v-\vartheta }\Vert y^{(\vartheta )}{\Vert }_{L^2(\mathrm{0,1})}.$

(35)

Lemma 3

Assume y ∈ H ^ϑ (0, 1), ϑ ≥ 0 and 0 ≤ ϱ < 1, so

$\begin{array}{ccc}{\Vert }_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}y{-}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\stackrel{\sim }{y}{\Vert }_{L^2(\mathrm{0,1})}& \le & \frac{1}{\mathrm{\Gamma }(2-\mathrm{\varrho })}\mathrm{c}{(M-1)}^{2v-\frac{1}{2}-\vartheta }\\ & \times & {(2^{\mathsf{k}-1})}^{v-\vartheta }{\Vert y^{(\vartheta )}\Vert }_{L^2(\mathrm{0,1})},\end{array}$

(36)

in which 1 ≤ v < ϑ.

Proof

We know the following relation holds for the Riemann–Liouville fractional integration and the Caputo fractional derivative

${}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}y(t){=}_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{n-\mathrm{\varrho }}{}_0^C{\mathfrak{D}}_t^{n}y(t).$

Now, using equations (8), (34), and (35), the above relation for n = 1, and (Richard, 2013)

${{\Vert y*z\Vert }_p\le \Vert y{\Vert }_1\Vert z\Vert }_p,$

we derive

$\begin{array}{ccc}& & {\Vert }_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}y(t)-{}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\stackrel{\sim }{y}(t){\Vert }_{L^2(\mathrm{0,1})}^2\\ & =& {\Vert }_0^{\mathcal{R}\mathcal{L}}{\mathfrak{I}}_t^{1-\mathrm{\varrho }}({}_0^C{\mathfrak{D}}_{t}y(t)-{}_0^C{\mathfrak{D}}_t\stackrel{\sim }{y}(t)){\Vert }_{L^2(\mathrm{0,1})}^2\\ & =& \Vert \frac{1}{t^{\mathrm{\varrho }}\mathrm{\Gamma }(1-\mathrm{\varrho })}*({}_0^C{\mathfrak{D}}_{t}y(t)-{}_0^C{\mathfrak{D}}_t\stackrel{\sim }{y}(t)){\Vert }_{L^2(\mathrm{0,1})}^2\\ & \le & {(\frac{1}{(1-\mathrm{\varrho })\mathrm{\Gamma }(1-\mathrm{\varrho })})}^2\Vert {}_0^C{\mathfrak{D}}_{t}y(t)-{}_0^C{\mathfrak{D}}_t\stackrel{\sim }{y}(t){\Vert }_{L^2(\mathrm{0,1})}^2\\ & \le & {(\frac{1}{\mathrm{\Gamma }(2-\mathrm{\varrho })})}^2\Vert y(t)-\stackrel{\sim }{y}(t){\Vert }_{H^v(\mathrm{0,1})}^2\\ & \le & {(\frac{1}{\mathrm{\Gamma }(2-\mathrm{\varrho })})}^2{\mathsf{c}}^2{(M-1)}^{4v-1-2\vartheta }{(2^{\mathsf{k}-1})}^{2v-2\vartheta }\Vert y^{(\vartheta )}{\Vert }_{L^2(\mathrm{0,1})}^2.\end{array}$

Taking the square root of the above equation, the proof is complete.

Theorem 5

Assume y ∈ H ^ϑ (0, 1), ϑ ≥ 0, M ≥ v, then

${\Vert }_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}y{-}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\stackrel{\sim }{y}{\Vert }_{L^2(\mathrm{0,1})}\le \mathfrak{C}\mathsf{c}{(M-1)}^{2v-\frac{1}{2}-\vartheta }{(2^{\mathsf{k}-1})}^{v-\vartheta }$

(37)

$\times {\Vert y^{(\vartheta )}\Vert }_{L^2(\mathrm{0,1})},$

where

${\Vert \xi (\mathrm{\varrho })\Vert }_{L^2(\mathrm{0,1})}\le \mathfrak{C}.$

Proof

Suppose that

$\begin{array}{ccc}& & {\Vert }_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}y{-}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\stackrel{\sim }{y}{\Vert }_{L^2(\mathrm{0,1})}\\ & =& {\Vert {\int }_0^1\xi (\mathrm{\varrho })({\Vert }_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}y(t)-{}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\stackrel{\sim }{y}(t))\mathrm{d}\mathrm{\varrho }\Vert }_{L^2(\mathrm{0,1})}\\ & \le & {\int }_0^1\Vert \xi (\mathrm{\varrho }){\Vert }_{L^2(\mathrm{0,1})}{\Vert }_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}y(t)-{}_0^C{\mathfrak{D}}_t^{\mathrm{\varrho }}\stackrel{\sim }{y}(t){\Vert }_{L^2(\mathrm{0,1})}\mathrm{d}\mathrm{\varrho }\\ & \le & {\int }_0^1\mathfrak{C}(\frac{1}{\mathrm{\Gamma }(2-\mathrm{\varrho })}\mathsf{c}{(M-1)}^{2v-\frac{1}{2}-\vartheta }{(2^{\mathsf{k}-1})}^{v-\vartheta }{\Vert y^{(\vartheta )}\Vert }_{L^2(\mathrm{0,1})})\mathrm{d}\mathrm{\varrho }\\ & \le & \mathfrak{C}\mathsf{c}{(M-1)}^{2v-\frac{1}{2}-\vartheta }{(2^{\mathsf{k}-1})}^{v-\vartheta }{\Vert y^{(\vartheta )}\Vert }_{L^2(\mathrm{0,1})}.\end{array}$

Algorithm implementation and numerical experiments

In this section, we propose the numerical results achieved by the present scheme. To confirm the theoretical prediction, three numerical experiments are carried out by considering the following items. The computations associated with the examples were performed using Mathematica 12.3 on a 2.67 GHz Corei5 personal computer with 4 GB of RAM.

Here, we consider some cases for the distribution functions as follows:

▶ Case (1): ξ(ϱ) = δ(ϱ−1),

▶ Case (2): ξ(ϱ) = δ(ϱ−0.99),

▶ Case (3): ξ(ϱ) = δ(ϱ−0.9),

▶ Case (4): ξ(ϱ) = δ(ϱ−0.8),

▶ Case (5): $\xi (\mathrm{\varrho })={\mathcal{N}}^{(\mathrm{0.8,0.05})}(\mathrm{\varrho })$ , and

▶ Case (6): $\xi (\mathrm{\varrho })={\mathcal{N}}^{(\mathrm{0.5,0.1})}(\mathrm{\varrho })$ ,

in which ${\mathcal{N}}^{(\mu ,\sigma )}(\mathrm{\varrho })=1/\sigma \sqrt{2\pi }\left(e^{-{(\mathrm{\varrho }-\mu )}^2/2{\sigma }^2}\right)$ represents a normal (Gaussian) distribution with mean μ and standard deviation σ.

Example 1

As the first example, we consider the following optimal control system with a distributed-order function described by

$\mathfrak{J}(\mathrm{x},\mathsf{u})=\frac{1}{2}{\int }_0^1\left[{\mathrm{x}}^2(t)+{\mathrm{u}}^2(t)\right]\mathrm{d}t,$

subject to

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathrm{x}(t)=\mathrm{u}(t)-\mathrm{x}(t),\hspace{1em}\mathrm{x}(0)=1.$

The optimal values of $\mathsf{x}(t)$ and $\mathsf{u}(t)$ for Case (1) are

$\mathsf{x}(t)=\theta \sinh (\sqrt{2}t)+\cosh (\sqrt{2}t),$

$\mathsf{u}(t)=(\theta +\sqrt{2})\sinh (\sqrt{2}t)+(\sqrt{2}\theta +1)\cosh (\sqrt{2}t),$

$\theta =-\frac{\cosh (\sqrt{2})+\sqrt{2}\sinh (\sqrt{2})}{\sqrt{2}\cosh (\sqrt{2})+\sinh (\sqrt{2})}\simeq -0.98,$

and

$\mathfrak{J}=0.1929092980931.$

Also, the other cases do not have optimal solutions. To approximate the value of $\mathfrak{J}$ for each of the cases, we apply the proposed method. The problem in all cases has been solved by different numerical schemes such as the spectral tau method (Youssri and Abd-Elhameed, 2018), the Legendre collocation method (Zaky, 2018), Bernstein wavelets (Rahimkhani and Ordokhani, 2021), the variational iteration method (Sahu and Saha Ray, 2018), the Adomian decomposition method (Zeid and Yousefi, 2016), the Chebyshev wavelet method (Sahu and Saha Ray, 2018), the Laguerre wavelet method (Sahu and Saha Ray, 2018), the Legendre wavelet method (Sahu and Saha Ray, 2018), and the cosine and sine wavelet method (Sahu and Saha Ray, 2018).

The numerical results of the performance index $\mathfrak{J}$ obtained by the present method and the aforementioned methods are reported in Table 1 and for all cases mentioned in this section. The graphs of absolute errors of state variable x(t) and control variable u(t) are plotted with k = 2, M = 7, for Case (1) in Figure 4. Moreover, the comparison of the exact and approximate solution of state variable x(t) and control variable u(t) with k = 1, M = 4 for different considered cases are shown in Figure 5. The approximate solutions for this example in Case (1) by the spectral tau method for N = M = 7 is plotted in Youssri and Abd-Elhameed (2018). From the figures in Youssri and Abd-Elhameed (2018) and Figure 5, it is clear that both the state and the control variables are in good agreement with the exact solutions.

Table 1.

The estimated values of $\mathfrak{J}$ in Example 1.

Method	Case (1)	Case (2)	Case (3)	Case (4)	Case (5)	Case (6)
Variational iteration method	0.192909	0.19153	0.17953	0.16711	—	—
Adomian decomposition method	0.192909	0.19155	0.17962	0.16740	—	—
Chebyshev wavelet method (k = 1, M = 4)	0.192921	0.19155	0.179641	0.167328	—	—
Laguerre wavelet method (k = 1, M = 4)	0.192604	0.19169	0.185942	0.181798	—	—
Legendre wavelet method (k = 1, M = 4)	0.192909	0.191531	0.179529	0.167078	—	—
Cosine and sine wavelet method (k = 1, M = 4)	0.196914	0.192572	0.180555	0.165906	—	—
Bernstein wavelet method (k = 2, M = 5)	0.183102	0.181852	0.169226	0.153519	0.144717	0.118702
Legendre collocation method	0.192909	—	0.179690	0.167347	0.165906	—
Present method (k = 1, M = 4)	0.19290927155	0.191530523	0.179524569	0.167057138	0.156200587	0.135489281
Present method (k = 1, M = 6)	0.19290929799	0.191530630	0.179528341	0.167075103	0.156208392	0.135360601
Present method (k = 2, M = 4)	0.19290929799	0.191530630	0.179528260	0.167075048	0.156216249	0.104684079
Present method (k = 2, M = 6)	0.19290929809	0.191530631	0.179530541	0.167075422	0.152703040	0.104641661
Exact value	0.19290929809	—	—	—	—	—

Figure 4.

The absolute errors of state variable x(t) and control variable u(t) with k = 2, M = 7, for Case (1), in Example 1.

Figure 5.

Comparison of the exact and approximate solution of state variable x(t) and control variable u(t) with k = 1, M = 4, for different considered cases in Example 1.

Example 2

We consider the following C-O FOCP

$\mathfrak{J}(\mathrm{x},\mathrm{u})={\int }_0^1\left[0.625{\mathrm{x}}^2(t)+0.5{\mathsf{u}}^2(t)+0.5\mathsf{u}(t)\mathrm{x}(t)\right]\mathrm{d}t,$

subject to

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathrm{x}(t)=\mathrm{u}(t)+0.5\mathrm{x}(t),\hspace{1em}\mathrm{x}(0)=1.$

The optimal values of $\mathrm{x}(t)$ , $\mathrm{u}(t)$ , and $\mathfrak{J}$ for Case (1) are

$\mathrm{x}(t)=\frac{\cosh (1-t)}{\cosh (1)},$

$\mathrm{u}(t)=-\frac{(\tanh (1-t)+0.5)\cosh (1-t)}{\cosh (1)},$

$\mathfrak{J}=\frac{e^2\sinh (2)}{{(1+e^2)}^2}\simeq 0.38079707797788315.$

This problem is solved by using some numerical methods such as the Chebyshev finite difference method (El-Kady, 2003), the Boubaker polynomials (Rabiei et al., 2017), the Bernoulli wavelet (Barikbin and Keshavarz, 2020), the Fibonacci wavelet (Sabermahani and Ordokhani, 2021a), the Bernstein wavelet method (Rahimkhani and Ordokhani, 2021), and so on. We compare the functional value $\mathfrak{J}$ obtained by the developed method and the mentioned methods for each case, in Table 2. Also, in Case (1), we report for k = 3, M = 5, the functional value is $\mathfrak{J}=0.38079707797788254$ and the absolute errors of the state variable x(t) and control variable u(t), in this case, are plotted in Figure 6. Figure 7 displays the comparison of the exact and approximate solution of the state variable x(t) and control variable u(t) with k = 1, M = 4, for all considered cases. Also, we consider the following relations

$\mathfrak{R}(t){=}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathrm{x}(t)-\mathrm{u}(t)-0.5\mathrm{x}(t),$

and use the residual norm in the following form

${\Vert \mathfrak{R}\Vert }^2={\int }_0^1{\mathfrak{R}}^2(t)\mathrm{d}t.$

Table 3 shows the residual norm for Cases (5) and (6) for various values of k and M.

Table 2.

The estimated values of $\mathfrak{J}$ in Example 2.

Method	Case (1)	Case (2)	Case (3)	Case (4)	Case (5)	Case (6)
Chebyshev finite difference method (m = 7)	0.380797	—	—	—	—	—
Boubaker polynomials (m = 5)	0.380797	—	—	—	—	—
Bernoulli wavelet (k = 1, M = 4)	0.380797	—	—	—	—	—
Fibonacci wavelet (k = 1, M = 4)	0.380797070	—	—	—	—	—
Fibonacci wavelet (k = 2, M = 4)	0.380797079	—	—	—	—	—
Bernstein wavelet method (k = 2, M = 5)	0.380797079	0.379459	0.367210	0.353386	0.364870	0.310670
Present method (k = 1, M = 4)	0.380797069467	0.37940725	0.36670442	0.35230693	0.35205564	0.30896981
Present method (k = 2, M = 4)	0.380797077976	0.37940728	0.36670514	0.35231049	0.36098546	0.27122067
Exact value	0.380797077978	—	—	—	—	—

Figure 6.

The absolute errors of state variable x(t) and control variable u(t) with k = 3, M = 5, for Case (1), in Example 2.

Figure 7.

Comparison of the exact and approximate solution of state variable x(t) and control variable u(t) with k = 1, M = 4, for different considered cases in Example 2.

Table 3.

The residual norm for Cases (5) and (6) for various values of k and M in Example 2.

	k = 1, M = 4	k = 2, M = 4
Case (5)	1.402 × 10−30	7.894 × 10−29
Case (6)	1.787 × 10−29	6.093 × 10−29

Example 3

We consider the following time-varying C-O FOCP

$\mathfrak{J}(\mathrm{x},\mathrm{u})=\frac{1}{2}{\int }_0^1\left[(\mathrm{x}(t)-t^a{)}^2+(\mathrm{u}(t)-t^a-\mathrm{\Gamma }(a+1){)}^2\right]\mathrm{d}t,$

subject to

${}_0^C{\mathfrak{D}}_t^{\xi (\mathrm{\varrho })}\mathrm{x}(t)=\mathrm{u}(t)-\mathrm{x}(t),\hspace{1em}\mathrm{x}(0)=0.$

The analytical solution of this example for ξ(ϱ) = δ(ϱ−a) is $\mathrm{x}(t)=t^a,\mathrm{u}(t)=t^a+\mathrm{\Gamma }(a+1)$ , with optimal functional value $\mathfrak{J}=0$ .

For Case (1) and a = 1, we implement the proposed scheme with k = 1, M = 2 to solve this problem. By expanding ${\mathsf{x}}^{\prime }(t)$ and $\mathsf{u}(t)$ using FWs, we get

${\mathsf{x}}^{\prime }(t)\simeq {\mathbf{X}}^T{\mathbf{\Psi }}_1^2(t),$

$\mathsf{u}(t)\simeq {\mathbf{U}}^T{\mathbf{\Psi }}_1^2(t),$

then, we obtain

$\mathsf{x}(t)\simeq {\mathbf{X}}^T\mathfrak{P}(1){\mathbf{\Psi }}_1^2(t).$

So, by substituting the above relations into the main problem, we find the minimum value of $\mathfrak{J}$ . To this aim, we derive

$\mathbf{X}={\left[\mathrm{1,0}\right]}^T,\mathbf{U}={\left[1,\frac{1}{\sqrt{3}}\right]}^T,\mathbf{S}={\left[\mathrm{0,0}\right]}^T.$

where S is a multiplier coefficient vector mentioned in equation (22). Therefore, we achieve $\mathsf{x}(t)=t,\mathsf{u}(t)=t+1$ and $\mathfrak{J}=0$ , which are the exact solutions. Furthermore, in this case, we compare the performance index $\mathfrak{J}$ obtained by our method (k = 1, M = 2), the Bernstein wavelets (Rahimkhani and Ordokhani, 2021) (k = 1, M = 2), the Chebyshev wavelet method (Sahu and Saha Ray, 2018) (k = 1, M = 4), the Laguerre wavelet method (Sahu and Saha Ray, 2018) (k = 1, M = 4), the Legendre wavelet method (Sahu and Saha Ray, 2018) (k = 1, M = 4), and the cosine and sine wavelet method (Sahu and Saha Ray, 2018) (k = 1, M = 4), in Table 4. Also, in Table 5, a comparison between the absolute error of $\mathfrak{J}$ is reported, for a = 1, obtained from the proposed method and Bernstein wavelets (Rahimkhani and Ordokhani, 2021) (k = 2, M = 4) for Cases (2)–(4). Table 6 presents a comparison between the absolute error of $\mathfrak{J}$ , for a = 1, obtained from the proposed method, Bernstein wavelets (Rahimkhani and Ordokhani, 2021) (k = 2, M = 4), the generalized fractional (GF) Chebyshev wavelet method (k = 2, M = 4) (Ghanbari and Razzaghi, 2022), and the Müntz–Legendre wavelet method (Razzaghi and Vo, 2022) for Cases (5) and (6). For more investigation, Figure 8 displays the graphs of $\mathsf{x}(t)$ and $\mathsf{u}(t)$ obtained by the present technique with k = 1, M = 4 and for ξ(ϱ) given in Cases (1)–(6).

Table 4.

The estimated values of $\mathfrak{J}$ , for Case (1) and a = 1, in Example 3.

Method	$\mathfrak{J}$
Cosine and sine wavelet method (k = 1, M = 4)	5.60618 × 10−3
Laguerre wavelet method (k = 1, M = 4)	4.64175 × 10−9
Chebyshev wavelet method (k = 1, M = 4)	7.62947 × 10−32
Legendre wavelet method (k = 1, M = 4)	5.86857 × 10−32
Bernstein wavelet method (k = 1, M = 2)	8.57528 × 10−32
Present method (k = 1, M = 2)	0
Exact value	0

Table 5.

The estimated values of $\mathfrak{J}$ , for Cases (2)–(4) and a = 1, in Example 3.

	Bernstein wavelet method (k = 2, M = 4)	Present method (k = 2, M = 4)
Case (2)	2.11834 × 10−8	1.53111 × 10−31
Case (3)	4.48006 × 10−6	1.92901 × 10−31
Case (4)	4.08630 × 10−5	4.19082 × 10−31

Table 6.

The estimated values of $\mathfrak{J}$ , for Cases (5) and (6) and a = 1, in Example 3.

	Bernstein wavelet method (k = 2, M = 4)	GF Chebyshev wavelet method (k = 1, M = 4, α = 1)	Müntz–Legendre wavelet method (k = 1, M = 4, α = 1)	Present method (k = 1, M = 4)
Case (5)	9.693 × 10−4	2.520 × 10−4	4.351 × 10−4	1.697 × 10−4
Case (6)	9.209 × 10−3	2.257 × 10−3	1.863 × 10−3	9.254 × 10−3

Figure 8.

Comparison of the exact and approximate solution of state variable x(t) and control variable u(t) with k = 1, M = 4, for different considered cases in Example 3.

Conclusion

A class of optimal control problems involving the distributed-order fractional derivative has been investigated. A computational scheme has been successfully employed for solving considered problems. The method is based on Fibonacci wavelets. Some useful error bounds have been established. The corresponding Riemann–Liouville operational matrix has been constructed, and using this, the operational matrix of the distributed-order fractional derivative is presented. These matrices together with the Gauss–Legendre numerical integration are used to discretize the considered optimal control problems. In general, the present scheme is a more accurate solution than some existing techniques. To verify the computational efficiency of the developed method, this method has been applied to several numerical experiments. It has been found that selecting an appropriate value for k and M is needed to derive a high accuracy and satisfactory convergence. The simulation results confirm this claim.

Considering the high effectiveness of wavelets in solving different types of problems, in the future, by defining fractional-order Fibonacci wavelets, we can have one more degree of freedom (α) to these functions in addition to K and M. Therefore, with the new tool, fractional Fibonacci wavelet, we can design effective numerical methods to solve different types of fractional problems arising in science and engineering.

Also, stability analysis of the suggested method for the numerical solution of the mentioned problems and applying this scheme for solving various kinds of fractional equations such as fractional partial differential equations, delay fractional optimal control problems, and fractional integro-differential equations are interesting problems for future works.