Ferromagnetic core coil hysteresis modeling using fractional derivatives

The modeling of a ferromagnetic core coil magnetic hysteresis has been considered. The measurement basis consisted of waveforms that have been recorded for various levels of the iron core saturation levels. The investigated models included classical cases as well as models including a nonlinear fractional coil. The possibilities of solutions for transient problems including such models have been recalled. The details of the estimation process have been described next, where each model evaluation made use of an original methodology dealing with periodic steady states. The influence of the model response on parameter changes has also been studied. Further on the parameter estimation procedure has been described, and the results for the various models have been presented.


Ferromagnetic core coil modeling
The modeling of distinct circuit elements is important due to the need for simulations of real phenomena of electrical engineering in power systems and electronic circuits. This, in no small part, is due to the need of foresight of the events that can occur in a studied system.
An element category, which is challenging when it comes to modeling, comprises of iron-core coils and transformers. This is largely caused by the need to reflect magnetic hystereses, which is especially difficult during transient states or nonsinusoidal responses with a significant contribution of higher harmonics.
The subject of accurate modeling of the more difficult cases of attempts to resemble magnetic hystereses is still an open topic [1,2]. Some of the most popular models are: the Jiles-Atherton model [3], the Preisach model [4] and the Chua model [5]. Over the years, there have been many propositions to modify these models in order to reflect the results of measurements made on iron-core coils more accurately [6,7]. Unfortunately, the application of a significant number of models yields at least one of the following problems: -a complex description, which leads to difficulties in implementations [8]; also-adaptability becomes a problem, e.g., when one wants to include the model in computations in a selected computational environment (e.g., in Matlab [9], Octave [10]), -a large number of equations, which directly leads to longer computation times [11], -the involvement of sophisticated functions, e.g., in models supporting anisotropy [12], -problems in numerical stability, due to some nonlinearities difficult to handle for typical algorithms [13].
The paper is focused on the introduction of models, where the above listed problems do not appear.
Often models are built with the B-H relationship in mind, which is applicable in electromagnetic field analyses. This paper concerns another category applied in modeling, which is a circuit approach, from which one can reconstruct a ψ-i relationship. This approach is useful when considering an element as a whole (for practical uses), when characteristics of inductive devices (e.g., transformers, fluxgate sensors) are investigated, instead of considering the physical properties of a material [14]. Also, phenomenological models can be applied, where although artificial variables are introduced-the computations still yield voltage and current waveforms. The circuit element approach has an advantage that it can be directly confronted with measurements performed on the modeled element.

Application of fractional calculus
In the study-fractional derivatives have been applied. They are a part of the mathematical field of fractional calculus [15,16], which includes definitions of both fractional derivatives and integrals.
Among the various definitions of the elements of fractional calculus (i.e., fractional derivatives and integrals [52][53][54]) the Caputo fractional derivative of order α ∈ [0, 1] has been considered: with t 0 being the initial time instance of an analysis, x (1) representing the first derivative of x(t) and (·) representing the gamma function:

Measurements
The studies in this paper make use of the results of measurements that have been accumulated throughout the last three years [48,50,55,56]. The measurements have been performed initially with a digital interference and event recorder RZ1 (manufactured by the Kared Company [57]). Over the successive investigation trials, an original measurement tool built on an Arduino platform [50] has been developed. Throughout the collection of recorded waveformsthe measurements of one selected coil have been chosen for further investigations. These comprise of periodic steady state voltage and current waveforms for various levels of the coil core saturation (exemplary cases in Fig. 1). More information on the measurements of the selected coil can be found in [55].
The flux linkage in a coil is defined as: where t 0 is the initial time instance of the analysis and u is the voltage across the coil. Because this quantity is not measured directly, the study applies the above  formula to obtain the flux linkage values from the voltage waveforms. The realization of the above is performed numerically, through Newton-Cotes formulae [58]. The flux linkage and current waveform pairs allow to reconstruct the hysteresis loops. An example of the obtained hysteresis loops is depicted in Fig. 2.

Investigated models
Profound experience earned in numerous previous studies [48,50,55,[59][60][61][62] has allowed to narrow down the class of the considered models. A model, which is often considered up-to-date and one of the easiest in its implementation, is the classical model, comprising of a parallel connection between a resistor and a nonlinear coil [63,64]. For simpler reference-this model is further on called the RL N model. Also, one can slightly modify the model by introducing a nonlinear resistance (this, in turn, is further on called the R N L N model) Both variants are presented in Fig. 3. The imperfection of the above-mentioned model variants in terms of the reflection of measurements (especially concerning the current waveforms [55,60]), during a large saturation condition of the core has forced an upgrade of the model. The resulting model consists of a resistor and a nonlinear fractional coil (this model is called the RL α N model). One can also consider its simpler version (which consists of a fractional, nonlinear coil alone-this is called the L α N model) and an extension of the model, which introduces a nonlinear resistance (called the R N L α N model). The nonlinear fractional coil is an element, which is described by the differential equation: where u is the voltage across that element, while ψ is an artificial variable (with the unit Wb · s α−1 ). Along with the above equation, the ψ-i relationship needs to be defined through a nonlinear equation. The computed ψ variable is only considered during a simulation process. The final flux linkage waveform for all of the models is defined through (3).
A further extension of the model can be made by including a fractional capacitor (Fig. 5). This is further on referred to as the R N L α N C β model. The fractional capacitor response is described through the fractional differential equation: with C β being an artificial parameter, with a unit of F · s β−1 .

Solving problems with fractional derivatives
With the conventional elements appearing in a circuit problem, the solution tools that are available generally fall under solvers for DAE (differential algebraic equations). When fractional elements appear, the problems can be described by the following general system of FDAE (fractional differential algebraic equations): where x(t) contains the variables that appear under fractional derivatives (and the orders are put in the vector α), while y(t) holds the remaining variables. D α t x(t) denotes the vector of fractional derivatives, i.e., the k-th entry is the fractional derivative of x k (t) with the order being the k-th entry of α (that is α k ). For such problems, the possibilities are much less numerous. For linear problems, one can apply some analytical solutions [66]; however, for a single nonlinear element appearing in the problem, these are no longer valid. Also, analytical solutions are only obtained for some specific source time functions (like unit step or sinusoidal functions [67]) and the most often considered system is one where all fractional derivatives are of the same order α [68] (the solution is then much simpler, but it is uncommon for such a system to show up in practice if more than one fractional element appears).
There are various numerical methods that have been described in the literature that deal with FDAE [these are often considered in various general forms [69][70][71], not necessarily as in (6)]. However, the availability of computational tools that allow to solve FDAE is very low. Let us consider a problem where the following elements appear: -autonomous voltage and current sources (described by any type of time functions) and controlled voltage or current sources, -linear and nonlinear resistors, -coils that are linear or nonlinear, which can also be regular coils or fractional coils, -capacitors that are linear or nonlinear, which can be regular capacitors or fractional capacitors (hence, leaving the possibility of a nonlinear fractional capacitor).
For certain implementation conveniences-the following form can be applied for the FDAE [72,73]: where x is the vector of variables under fractional derivatives (length of this vector: n x ), -y is the vector for the remaining variables (its length: n y ), -w joins the above two vectors: v(t) consists of the functions representing the time waveforms of the autonomous sources (the length of this vector is n v ), -M I is an n y × n y matrix, M II is an n y × n x matrix, M III is an n x ×n y matrix, M IV is an n x ×n x matrix and T is an n y × n v matrix, -the notation 0 k is applied in general for the description of a vector of k zeros.
Additionally, the vector F NL (w(t)) is introduced (with the length n NL ). In (7), it has been described more generally; however, it actually contains functions where each only depends on a single different variable. The only solvers that are available for such problems are the adaptive step size solvers [73,74] based on the method called SubIval [75,76] (this is an acronym for "the subinterval-based method for fractional derivative computations," which first appeared in [77]) and constant step size solvers based on fractional linear multistep methods [78] and product integration rule methods [79,80]. For periodic steady state problems (taking into account certain simplifications concerning the involvement of higher harmonics), one can apply the harmonic balance methodology [65], but this method is not available in a solver. When an analysis is required that includes one of the proposed models, it is best to apply one of the mentioned numerical solvers for transient problems. The estimation process alone in the case of this study does not require a transient solution to be computed in each model evaluation as it bases on periodic steady state solutions. The manner, in which a model evaluation is performed for selected parameter values, is described in Sect. 5. An example of how the w(t) vector is established and how the system described through (7) is formulated, when solving a circuit problem featuring the R N L α N C β model, is given in "Appendix A."

Model response in the estimation process
This section explains how a single evaluation of a model is performed, i.e., how an output waveform is produced for a provided input waveform. As it has been mentioned before-the measurement basis of the modeling process comprises of the voltage and current waveforms (more accurately-the recorded signal sample pairs). In order to reduce the time of the estimation process-for a single period of the 50 Hz voltage source-the final amount of samples taken into account is equal to n pts = 300 for the estimation process (in single model evaluations). For final observations of the results (comparisons between model and measurement waveforms), n obs = 2000 time axis nodes are taken into account. The voltage samples selected for the estimation process are expressed as the vector u meas , while the current samples are further on called i meas .
When having access to the voltage and current on the coil itself -one can inspect the response of the model for an assumed input of a waveform reconstructed from the voltage or current samples. Then, the other quantity serves as the output of the model, which, when compared with measurements, allows to ascertain the accuracy. Because in each of the studied models the voltage is common for all elements-it is selected as the input quantity.
Because of the nonlinearities appearing in the models, a single evaluation could prove demanding (and, most of all, relatively time consuming) when considering that evaluation as a time-dependent problem. Hence, the periodic steady state is considered, which allows for a particular approach. It involves procedures, where the following operations appear: (A) a time frequency analysis for the input waveform, resulting in a vector of odd harmonics up to h max (e.g., then the notation for the voltage harmonics is u), (B) a fractional differentiation or integration for the periodic steady state, which is in this approach performed on the complex valued representation; e.g., the fractional derivative of order α for a periodic waveform f (t) (expressed through a vector of harmonics f ) and resulting in a waveform g(t) (expressed as the vector of harmonics g) is computed through the Hadamard product [65]: where and: (C) for a set of samples in a vector f (expressing a waveform f (t)), for a nonlinear dependency f -g-the evaluation of a set of samples g (expressing a waveform g(t)), (D) a reconstruction of the output waveform (expressed as a vector g) from a vector of harmonics g.
The above operations are similar to what one can observed when solving a problem during the application of the harmonic balance methodology [65] mentioned in Sect. 4. For each of the elements of a model (connected in parallel), the voltage waveform is taken into account as the input and the current output for that element is evaluated. The final output of each model (depicted previously in Figs. 3, 4, 5) is then the sum of the currents. Figure 6 shows block diagrams for exemplary (most advanced) considered model elements.
What should be emphasized when the model performance is assessed during a parameter estimation procedure is that the task for a model is to accurately resemble hystereses for each core saturation level while using one set of parameters. Hence, such an evaluation, as described in this section, is performed for the waveforms of each considered core saturation level. More on how the model performance is then examined (and on the method of obtaining the parameters) can be found in Sect. 7.

Susceptibility on parameter changes
Before a parameter estimation procedure is conductedone can examine the influence of the particular parameters on the model output; where one can especially observe the resulting hysteresis (this can also be useful in an initial selection of the parameter values). The examples of this section are all given for the R N L α N C β model with the following initial parameter values: The above initial parameter values have been selected empirically. They allowed the model to yield a hysteresis shape similar to one that could appear in practice (with β and C β introducing only a slight influence).
The ψ-i relationship is described by the function: where the coefficient values: are applied. This function has been selected initially, for the trials of this section. It is actually one that has been used in previous studies [59]. The function used in the parameter estimation procedure (later on) has a more general form (more on this in Sect. 7). The nonlinear resistance, applied for the simulations of this section, is described by the following nonlinear u-i dependence: where initially: Similarly as in the case of the ψ-i relationship, this function has only been selected for the trials of this section, while a more general form is selected for the actual parameter estimation (Sect. 7). For a studied susceptibility on a selected parameter (or set of parameters concerning one nonlinear function), only the selected one (or ones) are modified, while the other ones remain as given above. The influence of the ψ-i function changes on the model output is presented in Fig. 7.
The effects are much like those that could be observed when modifying a nonlinear function in an RL N model, even though ψ-i does not directly represent the anhysteretic curve when α = 1.
The influence of the modification of the nonlinear resistor u-i function on the model output is presented in Fig. 8.
What has been observed is that the above modifications allow to influence the width of the hysteresis either for a smaller range of the flux and current levels or for the entire range.
Next, the influence on the fractional derivative order (denoted by α) applied in the nonlinear fractional coil has been studied. The results are depicted in Fig. 9.
A strong influence of the α parameter on the model output can be observed. What is interesting, in partic- ular, is a deformation that one could not obtain when modifying the previously studied features of the model. From these trials, one can determine that values at around 0.9 and above might be adequate and that this value should not be lower than 0.8 (even in this case the hysteresis shape does not look like one observed for ferromagnetic core coils). The modification of the fractional capacitor parameters is studied next. The influence of β is presented in Fig. 10, while the influence of C β is depicted in Fig. 11.
Some of the resulting features are like no others that could be imposed when modifying the previous parameters. One can also observe some influences similar to those that have been noticed earlier (e.g., concerning the hysteresis width). However, a distinct feature that can be introduced is what could be interpreted as an opposite slope of the anhysteretic curve, e.g., similar to what can be observed when taking into account capacitances in transformer windings [81].

Parameter estimation procedure
The computations of the study have been performed in the GNU Octave [10] environment. The introduction of the measurement waveforms, their processing, the possibility of model evaluations and the implementation of all other necessary operations all required individual computational scripts to be written.
Obviously, the model type determines the amount of the estimated parameters. What is also important, in the case of appearing nonlinear elements, is the form in which the nonlinear function is considered. These   The above has proven to work well [50,55], but an even better approach has been established in order to assure that the ψ-i relationship results in a monotonically increasing function. Instead of the ψ values-the estimated values can be given as ψ k and for each ψ k , k = 1, 2, . . . n ψ-i : For this study, ψ min = 10 −5 Wb · s α−1 has been assumed, along with n ψ-i = 15, i max = 1.5 A and i = i max n ψ-i −1 . The same concept can be applied for a nonlinear resistor. However, here the u values are assumed from the start, while the i R values are estimated. The number of u-i R pairs is n u-i R = 15, u max = 450 V and u = u max n u-i R −1 . Taking the above into account, as an example, a single set of parameters for the RL α N model consists of the ψ k values and R, α values, while for the R N L α N C β model a set of parameters consists of α, β, C β along with the ψ k and i k values.
The evaluation of each model involved the consideration of 6 sets of voltage and current waveform pairs and the computation of the output quantity (which is the current) according to what has been described in Sect. 5. The current output of the model is then expressed through samples that are put in a vector i. These are compared with the measurement results of i meas . The following value is computed for each considered set of the voltage-current sample pairs: where k is the number of the voltage-current waveform pair. In the estimation process, the f k values for all the considered sets are summed up. Then, the sum of these: is minimized. The estimation process for each model has made use of an algorithm, which consisted of series executions of the Nelder-Mead method (implemented in the nelder_mead_min function in the optimization package in Octave) followed by an application of the multidirectional search method (through the fmins function). Such hybrid optimization algorithms (combining two optimization methods) can be found in [82]. In order to possess an insight into the accuracy of the agreement of the waveforms one can observe the model and measurement waveforms and the resulting hystereses. Additionally, the following error value is defined: where i max is the maximum of absolute values of the vector i meas and, as previously in the case of f k , k is the number of the voltage-current waveform pair. Note that the above is evaluated for each set of the voltagecurrent sample pairs.

Modeling results
For each more advanced model, both strategies have been applied, where the initial parameters have been selected either through: -parameter tweaking, where the values have been determined through observations of the model output (hysteresis, as in Sect. 6), -a selection of parameters, which have been optimal for a less advanced model (e.g., the starting parameters for the R N L α N model have been directly determined by the optimal parameters for the RL α N model).
The estimation process for each model has been executed for both of the above starting parameter strategies, while the final results have been taken from the strategy that has produced the lower F value. A comparison of the objective function values for the considered models is presented in Table 1.
In a general comparison, one can observe: -a slight improvement when adding a nonlinear resistance in the classical model (i.e., resulting in the R N L N model),  Table 2 Comparison of errors (in percentages) for the considered models (values rounded to the ten thousandth) -an unexpected advantage of the L α N model in comparison with the R N L N model, -the most noticeable improvements (of the agreement between the measurements and simulation results) have been observed up to the application of the RL α N model, while consecutive upgrades to R N L α N and R N L α N C β models have produced only slight improvements.
In this study, an important fact is that the modelmeasurement resemblance is fitted for various levels of the core saturation (with F actually presenting a measure of the total accuracy) for a single set of parameters. In order to determine the ability to reflect the measurements for each model for each of the core saturation levels separately-the k values have been determined [defined by (17) in Sect. 7]. The results are presented in Table 2.
Additionally to the previous observations, what can be noticed is that for some improvements of the objective function value F there have been some reduced accuracies when it comes to the errors in the cases of distinct core saturation levels (e.g., for 6 when upgrading to the R N L N model). The classical models (without fractional elements) generally seem to fail in attempts to adequately resemble each considered set of measurements at the same time for one selected set of parameters. This can also be noticed when observing the model and measurement hysteresis comparisons for one of the classical models. As an example-the comparison for the R N L N model is depicted in Fig. 12. An additional flaw of the integer-order models is the inability to accurately resemble the hysteresis for when the core becomes more saturated (magnification in Fig. 12). This can be especially observed above the saturation knee.
As mentioned before-an unexpected improvement in the accuracy has been observed when applying the L α N model (consisting of just a nonlinear fractional coil-in this case the parameter α = 0.9664). An improvement can be observed also when upgrading it to an RL α N model (in this case α remained unchanged), where not only F is a smaller value but also all the k values indicate a better result. When applying the R N L α N model instead, a slight improvement is observed in almost all cases (except 6 , though that value is very similar to its equivalent in the RL α N model). A compilation of the hysteresis comparisons for the various iron core saturation levels for the R N L α N model is presented in Fig. 13. One could also notice that the value of the α parameter has not changed.
The addition of a fractional capacitor to the model has also produced an insignificant improvement, which is why the hysteresis comparisons are not depicted in this paper. The parameter estimation procedure has yielded the fractional derivative orders α = 0.9662 and β = 0.7331.
A very good resemblance can be observed when comparing the model results with the measurements, which can also be noticed when the current waveforms are compared (Fig. 14). As a reminder: the voltage waveforms are not compared as for each of the simulations the input for each model has been that voltage waveform.

Summary and conclusions
The investigations of this paper concerned the modeling of a ferromagnetic core coil response, where an emphasis has been put on the reflection of the hysteresis. A circuit approach has been applied, meaning that the element is treated as a whole and the ψ-i relationship is taken into account. The measurement basis consisted of sets of voltage and current waveforms that have been gathered throughout the authors' previous studies. The collections of recorded waveforms have allowed an insight to various levels of the core saturation levels (examples in Figs. 1, 2), where the flux linkage waveforms have been recreated from the voltage waveforms (this is described in Sect. 2, where (3) has been applied). From the above-six of the sets (each representing a different level of the core saturation) have been taken into account for further investigations.
The task for a model was to have the ability to accurately reflect the measurement waveforms (and, hence, the hystereses) for the various core saturation levels. Unlike other studies, where often only one or two scenarios are taken into account [3,59] the study aimed at an ability of a model to reflect all the cases of the core saturation levels for a single set of parameters.
Experience gained in previous studies has allowed to extract six different circuit element-based models that could be useful in this study. These included classical circuit elements and elements applying fractional calculus (their detailed description has been given in Sect. 3

, Figs. 3, 4 and 5):
-the classical model (often appearing in the literature) featuring a resistor and a nonlinear coil-this has been called the RL N model, -an extension of the classical model, replacing the linear resistance with one that is nonlinear (the R N L N model), -a simple (yet what has later turned out to be surprisingly accurate) model featuring a nonlinear fractional coil alone (the L α N model), -a model featuring a resistance and a nonlinear fractional coil (the RL α N model), -an extension of the above, which replaced the linear resistance with a nonlinear one (the R N L α N model), -a further extension, adding a fractional capacitor (the R N L α N C β model). The inclusion of fractional calculus-based elements adds a certain complexity to the mathematical form of the considered problems [which are then either considered in a general form given by (6) or, more specifically, by (7)]. The ability to solve such problems is addressed in Sect. 4.
The models generally all consisted of single elements connected in parallel, which allowed to apply a particular approach when each of the models has been evaluated in an estimation process. A single model For a better familiarity on how each feature of the considered model influences their response (especially on the shape of the resulting hysteresis), a susceptibility analysis is performed and presented in Sect. 6. Further on, this has also been helpful when selecting a set of initial parameters for each of the models. The final general forms of the nonlinear functions (for nonlinear resistances and nonlinear coil ψ-i relationships) are discussed in Sect. 7 along with details on the estimation procedure and the form of the objective function [described by (15) and (16)]. Additionally, the reflection of the hysteresis for each core saturation level separately has also been considered, where the error given by (17) has been defined.
The modeling results are presented in Sect. 8. Both the final values of the objective function and the error values (for each considered core saturation level) have been addressed (Tables 1, 2). From the analysis in that section, the following can be concluded: -the R N L N model has introduced some slight improvements-the objective function value has improved in general; for some core saturation levels, one can observe a significant improvement (given by the 1 , 2 and 3 errors); however, the accuracy for the remaining ones has been decreased (especially for the lowest core saturation level, which could be observed through the 6 error value), -for the classical models, one can generally observe two flaws: the first is that they are not able to reflect all the hysteresis saturation scenarios at once for one set of parameters (although a very good agreement can be achieved for a single hysteresisthis has been shown in previous studies [59,61]), while the second is the difficulty in properly resembling the hysteresis width above the saturation knee (example in Fig. 13), -as mentioned in Sect. 8-a very interesting result can be observed when applying a model comprising of a nonlinear fractional coil alone (the L α N model)-one can observe a significantly better result than for the models basing on classical circuit elements (both in terms of the k error values and the objective function F), -a still noticeable improvement can be observed when adding a linear resistor to the nonlinear fractional coil model (resulting in the RL α N model), -an insignificant improvement could be observed when the model is extended to the case of the R N L α N model-the error values have improved in five of the six cases (by an insignificant amount) and one case (described by 6 ) suffered from a slight loss in accuracy; the objective function also suggested an insignificant improvement, -the R N L α N C β model also introduced only slight improvements in the resemblance of the studied coil hystereses.
The RL α N model, in some terms, could be considered as the optimal model for the studied coil because of the accuracy improvement in comparison with less advanced models. It is worthwhile to point out that this model only requires one additional parameter in comparison with the classical RL N model (i.e., the parameter α). Further additions to the model (a nonlinear resistance or a fractional capacitor) have not introduced significant improvements. However, one can already notice a very good agreement for the R N L α N model, leaving not much space for further improvements. The hystereses (Fig. 13) and waveforms (Fig. 14) have been very accurately reflected, where in some cases the differences are almost indistinguishable. This suggests at a very useful application of this model in simulation practices on real objects (e.g., for ferroresonance analyses [55], even transient state simulations [48]).
The inclusion of the fractional capacitor has not introduced significant improvements for the reflection of the studied hystereses; however, the inclusion of this element allowed to introduce features that could not be achieved otherwise (discussed in Sect. 7, observed in Figs. 10 and 11). This can be useful when circumstances would lead to such features being observed (it can be pointed out that such cases are possible for real objects applying ferromagnetic cores [81]).
Future studies will be aimed at the design of models that could accurately resemble measurements not only obtained for various ferromagnetic core saturation levels, but also for scenarios, where very significant components of higher harmonics will be introduced on purpose.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Example of FDAE formulation
This appendix explains how a system in the form of (7) is formulated for a circuit containing an exemplary model that has been considered in this paper.
The R N L α N C β model has been selected for the example as it is the most advanced one that has been studied. The system will be formulated in a general manner, i.e., the voltage and all the selected currents (on nonlinear and inertial elements) along with node potentials will be included in the equations. The mentioned variables are marked as in Fig. 15.
Initially, only the equations will be formulated, which the model introduces, while further ones will be introduced depending on where the model is placed in a circuit. The x(t) vector consists of the variables under fractional derivatives, hence: where the derivative orders are given by: while: The source vector v(t) will be introduced later, when the placement of the model in the circuit will be determined. The potential difference introduces the equation: which determines a row in the M I and M II matrices. For the purpose of this example, it is assumed that the equations start from n start , hence: M I n start ,1 = 1, M I n start ,2 = −1, M II n start ,2 = −1.
The remaining equations that the model introduces are the nonlinear equations and the fractional differential equations. For the nonlinear equations, when implementations are considered [72,73], it is convenient to introduce an auxiliary vector i arg , whose entries denote the indices of the arguments of the functions appearing in F NL , e.g., for n NL = 3 and: the vector F NL takes the form: The nonlinear equations describe the u-i NL and ψ-i ψ relationships. For a u(i NL ) description: The first entry in F NL is then the f u function, while: According to the form given by (7), the nonlinear equations appear as the last n NL algebraic equations.
Assuming that the equations the model introduces appear starting from n NLstart -the left-hand side in (25) introduces: M II n NLstart ,2 = 1.
For an i NL (u) function: i NL = f i (u).
where f i then appears as the first entry in F NL , while: and: M I n NLstart ,3 = 1.
Similarly for the ψ-i relationship-assuming that the nonlinear equation is: the second entry in F NL is then the f ψ function and: with: M II n NLstart +1,1 = 1.
On the other hand, for a nonlinear function: f appears as the second entry in F NL and: with: M I n NLstart +1,4 = 1.
The remaining equations are the fractional differential equations: and: Assuming that already n rem. differential equations have been introduced, the following entries are added: M IV n rem. +1,2 = −1. and: M III n rem. +2,5 = − 1 C β .
Further additions to the matrices and the F NL vector are introduced when the remaining part of the circuit is determined. For a simple case of a circuit when only a voltage source, described by e(t), is connected to the terminals of the model, i.e., the V 1 and V 2 nodes, and directed at the V 1 node, no additional variables are introduced to x(t) and y(t), while v(t) becomes:
When all the necessary matrices, vectors and functions have been defined-the system can be solved through a solver basing on the SubIval method as mentioned in Sect. 4, e.g., the one that is available for GNU Octave and Matlab [73].