King algorithm: A novel optimization approach based on variable-order fractional calculus with application in chaotic financial systems

Chaos, Solitons and Fractals 132 (2020) 109569

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

King algorithm: A novel optimization approach based on variable-order fractional calculus with application in chaotic financial systems Samaneh Soradi-Zeid a, Hadi Jahanshahi b,∗, Amin Yousefpour c, Stelios Bekiros d,e a

Faculty of Industry and Mining (Khash), University of Sistan and Baluchestan, Zahedan, Iran Department of Mechanical Engineering, University of Manitoba, Winnipeg R3T 5V6, Canada School of Mechanical Engineering, College of Engineering, University of Tehran, Tehran 11155-4563, Iran d Department of Economics, European University Institute, Viadelle Fontanelle,18, I, Florence, 50014, Italy e Rimini Centre for Economic Analysis (RCEA), Wilfrid Laurier University, 75 University Ave W., Waterloo N2L3C5, Ontario, Canada b c

a r t i c l e

i n f o

Article history: Received 26 November 2019 Accepted 15 December 2019

2010 MSC: 26A33 65M70 42A15 49M25 Keywords: Variable order fractional calculus Variable order fractional system Variable order fractional optimal control problem King algorithm Volterra integral equation Direct trajectory optimization

a b s t r a c t In this study, a new optimization algorithm, called King, is introduced for solving variable order fractional optimal control problems (VO-FOCPs). The variable order fractional derivative is portrayed in the Caputo sense through the dynamics of the system as variable order fractional differential equation (VO-FDE). To this end, firstly, the VO-FOCP is converted into a system of VO-FDEs. Then, according to the fact that the VO-FDE is equivalent to a Volterra integral equation, the system of VO-FDEs is transformed into an equivalent system of variable order fractional integro-differential equations. In the next step, both the context of minimization of total error and a joint application of Banach’s fixed-point theorem are used to solve a nonlinear optimization problem. Actually, using the existing mechanism, the synchronization problem is recast to an optimization problem. Due to the discretization and its board range of arbitrary nodes, the proposed method provides a very adjustable framework for direct trajectory optimization. Finally, the proposed algorithm is verified using several common optimization functions and a chaotic financial system. Also, through simulation results, the proposed method is compared with some popular methods. Simulation results demonstrate the appropriate performance of the introduced method.

1. Introduction According to variability of arbitrary value for the order of fractional derivatives and integrals, we can consider an order that is not constant which generates an extension of the classical fractional calculus, namely variable order fractional calculus (VO-FC). More precisely, this topic is a development of fractional calculus where the orders of fractional operators are estimated by any arbitrary value which depends on the time. This subject applied through many physical applications evidently, when the memory effects are sensitive to the time and place of the event. Several ways exist that define variable order fractional derivatives but the Riemann-Liouville and the Caputo fractional derivatives of variable order have been used more than the other definitions. Many

∗

Corresponding author. E-mail address: [email protected] (H. Jahanshahi).

https://doi.org/10.1016/j.chaos.2019.109569 0960-0779/© 2019 Elsevier Ltd. All rights reserved.

© 2019 Elsevier Ltd. All rights reserved.

research studies on developing applications of VO-FC have been done to prove the diversity of scientific, engineering fields and physical models [22–24]. More specifically, fractional formulations used to exemplified the mechanics and dynamic systems [25,26,41], in this regard, we could consider chaotic financial and economic systems as suitable examples [14,20]. An extension of fractionality concepts in financial mathematics has also been developed in [21]. It must be considered that due to the high complexity of variable order fractional derivatives, the analytically handling equations described by this derivatives is extremely difficult and even impossible [7]. To overcome this challenge, practical numerical/approximate methods have been presenting to solve them. These processes turn into rapid increasing developments of numerical methods for VO-FDEs. Accordingly, some numerical methods for solving FDEs have been generalized to solve VO-FDEs have been studied in [8,26–40] A generalization of classic FOCPs when the dynamics system is described by VO-FDEs is considered as VO-FOCPs. It is worth

2

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

noting that the research on VO-FOCPs is completely new and numerical studies of these problems are still at an early stage of growth. The reason to formulate and solve VO-FOCPs has recently been answered affirmatively turn into significant increasing of these problems in control systems. Although many computational methods have been proposed for solving FOCPs, [6], only a small number of these methods have so far been generalized for solving VO-FOCPs [1–5]. In these references, the authors have tried to use the interpolate approximate basis functions to reformulated the VO-FOCPs as an equivalent system of algebraic equations that makes the problem significantly simpler. The overall aim of the current paper is to introducing a class of VO-FOCPs and a king algorithm based on pontryagin maximum principle (PMP). Variable transformation and fixed-point theorem are used to obtain a new equation with a smoother solution. Thus, we introduce the following VO- FOCP:

min J (u ) =



tf

t0

F (t , x(t ), u(t ))dt ,

(1)

(t ) = G(t , x(t ), u(t )), n − 1 < α (t ) ≤ n,

(2)

N 

1

f (t ) = lim α (t ) h→0 h



(−1 )

j

α (t )



j

j=0

f (t − jh ).

(4)

Definition 2. The Riemann–Liouville variable order fractional integral operator of order α (t) ≥ 0 for a function f(t) defined by: α (t )

t0 It

f (t ) =



1

(α (t ))

t

t0

(t − s )α (t )−1 f (s )ds,

(5)

where (. ) : (0, ∞ ) → R is known as the Euler-Gamma function. Definition 3. The left and right variable order fractional derivatives in Riemann–Liouville implication of f(t) for n − 1 < α (t ) ≤ n are defined by: RL α (t ) t0 Dt

f (t ) =

1 dn (n − α (t )) dt n

RL α (t ) t Dt f

f (t ) =

(−1 )n dn (n − α (t )) dt n



t

t0

 t

(t − s )n−α (t )−1 f (s )ds,

tf

(6)

(s − t )n−α (t )−1 f (s )ds,

(7)

respectively.

and the initial conditions:

x(i ) (t0 ) = xi , i = 0, 1, · · · , n − 1, u(t ) ∈ U, t ∈ [t0 , t f ],

GL α (t ) t0 Dt

and

with the VO-FDE: C α (t ) x t0 Dt

Definition 1. The Grunwald–Letnikov definition of fractional variable order derivative is defined as follows:

(3)

where n ∈ N, xi ∈ R, i = 0, 1, . . . , n − 1, t0 and tf are two positive constant, F and G are continuous functions, the set U ⊂ Rm represents the allowable inputs, which are considered to be continuous functions and Ct Dtα (t ) denotes the Caputo fractional derivative 0 of variable order α (t), which will be introduced in the next section. Here, in the suggested technique, the VO-FOCP (1)–(3) are reformulated into an optimization problem with discrete parameters in a way that the main computational cost of the VO-FOCP (1)– (3) comes from solving the discrete optimization problem. For this purpose, by using the PMP, which consists of VO-FDE (2) in conjunction with the performance index (1) by a set of indefinite Lagrange multipliers, the VO-FOCP in (1)–(3) is recast into a system of VO-FDEs, at first. Then, adjoining the application of Banach’s fixed-point theorem, this system is turned into an equal system of variable-order fractional integro-differential equations. In addition, we use the context of minimization of total error to access a nonlinear optimization problem. To improve the quality of numerical solutions we applied the necessary and sufficient condition for obtaining the optimal extremum to obtain the unknown coefficients. The total scheme of the paper is formed as follows: Section 2 includes the necessary definitions and mathematical preliminaries of the variable order fractional calculus. In Section 3, the proposed method is described. In Section 4, we described our numerical simulation for VO-FOCP (1)–(3). Some illustrative numerical examples which are solved by the proposed method are provided in Section 5. Finally conclusions are drawn in the last section. 2. Basic definitions and notations of variable order fractional calculus The issue of the dynamical problem is considered as a basic solution for variable-order calculus that has been proposed by Lorenzo and Hartley [42,44], Ingman et al. [43] and Coimbra [45] and recently developed in some literature [49,50,52]. Here we express some basic definitions and preliminaries of the variableorder fractional calculus. Moreover, we refer the interested readers for recent advances in variable order derivative and some of their properties to [46–48,51]. In order to better illustrate the notations, we define h = (t f − t0 )/N to be the uniform time step, N is a positive integer and n ∈ N.

Definition 4. The Caputo variable order fractional derivatives of f(t) for n − 1 < α (t ) ≤ n in the left and right concept are defined respectively by:

f (t ) =

1 (n − α (t ))

f (t ) =

(−1 )n (n − α (t ))

C α (t ) t0 Dt

and C α (t ) t Dt f



t

t0

 t

tf

(t − s )n−α (t )−1 f (n) (s )ds,

(8)

(s − t )n−α (t )−1 f (n) (s )ds.

(9)

The above definitions directly result in: C α (t ) t0 Dt

n α (t ) d f

f (t ) =

t0 It

 C α (t ) m t t0 Dt

=

dt n

(t ),

n − 1 < α (t ) < n,

(10)

(m+1 ) m−α (t ) , (m−α (t )+1 ) t

n − 1 < α (t ) ≤ n, n ≤ m ∈ N

0,

o.w. (11)

The linear property for these operation can be represented as:

Dα (t )



 λ f + μg (t ) = λDα (t ) f (t ) + μDα (t ) g(t ),

(12)

where λ and μ are constants. Several methods have been proposed recently to estimate variable order fractional operators. For instance, Moghaddam and Machado [53] have developed a finite difference method for approximating variable order fractional derivatives of arbitrary order. In [54], the authors described one approximation based on the Grunwald–Letnikov and Mittag-Leffler functions. Authors in [35,37] have presented the spline interpolation for their discretization technique to estimate variable order fractional integral operators. 3. Approximation method Consider the general form of VO-FDE as follows: C α (t ) x t0 Dt

(t ) = f (t , x(t )),

(13)

and the initial condition x(t0 ) = x0 , in which x(t ) ∈ X ⊂ is the state variable, t ∈ [t0 , tf ] is independent time variable, x0 is the initial state, 0 < α (t) ≤ 1 and f(.) is a continuous function. Without losing totality, we assume that t0 = 0 and t f = 1. Rn

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

Lemma 5. Let function f(.) is continuous. The problem (13) can be recast as the form of the following Volterra integral equation

x(t ) = x0 +



1

(α (t ))

t

0

f (τ , x(τ ))(t − τ )α (t )−1 dτ .

(14)

It means, every continuous solution of (14) is also a solution of original problem (13) and vice versa. If α (t ) = const, (14) would reduce to a weakly singular Volterra integral equations which can be studied in the general theory of integral equations with algebraic singularity, or in the frame of fractional calculus. Proof. The result follows straightly from the variable order fractional operators definitions. 

3

Definition 7. Based on the VO-FDE which is represented by (13), we defined the total error of the variable order fractional functional as follows:

E (t, x(t )) =



1

0

C0 Dtα (t ) x(t ) − f (t, x(t )) dt,

(19)

where E : X × [0, 1] → R is a continuous function. Theorem 8. For precise modeling of the variable-order fractional nonlinear system behavior C0 Dtα (t ) x(t ) by nonlinear function f(t, x(t)) with initial conditions x(0 ) = x0 , the necessary and sufficient conditions is that the following condition should be acquiescent:

E (t, x(t )) = 0.

(20)

In the next, we will design an effective scheme to stabilize the variable order fractional system (13).

Proof. Let h be in the following form:

Theorem 6. Providing that α (t): [0, 1] → (0, 1) is continuous. Also, presume that there exists a constant K > 0 in which f (t, x ) : [0, 1] × R → R be continuous and fulfills the condition

where C0 Dtα (t ) x(t ) − f (t , x(t )) is a non-negative continuous function

| f (t, x ) − f (t, y )| ≤ K |x − y|, ∀t ∈ [0, 1], x, y ∈ R. If αK∗ < 1, where α ∗ = inf{α (t ); 0 ≤ t ≤ 1}, then the initial value VOFDE (13) has only unique solution defined on x ∈ C[0, 1]. Proof. It was proved that problem (13) is reducible to the integral problem in (14). Using Lemma 5 and following Ref. [9], we introduce the following operator:

(T x )(t ) = x0 +



1

(α (t ))

t 0

f (τ , x(τ ))(t − τ )α (t )−1 dτ ,

(15)

in which T : C ([0, 1], R ) → C ([0, 1], R ). Obviously, T is a continuous operator. Since the initial value VO-FDE (13) is equivalent with the integral Eq. (14), and since this integral equation can be rewritten in the form (T x )(t ) = x(t ), the existence and uniqueness of the fixed point for the operator T prove this Theorem. In this way, we ought to employ Banach contraction principle. Let {xn } be a sequence such that xn → x in C ([0, 1], R ). Then for each t ∈ [0, 1], we have

|T (xn (t )) − T (x(t ))| ≤

≤

1



t



t

(t − τ )α (t )−1 (α (t )) 0 f ( τ , x ( τ )) − f ( τ , x ( τ )) dτ n 1

(t − τ )α (t )−1 sup (α (t )) 0 0≤τ ≤1 f (τ , xn (τ )) − f (τ , x(τ )) dτ .

≤

C Dα (t ) x (t ) 0 t

on [0, 1]. Finally, t ∈ (0, 1) and then

= f (t , x(t )), ∀t ∈ [0, 1], E (t, x(t )) = 0.



Now, we define the following calculus of variations problem for variable order fractional system (13):

E (t, x(t )) = min x

 0

1

C α (t ) 0 Dt x(t ) − f (t, x(t )) dt,

(17)

∗ Kt α x − y

∞

≤



t

0

K

α∗

E (t, x∗ (t )) = 0 ⇒



C α (t ) x 0 Dt

(t )



∗

= f (t, x∗ (t ))

(23)

Based on Theorem 8, in order to solve VO-FDE in the form of (13), we can solve an optimization problem in the form of (22) in which, x∗ (t) is the optimal reactance of the system indicated by (13). Accordingly, based on Lemma 5, Eq. (13) is reducible to the fractional integral Eq. (14). So, as discussed above, in optimal design of a system in the form of Eq. (14), by using the finite horizon rule, the objective is provided to minimizes the following functional: 1 0

 t 1 f (τ , x(τ ))(t − τ )α (t )−1 dτ dt. x(t ) − x0 − (α (t )) 0 (24)

Now, for solving the optimization problem (24), we divide the continuous interval t ∈ [0, 1] into N equal subintervals. Therefore, problem (24) reaches the following form:

min x

 t 1 f (τ , x(τ ))(t − τ )α (t )−1 dτ dt, x(t ) − x0 − (α (t )) ti−1 0

N   i=1

ti

(25) i N,

where t0 = 0, tN = 1 and ti = i = 1, 2, · · · , N − 1. Using the trapezoidal rule in each subinterval, the minimization problem (25) is modified as follows:

(t − τ )inft∈[0,1] α (t )−1 dτ

x − y ∞ .

(22)

in which x(0 ) = x0 . It is assumed that x∗ (t) is the optimal response of the above problem which represent the variable state. Therefore, based on Theorem 8, we have:

t

K x − y ∞ (supt∈[0,1] α (t ))

α∗

it is clear

C Dα (t ) x (t ) − f (t , x (t )) is continuous 0 t

1 that if 0 h(t )dt = 0, h ≡ 0 for each

(16)

1 (t − τ )α (t )−1 f (τ , x(τ )) |T ( x ) − T ( y )| ≤ (α (t )) 0 − f (τ , y(τ )) dτ ≤

1]. Hence,



Now, let x, y ∈ C ([0, 1], R ). By using the fact that sup{α (t ); 0 ≤ t ≤ 1} = 1, for every t ∈ [0, 1] we have:



(21)

and f(t, x(t)) is a continuous function. Therefore, h(t) is a continuous also, the variables C0 Dtα (t ) x(t ) and x are continuous on [0,

J=

Since f is a continuous function, we have:

T (xn ) − T (x ) ∞ → 0 as n → ∞.

h(t ) = C0 Dtα (t ) x(t ) − f (t , x(t )) L1 ,

N  1 1  x(ti−1 ) − x0 − (α (t )) 2N i−1 i=1  ti−1 × f (τ , x(τ ))(ti−1 − τ )α (ti−1 )−1 dτ

min (18)

Therefore, as a sequel of Banach fixed point theorem, we elicit that T has a fixed point that is a solution of problem (13). 



0

+ x(ti ) − x0 −

1

(α (ti ))



ti 0

f (τ , x(τ ))(ti − τ )α (ti )−1 dτ

(26)

4

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

The function f(t, x(t)) is approximate on interval [0, 1] using the piecewise linear function as follows:

f (τ , x(τ )) =

N 

f (tn , x(tn ))χ[ n−1 , n ] (τ ), N

n=1

(27)

N

where χAn (x ) is the characteristic function in which An , n = 1, 2, · · · , N, are a sequence of disjoint sets on [0, 1]. It is quite obvious that when N tends to infinity, the approximation (27) tends to f(τ , x(τ )) uniformly on [0, 1]. Then the minimization problem (26) is altered to the following discrete approximate form:

According to Lemma 9, we choose the following change of variables:

i−1  1 f (tn , x(tn )) pn,i x(ti−1 ) + x(ti ) − 2x0 − (α (t )) i−1 n=1 i  1 − f (tn , x(tn ))qn,i = ui + vi , (α (ti )) n=1 x(ti−1 ) + x(ti ) − 2x0 −

N  i−1  1 1  min x(ti−1 ) − x0 − f (t , x(tn )) 2N (α (ti−1 )) n=1 n i=1  ti−1 × (ti−1 − τ )α (ti−1 )−1 χ[ n−1 , n ] (τ )dτ N

0



+ x(ti ) − x0 −  ×

ti

0

1

i 

(α (ti )) n=1

n=1

−

(28)

N  i−1  1 1  min x(ti−1 ) − x0 − f (tn , x(tn )) pn,i 2N (α (ti−1 )) + x(ti ) − x0 − in which,

pn,i

(α (ti )) n=1

f (tn , x(tn ))qn,i ,

qn,i =

hα (ti )



1,

n=i

α (ti ) (i − n + 1 )α (ti ) − (i − n )α (ti ) , n < i

− (29)

1

i 

(α (ti )) n=1

f (tn , x(tn ))qn,i = ui − vi , (35)

Accordingly, we imply the main theorem of this section as follows.

n
(30)

(31)

⎧ N ⎪ ⎪min i=1 vi ⎪ ⎪ ⎪ ⎨s.t.

(32)

Theorem 10. The truncation error of the Volterra integral Eq. (14) by using the approximation (27) is of order O(hn+1−α (t ) ). Proof. Let x(t ) ∈ C n+1 (), n − 1 < α (t ) ≤ n, ti = ih ∈ , i = 1, 2, · · · , N, h = N1 , be the exact solution of Eq. (13). Then, the approximate solution of using the piecewise linear approximation (27) satisfies

x(t ) − x˜(t ) ∞ = O(h ).

(36)

Thus, we have:

C0 Dtαi (t ) x(t ) − C0 Dtαi (t ) x˜(t ) ∞ ≤ C0 Dtαi (t ) (x(t ) − x˜(t )) ∞ ≤ x(n) (t )  ti 1 − x˜(n ) (t ) ∞ (ti − τ )n−α (ti )−1 dτ (n − α (t )) 0

i

tin−α (ti )

x(n) (t ) − x˜(n) (t ) ∞ (n + 1 − α (ti )) = Cα (ti ) hn−α (ti ) x(n ) (t ) − x˜(n ) (t ) ∞ ≤ O(hn+1−α (ti ) ), =

(37)

in−α (ti ) where Cα (ti ) = . From this theorem, we can easily (n + 1 − α (ti )) conclude that | pn,i | ≤ Cα (ti ) . Similarly, this result holds for qn,i .  3.1. Application in fractional order financial system

u∗i ,

where I is a compact set. Then i = 1, 2, · · · , N, is the optimal solution of the following non-linear programming problem (NLP):

u∈I

i−1  1 f (tn , x(tn )) pn,i (α (ti−1 ))

ui , vi ≥ 0, i = 1, 2, · · · , N.

Lemma 9. Let pairs (v∗i , u∗i ), i = 1, 2, · · · , N, be the optimal solutions of the following LP:

min

i=1

s.t.

n=1

where h = N1 and (29) is the common demonstration of a discrete nonlinear optimization problem for VO-FDE (13). The exact solutions x∗ (t) of problem (29) are not known. Here, there exists an absolute value function (29) and also a function f(x) that was a linear function by approximation (27). So to significantly reduce the calculations, we used the linear approximation. Therefore, the minimization problem (29) turns into a linear programming problem (LP) with the following change of variables.

⎪ vi ≥ ui , vi ≥ −ui , ⎪ ⎪ ⎪ ⎪ ⎩ vi ≥ 0, ui ∈ I,

(34)

N 1  ( ui + vi ) 2N

x(ti−1 ) + x(ti ) − 2x0 −



n=i−1 hα (ti−1 ) 1, = α (ti−1 ) (i − n )α (ti−1 ) − (i − n − 1 )α (ti−1 ) ,

and

min

n=1

i 

f (tn , x(tn ))qn,i = ui − vi ,

where pn,i and qn,i will be obtained from (30) and (31), respectively. Then, the optimization problem (29) is converted to the following LP:

Since every bounded Riemann integrable function defined on a bounded interval is Lebesgue integrable and the two integrals are equal, so by calculating the integrals in problem (28) we have:

1

1

i = 1, 2, · · · , N,

f (tn , x(tn ))

(ti − τ )α (ti )−1 χ[ n−1 n ( τ )d τ . N ,N]



i 

(α (ti )) n=1

N

i=1

i−1  1 f (tn , x(tn )) pn,i (α (ti−1 ))

N  ui .

(33)

i=1

Proof. See [12] for more details.



A straightforward application of system (13) is extended by introducing a chaotic fractional-order financial system. Generally speaking, a government would like to stabilize its financial system in a certain finite time. Since there are some memory effects in the financial systems, the fractional calculus is introduced to the mentioned system. The fractional-order form of a financial system can

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

be written as follows:

of Eqs. (39)–(41) solves the VO-FOCP (1)–(3). See [10,11] for more details. 

⎧ α (t ) C x(t ) = z + (y − a )x + m1 w, ⎪ t0 Dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎨Ct Dtα (t ) y(t ) = 1 − by − x2 + u + m2 w, 0

In a special way, we describe the quadratic VO-FOCP:

(38)

⎪ C ⎪ z(t ) = −x − cz + m3 w, ⎪ t Dt ⎪ ⎪0 ⎪ ⎩C α (t ) u(t ) = −xyz t0 Dt α (t )

Theorem 11. (Necessary optimality conditions) Let (x(t), u(t)) be a minimizer of (1) under dynamic constraint (2) with the initial condition x(0 ) = x0 . Then, there exists a function λ(t) such that for all t ∈ [0, 1], the triplet (x(t), u(t), λ(t)) satisfies the following conditions: the Hamiltonian system

x(t ) = G(t , x(t ), u(t )),

RL α (t ) t D1

λ(t ) = ∂2 H (t , x(t ), u(t ), λ(t )),

(39)

the stationary condition

∂3 H (t, x(t ), u(t ), λ(t )) = 0,

(40)

and the transversally condition

λ ( 0 ) = 0,

x ( 0 ) = x0 ,

(41)

where ∂ j H, j = 1, 2, 3, denote the derivative of H(., ., .) with respect to its jth argument and H is the Hamiltonian, defined by

(42) Proof. To this end, we can convert (1)-(3) into an equivalent augmented performance index:



1 0







xT (t )Q (t )x(t ) + uT (t )R(t )u(t ) dt ,

0

(44)

(t ) = A(t )x(t ) + B(t )u(t ), t ∈ [0, 1],

(45)

and the initial conditions:

x ( 0 ) = x0 ,

(46)

where, B(t) = 0, Q(t) ≥ 0, Q (0 ) = 0 and R(t) > 0. Following Theorem 11, it can be easily shown that the necessary optimality conditions for VO-FOCP (44)–(46) lead to the following equations: C α (t ) x 0 Dt RL α (t ) t D1

(t ) = A(t )x(t ) + B(t )u(t ),

x ( 0 ) = x0 ,

λ(t ) = Q (t )x(t ) + λ(t )A(t ), λ(1 ) = 0,

(47)

u(t ) = −R−1 (t )B(t )λ(t ). We apply the approximation scheme in the last section to this system. So, according to Lemma 5, the first and second equations of (47) transfer to the following equivalent Volterra integral equations:

x(t ) = x0 +



1

(α (t ))

t



0



A(τ )x(τ ) + B(τ )u(τ ) (t − τ )α (t )−1 dτ ,



H (t, x, u, λ ) − λT C0 Dtα (t ) x dt

λ(t ) =



1

(α (t ))

t

1





Q (τ )x(τ ) + λ(τ )A(τ ) (τ − t )α (t )−1 dτ .

Based on Definition 7 and also carrying out the error function in the general form (24), we can obtain the numerical solution of VOFOCP (44)–(46) with placement system (48) by an equivalent optimization problem as follow:

 t   1 x ( t ) − x − A ( τ )x ( τ ) + B ( τ )u ( τ ) 0 (α (t )) 0 0 α (t )−1 (t − τ ) dτ  1  1 + λ(t ) − Q ( τ )x ( τ ) + λ ( τ )A ( τ ) (α (t )) t  α (t )−1 (τ − t ) dτ dt. (49) 

min

1

Accordingly, by approximating the integrals as before and using Lemma 9 due to the influence of the nonlinear term, we construct the optimization problem (49) to a smooth optimization problem. Finally, by acquiring the optimal solutions of this problem, we can recognize the approximate solutions of the primary problem (44)– (46). 5. Results and discussion

H (t, x(t ), u(t ), λ(t )) = F (t, x(t ), u(t )) + λ(t )G(t, x(t ), u(t )).

Jˆ(u, λ ) =

1

(48)

The VO-FOCP consists of finding the optimal control u(t) and the corresponding state variable x(t) which minimize functional system (1) and that applies the conditions (2), (3) with 0 < α (t) ≤ 1 and t ∈ [0, 1]. So, we get through the PMP conditions corresponding to VO-FOCPs based on the Caputo definition for the variable order fractional derivatives.

α (t )



with the following VO-FDE:

4. Numerical simulation of VO-FOCPs

0 Dt

1 2

min J (u ) = C α (t ) x 0 Dt

where x denotes the interest rate, y indicates the investment demand, z denotes the price index, u represents the investment incentive, a ≥ 0 is the saving amount, b ≥ 0 is the cost per investment, c ≥ 0 is the demand elasticity of commercial markets, w denotes the market confidence, m1 , m2 , m3 are the impact factors and 0 < α (t) < 1. When α (t ) = constant, a king algorithm based on Eq. (14) is employed to identify chaos and a robust controller is designed to stabilize the fractional-order chaotic system (38) in finite time [14– 16]. According to Theorems 6 and 8, this scheme can keep the original structure of this system as much as possible and can be applied to stabilizing other chaotic systems including dynamic economic systems. The authors in [18] investigate the stabilization of fractional order chaotic three-dimensional economic system with uncertainties based on adaptive sliding mode control. An adaptive terminal sliding mode control with neural network estimator for finite-time stabilization and synchronization of a four-dimensional fractional-order system has been presented in [17,19].

C

5

(43)

in which λ(t) denotes the adjoint or costate variable for the continuous VO-FDEs constraint (2). By computing the first variation of (43) and setting this variation to zero, we obtain the necessary conditions for the VO-FOCP (1)-(3). In general, the system

In this section, we present some numerical examples and solved them by means of the proposed method. The results are provided for various values of α (t) and the calculations are performed using the MATLAB software. The absolute error is defined as

E (x, x¯ ) = max

1≤ j≤N





x(t j ) − x¯ (t j ) ,

(50)

to implement the applicability and effectiveness of our method, in which x(tj ) and x¯ (t j ), j = 1, 2, · · · , N, are the exact and approximated solutions of our problem, respectively.

6

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

Fig. 1. Approximate values of x(t) at different constant values of α (t) for Example 1.

Table 1 Estimated values of J at different choices of α (t) and N for Example 1.

α (t)

[10]

N = 100

N = 150

N = 300

1 − 0.2t 1 − 0.1t 1

0.484268 0.476446 0.460500

0.484723 0.476351 0.468051

0.484513 0.475881 0.462023

0.464193 0.475564 0.460501

Table 2 Absolute errors of x(t) with different values of α (t) at t = 1 for Example 1.

1 2

 0

1





x2 (t ) + u2 (t ) dt ,

(51)

(t ) = tx(t ) + u(t ),

N = 150

N = 200

N = 400

2.457 × 10−4 1.132 × 10−5 3.769 × 10−6 2.964 × 10−6

2.121 × 10−4 1.038 × 10−5 3.256 × 10−6 3.017 × 10−7

1.170 × 10−4 3.634 × 10−5 2.299 × 10−6 2.103 × 10−6

1.136 × 10−4 3.464 × 10−5 2.131 × 10−6 1.805 × 10−7

0 < α (t ) ≤ 1,

N = 100

t 0.1 0.3 0.5 0.7 0.9

by subject to the following VO-FDE: C α (t ) x 0 Dt

N = 100

1 − 0.2t 1 − 0.1t 1 − 0.01t 1

Table 3 Absolute errors of u(t) for different values of N and t with α (t ) = 1 − 0.15t for Example 1.

Example 1. Minimize the performance index [10]:

J (u ) =

α (t)

(52)

with the initial condition x(0 ) = 1. Figs. 1 and 2 show the behavior of the numerical solutions of x(t) and u(t), respectively, with different constant values of α (t). In Table 1, we list the approximate values of the optimal cost of this problem with those obtained in [10] for different values of N and α (t). In addition, Tables 2 and 3 show the absolute errors of x(t) and u(t), respectively, for various choice of α (t), N and t. From these results, it can be concluded that the proposed method provides approximate solutions, accurately. Moreover, the accuracy of the numerical results has been improved by increasing the number of N.

N = 150 −5

6.435 × 10 2.153 × 10−5 2.340 × 10−6 4.135 × 10−7 3.041 × 10−7

N = 300 −5

4.164 × 10 5.614 × 10−6 7.045 × 10−6 3.896 × 10−7 3.012 × 10−8

1 0.75 + 0.2sin(10t ) 0.75 + 0.2sin(30t ) 0.75 + 0.2sin(50t )

t = 0.2

Example 2. Consider the following VO-FOCP [2]:

min J (u ) =

 0

1



1 2

(x(t ) − t 2 )2 + (u(t ) − te−t + et

5.921 × 10 8.327 × 10−10 1.097 × 10−9 1.946 × 10−10

2

 ) dt,

−t 2

(53)

subject to the dynamical system C α (t ) x 0 Dt

t = 0.4 −12

2.998 × 10−6 4.321 × 10−7 1.145 × 10−8 0.801 × 10−8 0.582 × 10−9

5.431 × 10 6.266 × 10−6 3.776 × 10−7 2.509 × 10−8 1.946 × 10−8

(t ) = ex(t ) + 2et u(t ),

0 < α (t ) ≤ 1,

Table 4 Absolute errors of x(t) for different values of α (t) and t with N = 300 for Example 2.

α (t)

N = 400 −6

t = 0.6 −12

2.655 × 10 7.13 × 10−10 9.743 × 10−9 3.402 × 10−10

t = 0.8 −14

2.632 × 10 2.942 × 10−11 6.655 × 10−10 3.720 × 10−11

t=1 −14

1.849 × 10 6.715 × 10−12 7.137 × 10−10 3.718 × 10−12

2.225 × 10−15 3.747 × 10−14 1.129 × 10−13 2.197 × 10−13

(54)

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569

7

Fig. 2. Approximate values of u(t) at different constant values of α (t) for Example 1.

Table 5 Absolute errors of u(t) for different values of α (t) and t with N = 300 for Example 2.

α (t)

t = 0.2

1 0.75 + 0.2sin(10t ) 0.75 + 0.2sin(30t ) 0.75 + 0.2sin(50t )

t = 0.4 −5

2.310 × 10 1.876 × 10−5 1.559 × 10−8 6.113 × 10−5

t = 0.6 −5

2.645 × 10 2.221 × 10−4 2.128 × 10−6 6.434 × 10−7

t = 0.8 −8

8.215 × 10 2.885 × 10−4 4.226 × 10−4 9.488 × 10−4

t=1 −7

7.054 × 10−8 3.112 × 10−5 2.987 × 10−4 6.658 × 10−7

1.106 × 10 6.409 × 10−5 3.564 × 10−4 3.981 × 10−8

Table 6 Absolute errors of x(t) for different values of α (t) and t with N = 300 for Example 3.

α (t)

t = 0.2

t = 0.4

0.6

t = 0.8

t=1

1 − 0.2t 1 − 0.1t 1

8.260 × 10−4 2.187 × 10−6 3.217 × 10−6

1.261 × 10−4 3.288 × 10−5 8.260 × 10−7

1.306 × 10−4 1.436 × 10−5 4.083 × 10−6

1.722 × 10−5 1.371 × 10−5 1.649 × 10−6

2.165 × 10−5 2.187 × 10−6 5.592 × 10−7

with the initial condition x(0 ) = 0. The exact solution of this prob1 2 lem for α (t ) = 1 is given by x(t ) = t 2 and u(t ) = te−t − et −t . The 2 absolute errors of our method are reported in Tables 4 and 5 for x(t) and u(t), respectively, with different values of α (t) which implies that the proposed method in this study has high accuracy. Example 3. Minimize the performance Eq. (51) subject to the following dynamic: C α (t ) 0 Dt

= −x(t ) + u(t ),

0 < α (t ) ≤ 1,

index

J

given

Table 7 Absolute errors of u(t) with different values of α (t) at t = 0.5 for Example 3.

α (t)

N = 100

N = 150

N = 300

N = 400

1 − 0.2t 1 − 0.1t 1

2.941 × 10−5 3.438 × 10−6 8.858 × 10−7

1.320 × 10−6 2.139 × 10−6 5.878 × 10−7

4.341 × 10−7 4.926 × 10−7 1.678 × 10−8

2.517 × 10−6 2.517 × 10−7 1.713 × 10−9

by

(55)

with the initial condition x(0 ) = 1. The exact solution of this√problems is√obtained for α (t ) = 1 and (t ) = cosh( 2√t ) + √ performed √ by x√ θ sinh( 2t ) and u(t ) = (1 + 2θ )cosh( 2t ) + ( 2 + θ )sinh( 2t ) where θ ≈ −0.98. The absolute errors of these solutions for variable choices of α (t) are reported in Tables 6 and 7. For evaluation the accuracy of the obtained results, we compared our results with the corresponding values [1,10] in Table 8.

Table 8 Estimated values of J at different choices of α (t) for Example 3.

α (t)

[1]

[10]

N = 100

N = 200

N = 400

1 − 0.2t 1 − 0.1t 1

0.16711 0.17953 0.192909

0.167347 0.179690 0.192909

0.167127 0.179548 0.192909

0.167119 0.179540 0.192909

0.167112 0.179534 0.192909

8

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569 Table 9 Absolute errors of x(t) with different values of α (t) for Example 4. N = 300 t 0.1 0.3 0.5 0.7 0.9

N = 400

α 1 (t)

α 2 (t)

3.5718 × 10−11 5.2662 × 10−12 7.6263 × 10−13 1.7608 × 10−11 5.9334 × 10−12

3.4924 × 10−11 3.8532 × 10−12 6.7655 × 10−12 1.3111 × 10−11 3.7747 × 10−12

α 1 (t)

3.7743 × 10−12 5.0576 × 10−13 7.3242 × 10−13 1.0929 × 10−12 5.2662 × 10−13

Heydari [3]

α 2 (t)

2.8871 × 10−12 4.3855 × 10−13 6.5035 × 10−13 1.2079 × 10−12 3.5924 × 10−13

α 1 (t)

3.6056 × 10−13 5.3421 × 10−13 7.6626 × 10−13 1.1523 × 10−12 5.3706 × 10−13

α 2 (t)

2.3954 × 10−13 3.9430 × 10−13 6.7547 × 10−13 1.3163 × 10−12 3.6390 × 10−13

Table 10 Absolute errors of u(t) with different values of α (t) for Example 4. N = 300 t 0.1 0.3 0.5 0.7 0.9

N = 400

α 1 (t)

α 2 (t)

6.3797 × 10−6 3.7091 × 10−6 2.9653 × 10−7 1.9719 × 10−7 8.8668 × 10−7

9.8508 × 10−7 5.0843 × 10−7 3.2174 × 10−7 2.16535 × 10−7 2.1729 × 10−6

α 1 (t)

8.3821 × 10−7 3.1599 × 10−7 2.7350 × 10−7 1.6495 × 10−7 8.2605 × 10−8

Heydari [3]

α 2 (t)

7.6263 × 10−7 5.2577 × 10−7 3.0930 × 10−7 2.0579 × 10−7 1.0610 × 10−7

α 1 (t)

6.8739 × 10−7 3.9493 × 10−7 2.5712 × 10−7 1.7821 × 10−7 9.4306 × 10−8

α 2 (t)

7.9052 × 10−7 4.60 0 0 × 10−7 3.0082 × 10−7 2.0901 × 10−7 1.1078 × 10−7

Table 11 Absolute errors for the optimal cost of J at different choices of α (t) and N for Example 4.

α (t)

Heydari [3]

N = 200

N = 300

N = 400

α 1 (t)

7.0647 × 10−24 6.4024 × 10−24 5.4883 × 10−24 4.3010 × 10−24

9.0062 × 10−23 2.7905 × 10−22 1.7615 × 10−23 2.4869 × 10−22

8.9267 × 10−23 3.4456 × 10−23 5.9122 × 10−24 5.7176 × 10−23

7.0934 × 10−24 8.8585 × 10−24 5.7721 × 10−24 5.0089 × 10−24

0.25 + 0.2sin(2π t ) α 2 (t) 0.25 + 0.5t 5

Example 4. Consider the following VO-FOCP [3]:

min J (u ) =





1 0



(x(t ) − cos(t ))2

+ u(t ) + cos(t ) − t −α (t ) 2,1−α (t ) (−t 2 )

Declaration of Competing Interest

 2

dt

(56)

subject to the following VO-FDE: C α (t ) x 0 Dt

(t ) = x(t ) + u(t ),

0 < α (t ) ≤ 1,

(57)

with the initial condition x(0 ) = 1. Here, ξ ,ζ (z) demonstrates the Mittag-Leffler function [13]. The optimal solution of this problem are presented by x(t ) = cos(t ) and u(t ) = t −α (t ) 2,1−α (t ) (−t 2 ) − cos(t ). The absolute errors for different values of N and t are represented in Tables 9–11. In these tables we consider α1 (t ) = 0.25 + 0.2sin(π t ) and α2 (t ) = 0.25 + 0.2t 5 . We corroborate that these accurate results have been acquired by applying only a small number of N and the accuracy can be enhanced by increasing the number of N. 6. Conclusion A new computational approach based on a king algorithm is constructed for numerically approximating the solutions of a class of VO-FOCPs. Getting through this algorithm, firstly the VOFOCP is transformed into a system of VO-FDEs. Then, according to the fact that the initial value problem of a VO-FDE is equivalent to a Volterra integral equation, this system of VO-FDEs is reformulated into an equivalent system of variable order fractional integro-differential equations. In the last step of this method, we used the context of Banach’s fixed-point theorem with a simultaneous application of minimization of total error in order to reduce this system into a smooth optimization problem. Lastly, through several examples it was confirmed that the results of this method are in perfect agreement with those obtained in the literature.

This statement is to certify that no conflict of interest exits in the submission of this manuscript. Also, the manuscript is approved for publication by all authors. I would like to declare that the work described was an original research that has not been published previously, and not under consideration for publication elsewhere. All the authors listed have approved the manuscript that is enclosed. CRediT authorship contribution statement Samaneh Soradi-Zeid: Conceptualization, Writing - original draft, Project administration. Hadi Jahanshahi: Conceptualization, Writing - original draft, Project administration. Amin Yousefpour: Conceptualization, Writing - original draft, Project administration. Stelios Bekiros: Conceptualization, Writing - original draft, Project administration. References [1] Zahra W K, Hikal M M. Non standard finite difference method for solving variable order fractional optimal control problems. J Vib Control 2017;23(6):948–58. [2] Heydari MH, Avazzadeh Z. A new wavelet method for variable-order fractional optimal control problems. Asian J Control 2018;20(5):1804–17. [3] Heydari MH. A new direct method based on the Chebyshev cardinal functions for variable-order fractional optimal control problems. J Frankl. Inst. 2018;355(12):4970–95. [4] Heydari MH, Avazzadeh Z. A computational method for solving two-dimensional nonlinear variable-order fractional optimal control problems. Asian J Control 2018:1–15. [5] Hassani H, Avazzadeh Z. Transcendental Bernstein series for solving nonlinear variable order fractional optimal control problems. Appl Math Comput 2019;362:124563. [6] Zeid SS, Effati S, Kamyad AV. Approximation methods for solving fractional optimal control problems. Comput Appl Math 2018;37(1):158–82. [7] Zeid SS. Approximation methods for solving fractional equations. Chaos, Solitons Fractals 2019;125:171–93.

S. Soradi-Zeid, H. Jahanshahi and A. Yousefpour et al. / Chaos, Solitons and Fractals 132 (2020) 109569 [8] Solis-Perez JE, Gomez-Aguilar JF, Atangana A. Novel numerical method for solving variable-order fractional differential equations with power, exponential and Mittag-Leffler laws. Chaos, Solitons Fractals 2018;114:175–85. [9] Zeid SS, Kamyad AV, Effati S, Rakhshan SA, Hosseinpour S. Numerical solutions for solving a class of fractional optimal control problems via fixed-point approach. SeMA J 2017;74(4):585–603. [10] Zaky MA. A Legendre collocation method for distributed-order fractional optimal control problems. Nonlinear Dyn 2018;91(4):2667–81. [11] Zaky MA, Machado JT. On the formulation and numerical simulation of distributed-order fractional optimal control problems. Commun Nonlinear Sci Numer Simul 2017;52:177–89. [12] Zeid SS, Kamyad AV. On generalized high order derivatives of nonsmooth functions. Am J Comput Math 2014;4(04):317. [13] Oldham K, Spanier J. The fractional calculus theory and applications of differentiation and integration to arbitrary order (Vol. 111). Elsevier; 1974. [14] Xin B, Zhang J. Finite-time stabilizing a fractional-order chaotic financial system with market confidence. Nonlinear Dyn 2014;79:1399–409. [15] Xin B, Li Y. 0 − 1 Test for chaos in a fractional order financial system with investment incentive. Abst Appl Anal 2013. [16] Jahanshahi H, Yousefpour A, Weic Z, Alcaraz R, Bekiros S. A financial hyper chaotic system with coexisting attractors: dynamic investigation, entropy analysis, control and synchronization. Chaos, Solitons Fractals 2019;126:66–77. [17] Yousefpour A, Jahanshahi H, Munoz-Pachecoc JM, Bekiros S, Wei Z. A fractional-order hyper-chaotic economic system with transient chaos. Chaos, Solitons Fractals 2020;130:1–12. [18] Wang Z, Huang X, Shen H. Control of an uncertain fractional order economic system via adaptive sliding mode. Neurocomputing 2012;83:83–8. [19] Jahanshahi H, Shahriari-Kahkeshi M, Alcaraz R, Wang X, Singh VP, Pham VT. Entropy analysis and neural network-based adaptive control of a non-equilibrium four-dimensional chaotic system with hidden attractors. Entropy 2019;21:156. [20] Chen W. Nonlinear dynamics and chaos in a fractional-order financial system. Chaos, Solitons Fractals 2008;36:1305–14. [21] Laskin N. Fractional market dynamics. Physica A 20 0 0;287:482–92. [22] Sun HG, Chen W, Wei H, Chen YQ. A comparative study of constantorder and variable-order fractional models in characterizing memory property of systems. Eur Phys J Special Top 2011;193(1):185. [23] Sun H, Chen Y, Chen W. Random-order fractional differential equation models. Signal Process 2011;91(3):525–30. [24] Sweilam NH, AL-Mekhlafi SM. Numerical study for multi-strain tuberculosis (TB) model of variable-order fractional derivatives. J Adv Res 2016;7(2):271–83. [25] Sheng H, Sun H, Chen Y, Qiu T. Synthesis of multi-fractional gaussian noises based on variable-order fractional operators. Signal Process 2011;91(7):1645–50. [26] Orosco J, Coimbra CFM. On the control and stability of variable-order mechanical systems. Nonlinear Dyn 2016;86(1):695–710. [27] Moghaddam B P, Machado J A T. A stable three-level explicit spline finite difference scheme for a class of nonlinear time variable order fractional partial differential equations. Comput Math Appl 2017;73(6):1262–9. [28] Moghaddam BP, Machado JAT. SM-Algorithms for approximating the variable-order fractional derivative of high order. Fundam Inf 2017;151(1–4):293–311. [29] Bhrawy AH, Zaky MA. Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dyn 2015;80(1–2):101–16. [30] Bhrawy AH, Zaky MA. Numerical algorithm for the variable-order caputo fractional functional differential equation. Nonlinear Dyn 2016;85(3):1815–23. [31] Zhuang P, Liu F, Anh V, Turner I. Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term. SIAM J Numer Anal 2009;47(3):1760–81.

9

[32] Lin R, Liu F, Anh V, Turner I. Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation. Appl Math Comput 2009;212(2):435–45. [33] Sun H, Chen W, Li C, Chen Y. Finite difference schemes for variable order time fractional diffusion equation. Int J Bifurc Chaos 2012;22(04):1250085. [34] Shen S, Liu F, Chen J, Turner I, Anh V. Numerical techniques for the variable order time fractional diffusion equation. Appl Math Comput 2012;218(22):10861–70. [35] Yaghoobi S, Moghaddam BP, Ivaz K. An efficient cubic spline approximation for variable-order fractional differential equations with time delay. Nonlinear Dyn 2017;87(2):815–26. [36] Moghaddam BP, Yaghoobi S, Machado JT. An extended predictor corrector algorithm for variable-order fractional delay differential equations. J Comput Nonlinear Dyn 2016;11(6):061001. [37] Keshi FK, Moghaddam BP, Aghili A. A numerical approach for solving a class of variable-order fractional functional integral equations. Comput Appl Math 2018;37(4):4821–34. [38] Moghaddam BP, Machado JAT, Behforooz H. An integro quadratic spline approach for a class of variable-order fractional initial value problems. Chaos, Solitons Fractals 2017;102:354–60. [39] Samko S. Fractional integration and differentiation of variable order: an overview. Nonlinear Dyn 2013;71(4):653–62. [40] Zhang H, Liu F, Phanikumar MS, Meerschaert MM. A novel numerical method for the time variable fractional order mobile-immobile advection-dispersion model. Comput Math Appl 2013;66(5):693–701. [41] Morales-Delgado VF, Gomez-Aguilar JF, Taneco-Hernandez MA, Escobar-Jimenez RF. A novel fractional derivative with variable-and constant-order applied to a mass-spring-damper system. Eur Phys J Plus 2018;133(2):78. [42] Lorenzo CF, Hartley TT. Initialization, conceptualization, and application in the generalized (fractional) calculus. Crit Rev Biomed Eng 2007;35(6):1–106. [43] Ingman D, Suzdalnitsky J, Zeifman M. Constitutive dynamic-order model for nonlinear contact phenomena. J Appl Mech 20 0 0;67(2):383–90. [44] Lorenzo CF, Hartley TT. Variable order and distributed order fractional operators. Nonlinear Dyn 2002;29(1–4):57–98. [45] Coimbra CF. Mechanics with variable-order differential operators. Ann Phys 2003;12(11–12):692–703. [46] Sierociuk D, Malesza W, Macias M. Numerical schemes for initialized constant and variable fractional-order derivatives: matrix approach and its analog verification. J Vib Control 2016;22(8):2032–44. [47] Valrio D, Da-Costa JS. Variable-order fractional derivatives and their numerical approximations. Signal Process 2011;91(3):470–83. [48] Zhao X, Sun ZZ, Karniadakis GE. Second-order approximations for variable order fractional derivatives: algorithms and applications. J Comput Phys 2015;293:184–200. [49] Sierociuk D, Malesza W, Macias M. Numerical schemes for initialized constant and variable fractional-order derivatives: matrix approach and its analog verification. J Vib Control 2016;22(8):2032–44. [50] Xu Y, He Z. Existence and uniqueness results for cauchy problem of variable-order fractional differential equations. J Appl Math Comput 2013;43(1–2):295–306. [51] Lifshits M, Linde W. Fractional integration operators of variable order: continuity and compactness properties. Math Nachr 2014;287(8–9):980–10 0 0. [52] Sierociuk D, Malesza W, Macias M. On the recursive fractional variable-order derivative: equivalent switching strategy, duality, and analog modeling. circuits. Syst Signal Process 2015;34(4):1077–113. [53] Moghaddam BP, Machado JAT. Extended algorithms for approximating variable order fractional derivatives with applications. J Sci Comput 2017;71(3):1351–74. [54] Ortigueira MD, Valerio D, Machado JT. Variable order fractional systems. Commun Nonlinear Sci Numer Simul 2019;71:231–43.