Neural minimization methods (NMM) for solving variable order fractional delay differential equations (FDDEs) with simulated annealing (SA)

To enrich any model and its dynamics introduction of delay is useful, that models a precise description of real-life phenomena. Differential equations in which current time derivatives count on the solution and its derivatives at a prior time are known as delay differential equations (DDEs). In this study, we are introducing new techniques for finding the numerical solution of fractional delay differential equations (FDDEs) based on the application of neural minimization (NM) by utilizing Chebyshev simulated annealing neural network (ChSANN) and Legendre simulated annealing neural network (LSANN). The main purpose of using Chebyshev and Legendre polynomials, along with simulated annealing (SA), is to reduce mean square error (MSE) that leads to more accurate numerical approximations. This study provides the application of ChSANN and LSANN for solving DDEs and FDDEs. Proposed schemes can be effortlessly executed by using Mathematica or MATLAB software to get explicit solutions. Computational outcomes are depicted, for various numerical experiments, numerically and graphically with error analysis to demonstrate the accuracy and efficiency of the methods.


Introduction
In the old days, fractional calculus was only used by pure mathematicians due to its imperceptible applications at that time. When mathematicians were trying to implement fractional calculus in modeling of physical phenomena then applicability of this marvelous tool was successfully revealed. Therefore, in recent years, fractional derivatives have been used in many phenomena in electromagnetic theory, fluid mechanics, viscoelasticity, circuit theory, control theory, biology, atmospheric physics, etc. Many real-world problems can be accurately modeled by fractional differential equations (FDEs) such as damping laws, fluid mechanics, rheology, physics, mathematical biology, diffusion processes, electrochemistry, and so on. Nowadays, there has been a tremendous increase of interest in theory of FDDEs, due to the fact that it can express the system of many dynamics population with more accuracy as demonstrated by past research in science and engineering. In the systems of real world, delays can be recognized and implemented everywhere and due to this advancement in modeling it has captured a lot of attention of the scientific community towards it.

Literature review
Many researchers were inspired by the subject of fractional order and addressed these phenomena in different scenarios. For example, Lie et al. [1] used phase portraits, time domain waveform and bifurcation to examine the behavior of nonlinear system of fractional order. Modified Kalman filter was proposed by [2] to deal with fractional order system. To achieve the suitable generalized projective synchronization (GPS) of incommensurated fractional order system, [3] has introduced fuzzy approach. Coronel-Escamilla et al. [4] used Euler-Lagrange and Hamilton formalisms for describing the the fractional modeling and control of an industrial selective compliant assembly robot arm(SCARA).
Many researchers have made their efforts to investigate DEs and FDEs such as Zhang et al. [5] studied generalized Burgers equation and a generalized Kupershmidt equation through liegroup analysis method for similarity reductions and exact solutions, Zhang and Zhou [6] also carried out analysis of Drinfeld-Sokolov-Wilson system through symmetry analysis method, Yang et al. [7,8] implemented travelling wave solution to local fractional two-dimensional Burgers-type equations and Boussinesq equation in fractal domain. Atangana and Gómez-Aguilar [9] presented exact solution and semi group principle for evolutions equation by using three different definitions of fractional derivatives. Lie et al. [10] calculated the Haar wavelet operational matrix together with Block phase function to find the solution of FDEs. Li et al. [11] have also made a successful attempt to approximate the same problem by using Chebyshev wavelet operational matrix method. Adomian decomposition method and variational iteration method were implemented in [12][13][14][15][16] to solve a variety of FDEs.While differential transform method and power series method are also noteworthy for the solution of FDEs [17][18][19][20][21][22].
In recent years the approach of neural architecture has been used to solve Des and FDEs. In [23] Aarts and Veer implemented the multilayer neural algorithm for the solution of partial differential equations (PDEs) with evolutionary algorithm for training of weights. Feed forward NN technique, with the blend of piece splines of Lagrange polynomial, was proposed by [24]. Same approach was applied by [25] in which genetic algorithm was used as an evolutionary algorithm for training of network of nonbed-catalytic gas reactor system. For solving PDEs [26] has put an effort by implementing NN together with Broden-Flecher-Goldfarb-Shanno algorithm. Nelder Meade optimization procedure with hybrid neural network (HNN) was adopted by [27] for the numerical simulation of higher order DEs. Levenberg-Marquardt algorithm with ANN and Mittag-Leffler kernel was implemented by [28] to solve FDEs.
Those systems that are governed by their past are modeled in the form of delay differential equations.DDEs are proved useful in control systems [29], lasers, traffic models [30], metal cutting, epidemiology, neuroscience, population dynamics [31], and chemical kinetics [32] condition. Due to the infinite dimensionality of delay systems it is very challenging to analytically analyze DDEs, therefore numerical simulations of DDEs play a key role for study of such systems. A noteworthy study of FDDEs through neural network can be visualized in [33] Existence and uniqueness theorems on FDDEs are discussed in [34][35][36] This study will generate an approximate solution for solving the DDEs and FDDEs by using ChSANN and LSANN, which were first developed by Khan et al. [37,38] to solve Lane Emden equations and fractional differential equations on a discrete domain. In the following paper we have developed an approach to solve the higher order DDEs and FDDEs on a continuous domain. This paper is organized as follows: first section of paper describes introduction and literature review while second section concerns with details of the methodologies with well explained algorithm and implication procedure. Error analysis procedure is explained in third section where as fourth, fifth and sixth sections describe the numerical experiments along with results and their discussions.

Methodology
The proposed methodologies are based on functional link neural network with optimization through thermal minimization. In this study, the Caputo definition will be used for working out the fractional derivative in the subsequent procedure. These definitions of commonly used fractional differential operators are discussed in [39].
Initially introduced by Pao [40], ChSANN and LSANN are the revised version of functional link artificial NN coupled with optimization strategy for learning. Functional link architecture of NN was design to build the connection between the linearity in a single layer NN and the computationally challenging multilayer NN.
Inspired by the physical process of annealing, SA is basically a kind of combinatorial optimization process. This process is based on two steps: first to perturb and then to measure quality of the solution. MSE is basically the fitness function denoted by E r , that can be minimized by the use of SA.

Algorithm
Due to the structural similarity in ChSANN and LSANN, the steps of algorithm are described in combination for both methods. Accuracy of results depends on the selection of base polynomial.
Step 1: Initialize the network by applying number of Chebyshev polynomials or Legendre polynomials (to independent variable x) k = 0 to n.
Step 2: Provide in each polynomial a network adaptive coefficient (NAC).
Step 3: Calculate the summation of product of NAC and Chebyshev polynomial or Legendre polynomial and store the value in ϕ or ψ respectively.
Step 4: Activate ϕ and ψ, by first three terms of Taylor Series expansion of tanh function.
Step 5: As given by Lagaris and Fotiadias [41] trial solution will be generated by the help of initial conditions and activated ϕ (in case of LSANN ψ) Step 6: Calculate the value of delay trail solution by repeating step 1 to 5 with delay independent variable.
Step 7: Calculate the MSE of DDEs or FDDEs by discretizing the domain in β number of points.
Step 8: Set Tolerance for accepting minimized value of MSE.
Step 9: Minimize the MSE by thermal minimizing methodology with the following settings from Mathematica 11.0.
➢ Tolerance for accepting constraint violations! 0.001 Step 10: If the value of MSE falls in the range of pre-defined criteria, then substitute the value of NAC in trial solution to get the output otherwise go to step 1 change the value of n and repeat whole procedure till the acceptable MSE is obtained.
Pictorial presentation of above algorithm can be observed in Fig 1

Employment on delay differential equation
Now we apply ChSANN and LSANN on DDEs of the following type, With initial conditions as follows, For implementation of ChSANN and LSANN Eq (1) can be written as For trial and delay trial solution of above differential equation consider, where T k is Chebyshev polynomial with the following recursive formula defined as Here T 0 = 1 and T 1 = x are the fundamental values of Chebyshev polynomialsand where L j is Legendre polynomial with the following recursive formula defined as: where, L 0 (x) = 1 and L 1 (x) = x are the fundamental values of Legendre polynomials.
Here we are using first three terms of Taylor series expansion of tanh function to activate ϕ and ψ. As defined by Lagaris and Fotiadis [41], the trial solution of Eq (1) can be written as, Where N is the activated ϕ or ψ, depends on the method, while the delay trial solution can be given by replacing the x by x−τx.
The MSE of the Eq (1) will be calculated from the following: here β represents the number of trial points while Eq (8) will be the fitness function for learning of NAC. For implementation on FDDEs only the method of computing the fractional derivative of trial solution will vary, which will be taken in this study according to Caputo definition as follows.

Definition
According to [39] Caputo operator for λ>0 can be defined as: Mathematica 11.0 is the minimization implementation tool in this study but details can be seen from [42].

Error analysis
The error analysis of numerical experiments for ChSANN and LSANN methodologies can be observed by following procedure. By substituting the values of NAC, after learning from SA algorithm, into the trial solution it will become ChSANN or LSANN solution that can be further substituted into Eq (9) for analyzing the accuracy of method on the domain of [0,1].
EðxÞ ¼ jf n ðxÞ À Fðf ðx À txÞ; f 0 ðxÞ; f @ ðxÞ; f ‴ ðxÞ; . . . ; f ðnÀ 1Þ ðxÞÞ À gðxÞj ffi 0 While f(x) is the obtained approximated continuous solution by ChSANN or LSANN. E r (x i ) tends to 0 as the value of MSE obtained by ChSANN and LSANN is in the predefined range. Convergence of solution is totally dependent on the learning methodology of respected NN architecture which is SA in the current case.  Numerical experiments Experiment 1. Consider 2 nd order DDE along with the initial conditions as: The exact solution when α = 1 is given as: In this experiment we employed proposed methodologies on above second order linear DDE on domain of [0,1].Both the methods were employed by dividing the domain with 10 equidistant training points and 6 NAC. For ChSANN and LSANN at α = 2, the MSE at defined conditions is found to be 1.89855 ×10 −11 and 1.32344 ×10 −14 respectively. (Fig 2) depicts the comparison of both methods with true solution at continuous domain of [0,1], while (Fig 3) displays the error analysis for both methods at α = 2. For the above experiment the trial and delay trail solutions are found to be following: Where, N and M are the structural outputs of both NNs. Table 1 displays the final values of NAC after training by SA algorithm and Figs 4 and 5 is displaying the data for 100 independent runs by altering the scale for random jumps. Experiment 2. Consider 3 rd order nonlinear DDE along with the initial condition as: The exact solution when α = 3 is given below:  values of weights after training by SA algorithm while Figs 8 and 9 are displaying the results for 100 independent runs for elapsed time in seconds, fitness and number of iterations. Following are trial and delay trail solution for the current experiment:

N:
While M and N are structural NN outputs of both the methods for in progress experiment. Example 3. Consider FDDE along with the initial conditions as: The exact solution at α = 1 is given by: yðxÞ ¼ e x ChSANN and LSANN have been employed successfully on above FDDE with 10 training points and 6 NAC. (Fig 10) depicts the comparison of ChSANN and LSANN solutions with true values at α = 1 while values of final NAC after learning by SA algorithm can be visualized in Table 3. Both the methods were executed for different fractional values of α for which results can be visualized in Tables 4-6. Error analysis for all fractional values can be seen in (Fig 11) and Figs 12 and 13 is displaying the results for 100 independent runs for both the proposed methods  Trial and delay trial solutions are taken to be: Example 4. Consider Non-Linear FDDE along with the initial condition as: The exact solution at α = 2 is given by:

yðxÞ ¼ CosðxÞ
ChSANN and LSANN are implemented on the above experiment of nonlinear FDDE. (Fig  14) shows the comparison of both methods with true values at α = 2. While (Fig 15) represents the error analysis for different fractional values of α. Trial solution was taken in the same manner as in previous experiments. ChSANN and LSANN solutions with obtained MSE at

Discussion
The above work is concerned with the successful implementation of ChSANN and LSANN on higher order DDEs and FDDEs. Some benchmark examples were considered for experimental cases and validity of implementation has been judged by standard error analysis procedure, data analysis of 100 number of runs of algorithm and comparison with other methods.

Error analysis
For test experiment 1, MSE for ChSANN and LSANN for α = 2 were found to be 1.89855 ×10 −11 and 1.32344 ×10 −14 that gave the minimum error for ChSANN is 2.8 ×10 −5 and for LSANN is 2.6 ×10 −7 that can be easily visualized from (Fig 3). It shows that accuracy of both the methods is inversely proportional to the value of MSE and it can also be noticed that change of polynomials in both methods has strongly influenced the learning of NAC by SA algorithm that can be witnessed from  Analysis from (Fig 7) exhibiting the better performance of LSANN with less MSE as described above. Moreover it can also be visualized that this structure of neural network can be easily implemented on higher order nonlinear differential equations with ease.
In experiment no 3, both the proposed architectures were employed on linear delay fractional differential equation. MSE for ChSANN and LSANN at α = 1 were found to be 4.63009 ×10 −11 and 3.88356×10 −11 respectively that gave excellent results that can be seen in (Fig 10).  Further both the methods were also executed for α = 0.7,0.8 and 0.9. Obtained MSE for all the fractional values obtained by both the methods were in the same range so the accuracy of results for all the fractional values is approximately similar that can be visualized in (Fig 11). Experiment no 4 is a case of nonlinear FDDE. For α = 2 values of MSE with 6 NAC were found to be 7.82099 ×10 −6 and 1.050097 ×10 −5 . For fractional values of α ChSANN is showing better results than LSANN at α = 1.5 and α = 1.9, with better values of MSE for ChSANN and for α = 1.7 both the methods are exhibiting the similar accuracy. Accuracy at fractional values can be visualized from (Fig 15).
Data analysis for 100 numbers of independent runs. For each test experiment, algorithms of proposed techniques were executed 100 times by altering the scale of random jumps to assess the precision, performance and reliability. Results of obtained data can be visualized in form of figures which shows that for test experiment 1 fitness function is revolving between 10 −4 to 10 −14 and 10 −6 to 10 −14 for CHSANN and LSANN respectively, Elapsed time in second is found to be within three seconds for both the methods while number of iterations were between 600-1200 and 400-1000 for CHSANN and LSANN respectively. Results of 100 independent runs for test experiments 2-4 can be visualized in Figs 7 and 8, Figs 11 and 12 and Figs 16 and 17 respectively, which demonstrate a similar trend except for the case of nonlinear models for which the maximum elapsed time is found to be 20 seconds due to computational complexity.
Comparison with other methods. We compared the proposed techniques in terms of accuracy, elapsed time, ease of calculation and error prediction with the methods presented in [43], [44] and [45]. Methods in [43][44][45] have been implemented on similar type of problems as in current study. Test example number 3 by Radial basis method presented in [43] and test example no 4 presented in [45] is similar to test experiment no 2 presented in following paper, Test example number 5 presented in [44] is similar to test experiment number 4 presented above and test experiment number 2 presented in [45] is similar to test experiment no 1 by proposed methods.
Following key points of comparison can be noticed.
➢ Radial basis method, method in [44] and method in [45] are providing results on collocation points while proposed schemes are providing a continuous solution. ➢ On the other hand method in [43] is taking 10 to 85 seconds for solving a linear problem while proposed techniques are consuming 6 to 12 seconds and 3 to 5 seconds (Figs 8 and 9) for solving nonlinear problem by ChSANN and LSANN respectively. ➢ Computational complexity of method presented in [43][44][45] is very large due to solving a large system of nonlinear equations while proposed techniques are too simple in terms of implementation that can be observed through the computational time difference.
➢ There is no way to predict the accuracy in method proposed in [43][44][45] when there is no exact solution present while proposed schemes can predict accuracy of solution through fitness function. In terms of accuracy methods in [43][44][45] is giving more accurate results than proposed schemes but limitations of [43][44][45] is making proposed schemes more powerful.

Conclusion
In above study we have developed two methods ChSANN and LSANN for simulation of fractional delay differential equation. After analyzing procedure and numerical experiments following points can be concluded.
➢ Proposed methods can be successfully implemented on linear and nonlinear FDDEs with ease of calculation.
➢ Accuracy of method can be increased by improving the learning methodology of NAC.
➢ Accuracy of both the methods is inversely proportional to MSE.
➢ Both the methods can easily handle the nonlinear terms.
➢ Accuracy prediction can be obtained for fractional values of derivatives by observing the MSE values.
In future the proposed schemes can be further developed for accuracy by refining the learning methodology of NAC and by improving the neural architecture. However, it can also be successfully implemented on partial differential equations with some alterations in methodology.