Fractional rheology-informed neural networks for data-driven identiﬁcation of viscoelastic constitutive models

Developing constitutive models that can describe a complex ﬂuid’s response to an applied stimulus has been one of the critical pursuits of rheologists. The complexity of the models typically goes hand-in-hand with that of the observed behaviors and can quickly become prohibitive depending on the choice of materials and/or ﬂow protocols. Therefore, reducing the number of ﬁtting parameters by seeking compact representations of those constitutive models can obviate extra experimentation to conﬁne the parameter space. To this end, fractional derivatives in which the differential response of matter accepts non-integer orders have shown promise. Here, we develop neural networks that are informed by a series of different fractional constitutive models. These fractional rheology-informed neural networks (RhINNs) are then used to recover the relevant parameters (fractional derivative orders) of three fractional viscoelastic constitutive models, i.e., fractional Maxwell, Kelvin-Voigt, and Zener models. We ﬁnd that for all three studied models, RhINNs recover the observed behavior accurately, although in some cases, the fractional derivative order is recovered with signiﬁcant deviations from what is known as ground truth. This suggests that extra fractional elements are redundant when the material response is relatively simple. Therefore, choosing a fractional constitutive model for a given material response is contingent upon the response complexity, as fractional elements embody a wide range of transient material behaviors.


Introduction
The quest for finding appropriate closed-form constitutive models to represent certain material behaviors has witnessed decades of rigorous research, which was initially galvanized mainly by those observed behaviors (Bingham 1916;Morrison 2001;Bird et al. 1987).Viscous, elastic, plastic, and thixotropic behaviors are commonly capable of categorizing material responses to an imposed actuation.In particular, viscoelastic materials, in which the stress response decays over some characteristic time, inherit the viscous response of fluids and the elastic one from solids.The time for the decay to occur is defined as the relaxation time of the fluid, i.e., the Donya Dabiri, Milad Saadat, and Deepak Mangal contributed equally to this work.B Safa Jamali s.jamali@northeastern.edu 1 Department of Mechanical and Industrial Engineering, Northeastern University, Boston 02115, Massachusetts, USA time needed for the material to exhibit an elastic response and relax, which leaves the material with a steady-state viscous response.Canonically, the simplest constitutive model that captures linear viscoelasticity is the Maxwell model1 (Morrison 2001), whose (scalar) shear stress response (σ 21 = σ ) under the imposition of a simple shear is as follows: where η and G are the viscosity (in Pa•s) and elastic modulus (in Pa), respectively, the ratio of which (η/G) may be thought of as a relaxation time, t is time (in s), and ˙ (t), in s −1 , is the imposed deformation rate 2 .In fact, the Maxwell model is derived by assembling a spring (representative of the elastic response) and a dashpot (the viscous response) in series, superimposing the deformation ( (t)) of these two elements, and solving for σ (t).Polymer melts with a small and narrowly distributed molecular weight under small deformations are adequately modeled by the Maxwell constitutive equation.In another canonical model, the Kelvin-Voigt assembly considers the spring-dashpot elements stacked in parallel, which is suitable for creep (stress-controlled) experiments, for instance.While these simple models are based on a single characteristic relaxation time, in real-world viscoelastic fluids, a wide range of relaxation times are commonly present.
Classically, this can be addressed by assembling a number of spring-dashpot (Maxwell) branches in parallel and thus capturing a range of relaxation times, also known as the generalized Maxwell model.On the other hand, in many rheologically relevant materials such as gels and worm-like micellar solutions and a range of polymeric systems, the material exhibits a power-law response, i.e., (t) ∝ t a .Such power-law responses are not well described by canonical rheological models such as the Maxwell or Kelvin-Voigt models (Tschoegl 2012;Jaishankar and McKinley 2013).
As mentioned above, the viscoelastic response inherits the elastic response of a solid (σ (t) ∝ ) and the viscous behavior of a Newtonian fluid (σ (t) ∝ ∂ ∂t ).Scott-Blair (Scott-Blair 1947) proposed a compact way to describe viscoelasticity by borrowing the concept of fractional derivatives (Stankiewicz 2018;Jaishankar and McKinley 2013): where E and τ are the elastic modulus (in Pa) and relaxation time (in s), respectively, and 0 ≤ α ≤ 1 is the fractional derivative order, effectively creating an element that interpolates between the constitutive responses of a spring and a dashpot; see Fig. 1.Such elements are usually referred to as Scott-Blair or spring-pot elements.In this formalism, the product of E and τ α (in Pa•s α ) may be thought of as a quasiproperty, V, whose unit depends on the fractional derivative.
The seemingly unphysical perception of this quasi-property (and its non-integer dimension in time) may be justified by the fact that V represents the firmness3 of a material (Scott-Blair and Coppen 1942).
Fruitful efforts have been made to approximate (and model) the observed behavior of different fluids using the concept of fractional derivatives.Highly anomalous butyl rubber (Scott-Blair 1947), soft particulate gels (Bantawa et al.  (t)  dt α , that compactly describes a material that embodies both solid (α = 0) and viscous (α = 1) responses 2022), cheese (Bonfanti et al. 2020;Faber et al. 2017), Xanthan gum (Jaishankar and McKinley 2014), and biological samples such as highly viscoelastic bovine serum albumin (BSA) (Jaishankar and McKinley 2013) have been successfully modeled with respective fractional models, while the classical counterparts of these models often require an exceedingly large number of the canonical Maxwell (or Kelvin-Voigt) elements and thus circumscribing the applicability of classical approaches to model viscoelastic samples with power-law responses (Tschoegl 2012;Jaishankar and McKinley 2013).Consistent and approachable methods, therefore, are called for to identify fractional models apposite to such viscoelastic samples.The issue, however, is that solution of fractional constitutive models, compared to classical ordinary differential equations that are used in rheological constitutive models, is far from trivial.As such, platforms that can automatically solve for different forms of fractional models may pave the way for further development of this class of compact models.
In a forward problem, a constitutive equation, e.g., a fractional viscoelastic model (FVEM), is at hand with all the initial conditions and parametric information, i.e., the fractional derivative order (α) and quasi-properties, and the task is to solve for a rheological property such as σ (t).Several methods have been proposed and employed in the literature to solve FVEMs, such as the Laplace transform (Mainardi 2010;Feng et al. 2020;Fang et al. 2020), finite-difference (Shen 2020;Lin et al. 2007;Diethelm et al. 2005), finite-element (Alotta et al. 2018), and the Adomian decomposition method (Tripathi et al. 2010), among others.However, in inverse problems, measurable observations are gathered from a sample with missing information about the constitutive equation that governs the rheological process, and the goal is to recover those missing parameters by leveraging the collected data.The algorithms mentioned above, however, often struggle in the face of inverse problems and ill-posedness (Jagtap et al. 2022), meaning that a solution may be either nonexistent or not unique.The parameter recovery of generalized Maxwell models (by stacking finite spring-dashpot units) can quickly become prohibitive as the number of units increases.Data-driven methodologies, on the other hand, can tolerate ill-posedness and thus are suitable choices when tackling inverse (and forward) problems (Chen et al. 2020) to provide unified platforms for the identification and discovery of FVEMs.
Machine learning platforms have shown promise to obtain forward and inverse solutions of differential equations alike (Raissi et al. 2019;Karniadakis et al. 2021;Raissi et al. 2020;Raissi 2018;Karnakov et al. 2022;Li et al. 2022).By inducing physical intuition into the training process, physics-informed neural networks (PINNs) were developed for the forward and inverse analysis of many scientific problems (Cai et al. 2022;Lawal et al. 2022;Ishitsuka and Lin 2023;Asrav and Aydin 2023;Tang et al. 2023).By implicitly (or explicitly) informing the neural net (NN) of a rheological constitutive model, several rheology-informed neural network (RhINN) platforms were recently developed (Mahmoudabadbozchelou and Jamali 2021;Mahmoudabadbozchelou et al. 2021Mahmoudabadbozchelou et al. , 2022b, a;, a;Saadat et al. 2022).Other efficient and promising data-driven algorithms and tools have also been developed to recover parameters in rheological problems of interest (Freund and Ewoldt 2015;Armstrong et al. 2017;Singh et al. 2019;Thakur et al. 2022;Reyes et al. 2021).Generally, RhINNs can be used to detect, construct, and solve a variety of rheological constitutive models in the form of ordinary differential equations.Nonetheless, as these platforms have not been developed and tested for fractional derivatives, this work specifically focuses on the identification of FVEMs using RhINNs in an inverse implementation.
Here, as a subset of all the possible combinations of fractional viscoelastic models, three configurations of FVEMs, namely, the fractional Maxwell model (FMM), fractional Kelvin-Voigt model (FKVM), and fractional Zener model (FZM), which is a fractional Maxwell unit connected to a spring-pot in parallel are selected.The goal is to recover the fractional derivative order of each spring-pot element, assuming that other quasi-properties are known from previous experiments or rheometry.This task is far from trivial, as inverse-problem solvers for fractional constitutive equations are yet to be popularized, and the pressing need for material discovery propels the numerical and data-driven methods to confine the parameter space of inverse problems and eventually recover those parameters.While these cases by no means encompass all the possible combinations and intricacies of FVEMs, the selected configurations nonetheless hand valuable information about the applicability of RhINNs in solving fractional inverse problems.The ultimate goal is to develop a generalizable platform for the data-driven solution of frac-tional constitutive models without compromising accuracy as the complexity of the underlying response/model increases.

Problem setup and methodology
In "Fractional Maxwell model," "Fractional Kelvin-Voigt model," and "Fractional Zener model," the fractional Maxwell, Kelvin-Voigt, and Zener models are introduced, respectively, in terms of their fractional differential equations (FDEs) and analytical solutions.These three models are also mechanistically depicted in Fig. 2.Then, in "Rheology-informed neural networks and fractional derivatives," the concepts behind rheology-informed neural networks (RhINNs) and the numerical procedure to solve FDEs are discussed.
In this work, the analysis is performed by studying the synthetically generated data using the analytical solutions (as opposed to seeking digitized experimental measurements from the literature) because by doing so, the ground truth fractional derivative orders can be tuned to evaluate the robustness in the recovery of fractional derivative orders.
The solution of Eq. 3 upon imposing a step strain ( (t) = 0 H (t), where H (t) is the Heaviside step function) can be derived analytically, and the relaxation modulus, G(t), in Pa, is expressed as follows: 123 where M κ,μ is the generalized Mittag-Leffler function defined as follows (Schiessel et al. 1995): where (.) is the gamma function.

Fractional Kelvin-Voigt model
For the case of the fractional Kelvin-Voigt model (FKVM), two Scott-Blair elements (with parameters similar to that in "Fractional Maxwell model") are stacked in parallel; see Fig. 2b.The deformation, similar to the classical Kelvin-Voigt model, is identical in each branch, and the stresses are additive.By simply adding the stress response of two Scott-Blair elements (Eq.2) and setting 1 = 2 = (t), one may derive the stress response of FKVM (Schiessel et al. 1995): where E and τ continue to be defined as in "Fractional Maxwell model".The creep compliance, J (t), in Pa −1 , is then defined as follows: where σ 0 is the step-stress in a creep experiment.By assuming a unit step-stress (σ (t) = σ 0 = 1), the solution of Eq. 6 ( (t)) and the exact solution (Eq.7) will become identical.

Fractional Zener model
The classical Zener model consists of a Maxwell model in parallel with a spring.In the fractional Zener model (FZM), three fractional elements with parameters (E 1 , τ 1 , α), (E 2 , τ 2 , β), and (E 3 , τ 3 , γ ) are arranged, as displayed in Fig. 2c.The total mechanical response of the model with constraints 0 ≤ β < α ≤ 1 and 0 ≤ γ ≤ 1, and by , and E = E 3 (τ 3 /τ ) γ is given as follows (Schiessel et al. 1995): The relaxation modulus, G(t), of FZM under a stressrelaxation test with constraints mentioned above is given by the following: Similar to FMM, a stress relaxation test under a unit step strain is assumed here, i.e., (t) = H (0 + ) = 1, and G(t) = σ (t).Thus, the solution of Eq. 8 and the analytical solution (Eq.9) can be compared side-by-side.
It is worth mentioning that the range of material response using the analytical solutions, i.e., Eqs. 4, 7, and 9, highly depends on the fractional derivative orders (and quasiproperties, which are assumed given here; see the outset of Fig. 3 Schematic illustration of the RhINN architecture.In this figure, the material response is assumed to be the relaxation modulus, G(t), which is applicable to the fractional Maxwell and fractional Zener models.For the fractional Kelvin-Voigt model, this output is replaced with the creep compliance, J (t).The fitting parameters (Params) in RhINNs, which are the fractional derivative orders, are defined trainbale "Results and discussion").Consequently, to keep the discussion concise, readers are encouraged to refer to the material response plots provided in "Results and discussion" for more quantitative information relevant to the respective data sets.

Rheology-informed neural networks and fractional derivatives
For all three cases introduced above, the only input to the neural network (NN) is time, t, while the NN output is the material response, which is σ (t) (or G(t)) for the FMM and FZM cases, and (t) (or J (t)) for the FKVM; see Fig. 3.The neural network is composed of several hidden layers, each containing a predefined number of neurons.The trainable variables, i.e., the weights and biases of the neurons plus the fitting parameters (Params in Fig. 3), which are the fractional derivative orders in each case, are learned by minimizing a composite loss function: where φ d is the loss due to the discrepancy between the predicted material response (G p4 ) and the ground truth material response (G gt ), defined as the following mean-squared error (MSE): Moreover, the residual (Res in Fig. 3) between the two sides of each FDE (φ f ) can be calculated, which is also a mean squared error.As an example, this residual for the case of FMM will be as follows: where n r is the number of residual points that are used to define artificial input arrays in time.In writing the second equality in Eq. 12, a unit step-strain function is assumed for i , where the α th fractional derivative of a constant function y(t) = 1 is known to be t −α / (1−α) (Garrappa et al. 2019).
The objective in an inverse problem, as stated in "Introduction," is to recover the fitting parameters of an embedded constitutive equation, which are the three FDEs introduced in "Fractional Maxwell model," "Fractional Kelvin-Voigt model," and "Fractional Zener model."In a RhINN implementation, these FDEs are interrogated at discrete points in time (n r in Eq. 12).To calculate the residual loss component (φ f ), it is thus necessary to adhere to a definition of fractional derivatives best suited for the problem at hand.Among the existing definitions (or senses) of fractional derivatives (Baleanu 2016), the Riemann-Liouville (RL) and Caputo definitions seem to receive more attention (Jiang and Zhang 2020).Despite the equivalency of these two senses for certain viscoelastic cases (e.g., the fractional Zener model (Bagley 2007)), the Caputo solution is the suitable choice when (integer-order derivative) information about the boundary and/or initial conditions is available (Mainardi 2010;Li et 2011), which is the case here: the zeroth-order derivative of G(t) and J (t) curves is given, and thus, the initial conditions are accessible5 .Thus, the Caputo sense is employed in this work to handle the fractional derivative of the outputs (σ (t) or (t)) with respect to the input, t.The α th order (0 ≤ α ≤ 1) fractional derivative of f (t) w.r.t.time in the Caputo sense is defined as follows (Baleanu et al. 2016): In writing Eq. 13, the lower limit of integration is set to t = 0.
It is worth mentioning that fractional derivative orders higher than unity are calculable using the Caputo sense; however, fractional derivatives are employed for viscoelastic materials here, and thus, the derivative order is deemed to remain within zero and one.Equation 13is not readily implementable on NNs, and thus, one will need to approximate the integral term.In this work, a well-established finite-difference algorithm (Diethelm et al. 2005) is used to calculate the Caputo derivative of f (t): where h is the (uniform) step size (h = t/n r ), f 0 is the value of f (t) at t = 0, and a n,n r is the quadrature weights derived from a product trapezoidal rule: which leads to an error of order h 2−α .The hereditary nature of the Caputo fractional derivative lies in the fact that the derivative at each point depends on the past information of f (t) in time.
In this study, n r = 100 artificial points are uniformly distributed between 0.001 and 5s in time to calculate the residual loss, Eq. 12, and train RhINNs.Further increasing n r resulted in negligible improvement in parameter recovery.One justification for selecting a small training time window is that the chosen numerical approach for the fractional derivative (Eq.14) becomes unstable for non-uniform grids.Therefore, by increasing the time window, fewer points are available at the beginning of time, which can compromise the parameter recovery.The exact number of points is used to generate the ground truth data (G(t) or J (t)), which are then fed into the RhINNs.
One crucial step in solving inverse problems using RhINNs is the selection of NN hyperparameters, i.e., the NN parameters that are selected (and optimized) before the training process.For all three FDEs, four hidden layers, each containing ten neurons with a tanh activation function, are employed.Insignificant variation in the achieved convergence (and the execution time) was observed for other values of neuron or layer count.The trainable fractional derivative orders to be recovered are constrained to remain between 0.01 and 0.99 for numerical stability and also to respect the assumptions made to derive Eq. 14.A piecewise learning rate was employed that gradually decays in time to ensure gradient stability and also to find a compromise between the accuracy and the code execution time.Rather than a hard error threshold, the training process is halted once the recovered parameters (or the composite loss, φ) plateaus and does not change.Each FDE inherits a particular degree of complexity, and a hard loss threshold that can effectively terminate the training process of one FDE might not apply to the other two.The initial conditions for the fitting parameters are also tabulated in Table 1.The neural network can tolerate initial conditions different than the ones tabulated below, but for unphysical initial conditions (e.g., α < β in the FMM cases), there may be convergence issues.The loss minimization task is handled by TensorFlow's builtin tf.keras.optimizers.Adam optimizer.The entire optimization task can also be performed using SPSVERBc1's L-BFGS-B method in our code, but since this method is more computationally expensive than Adam, it is more convenient to perform Adam and use L-BFGS-B if need be.

Results and discussion
For all the FDE cases studied here, the compound elastic modulus (E) and relaxation (τ ) are assumed given and were set to their respective values used when generating the exact data.The inevitability of this bold assumption lies in the fact that the (unique) recovery of fractional derivative orders, along with all the properties of the Scott-Blair elements, was infeasible for some cases.Nonetheless, assuming missing values of E and τ , the compound elastic modulus is accessible by taking the initial slope of the stress (or the relaxation modulus in a step-strain stress relaxation test) response in time.The relaxation time, in a similar fashion, is derivable by extracting the time needed for the material to relax the elastic response, which is the time range before the stress (or G(t)) plateau.Admittedly, these pieces of information may demand rounds of experiments.However, it is safe to assume that these parameters can be at least estimated with an acceptable level of confidence.Also, it is worth mentioning that our preliminary work suggests that all parameters can be recovered through RhINN solutions in the majority of the cases studied, but further work on that front is required to provide a robust and generalizable framework.As such, here and in this work, we limit the study to the recovery of fractional derivative orders.These fractional derivative orders are not readily accessible from rheometry, and RhINN is interrogated to precisely pinpoint these derivative orders.

Fractional Maxwell model
For the fractional Maxwell model, as mentioned in "Fractional Maxwell model," RhINN is interrogated to recover the fractional derivative orders of the elements, i.e., α and β, as shown in Fig. 2a, under a stress relaxation test with a unit step-strain ( (t) = H (t)).The neural network is trained on a broad range of α and β values with 0 ≤ β < α ≤ 1.Moreover, the other two model parameters, i.e., E = 1Pa and τ = 20s, are assumed given; see the explanation in "Results and discussion."In Fig. 4, the total relative absolute error (RAE) (in %) between the RhINN predictions and the ground truth values for a wide range of α and β values is displayed.In each tile, the top and bottom numbers represent the recovered α and β, respectively, which can be cross-checked with their corresponding exact values given for each row and column.The maximum total RAE for the entire parameter space of α and β is less than 10%, which indicates that RhINN is capable of recovering the fractional derivative order of FMMs.
From a rheologist's vantage point, having a precise fractional derivative order recovery might seem secondary to having a measurable rheological property, which is the relaxation modulus, G(t), in the case of a stress relaxation test.That is, the ultimate goal in constitutive modeling is enabling accurate prediction of material response, and as such, the Fig. 4 Total relative absolute error (RAE) in fractional derivative order values (α and β) compared to their respective exact values for the fractional Maxwell model under a stress relaxation test with a unit stepstrain function.The numbers in each tile, from top to bottom, represent the recovered α and β values, respectively, with their exact values given by their corresponding row and column.The upper triangle is inaccessible since in deriving Eq. 3, it is assumed that 0 ≤ β < α ≤ 1 absolute value of the parameters involved in a model may not necessarily mean as much.Since the ground truth formula (Eq.4) to generate the data is available, it is possible to calculate the recovered relaxation modulus (G(t)) by inserting the recovered α and β values along with the given values of E and τ and then comparing the recovered G(t) with the ground truth one.Such material responses are plotted in Fig. 5 for four combinations of α and β, with the RhINN predictions and the exact G(t) values shown in solid and dashed lines, respectively.The plots in Fig. 5 indicate that even for the case of the largest RAE reported, recovered relaxation spectra is closely tracking the ground truth solution.Recalling the definition of Scott-Blair elements (Eq.2), the lower value of the fractional derivative order (β) is a manifestation of the high-frequency material response, which is the elastic response in viscoelastic materials.Therefore, changing β should be observable in the initial G(t) response in Fig. 5. Here, The relaxation modulus is lower for the smaller values of β since G(t) is proportional to t −β , which decreases as β drops; see Eq. 4.

Fractional Kelvin-Voigt model
RhINNs are then applied to recover fractional derivative orders of the elements in the fractional Kelvin-Voigt model (FKVM) under a creep test using a unit step-stress function for various (α, β) combinations with constraint 0 ≤ β < α ≤ 1.In each case, the other model parameters, E = 1Pa Fig. 5 The relaxation modulus, G(t) of the fractional Maxwell model.The RhINN (solid line) response is calculated inserting the recovered sets of α and β into Eq.4, plotting the material response for a time window equal to that of the training data (0.001 and 5s), and is compared with the ground truth (dashed line) relaxation modulus.The elastic modulus (E = 1Pa) and relaxation time (τ = 20s) are assumed given and τ = 20s, are assumed given.The total RAE between RhINN predictions and the ground truth values for different fractional derivative orders are shown in Fig. 6.The RAE values are mostly within a 20% band, with the exception of some scenarios where the overall relative absolute inaccuracy reaches up to 50%.Note that in those examples, the small value of β results in a large RAE.The higher RAE scenarios were further investigated by examining the recovery of each particular fractional derivative order (as represented by the values in the tiles).It was found out that the higher order, α, is consistently recovered more accurately than β.One possible reason for higher deviation in β recovery could be that RhINN was trained over a short time window, and the higher fractional derivative order (or more viscous element) dominates the material response in the Kelvin-Voigt model within short timescales.
As previously presented for the case of FMMs, and to support our interpretation, the material response is also plotted using the predicted and ground truth fractional derivative orders.In Fig. 7, the creep compliance (J (t)) is depicted for cases where there are significant differences between the predicted and ground truth fractional derivative orders.In this figure, the recovered derivative orders are inserted into Eq.7 over a larger time window, i.e., 1 × 10 −2 ≤ t/τ ≤ 1 × 10 4 to monitor the out-of-sample response of J (t) using the recovered parameters.Across the range of derivative orders, the predicted and the ground truth responses are in good agreement at shorter times, e.g., t/τ < 10.Despite the high RAE in recovering α and β at smaller values of α (the upper left cases in Fig. 6), the J (t) response nonetheless commits a much smaller error, indicating that the combination of αandβ determines the accuracy of J (t) and not their individual values.However, disparity between the predicted and ground truth J (t) responses arises at longer times due to the relatively substantial error in recovering β.This can also be alleviated by increasing the training time window, which was set to be within 0.001 and 5 s throughout this study.Fig. 7 The creep compliance, J (t), response of the fractional Kelvin-Voigt model using the recovered (solid) and exact (dashed) fractional derivative orders.Other parameters, E = 1Pa and τ = 20s, are assumed known.The recovered derivative orders are inserted into Eq.7 over a larger time window, i.e., 1 × 10 −2 ≤ t/τ ≤ 1 × 10 4 to assess the out-of-sample response of J (t) using the recovered parameters

Fractional Zener model
Finally, RhINN is used in this case to recover the fractional derivative orders (α, and γ ) of the three Scott-Blair elements in the fractional Zener model (FZM) under a stress relaxation test using a unit step-strain function for a wide range of derivative order values with constraints 0 ≤ β < α ≤ 1 and 0 ≤ γ ≤ 1.Other model parameters (E 0 = E = 1Pa and τ = 20s) remained constant, as before.The NN hyperparameters are also similar to those used in "Fractional Maxwell model" and "Fractional Kelvin-Voigt model." Figure 8 summarizes the total relative absolute error in recovering the three fractional derivative orders.The RAE values are generally within a 30% band for the derivative orders, with the exception of a few scenarios where the relative absolute error hits 100%.In Fig. 8 b and c, these higher-error cases occur when γ is equal to one of the other two derivative orders (β and α, respectively).On the other hand, it should be noted that some recovered fractional derivative orders match their respective exact values.Thus, the true test of the RhINN order recovery can be done by comparing the relaxation moduli using the recovered parameters and benchmarking against the ground truth solution.
In Fig. 9, the recovered derivative orders are inserted into Eq.9 over a larger time window, i.e., 1 × 10 −3 ≤ t/τ ≤ 1 × 10 4 to assess the out-of-sample response of the material response using the recovered derivative orders.The relaxation modulus, G(t), using both the recovered and exact fractional derivative orders, is shown in Fig. 9 for the cases in which the highest values of RAE were observed.Despite the significant relative absolute inaccuracy in the recovered fractional derivative orders, the predicted relaxation moduli closely track the exact values for all scenarios at shorter times t/τ < 1.At longer times (unseen by the RhINN), however, a significant discrepancy between the predicted and actual responses is observed for these cases.Once again, this can be explained by recalling the fact that the training time window (0.001 and 5s) was not large enough to encompass all the intricacies of the material response at longer times.The selection of the time window was also constrained by the computational expenses of fractional derivatives.We remind the reader that in this particular algorithm, time has to strictly be spaced linearly, and since small time intervals are required at short times, this prevents one from a comprehensive and long time training process.Therefore, employing more versatile and efficient fractional derivative algorithms compatible with nonlinear step sizes should be at the core of future efforts on this topic.It should be stressed, though, that out-of-sample prediction using RhINNs (or any other numerical toolkit) is accurate only up to a threshold, and the numerical solver usually cannot stray too far from the data that were given during the training stage.9 The relaxation modulus, G(t), of the fractional Zener model using the recovered (solid) and exact (dashed) fractional derivative orders.Other parameters (E 0 = E = 1Pa and τ = 20s) are assumed known and constant.The recovered derivative orders are inserted into Eq.9 over a larger time window, i.e., 1 × 10 −3 ≤ t/τ ≤ 1 × 104 to monitor the out-of-sample response of G(t) using the recovered parameters

Conclusion
In this work, RhINNs were developed and employed to recover the fractional derivative orders of the elements in fractional viscoelastic Maxwell, Kelvin-Voigt, and Zener models based on single sets of (synthetically generated) data.Moreover, the derivative orders were used to generate the stress relaxation response of the fractional Maxwell and Zener models upon the imposition of a unit step-strain and the creep compliance response of the fractional Kelvin-Voigt model under a creep test using a unit step-stress function.A systematic study on the recovery of the derivative orders over their permitted parameter space (between zero and unity) was performed.With various combinations of the fractional derivative orders (and their respective constraints), multiple scenarios were examined in each of the three fractional models.Overall, in most situations, RhINNs are found to retrieve the fractional derivative orders of all the components in each of the models in the same range as the ground truth values used in the generation of the data.In other instances, the total relative absolute error (RAE) for some of the derivative orders is found to be substantial.However, the relaxation modulus (or creep compliance) using the recovered derivative orders closely mimics the input experiments with minimal deviations.This is mainly because the solution of these models to an applied flow protocol is not unique within the time window under question, meaning that several combinations of the derivative orders can replicate the same overall behavior.This is more noticeable when more complex models (with more spring-pot elements) are considered.As such, for the derivative order recovery of the fractional Zener model, the recovered parameters do not necessarily follow the ones used to generate the data; nonetheless, the retrieved order values do indeed result in an accurate recovery of the behavior.Also, in this work, the neural network was trained over a short time window in terms of the data provided, and in some cases, those time windows are likely not long enough for all model components to respond appropriately to an applied stimulus.Moreover, due to the hereditary nature of fractional derivatives and the need to sum all the past information to get the fractional derivative at each point, the RhINN convergence was found to be time-consuming (typically within tens of minutes, or an hour), which leaves room for further improvement in the computation time and efficiency of fractional problems using RhINNs.While this work presents a robust platform for the inverse data-driven solution of fractional constitutive models, further investigation is warranted on recovering all model parameters (including elastic moduli and relaxation times) and other flow protocols.Also, training the RhINNs over a longer time window could improve the efficiency of fractional derivative recovery.

Fig. 2
Fig. 2 Schematic description of the fractional a Maxwell, b Kelvin-Voigt, and c Zener models using spring-pot elements.The objective of this work is to recover the fractional derivative orders, i.e., α and β in fractional Maxwell and Kelvin-Voigt models and α, β, and γ in the fractional Zener model

Fig. 6
Fig. 6 Total relative absolute error (RAE) in the fractional derivative orders (α and β) with respect to their exact values in the fractional Kelvin-Voigt model under creep compliance for various combinations of fractional derivative orders with constraint 0 ≤ β < α ≤ 1.The numbers in each tile, from top to bottom, represent the recovered α and β values, respectively, with their exact values given by their corresponding row and column

Fig.
Fig.9The relaxation modulus, G(t), of the fractional Zener model using the recovered (solid) and exact (dashed) fractional derivative orders.Other parameters (E 0 = E = 1Pa and τ = 20s) are assumed known and constant.The recovered derivative orders are inserted into Eq.9 over a larger time window, i.e., 1 × 10 −3 ≤ t/τ ≤ 1 × 104 to monitor the out-of-sample response of G(t) using the recovered parameters

Fig. 8
Fig. 8 Total relative absolute error (RAE) in the recovered fractional derivative orders (α, β, and γ ) with respect to their exact values for fractional Zener model under a stress relaxation test using a unit step-strain function for various combinations of derivative orders with

Table 1
The initial conditions (ICs) for the fractional derivative orders (fitting parameters) The third Scott-Blair element is present only for the fractional Zener model (FZM).The NN can withstand initial conditions other than the ones tabulated here; however, for unphysical ICs (e.g., α < β for the FMM case), convergence issues may occur