Model-order reduction of lumped parameter systems via fractional calculus
This study investigates the use of fractional order differential models to simulate the dynamic response of non-homogeneous discrete systems and to achieve efficient and accurate model order reduction. The traditional integer order approach to the simulation of non-homogeneous systems dictates the use of numerical solutions and often imposes stringent compromises between accuracy and computational performance. Fractional calculus provides an alternative approach where complex dynamical systems can be modeled with compact fractional equations that not only can still guarantee analytical solutions, but can also enable high levels of order reduction without compromising on accuracy. Different approaches are explored in order to transform the integer order model into a reduced order fractional model able to match the dynamic response of the initial system. Analytical and numerical results show that, under certain conditions, an exact match is possible and the resulting fractional differential models have both a complex and frequency-dependent order of the differential operator. The implications of this type of approach for both model order reduction and model synthesis are discussed.
Numerical accuracy and computational efficiency have been long-standing challenges in the simulation of dynamical systems. Numerical solutions of complex continuous systems with non-trivial boundary and loading conditions typically require a discretization process to obtain lumped parameter models. The more complex the property spatial distribution (e.g. external loads, material or geometric parameters, boundary conditions), the higher the level of discretization needed to achieve a satisfactory representation of the original continuum system. The increased discretization directly impacts the computational performance and, for large systems, limits the level of achievable accuracy. This issue is even more accentuated when dealing with active control or real-time prediction of the dynamic response of systems under operating conditions for which fast state estimation is a key requirement.
Over the last several decades, these challenges have motivated the rapid growth and development of methodologies for the synthesis of computationally efficient approaches able to reduce the overall size of the models (i.e. of the total number of Degrees Of Freedom - DOF) while maintaining high numerical accuracy and fidelity to the actual dynamics. These techniques, referred to as model order reduction, are typically pursued when the dynamic response is sought only at selected locations (the so-called active DOFs) such that the DOFs associated with the remaining locations can be omitted. The reduction procedure is not trivial because it must account for the coupling between active and omitted DOFs in order to not change the underlying dynamics of the system.
Many sources in the literature [1, 2, 3] provide an extensive review of reduction techniques across a variety of disciplines. In the following, we concentrate only on applications to structural dynamics since this is the emphasis of our study. One widely used technique is Guyan reduction. This approach is also known as static reduction because it does not account for the system inertia and is therefore limited to statics. Reduction techniques for dynamic systems (therefore accounting for the system’s inertia) are often based on mode superposition or component mode synthesis . Perhaps one of the most widely used component mode synthesis technique, also known as the Craig-Bampton method [4, 5]. The Craig-Bampton method divides the system into several substructures which are required to be compatible along their shared boundaries. Assuming these boundaries are held fixed, the Craig-Bampton method is able to combine the motion of these boundary points with the displacement modes of the substructures (known as constraint modes). The dynamics of the system can be reduced to a set of both fixed-interface and constraint modes .
Many existing order reduction techniques can only provide an approximation of the local response. For instance, the accuracy of the approximation of the Craig-Bampton method is strictly dependent on the number of modes retained in the modal basis. The truncation of the basis should be assessed with respect to either the modal densities associated with the omitted modes or the dynamic content that should be transferred to the active degrees. In a similar way, enlarging the modal basis (i.e. extending the truncation order) comes at the expenses of computational performance.
To address these limitations, we explored a reduction order technique based on fractional calculus. We will show that while fractional models contain less DOFs than the original system, often times the dynamic response can be matched exactly at the active degrees. This is an important advantage of our fractional order reduction over reduction order techniques typically used. In addition to the order reduction capabilities, we anticipate that the proposed approach has possible applications to the system identification based only on measured or experimental data. As an example of these capabilities, we will show the application of the fractional models to perform broadband system identification of discrete parameter systems.
While the mathematics of fractional calculus has been extensively studied in the past century, applications are relatively recent. In particular, fractional calculus has seen applications in engineering areas such as controls, visco- and thermo-elasticity [7, 8, 9], and wave propagation in complex media [10, 11, 12]. The reader is referred to [13, 14, 15] for detailed reviews on the fundamentals of fractional calculus. In §6 we define some basic fractional calculus quantities. We anticipate that the methodology discussed below produces fractional differential models of complex and frequency-dependent order. Love , Ortigueira , Ross , and Andriambololona  among others have developed the mathematics of complex fractional derivatives as well as potential uses. Adams  and Neamaty  have developed solutions to certain complex fractional order differential equations. Other authors including Atanackovic , Makris , and Park  have successfully applied complex fractional calculus to viscoelasticity. While the mathematics of complex fractional calculus has been explored and developed, its connection to the actual physical processes being represented can be difficult to grasp. Perhaps Makris has one of the clearest interpretations of a complex order derivative: ”…one may interpret the complex derivative of an arbitrary function as the superposition of complex derivatives of harmonic functions. Evidently, complex-order derivatives modulate the phase and amplitude of harmonic components of a time-dependent function in a more complicated way than real-order derivatives. An important difference between real-valued and complex-valued time derivatives is that phase modulation in the latter case is frequency dependent whereas in the former is not” . As pointed out here by Makris, complex fractional derivatives can produce functions whose amplitude and phase are both frequency-dependent; this attribute plays a key role in the development of the fractional models. Valerio [25, 26] gives a thorough review of complex and variable-order fractional derivatives through the use of Laplace transforms and transfer functions.
The remainder of the paper is structured as follows: in the next section, we present how to obtain an undamped fractional single degree of freedom model having the same dynamic response of an integer damped single degree of freedom. While this does not technically qualify as order reduction, it illustrates the basic methodology that will be used throughout this work. Then, we present the procedure to reduce a multiple integer-order degree of freedom system to a single fractional-order degree of freedom system. This represents the fundamental step of the order reduction methodology. Next, we extend the methodology to reduce a M-degree of freedom integer-order system to a N-degree of freedom fractional-order model, where and . Lastly, we briefly discuss how the same methodology can be used in the frame of system identification and we show how to synthesize accurate fractional dynamic models based only on the knowledge of numerically obtained or experimentally measured data.
An important remark should be made concerning the terminology used in this paper. It is well-known that, for integer order systems, the overall order of the system is strictly related to the number of degrees of freedom because each individual degree is assumed to behave as a second order system. Therefore, order and dimension are typically considered as equivalent concepts and used interchangeably. On the contrary, due to the infinite dimensional character of a fractional derivative, the connection between the overall order of a system and the number of its physical degrees of freedom is somewhat more ambiguous. Such a discussion goes well beyond the scope of this paper and we simply highlight that, in the remainder, the term degrees of freedom will refer to the number of discrete masses while the term order will refer to the order of the individual differential equations.
2 Reduction to SDOF fractional systems
This section discusses the fundamental approach to obtain a fractional undamped single degree of freedom (SDOF) model exhibiting a dynamic behavior equivalent to an integer damped single degree of freedom system. To facilitate the understanding of the procedure, we first discuss the conversion from an integer-order damped SDOF (I-SDOF) to a fractional-order SDOF (F-SDOF). Then, we will extend the methodology to the case of an integer-order multiple degree of freedom (I-MDOF) to be reduced to a F-SDOF.
2.1 The fractional order SDOF model
Consider the fractional SDOF oscillator in Figure 1b. The equation of motion of the fractional model is given by:
where is the displacement from the equilibrium condition, is the mass of the fractional model, is the stiffness, is the time dependent load acting on the mass, and is the frequency-dependent order of the fractional derivative. In Equation (1), the fractional derivative should be intended in the Caputo form since this fractional derivative best represents an initial value problem of a physical system. However, it should be noted that the Riemann-Liouville derivative is equivalent to the Caputo derivative if all initial values are zero (which is also the general case considered in this paper). When 0 1, Equation (1) is known as the fractional relaxation equation [13, 14], while when 1 2, Equation (1) is known as the fractional oscillation equation [13, 14]. Our model will use the form of Equation (1), but will not necessarily be restricted to 0 2. Also, our model will use a complex value of the fractional order = , where and are the real and imaginary part of , respectively.
Taking the Laplace transform of Equation (1) with zero initial conditions, we obtain:
where is the Laplace variable. The transfer function of the F-SDOF is:
Substituting and , Equation (3) becomes:
Through some algebraic manipulation, Equation (4) can be shown to be equivalent to:
In terms of sines and cosines, Equation (6) can be written as:
Equation (5) can now be rewritten in terms of its real and imaginary parts:
As a result, the magnitude and phase of the transfer function are:
If the forcing function is given by a harmonic load , then the steady state response of the F-SDOF is given by:
2.2 Exact conversion from a damped I-SDOF to a F-SDOF
The general process to obtain the fractional derivative value from a non-homogeneous integer order system is first illustrated by converting an I-SDOF into a F-SDOF having an equivalent dynamic response.
Consider the differential equation of a classical damped oscillator (Figure 1a):
where is the mass, is the damping coefficient, is the spring stiffness, is the displacement of the mass from equilibrium, and is the external force.
Assuming the forcing to be the same in the fractional model as in the integer model (that is, = ), we can set the Laplace transform of Equation (14) to be equal to the Laplace transform of the fractional model. Letting the displacements of the I-SDOF and the F-SDOF be equivalent (that is, = ), we obtain a polynomial equation of fractional order:
Since we are converting a SDOF to another SDOF, we set the coefficients of the fractional model equal to their integer model counterparts; that is, = and = . Substituting this into Equation (15) and solving for results in:
Substituting s = i into Equation (16) yields:
where is the forcing frequency and Ln is the complex logarithm function. Recall that the complex logarithm function Ln in Equation (17) is defined as:
where is the modulus () and Arg is the principle argument (Arg).
Equation (17) allows an important observation: the order of the fractional model is both a complex and frequency-dependent quantity. Note that both characteristics are due to the fact that the damping term of the integer order model is now included into the fractional derivative; this is a well-known property of fractional oscillators [14, 27, 28, 29, 30]. As expected, setting would return an integer order derivative.
The above results suggest that if in Equation (1) is chosen according to Equation (17), then the steady state response of the two systems are exactly equivalent. To confirm this, we calculated numerically the response of the two systems. The steady state response of the I-SDOF subject to a force is given by  while the steady state of the fractional model is given by Equation (13). Consider an I-SDOF where = 2, = 1, and = 10 (unit system is arbitrary, although we shall use units of seconds for time so to conveniently refer to frequency in terms of rad/s or Hz). As already discussed, we then let = = 2 and = = 10 in the fractional model. The value of can then be calculated using Equation 17 for a specified forcing frequency. The steady state response of the fractional model can be obtained using the magnitude and phase as given by Equations 11 and 12. As an example, for a harmonic excitation of frequency = 10 rad/s, the value of is calculated as 1.9903 - 0.0151i. Figure 2a depicts the steady state response of both the integer order model and its equivalent fractional model, each obtained from their associated transfer functions, for = 10 rad/s. Both models return an equivalent steady state response, as expected due to the exact correspondence of the analytical transfer functions. Figure 2b depicts the magnitude and phase of the I-SDOF and F-SDOF models across a range of frequencies. Both sets of magnitude and phase show that the transfer function of the F-SDOF model matches exactly the transfer function of the corresponding I-SDOF. It should be pointed out that the phase of the fractional models is actually the negative value of the phase of the integer models as evidenced by the previously derived equations of the steady state response.
Figure 2c shows a plot of the value of as a function of . The natural frequency of both the I-SDOF and the F-SDOF in this numerical example is . Figure 2c shows that for excitation frequencies above the critical frequency, the real part of tends towards 2 while the imaginary part of tends toward 0. Thus, when the forcing frequency is well above the resonance frequency (i.e. the high-frequency asymptotic limit), the fractional model converges to an integer order model. On the contrary, for frequencies below the critical frequency, the model order is fractional and as expected due to the presence of viscous damping in the initial integer order model.
2.3 Model order reduction from I-MDOF to F-SDOF
Considering the general procedure discussed above, we can discuss the case of model order reduction from a multiple integer-order degree of freedom (I-MDOF) system to a fractional SDOF. This is an extreme form of reduction where the initial MDOF system is collapsed to a single DOF model. This study case is chosen to show the flexibility of the fractional model approach and the ability to yield highly accurate, or even exact, solutions over an extended frequency range.
Assume that the response of one of the degrees of freedom in the I-MDOF (Figure 3a) model is the active degree of interest and that all the remaining DOFs are omitted. We seek an equivalent F-SDOF representation (Figure 3b) of the system such that its response matches exactly or approximately the response of the corresponding active degree of freedom of the I-MDOF model. Recall that the transfer function of the fractional SDOF is given by Equation 3. In order to obtain the order of the F-SDOF, we first convert the I-MDOF to a state-space form. From the state-space form, the transfer function of any of the nodes in the integer MDOF can be obtained. The process of obtaining the transfer function from state-space is well-established [32, 33] and it is also briefly reviewed in §7.
Let be the transfer function for the active degree of freedom in the I-MDOF model, obtained from the state-space form. Assume that the corresponding fractional model has a mass of and a stiffness of , as reflected in Equation 1 and Figure 3b. By equating the transfer function of the fractional SDOF to the transfer function of the degree of interest in the I-MDOF, we obtain the fractional order :
where is written as . While we obtained by converting from state space to a transfer function, it should be noted that Equation 19 still holds true if the transfer function was obtained via another method. This latter characteristic will be the basis for the model synthesis approach discussed in §4. Equation 19 can be made a function of frequency by substituting . Note that Equation 19 provides a fractional order which guarantees an exact match of the transfer functions of the two systems. Therefore, as far as the individual functions and are exact, the response of the equivalent F-SDOF is an exact match of the initial I-MDOF.
Contrarily to the I-SDOF case, here, the equivalent mass and stiffness of the F-SDOF need to be determined. In fact, while in §2.2 the selection of the equivalent parameters and was straightforward given the existence of only one set of parameters, in the current configuration, multiple choices can be made. Among the possible approaches, we decided to set so that the total mass of the I-MDOF matches that of the F-SDOF:
where is the mass of the -th degree of freedom of the I-MDOF. Since the I-MDOF model consists of springs in series, we set equal to the equivalent stiffness of springs in series:
As in the previous case, we numerically test the methodology. We consider an I-MDOF with and non-uniform coefficients. Specifically, we consider = = 1, = = 2, = = 1, = = 2, = =1, and = = 2. The natural frequencies are all within the range 0.36 rad/s to 2.22 rad/s. Let the dynamic response of the first mass ( in Figure 3a) be the active degree and therefore the one whose response is reduced to a fractional SDOF. Using Equations 20 and 21, the parameters of the fractional model are taken as = 6 and = . Assuming that we have obtained the transfer function for the first mass in the integer MDOF model, Equation 19 can be applied to find for a desired forcing frequency. The steady state response of the fractional model can be obtained using the magnitude and phase as given by Equations 11 and 12. For a harmonic forcing function of frequency = 1 rad/s, the value of is calculated as 1.3807 + 0.7731i. Figure 4a shows the steady state response of both the first mass of the integer order model and its equivalent fractional model, each obtained from their associated transfer functions. As expected, the response of the equivalent fractional model is identical to the response of the first mass of the I-MDOF. Figure 4b depicts the match between the magnitude and phase of the first mass of the I-MDOF and the magnitude and phase of the F-SDOF.
Recall that the transfer function is a function of the forcing frequency ; thus Equation 19 is a function of . Figure 4c shows a plot of the value of as a function of for our example system. The trend of is highly dependent on the methodology used for determining the variables and . Following our suggested methodology (Equations 20 and 21) results in a trend of that closely follows that observed in §2.3. In the asymptotic limit, the fractional system tends to a second order system () with no phase modulation (); that is, no frequency dependence of the phase. Thus, in the asymptotic limit, the phase of the response becomes nearly constant and no longer changes with frequency. In the low frequency limit, we observe a fractional oscillator behavior () with phase modulation (). Note that the original system is of order eight with resonances clustered in a narrow frequency range. For this reason, in the frequency range of the local resonances, the equivalent fractional system can exhibit order . As highlighted above, the fractional order is affected by the choice made to define the equivalent parameters. Results also confirm that non-monotonic changes in the fractional order should be expected in between resonances, suggesting that the evolution from one local resonance to the next (typically modeled in conventional dynamic theory as a SDOF second order oscillator) occurs via a dynamic behavior that is locally fractional.
We conclude this section noting the two main advantages of the technique presented above: 1) the approach allows a remarkable reduction in order without any loss in the characteristic features of the original dynamics, and 2) the resulting F-SDOF can pave the way to the use of analytical solutions (e.g. Mittag-Leffler function ) for the simulation of complex dynamic systems that typically allow only numerical approaches.
3 Reduction to NDOF fractional models
In the previous section, we successfully reduced an I-MDOF to a F-SDOF while exactly retaining the dynamic response of the degree of interest in the initial integer order model. We extend this approach in order to make the reduction procedure more general and account for the reduction of an I-MDOF to a fractional multiple degree of freedom system (F-NDOF) with a number of DOFs , such as the models depicted in Figure 5. This approach allows reducing the original I-MDOF model while retaining multiple active degrees of freedom.
3.1 NDOF fractional Model
We consider how to obtain the transfer functions for a F-NDOF from the equations of motion. The equations of motion of each degree in the F-NDOF are actually very similar to those of a I-MDOF, except for the fact that the inertia terms () are replaced by fractional terms (). Recall that the selected fractional models in our approach do not include any explicit damping term. Consider a F-2DOF as in Figure 5b. The equations of motion are given by:
We could arrange these equations in a state-space form and then obtain the transfer functions in an analogous manner to what was done for the I-MDOF. However, rather than using a fractional state-space to transfer function method, we opt to calculate the Laplace transform of Equations 22 and 23 and then apply Cramer’s rule. Note that this approach could also be applied to obtain the transfer functions of the I-MDOF. Applying the Laplace transform yields:
In matrix form:
Using Cramer’s rule, and are given by:
The desired transfer functions can then be obtained using Equations 27 and 28. As an example, if and we can define and . The transfer functions of the output displacements relative to the input force on the first node are:
The substitution is made to convert to the frequency domain.
3.2 Order reduction from I-MDOF to F-NDOF
Leveraging the transfer function representation of the F-NDOF system provided in the previous paragraph, we can illustrate the reduction procedure from I-4DOF to F-2DOF. While, in general, both the I-MDOF’s and the F-NDOF’s degrees can have a forcing function acting on it, this example shall only include an external force on the first node in both the I-4DOF and in the F-2DOF. In the reduction process, let the 1st and 3rd masses in the I-4DOF be the active degrees. Therefore, the reduction technique shall yield steady state responses of the F-2DOF’s degrees which are equivalent to the steady state responses of the active degrees in the I-4DOF. Assume we have obtained the transfer functions (that is, the transfer function between the displacement and the force) of active degrees in the integer order model by using state space to transfer function techniques. We refer to these transfer functions as and , respectively.
In order to match the responses of the two systems at the selected DOFs, we impose = and = where and are given by Equations 32 and 33, respectively. It is tempting to arrange Equation 25 (with =0) into a form of and match it to and solve for . However, this will not lead to the desired matching of both responses of the active DOFs. To understand why, consider the mass and stiffness parameters in the F-2DOF. In the above procedure, it was assumed that the values of , , , and were determined by lumping masses and stiffnesses from the integer model. If we set equal to and equal to , we obtain a system of two equations in one unknown (), resulting in an over-determined system. Unless , , , and are properly chosen, the resulting fractional order will not be able to guarantee the exact matching of the reduced DOFs. Although different procedures to identify the parameters can be selected, we propose the following method. Lump the masses accordingly so that the total mass of the fractional model is equal to the total mass of the integer model. In our example, we set = and = . Next, set a relationship among the stiffness values of the F-NDOF. Recall Equation 21, in which we obtained a parameter by finding the equivalent stiffness of springs in series. We define = and = , where is an unknown coupling stiffness parameter. Note that = . Equating the transfer functions of the fractional and integer order models, we obtain a set of two nonlinear equations in two unknowns ( and ). To solve for the complex order and the stiffness coupling parameter (which will also be a complex quantity), a nonlinear numerical solver can be used.
As previously done, the procedure is further illustrated using a numerical example. Assume = = 1, = = 2, = = 1, = = 2, = = 1, and = = 2. The equivalent mass coefficients in the F-NDOF are = = 3. The magnitudes and phases over a frequency range are given in Figures 6b and 7b. Figure 6b depicts the magnitude and phase of the first degree (mass) in the I-4DOF and the first degree of the F-2DOF while Figure 7b show the magnitude and phase of the third degree in the I-4DOF and the second degree of the F-2NDOF. Results illustrate that the dynamic response of the F-NDOF was able to exactly match the response of the active degrees of the I-MDOF.
Furthermore, we verify that the steady state responses of the fractional and integer models are equivalent for a specific forcing frequency. For a harmonic forcing function of frequency = 1 rad/s, we find = 1.5834 + 0.4983i and = 1.5896 + 1.0440i. The numerical steady state results are plotted in Figures 6a and 7a. Figure 6a shows the steady state response of the first mass in the I-4DOF and the first mass of the F-2DOF. Figure 7a depicts the steady state response of the third mass in the I-4DOF and the second mass of the F-2DOF. The perfect overlap between the two plots shows that the dynamic response of the F-NDOF exactly match the response of the active degrees of the I-MDOF.
Finally, Figure 6c shows a plot of the value of as a function of and Figure 7c plots the value of as a function of . Once again, we observe the familiar trend in Figure 6c where the while in the asymptotic limit. A similar trend is also exhibited by the coupling parameter . At frequencies well past the resonance frequencies, becomes purely real. For ranges in which the coupling parameter is complex, is representative of the stiffness in the I-4DOF while (along with the complex order ) contributed to the damping. For high frequencies, where both and are integers, damping no longer plays a significant role in the steady state dynamics.
4 System identification from numerical data
The formulation developed in the previous sections assumed that the initial I-MDOF system was provided and that its fractional form was sought. However, there are many situations of practical interest in the analysis of dynamical systems where measured experimental data is available and either the corresponding mathematical model or its coefficients are unknown; this class of problems is typically referred to as system identification. Over the years, a variety of system identification techniques have been proposed [33, 34]. One of the most common approaches for vibration problems relies on matching second-order systems to individual resonances, therefore approximating the response of the system at resonance as a second order oscillator. This approach is also at the basis of the conventional half-power bandwidth method  for damping estimation.
This approach has some important limitations. When the resonance frequencies of the MDOF system are too closely spaced, the local resonance is not well approximated by the single DOF oscillator. Also, as shown in the previous numerical results, in between resonances, the behavior of the system is typically fractional due to the coupling between two or more DOFs. This also explains why, when comparing numerical and experimental data, the largest discrepancies are often observed at frequencies off-resonance (regardless of the accuracy at resonance).
The category of system identification methods requires the selection of a dynamic model that is matched to the experimental data by properly tuning the model parameters. In this approach, the structure of the mathematical model (typically based on differential operators) is selected a priori without any detailed insight in the true physical nature of the system. We have shown above that different operating regimes would require different models to achieve an accurate representation. Fractional models offer a much more general approach to the formulation of the equations of motion because they are capable of capturing in a single mathematical model a variety of physical mechanisms that would otherwise require multiple integer order models. The well-known change in the dynamic behavior of a system when transitioning from the low to the high frequency regime is a classical example of this phenomenon. In addition, it should be considered that typically the most appropriate dynamical model to describe a complex system is not known a priori; therefore, the use of a fractional model would allow a general approach in which the system identification process is allowed to converge to the most appropriate form of the governing equations. Therefore, in some respect this approach would result in a model identification method.
The above considerations suggest that fractional order models might provide a powerful methodology for the dynamic characterization of complex systems from measured data. Here below, we illustrate the approach first for an equivalent F-SDOF and then for a F-NDOF system. Note that, in the following we will refer to measured data as the data used to synthesize the dynamical models. In practice, we generated this reference data numerically.
4.1 Fractional SDOF
According to the general approach followed above, we first consider how to determine the value of across a desired frequency range in a fractional SDOF. Assume that numerical (or experimental) data describing the dynamics of a system has been obtained at a single location or DOF. Source data will be presented in the form of Bode plots. Let the magnitude of the measured Bode plot at the desired frequency be and the phase be . Furthermore, assume that the total mass and stiffness of the system have been measured. Note that an estimation of the system’s stiffness could be obtained from the amplitude of the transfer function in the low frequency range (i.e. from the quasi-static limit). The mass and stiffness of the F-SDOF system are and , respectively. Using these quantities in Equation 11 and 12 we can solve for and at a specific frequency. Rearranging these equations yields:
Using the identity 1+tan() = sec(), Equation 36 can be written as:
Taking the positive root:
Once and are obtained for a certain frequency , the nonlinear Equations 9 and 10 can be numerically solved to find the coefficients and in = . This procedure can then be repeated over a range of frequencies to obtain the variable order for the corresponding F-SDOF model.
To illustrate the procedure, assume that we have numerically (or experimentally) obtained the Bode plot given in Figure 8a. We wish to accurately represent this data using a F-SDOF model. Assume that the mass and stiffness of the system are known. In this example, we choose = 12 and = 0.4167. Equations 38 and 39 are solved to find and for 100 different frequencies between 0.01 and 100 rad/s. For each frequency, the value of can be found using Equations 9 and 10. Figure 9b shows the plot of for the selected example. To verify that the magnitude and phase of the F-SDOF match that of the numerical measured data, Equations 11 and 12 are used to plot the magnitude and phase of the obtained F-SDOF. Figure 9a shows that the magnitude and phase of the F-SDOF match exactly the initial data. In order to assess the effectiveness of the fractional model approach, Figure 9a also reports the magnitude and phase of an I-SDOF model created from the same input data. The parameters (mass, stiffness, damping) of the I-SDOF model were obtained by matching the transfer function of a second-order system at a selected resonance frequency ( = 0.3 rad/s). Note that this approach follows the traditional method used to extract physical parameters (e.g. damping) from experimental data. In fact, for sufficiently spaced resonances, the peaks in the magnitude of the transfer function are fit locally by single DOF systems. While the obtained I-SDOF provides a good match of the data near the resonance, larger discrepancies are observable off-resonance.
4.2 Fractional NDOF
Let us consider how to determine the value of at a given frequency in a fractional NDOF. Recall in the fractional SDOF, we used data from one Bode plot (magnitude and phase) to find . This concept can be extended further. Consider that we have measured Bode plots obtained at different locations or degrees of interest and we want to represent the data by a fractional NDOF model. Furthermore, assume that the total mass and stiffness of the system have been measured. The masses in the fractional NDOF are determined by appropriately discretizing the total mass of the system into blocks. The stiffnesses of the F-NDOF are not yet known. However, the model will use stiffness coupling parameters to establish relationships between the stiffnesses in the F-NDOF. There will be a total of unknowns; that is, the fractional order and stiffness coupling parameters. If the analysis is conducted at a specific frequency, there are known quantities from the experimental data; that is, the values of the magnitude and phase of each measured Bode plot. Considering the real and imaginary terms of and the stiffness coupling parameters separately, there are unknowns. Furthermore, there are equations which are obtained by equating both the measured magnitude and phase of the transfer functions from the fractional model to the magnitude and phase of the available Bode plots. The set of nonlinear equations can be solved numerically.
The procedure is illustrated on a F-2DOF model using two measured Bode plots (see Figures 8a and 8b) obtained at different locations. Assume that the total mass and stiffness are set to = 12 and = 0.4167 and that the corresponding parameters of the F-2DOF are = = 6. The unknown coupling parameter is set up such that = and = . The analysis is performed for 100 different frequencies between the values of 0.01 and 100 rad/s. The results are presented in terms of the Bode plots of the transfer functions (Figure 10) which clearly show that the F-2DOf model very closely captures the original dynamic behavior. Following an approach equivalent to the previous paragraph, we performed a system identification using also an integer order MDOF model. Figure 10 shows the magnitude and phase of an I-2DOF model created from the same input data. The model parameters were obtained by matching the transfer functions of the second-order system to the magnitude and phase of the two sets of Bode plots at a selected resonance frequency ( = 0.3 rad/s). As in the previous case, we observe that at resonance the response of the integer order system is able to capture the overall dynamics of the measured data, although without providing an exact fit. For frequencies off-resonance the discrepancies become more evident and the model fails in capturing the effective trend of the data. The order of the corresponding F-2DOF model is plotted in Figure 11a while is given in Figure 11b.
This study investigated the use of fractional order differential models to simulate the dynamic response of non-homogeneous discrete systems and to achieve efficient and accurate model order reduction. A limitation of many existing reduction order techniques, especially when applied to complex non-homogeneous systems, lies in their ability to accurately capture the local and global broadband dynamic response. We showed that the use of fractional differential equations having complex and variable order provides a viable and powerful alternative to conventional integer order models often guaranteeing, at least for the range of problems discussed in the present study, an exact match of the dynamic response. The frequency-dependent properties of the fractional operator enable frequency-dependent modulation of the phase and amplitude which are at the basis of the broad spectrum of problems that can be addressed with this type of models. Analytical and numerical results showed that the reduced fractional models can accurately (often times exactly) represent the dynamics of the active degrees of the initial system across a wide frequency spectrum. This is a key advantage of fractional models over traditional reduction techniques which typically yield narrowband performance and can guarantee only local accuracy. The proposed fractional approach was also tested in the frame of a system identification application. We argue that the use of fractional models does not only largely increase the accuracy of the parameter identification but it effectively results in a model identification approach. Numerical results showed that the fractional models were able to very accurately reflect the dynamics of the measured data, even at frequencies far off-resonance. This latter attribute is typically not achievable with conventional system identification techniques.
In conclusion, this study laid the groundwork for model order reduction and system identification techniques based on complex, variable order fractional models. The methodology was shown to accurately represent the dynamics of non-homogeneous systems over a wide frequency spectrum.
6 Appendix A: Basic Derivatives of Fractional Calculus
Several different definitions of a fractional derivative are available in the literature. Here below, we provide the two definitions that are related to this work: the Riemann-Liouville and the Caputo derivatives.
The of the Riemann-Liouville fractional derivative of order , for  is:
where is the Gamma function and =. On the other hand, the Caputo fractional definition is:
The Laplace Transform of Riemann-Liouville fractional derivative is:
while the Laplace Transform of Caputo fractional derivative is:
From Equations 42 and 43, it is evident that the Laplace transform of the Caputo derivative uses the same initial values that a typical integer order problem does (first derivative, second derivative, etc). The initial values of the Riemann-Liouville definition are actually non-integer order derivative values of the function at . The physical meaning of the necessary initial conditions using the Riemann-Liouville definition is an open question. On the other hand, the Caputo derivative lends itself to initial values which have a well-defined physical interpretation (initial position, velocity, acceleration, etc). Therefore, we shall use the Captuo fractional derivative in our analysis.
7 Appendix B: Brief Review of State Space to Transfer Function
A discrete MDOF system is represented by a set of second-order ordinary differential equations (ODE). To simplify their numerical solution, it is common to convert the original system into a system of first-order ODEs. In dynamics, this approach is often referred to as the state-space form. This can be represented by the vector equation:
where the input is and the output is:
The matrix A is the system matrix, B the input vector, C is the output vector, and D is a scalar called the direct transmission term. To illustrate, consider a 2-DOF mass-spring-damper system where there is forcing on the first block and the interested output is the displacement of the first mass. If and are the displacements of the two blocks, then the state-space is:
and = 0 (as is usually the case). The transfer function of the desired degree is [in this case, since the first entry in C is 1] and is given by:
where is the identity matrix of the appropriate size.
-  Besselink B, Tabak U, Lutowska A, van de Wouw N, Nijmeijer H, Rixen DJ, et al. 2013 A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control. Journal of Sound and Vibration 332(19), 4403 – 4422.
-  de Klerk D, Rixen DJ, Voormeeren SN. 2008 General framework for dynamic substructuring: history, review, and classification of techniques. AIAA Journal 46(5), 1169 – 1181.
-  Schilders WH, van der Vorst HA, Rommes J. 2008 Model Order Reduction: Theory, Research Aspects and Applications. Berlin: Springer.
-  Craig J RR, Bampton MCC. 1968 Coupling of substructures for dynamic analyses. AIAA Journal 6(7), 1313 – 1319.
-  Craig RR, Kurdila AJ. 2006 Fundamentals of Structural Dynamics, 2nd ed. Hoboken, NJ: John Wiley & Sons.
-  Kuether RJ, Allen MS, Hollkamp JJ. 2016 Modal substructuring of geometrically nonlinear finite-element models. AIAA Journal 54(2), 691 – 702.
-  Torvik PJ, Bagley RL. 1984 On the appearance of the fractional derivative in the behavior of real materials. Journal of Applied Mechanics, Transactions ASME 51(2), 294 – 298.
-  Wharmby AW, Bagley RL. 2013 Generalization of a theoretical basis for the application of fractional calculus to viscoelasticity. Journal of Rheology 57(5), 1429 – 1440.
-  Narahari Achar BN, Hanneken JW. 2009 Microscopic formulation of fractional calculus theory of viscoelasticity based on lattice dynamics. Physica Scripta T136, 014011.
-  Fellah ZEA, Fellah M, Depollier C. 2006 Transient wave propagation in inhomogeneous porous materials: Application of fractional derivatives. Signal Processing 86(10), 2658 – 2667.
-  Casasanta G, Garra R. 2012 Fractional calculus approach to the acoustic wave propagation with space-dependent sound speed. Signal, Image and Video Processing 6(3), 389 – 92.
-  Tarasov VE. 2016 Acoustic waves in fractal media: Non-integer dimensional spaces approach. Wave Motion 63, 18 – 22.
-  Podlubny, I. 1999 Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. New York, NY: Academic Press.
-  Herrmann R. 2014 Fractional Calculus: An Introduction for Physicists, 2nd ed. New Jersey: World Scientific Publishing Company.
-  Diethelm K. 2004 The Analysis of Fractional Differential Equations. New York, NY: Springer.
-  Love ER. 1971 Fractional Derivatives of Imaginary Order. J London Math Soc. 3(2), 241–259.
-  Ortigueira MD, Rodriguez-Germa L, Trujillo JJ. 2011 Complex Grunwald-Letnikov, Liouville, Riemann-Liouville, and Caputo derivatives for analytic functions. Communications in Nonlinear Science and Numerical Simulation 16(11), 4174 – 4182.
-  Ross B, Northover F. 1978 A Use of a derivative of complex order in the fractional calculus. Indian J Pure Appl Math. 9(4), 400–406.
-  Andriambololona R, Tokiniaina R, Rakotoson H. 2012 Definitions of Complex Order Integrals and Complex Order Derivatives Using Operator Approach. International Journal of Latest Research in Science and Technology 1(4), 317–323.
-  Adams JL, Hartley TT, Adams LI. 2010 A solution to the fundamental linear complex-order differential equation. Advances in Engineering Software 41(1), 70 – 74.
-  Neamaty A, Yadollahzadeh M, Darzi R. 2015 On Fractional Differential Equation with Complex Order. Progress in Fractional Differentiation and Applications 1(3), 223–227.
-  Atanackovi TM, Konjik S, Pilipovi S, Zorica D. 2016 Complex order fractional derivatives in viscoelasticity. Mechanics of Time-Dependent Materials 20(2), 175 – 195.
-  Makris N, Constantinou MC. 1993 Models of viscoelasticity with complex-order derivatives. Journal of Engineering Mechanics 119(7), 1453 – 1464.
-  Park SW. 2001 Analytical modeling of viscoelastic dampers for structural and vibration control. International Journal of Solids and Structures 38(44-45), 8065 – 8092.
-  Valerio D, Sa Da Costa J. 2009 Fractional Derivatives and Their Numerical Approximations II – Complex Orders. in: Symposium on Fractional Signals and Systems, Caparica.
-  Valerio D, Sa Da Costa J. 2011 Variable-order fractional derivatives and their numerical approximations. Signal Processing 91(3), 470 – 483.
-  Narahari Achar BN, Hanneken JW, Clarke T. 2004 Damping characteristics of a fractional oscillator. Physica A 339(3-4), 311 – 19.
-  Ryabov YE, Puzenko A. 2002 Damped oscillations in view of the fractional oscillator equation. Physical Review B (Condensed Matter and Materials Physics) 66(18), 184201.
-  Stanislavsky AA. 2004 Fractional oscillator. Physical Review E (Statistical, Nonlinear, and Soft Matter Physics) 70(5), 51103.
-  Tofighi A. 2003 The intrinsic damping of the fractional oscillator. Physica A 329(1-2), 29 – 34.
-  Rao SS. 2011 Mechanical Vibrations, 5th ed. Indianapolis, IN: Prentice Hall.
-  Franklin GF, Powell JD, Emami-Naeini A. 2015 Feedback Control of Dynamic Systems, 7th ed. Indianapolis, IN: Pearson.
-  Juang JN, Phan MQ. 2001 Identification and Control of Mechanical Systems. Cambridge, UK: Cambridge University Press.
-  Isermann R, Munchhof M. 2011 Identification of Dynamic Systems. New York, NY: Springer.