Numerical simulations for fractional variable-order equations⁎

Proceedings of the 3rd IFAC Conference on Proceedings of the 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Proceedings of the 3rd IFAC Conference on Control Advances in Proportional-Integral-Derivative Control Available online at www.sciencedirect.com Ghent, Belgium, May 9-11, 2018 Advances in Proportional-Integral-Derivative Proceedings of the 3rd IFAC Conference on Control Ghent, Belgium, May 9-11, 2018 Ghent, Belgium, May 9-11, 2018 Advances in Proportional-Integral-Derivative Control Ghent, Belgium, May 9-11, 2018

ScienceDirect

IFAC PapersOnLine 51-4 (2018) 853–858

Numerical simulations for fractional Numerical simulations for fractional Numerical simulations for fractional  variable-order equations Numerical simulations for fractional variable-order equations variable-order equations  variable-order equations Dorota Mozyrska ∗∗ Piotr Oziablo ∗∗ ∗∗

Dorota Mozyrska ∗ Piotr Oziablo ∗∗ Dorota Mozyrska Piotr Oziablo ∗ Dorota Mozyrska Piotr Oziablo ∗∗of Technology, Faculty of Computer Faculty of Computer Science, Science, Bialystok Bialystok University University of Technology, Faculty of Computer Science, Bialystok of Technology, Wiejska 45A, Bialystok, Poland (e-mail:University [email protected]) ∗ Wiejska 45A, Bialystok, Poland (e-mail: [email protected]) ∗∗ Faculty Computer Science, Bialystok of 45A, Bialystok, Poland (e-mail:University [email protected]) Faculty of of Computer Science, Bialystok University of Technology, Technology, ∗∗Wiejska ∗∗ Faculty of Computer Science, Bialystok University of Technology, Faculty of Computer Science, Bialystok University of Technology, Wiejska 45A, Bia lystok, Poland (e-mail: [email protected]) Wiejska 45A, Bia l ystok, Poland (e-mail: [email protected]) Wiejska 45A, Bialystok, Poland (e-mail: [email protected]) ∗∗ Faculty of Computer Bialystok University of Technology, Wiejska 45A, Bialystok, Science, Poland (e-mail: [email protected]) Wiejska 45A, Bialystok, Poland (e-mail: [email protected]) Abstract: Abstract: The The objective objective of of the the paper paper is is to to present present the the method method of of fitting fitting finding finding constant constant Abstract: The objective of the paper is to presentofthe method of fitting finding variable-, constant λ coefficient and a parameter of an order function the processes described λ coefficient and a parameter of an order function of the processes described by by variable-, Abstract: objective of the is Gr¨ to present ofAsfitting finding constant λ coefficientThe andbackward a parameter of paper anoforder function ofthe themethod processes described by variable-, fractional-order difference the u nwald-Letnikov-type. a qualitative criterion fractional-order backward difference of the u As aa qualitative criterion fractional-order difference oforder the Gr¨ Gr¨ unwald-Letnikov-type. nwald-Letnikov-type. As22described qualitative criterion λ coefficient andbackward a parameter of an function of the processes by variable-, of the estimation the Coefficient of Determination (which we mark as R ) and the Mean of the estimation the Coefficient of Determination (which we mark as R2 ) and the Mean Square Square fractional-order backward difference of the Gr¨ uwere nwald-Letnikov-type. As)aand qualitative of the are estimation thethe Coefficient ofexperiments Determination (which we mark as R the Meancriterion Square Error used. All numerical done with MATLAB. Error are used. All the numerical experiments were done with MATLAB. 2 of the are estimation thethe Coefficient Determination (which mark as R ) and the Mean Square Error used. All numericalofexperiments were done we with MATLAB. ©Error 2018,are IFAC (International Federationexperiments of Automatic were Control) Hosting Elsevier Ltd. All rights reserved. used. All theequations, numerical done with by MATLAB. Keywords: Difference Dynamic systems, Eigenfunction, Fractional Keywords: Difference equations, Dynamic systems, Eigenfunction, Fractional variable-order variable-order Keywords: Difference equations, Dynamic systems, Eigenfunction, Fractional variable-order ∗ ∗ ∗

Keywords: Difference equations, Dynamic systems, Eigenfunction, Fractional variable-order 1. INTRODUCTION fractional order parameters was inspired by Almeida & fractional order parameters was inspired by Almeida & 1. 1. INTRODUCTION INTRODUCTION fractional order parameters waspaper inspired byfirst Almeida & Bastos & Monteiro (2018). The is the attempt Bastos & Monteiro (2018). The paper is the first attempt 1. INTRODUCTION fractional order parameters was inspired by Almeida & Bastos & Monteiro (2018). The paper is the first attempt to wider research of fitting values of fractional-, variableFractional order calculus became an useful tool in mod- to wider research of fitting values of fractional-, variableFractional order calculus became an useful tool in modBastos & Monteiro (2018). The paper is the first attempt to wider research of fitting values of fractional-, variableorder function. The method is planned to be used in fracFractional order calculus became an useful tool in modelling, which could be successfully used in both: continuous- order function. The method is planned to be used in fracelling, which could be used in both: continuouswider research of method fitting The values of fractional-, order function. The is planned to be used in state fractional PID-like controllers. important goal isvariableto elling, which couldcalculus be successfully successfully in involve both: continuousFractional order became an useful tool in mod- to and discrete-time systems. Modelsused which fractional tional PID-like controllers. important goal is state and discrete-time systems. Models which involve fractional PID-like controllers. The important goal isasto to state order function. The methodThe iscan planned to be used inbetter fracwhat type of order functions be considered and discrete-time systems. Models which fractional elling, whichand could be successfully used in involve both: continuousderivatives fractional difference operators can very tional what type of order functions can be considered as better derivatives and fractional difference operators can very tional PID-like controllers. The important goalmethods isastobetter state what type of order functions can be considered simulations of measurements. Moreover, the derivatives andreal fractional difference operators can than very in and discrete-time systems. Models which involve fractional often describe phenomena in more accurate way in simulations of measurements. Moreover, the methods of of often describe real phenomena in more accurate way than what type of order functions can be considered as better in simulations of measurements. Moreover, the methods of often describe real phenomena in more accurate way than derivatives and fractional difference operators can very it is possible with integer order models. The theory and fitting fitting parameters parameters of of equations equations and and order-functions order-functions are are the the it is possible with integer order models. The theory and in simulations of measurements. Moreover, the methods of fitting parameters of modelling equations and order-functions are the issue in possible of controllers in closed-loop it is possible integercalculus order The theory and core often describe real phenomena in models. more way than applications ofwith fractional can be beaccurate reviewed in (Hilcore issue in modelling of in applications of fractional calculus can in (Hilcore issue in possible possible of controllers controllers in closed-loop closed-loop parameters of modelling equations and order-functions are the systems. applications ofwith fractional calculus can be&reviewed reviewed in 2015, (Hilit is 2000; possible integer order models. The theory and fitting fer, Kaczorek, 2009; Mozyrska Wyrwas, systems. fer, 2000; Kaczorek, 2009; Mozyrska Wyrwas, core issue in possible modelling of controllers in closed-loop fer, 2000; Kaczorek, Mozyrska & Wyrwas, 2015, applications of fractional calculus can be& reviewed in 2015, (Hil- systems. 2017; Ostalczyk, 2016; 2009; Podlubny, 1999). For variable-order 2017; Ostalczyk, 2016; Podlubny, 1999). For variable-order 2. 2017; Ostalczyk, 2016; Podlubny, 1999). For variable-order fer, 2000; Kaczorek, 2009; Mozyrska & Wyrwas, 2015, applications the reader can see more in Mozyrska & Os- systems. 2. PRELIMINARIES PRELIMINARIES applications the reader can see more in Mozyrska & Os2. PRELIMINARIES applications theIn reader can see more inFor Mozyrska & Os2017; Ostalczyk, 2016; Podlubny, 1999). variable-order talczyk (2017). the paper we investigate discrete-time talczyk (2017). In the paper we investigate discrete-time & 2. PRELIMINARIES talczyk (2017). Inreader the paper we investigate discrete-time applications thevariable-orders. can see more & Os- Definition Definition 1. 1. (Mozyrska (Mozyrska & Ostalczyk Ostalczyk (2017)) (2017)) operators with One of in theMozyrska most important Definition 1. (Mozyrska & Ostalczyk operators with variable-orders. One of the most important For k, l ∈ N and a given order function(2017)) ν(·) we define the operators with variable-orders. One of the most important talczykin(2017). In the paper awe investigate discrete-time issues modeling is finding function which fits the exFor k, l ∈ N and a given order function ν(·) we a function which fits the exDefinition 1. (Mozyrska & Ostalczyk (2017)) For k, l ∈ N and a given order function ν(·) two we define define the the issues in modeling is finding oblivion function, as a discrete function issues in modeling is finding operators with variable-orders. of Surez the most important a One function which fits the ex- oblivion function, as a discrete function of of two variables, variables, perimental data, Garcia & Aguirre & (2008). In the [ν(l)] For k, l ∈ N and a given order function ν(·) we define the oblivion function, as a discrete function of two variables, perimental data, Garcia & Aguirre & Surez (2008). In the (k) given for k > 0 as by its values a [ν(l)] perimental data, Garcia & Aguirre & Surez (2008). In the issues we in modeling isthe finding a function which fits the ex- by its values a[ν(l)] (k) given for k > 0 as paper deal with problem of identifying parameters (k)ν(l) given for kfunction > 0· [ν(l) as of−two by its values a paper we deal with the problem of identifying parameters oblivion function, as a discrete variables, [ν(l) − 1] · · k + 1] paper we deal with problem of fractional identifying parameters perimental data, Garcia & Aguirre & Surez (2008). In the [ν(l)] k ν(l) [ν(l) − 1] · · · [ν(l) − k + 1] of a variable-order inthe discrete-time equations by (k) = (−1) aits [ν(l)] k ν(l) (k) given as − k + 1] ,, (1) values a[ν(l)] of a discrete-time equations by [ν(l)for − k1]> ·k! · 0· [ν(l) = (−1) of a variable-order variable-order inthe discrete-time fractional equations by by a [ν(l)] (k) k we deal within problem of fractional identifying parameters apaper variety of techniques, mainly based on MATLAB rou(k) = (−1) a , (1) (1) k! a variety of techniques, mainly based on MATLAB rouν(l) [ν(l) − 1] · · · [ν(l) − k + 1] [ν(l)] k atines. of techniques, mainly based on MATLAB k! of variety a variable-order in discrete-time fractional equations by anda[ν(l)] Variable-order derivatives and differences are rounew (0) = 1. a [ν(l)] (k) = (−1) , (1) (0) = 1. a tines. Variable-order derivatives and are new k! tines. Variable-order derivatives and differences differences are rounew and a variety of techniques, mainly based on MATLAB directions in the research of fractional systems. Moreover, and a[ν(l)] (0) = 1. directions in the research of fractional systems. Moreover, Formula 11 is directions in research of fractional systems. Moreover, tines.there Variable-order derivatives and are (0) in = Definition 1. a[ν(l)](1) now arethe research activities that differences are focused onnew de- and Formula (1) in Definition is equivalent equivalent to to the the following following now there are research activities that are focused on deFormula (1) in Definition 1 is equivalent to the following recurrence with respect to k ∈ N now there are research activities that are focused on dedirections in the research of fractional systems. Moreover, veloping new analysis and closed-loop closed-loop system synthesis synthesis recurrence with respect to k ∈ N veloping new analysis and system recurrence with respect to k ∈ N [ν(l)] (1) in Definition 1 is equivalent to the following veloping new analysis system now therefor are research and activities that are focused on de- Formula closed-loop synthesis methods fractional-order controllers being an extension a (0) = methods for fractional-order controllers being an a[ν(l)] = 11 ,, [ν(l)] (0) with  recurrence methods for fractional-order controllers being an extension extension veloping new analysis and An closed-loop system synthesis of classical control theory. appropriate choice of the a (0) = 1 , respect to  k ∈ Nν(l) + 1  (2) of classical control theory. An appropriate choice of the [ν(l)] [ν(l)]   for k  1 . (2) ν(l) + 1 of classical theory. An of methods for control fractional-order controllers beingchoice extension (k) = 1a[ν(l)] a order functions in fractional PIDappropriate controllers isanstill onethe of (0) , (k − 1) 1 − (2) a[ν(l)] order functions in fractional PID controllers is still one of [ν(l)] (k) = a[ν(l)] (k − 1) 1 − ν(l)k+ 1  for k  1 . order functions in fractional PIDappropriate controllers choice is still of onethe of (k) = a (k − 1) 1 − ν(l)k+ 1 for k  1 . (2) a of classical control theory. An the open problems. the open problems. k (2017)) Definition 2. (Mozyrska & Ostalczyk the open problems. = (Mozyrska a[ν(l)] (k − 1) 1− for k  1 . a[ν(l)] (k)2. order functions in fractional PID controllers is still one of Definition & Ostalczyk (2017)) In this work we calculate values of scalar eigenvalue funcDefinition 2. (Mozyrska & Ostalczyk (2017)) The Gr¨ unwald-Letnikov variable-, kfractional-order backthethis open problems. In work we calculate values of scalar eigenvalue funcGr¨ u variable-, fractional-order backIn thisforwork we calculate values of scalar eigenvalue tions fractional variable-order initial value problemfuncand The Definition 2. (Mozyrska & Ostalczyk (2017)) The Gr¨ unwald-Letnikov nwald-Letnikov variable-, fractional-order backward difference (GL-VFOBD) with an order function νν :: tions for fractional variable-order initial value problem and ward difference (GL-VFOBD) with an order function tions fractional variable-order initial value problem and The In thisfor calculate of scalar eigenvalue functhen wework are we fitting back values the noised eigenvalue function. Gr¨ u nwald-Letnikov variable-, fractional-order backward difference (GL-VFOBD) with an order function ν : Z → R ∪ {0} of function x(·) is defined as a finite sum + then we are fitting back the eigenvalue function. {0} of function x(·) is defined as finite → +∪ then we fractional areprocess fittingvariable-order theonnoised noised eigenvalue function. tionsfitting for initial valuethat problem and Z The isback based assumption the class Z →R Rdifference function defined as aa function finite sum sumν : ward with an order kx(·) is +∪ {0} of(GL-VFOBD)  The fitting process is based on assumption that the class k The fitting isback based onnoised assumption that functions the class Z → R+∪ {0} then we function areprocess fitting the eigenvalue function. of order is known. There are considered  [ν(k)] of  is defined − asi)a finite sum kx(·) ∆[ν(k)] function (k) = a of order function is known. There are considered functions  [ν(k)] x [ν(k)] (i)x(k ∆ x (k) = a of order function is known. There are considered functions The fitting process is based on assumption that the class [ν(k)] [ν(k)] (i)x(k − i) different class of monotonicity. Moreover, the constant λ k ∆ x (k) = a (i)x(k − i)   i=0  of different class ofis monotonicity. Moreover, the λ of different class monotonicity. Moreover, the constant constant λ  order function known. Therewith are considered functions i=0 [ν(k)] [ν(k)] coefficient of the of linear equation variable-, fractional∆ x (k) = a (i)x(k − i)    i=0 coefficient of the linear equation variable-, fractionalcoefficient of the linear equation with with variable-, fractional x(k)  of different class of monotonicity. Moreover, the constant λ order operator and the parameter of order function should x(k) i=0 order operator and the parameter of order function should x(k)   (3)    [ν(k)] order operator and thefitting parameter of order function should coefficient of the linear equation with variable-, fractionalx(k  be estimated during procedure. The idea of usx(k − − 1) 1) (3) [ν(k)] x(k)  . be estimated during fitting procedure. The idea of usx(k − 1) · · · = (1) · · · a (k) 1 a (3) [ν(k)] [ν(k)]   be estimated procedure. idea of usorder operator during and thefitting parameter of orderThe function should   ing particular MATLAB routine (lscurvefit) to estimate . = 1 a[ν(k)] (1) · · · a[ν(k)] (k)  · · ·    x(k − 1) ing particular MATLAB routine (lscurvefit) to estimate · · · . = x(1) (1) · · · a (k) 1 a (3)     [ν(k)] ing particular during MATLAB routine (lscurvefit) estimate be estimated fitting procedure. The to idea of usx(1)    [ν(k)] x(1) · · · . = x(0)   (1) · · · a (k) 1 a  work was supported by Polish founds of Nationalto Science CeningThe particular MATLAB routine (lscurvefit) estimate   x(0) The work was supported by Polish founds of National Science Cenx(1)  x(0)  ter, granted on the basis ofby decision The work was supported Polish DEC-2016/23/B/ST7/03686. founds of National Science Center, granted on the basis of decision DEC-2016/23/B/ST7/03686. x(0)  ter, granted on the basis ofby decision The work was supported Polish DEC-2016/23/B/ST7/03686. founds of National Science Cen-

ter, granted theIFAC basis (International of decision DEC-2016/23/B/ST7/03686. 2405-8963 © 2018, IFAC Federation of Automatic Control) Copyright © on 2018 853 Hosting by Elsevier Ltd. All rights reserved. Copyright © under 2018 IFAC 853 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 853 10.1016/j.ifacol.2018.06.122 Copyright © 2018 IFAC 853

IFAC PID 2018 854 Ghent, Belgium, May 9-11, 2018

Dorota Mozyrska et al. / IFAC PapersOnLine 51-4 (2018) 853–858

In the paper we consider the basic form of linear equation with variable-order   ∆[ν(k)] x (k) = λx(k − 1) + u(k − 1) , k ≥ 1

(4)

with initial condition x(0) = x0 , constant coefficient λ ∈ R. The function u(·) is given input signal. Then we can solve equation (4) by recurrence: x(k) = −

k  i=1

a[ν(k)] (i)x(k − i) + u(k − 1) , k ≥ 1 .

(5)

If in equation (4) we use constant order function equals zero, then we receive classical linear recurrence of first order. In other cases we include memory into the system. Let us define the following matrix, see for example Mozyrska & Ostalczyk (2017)

data can be described as yk = x(k) + , where  ∼ N (0, σ) is added random value. The simplest criterion to fit λ would be in this case the minimisation of the following function Sk (λ) := EkT Ek , (10) where   −1   λy(k − 1) Ek = y(k) − A[ν(k)] + u(k) . (11) y(0) To simplify the formulation of minimiser we slightly changed formula (11), which enables easier calculations and (based on the experimental results) should not have a big impact on the final results. The updated formula which we use is the following 

[ν(k)]

−1 

  −1 y(k − 1) u(k) . + A[ν(k)] y(0) (12)

 1 a[ν(k)] (1) · · · a[ν(k)] (k) 0 1 · · · a[ν(k−1)] (k − 1)   A[ν(k)] =  .  . (6) .. ..  ..  . . 0 0 ··· 1 Instead of working with recurrence, we can use a matrix form of defined matrices A[ν(k)] . Moreover, let us use the following notation:   x(k) x(k − 1)   ..  x(k) =  (7) .    x(1)  x(0) and similarly for input signal   u(k) u(k − 1)   .. . (8) u(k) =  .    u(1)  u(0)

Ek = y(k) − λ A

and it looks like the series of algebraic solutions   −1   λx(k − 1) [ν(k)] + u(k) , k ≥ 1 . (9) x(k) = A x(0)

Then, theoretical values y ˆ(k) we have from    −1  λ0 y(k − 1) [ν(k)] y ˆ(k) = A + u(k) . y(0)



Then, equation (4) can be written in a matrix form   λx(k − 1) + u(k) , k ≥ 1 A[ν(k)] x(k) = x(0)

3. ANALYTICAL APPROXIMATION During the first stage of our research we are trying to find constant coefficient λ in an analytical way. Our GLVFOBD calculations are based on (9) matrix definition. We assume that the order function is known and because of that we also know the matrices given by formula (6). As a qualitative criterion of the fitting process we are using Mean Squared Error (which in our case would be a function of λ which we mark as Sk (λ)) and the Coefficient of Determination which we mark as R2 . We are doing numerical experiment by taking the values of exact solution to (9). Then, we add random values with mean equal to 0 and standard deviation equal to 20% of the mean of the given values of function x(k). This simulates real time measurement data. In this case the simulated 854

The similar situation has been considered in the paper, where we investigated approximations of λ for equations without u(·) function, in Oziablo (2018). Let us use simplified Ek , then we can calculate the first derivative of Sk (λ) in the following way (13) Sk (λ) = EkT Ek + EkT Ek . From that the critical point of Sk (λ) is for Sk (λ) = 0, we receive λ0 which minimises Sk for the set of k + 1 values of measurements. It is given by the formula     T  1 yk−1 T yk A y0 (AT )yk + yk−1 λ0 = y0 2d   (14)  T  T yk−1 T T , A Au − u A A − yk−1 y0 k k y0 where

d=



T y0 yk−1





y A A k−1 y0 T



(15)

 −1 and A = A[ν(k)] . We also use the abbreviations that yk := y(redk). (16)

In our research we used four different order functions: 0.5 , ν4 (k) = ν1 (k) = e−0.2k , ν2 (k) = 1−e−0.2k , ν3 (k) = 1− k+1 2 sin (0.5k). The results of analytical lambda calculations are shown in Figures 1, 2, 3, 4. The invocation results summary are shown in Table 1. Order function ν1 (k) ν2 (k) ν3 (k) ν4 (k)

Estimated λ 0.26897 0.26458 0.24310 0.25784

R2 0.95305 0.65547 0.37761 0.83715

Sk (λ) 0.023426 7.816600 499.2796 0.027874

Table 1. Analytical λ estimation results. The analytical method of λ estimation returned the best result for the order function ν1 (k). Even when the estimated λ value 0.26897 significantly differed from the λ of

IFAC PID 2018 Ghent, Belgium, May 9-11, 2018

Dorota Mozyrska et al. / IFAC PapersOnLine 51-4 (2018) 853–858

855

Fig. 1. Analytical estimation of λ for order function ν1 (k) = e−0.2k .

Fig. 3. Analytical estimation of λ for order function 0.5 ν3 (k) = 1 − k+1 .

Fig. 2. Analytical estimation of λ for order function ν2 (k) = 1 − e−0.2k (Levenberg-Marquardt).

Fig. 4. Analytical estimation of λ for order function ν4 (k) = sin2 (0.5k).

the original, noised function (which was 0.3) Coefficient of Determination (R2 ) higher than 0.9 proofs that the function was properly fitted. Also for the order function ν4 (k) the fitting it gives a good result (R2 equal to 0.83715), even though the estimated λ value is equal to 0.25784 (in comparison to original 0.3). For the order function ν2 (k) the values of R2 is equal to 0.65547 and we get the satisfactory fit (R2 value between 0,6 - 0,8 usually is considered as satisfactory). What is interesting, that the estimation error of λ was in this case lower than e.g. for ν3 (k) order function. The worst fitting result we get for order function ν3 (k). Then, R2 equals to 0.37761, means that the fitting is not sufficient (R2 value between 0 - 0.5 usually is considered as not sufficient).

To fit our simulated experimental data with equation given by (9) we use lscurvefit MATLAB routine. The routine takes as an argument the function, which in our case calculates GL-VFOBD values, the set of experimental data and the initial value of searched parameter (in our case the initial value of λ). As the result the routine returns the parameter value for which the Mean Squared Error between function values (which calculates in our case GLVFOBD) and experimental data is the lowest. The tests were done for the same set of four order functions as in the previous paragraph. Also the input data were exactly the same to make it easier to compare the results with the analytical solution described in the previous paragraph.

4. APPROXIMATION OF UNKNOWN λ BY MATLAB ROUTINE During the second part of our research we are trying to find a constant coefficient λ using MATLAB routines. We assume that the order function is known. 855

We checked the results for two different searching algorithms supported by lscurvefit routine which are TrustRegion-Reflective Least Squares and Levenberg-Marquardt Method. In all the cases both methods gave exactly the same results. The results of the invocations are shown in Figures 5, 6, 7, 8 (with marked Coefficient of Determination and Mean Squared Error values).

IFAC PID 2018 Ghent, Belgium, May 9-11, 2018 856

Dorota Mozyrska et al. / IFAC PapersOnLine 51-4 (2018) 853–858

Fig. 5. Estimation (lsqcurvefit) of λ for order function ν1 (k) = e−0.2k .

Fig. 7. Estimation (lsqcurvefit) of λ for order function 0.5 ν3 (k) = 1 − k+1 .

Fig. 6. Estimation (lsqcurvefit) of λ for order function ν2 (k) = 1 − e−0.2k (Levenberg-Marquardt).

Fig. 8. Estimation (lsqcurvefit) of λ for order function ν4 (k) = sin2 (0.5k).

The invocation results summary is shown in Table 2.

lower even than for the order function ν1 (k) for which the estimation error was relatively high λ. The possible reason of getting such results is the noise added to the original function (with λ = 0.3) which may increase the calculated Coefficient of Determination.

Order function ν1 (k) ν2 (k) ν3 (k) ν4 (k)

Estimated λ 0.32682 0.30473 0.30404 0.30972

R2 0.99206 0.99286 0.98801 0.99601

Sk (λ) 0.0039607 0.1620900 9.6219000 0.0006827

Table 2. λ estimation with lsqcurvefit routine.

5. APPROXIMATION OF UNKNOWN λ AND ORDER FUNCTION VALUES

Using lscurvefit MATLAB routine we received very good results. For all the order functions Coefficient of Determination (R2 ) significantly exceeded 0.9 (in most of the cases it was higher than 0.99) which means very good fitting. Also estimated λ values were in most of the cases very close to the original λ value of the noised signal. The lscurvefit MATLAB routine returned the best result for the order function ν4 (k) (R2 equal to 0.99601). What is interesting is that while the Coefficient of Determination was the highest for ν4 (k), the estimation error of λ was the lowest for ν3 (k). At the same time for the order function ν3 (k) the value of R2 was the lowest one (but still significantly above 0.9),

In the last step of our research we tried to fit the constant coefficient λ and the parameter p of the order function assuming that the general class of the order function is known. For the tests we took given in previous section formulas of the order functions (where p is searched parameter). As in the previous steps we used lscurvefit, but in this case the routine was provided with the initial values of λ and p and as a result it returned the optimal values of these two parameters which minimize Mean Squared Error between experimental data and GL-VFOBD values. For the order functions ν1 (k) and ν2 (k) the original p parameter value was -0.2. For the order functions ν3 (k)

856

IFAC PID 2018 Ghent, Belgium, May 9-11, 2018

Dorota Mozyrska et al. / IFAC PapersOnLine 51-4 (2018) 853–858

857

and ν4 (k) as the original p parameter value we took 0.5. For lscurvefit routine in most of the cases both (LevenbergMarquardt and Trust Region) algorithms gave the same results. Just for the order function ν2 (k) = 1 − e−0.2k the results were different. The results of lscurvefit routine invocations are shown in Figures 9, 10, 11, 12, 13.

Fig. 11. Estimation (lsqcurvefit) of λ and p parameter for order function ν2 (k) = 1 − e−0.2k (Trust Region).

Fig. 9. Estimation (lsqcurvefit) of λ and p parameter for order function ν1 (k) = e−0.2k .

Fig. 12. Estimation (lsqcurvefit) of λ and p parameter for 0.5 order function ν3 (k) = 1 − k+1 .

Fig. 10. Estimation (lsqcurvefit) of λ and p parameter for order function ν2 (k) = 1 − e−0.2k (LevenbergMarquardt). The invocation results summary is shown in Table 3 (LM - Levenberg-Marquardt algorithm, TR - Trust Region algorithm). Using lscurvefit MATLAB routine to estimate Order function ν1 (k) ν2 (k)(LM ) ν2 (k)(T R) ν3 (k) ν4 (k)

Est. λ 0.2773 -1.9979 0.2592 0.2643 0.9535

Est. p -0.1808 -22.428 -0.2820 0.0336 0

R2 0.9954 0 0.9943 0.9940 0.7108

Sk (λ, p) 0.0022 32.5162 0.1188 4.8367 0.0516

Table 3. λ and p parameter estimation with lsqcurvefit routine. 857

λ and p parameters gave very good results for order functions ν1 (k) and ν3 (k). In both cases we had Coefficient of Determination (R2 ) higher than 0.99 which means a very good fitting. For the order function ν2 (k) lscurvefit routine returned good results only if we used Trust Region algorithm (R2 higher than 0.99). Fitting did not work at all for order function ν2 (k) using lscurvefit routine with Levenberg-Marquardt algorithm (R2 equal to 0). For order function ν4 (k) fitting results were not as good as for the rest of the functions but still one can considered as satisfactory (R2 value between 0.6 and 0.8). What is worth to be noticed, it is that estimated values of λ and p parameters were usually much different than the original parameter values, even while the Coefficient of Determination proved that the fitting process was correct and the fitting error was very low.

IFAC PID 2018 858 Ghent, Belgium, May 9-11, 2018

Dorota Mozyrska et al. / IFAC PapersOnLine 51-4 (2018) 853–858

ference on Mathematical Modelling, Vienna, Austria, February 21-23, 2018, p.85–86. Podlubny, I. (1999). Fractional Differential Equations. Mathematics in Sciences and Engineering. 198. Academic Press, San Diego.

Fig. 13. Estimation (lsqcurvefit) of λ and p parameter for order function ν4 (k) = sin2 (0.5k). 6. CONCLUSION Both λ and parameter of order function can be successfully estimated using MATLAB. Function fitting for most of the examples was very high, taking as a qualitative criterion the Coefficient of Determination R2 . Also presented analytical method to find λ gave good results in most of the examples, but in this cases the fitting quality depends on the order function. REFERENCES Almeida, R. & Bastos, N.R.O. & Monteiro, M.T.T. (2018). A fractional Malthusian growth model with variable order using an optimization approach. Published online in International Academic Press. Garcia, R. & Aguirre, B. & Surez, R.. (2008). Stabilization of Linear Sampled-Data Systems by a Time-Delay Feedback Control. Mathematical Problems in Engineering. Vol. 2008, Art. ID 270518, 15 pages, ISSN: 1024-123X. DOI: 10.1155/2008/270518. Hilfer, R. (2000). Applications of Fractional Calculus in Physics. World Scientific Publishing Company, Singapore. Kaczorek, T. (2009). Fractional positive linear systems. Kybernetes, 38(7/8):1059-1078. D. Mozyrska and M. Wyrwas. Stability of discrete fractional linear systems with positive orders. IFACPapersOnLine 50 (1), 8115-8120 (2017). Mozyrska, D. & Wyrwas, M. (2015). The Z-transform method and delta type fractional difference operators. Discrete Dynamics in Nature and Society, 2015:25 pages. Mozyrska, D. & Ostalczyk, P. (2017). Generalized Fractional-Order Discrete-Time Integrator. Complexity, Volume 2017, Article ID 3452409, 11 pages. Ostalczyk, P. (2016). Discrete fractional calculus. Applications in control and image processing. World Scientific Publishing Co Pte Ltd, vol. Series in Computer Vision - Vol. 4. Oziablo, P. (2018). Numerical simulations for fractional variable-order difference eigenfunctions. MATHMOD 2018 Extended Abstract Volume, 9th Vienna Con858