Application of Selected Numerical Methods to Model the Fractional–Order System Behavior of Nonlaminated Magnetic Actuators

—The electromagnetic dynamics of nonlaminated magnetic actuators are highly inﬂuenced by eddy currents and minor perturbations like core saturation, hysteresis as well as fringing and leakage ﬂuxes. In the literature, analytical high– ﬁdelity models describing these phenomena are known, which lead to complex reluctance networks or transcendental system descriptions with fractional–order characteristics. Therefore, they are not suitable for a direct implementation within the actuator control. Previously, we provided appropriate analytical rational approximations that allow a digital real–time implementation of these models on a microcontroller. However, the inclusion of the minor perturbations, if possible, leads to impractical model orders requiring simpliﬁcations, which compromise the model accuracy. This article studies numerical methods to reduce high model orders or directly approximate the transcendental systems or empirical measurement data. The greater degree of freedom allows for a possible higher model accuracy with sufﬁciently low orders. We review and improve existing approaches like Levy’s method and Vector Fitting and apply them to the frequency response of the underlying fractional–order system. Furthermore we propose an order reduction algorithm based on a pole–zero– cancellation with tracking error compensation. Using measurement data, a comparison shows that the numerical approaches match or excel our previously studied analytical approximation.


I. INTRODUCTION
The physical phenomena inside solid iron cores, like the magnetic skin effect, can be described with the diffusion equation, which naturally leads to a system description with half-order derivatives [1]. These so-called fractional-order systems [2] are characterized by the occurrence of the Laplace variable with a fractional exponent . In case of magnetic actuators, with = 1/2 or 1/4 usually forms the argument of transcendental functions like hyperbolic or Bessel functions. Because of the mathematical complexity, the established actuator models [3] are of no great use for actual control system applications like the model-based flux control proposed in [4].
For these reasons, we previously developed various approximation approaches, based on continued fractions as well as M. Hecht is with the Department of Information Technology and Electrical Engineering, ETH Zürich, 8092 Zürich, Switzerland. (email: hechtm@ethz.ch) R. Seifert and W. Hofmann are with the Chair of Electrical Machines and Drives, Technische Universität Dresden, 01069 Dresden, Germany (email: robert.seifert@tu-dresden.de, wilfried.hofmann@tu-dresden.de).
Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.36227/techrxiv.14899689 This manuscript is not peer-reviewed and may be subject to changes. This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.
Padé-and Matsuda-approximations, which lead to a rational function description eligible for digital implementation [4]. All proposed methods are analytical and keep the structure of the solution of the underlying diffusion equation. However, the hence imposed restrictions either necessitate high orders of the rational model function of at least = 25 or result in a trade-off between accuracy and bandwidth. Model refinements like reluctance networks [5], to model fringing and leakage fluxes, further increase the order . Nonlinear perturbations like hysteresis and saturation [6] cannot be included directly and may only be considered by additional high-order models.
At the cost of possibly unphysical solutions, numerical methods can overcome these drawbacks as the system's zeros and poles can be chosen arbitrarily as they are not bonded to any physical law or mathematical description. Only the stability of the system is to be ensured. The numerical fitting of transfer functions on empirical data is well studied, with the Vector fitting of Gustavsen et al. [7] as well as Levy's Method [8] and its various improvements (overview in [9]) as most prominent examples. However, these methods are derived for simple example systems and usually applied to systems with integer-order behavior. We only know a few applications to transcendental and fractional-order systems describing nonlaminated electromagnetic actuators. Both, Zhou et al. [10] and Zhong et al. [11] used Levy's Method as a means to an end to verify their measurement results for a nonlaminated magnetic bearing. But they also permit fractional exponents for the fitted model, which cannot be discretized and are therefore not constructive. Alternatively, Herzog et al. [12] as well as [10] propose to fit a low-order integer system ( ≤ 4), which might not be sufficiently accurate. As long as biquad cascades (second order sections systems) are used for the later digital implementation of the model in the control system, much higher model orders ≤ 50 are feasible with no apparent drawbacks [13].
This article studies the application of Levy's Method and Vector fitting considering various improvements. They are used to determine an accurate integer-order rational representation of a nonlaminated magnetic thrust bearing, by means of the transfer function of the controlled coil (with the resistance Cu and dead time σ ) The so-called effective inductance eff ( ), or eddy-inductance [1], can be roughly described as a half-order low-pass and therefore serves as the adequate initial example for our case. Furthermore, we use a pole-zero-cancellation with a newly proposed tracking error compensation to reduce the excessive orders resulting from the high-fidelity analytical model in  [4], to match the model orders of the numerical approaches. Additionally, it can be applied to the numerical approximations to cancel insignificant pole-zero pairs. Finally, we conclude with a comparison of all approximation methods, numerical and analytical, with the measured frequency response and a FE-analysis of the system in (1).
II. PRELIMINARIES The nonlaminated magnetic actuator employed as an example, is the magnetic thrust bearing from [4], but we are only interested in its approximated system description and the resulting frequency response. The effective inductance eff ( ) of the bearing stems from the reciprocal of the effective reluctance R eff of its magnetic circuit where is the number of coil turns. It consists of the sum of various transcendental functions T ( ), like the tangens hyberbolicus or modified bessel functions. Each of them resembles a half-or quarter-order PD-element, so the overall reluctance can be described approximately with the equivalent implicit fractional-order system (parameters in Tab. I) To our knowledge, this kind of system behavior was not approximated before with the numerical methods studied in this article -at least not in the context of magnetic actuators for reasonable high system orders . The mentioned exceptions [10], [12], which indeed apply Levy's Method in a similar case, limit to 3. . . 4 for no given reason. We presume coefficient-quantization errors [14] during the calculation, which affect the system stability, most likely prevented higher orders. However, they can be avoided by the consequent use of high-digit variable-precision arithmetic (VPA) with at least 80. . . 180 digits [13] (depending on ) until the final discrete poles and zeros for the digital implementation are determined.

III. LEVY'S METHOD
Levy's Method approximates a complex curve by minimizing the sum of squared errors at every sample point = 2 of a given dataset or function (j ) in the frequency-domain with = j = j 2 within the range [ min , max ]. We seek a rational frequency response function where 0 , . . . , , 1 , . . . , are the coefficients of the numerator (j ) and denominator polynomial (j ), respectively. By multiplying the error = (j ) − (j ) with the denominator (j ) at every sample point, one determines the biased error which is set to zero to define the minimizing function: (j ) · 1 + 1 ·(j ) + 2 ·(j ) 2 + · · · + ·(j ) − 0 + 1 ·(j ) + 2 ·(j ) 2 + · · · + ·(j ) = 0 . (6) It follows the system description for every sample point = j with the sampled system vector As an ordinary least squares problem, one can solve for by minimizing the sum over all quadratic errors min − 2 , where · 2 denotes the Euclidean vector norm. In the original approach [8], Levy solves that problem by summing the squared absolute errors across all sample frequencies and setting the partial derivation with respect to each of the polynomial coefficient to zero. A more elegant way is chosen in [9]: since it is a classic least squares problem, its solution must fulfill the normal equations as proofed in [15]. There are several ways to solve the resulting SLE, whereas the QR-decomposition [15] is favorable as it is fast and less prone to quantization errors. In addition to that, we remind, that especially for wider frequency ranges, the use of VPA or the alternative partitioning of the original function into sections along the frequency axis [16] is essential.

Weighting and Bias Correction
Since the error is multiplied by the denominator in (5), not all errors are equally weighted and one may not obtain the best least squares solution possible. This bias can be overcome through an iteration procedure,where the biased error (j ) is divided by the denominator from the previous iteration [17]: leading to the weighting where [ ] indicates the current iteration step. This procedure, often referred as the Sanathanan-Koerner (SK) iteration [18], is no true weighting in the strict sense as it solely corrects the induced bias. That is why, we adapt Noda's favorable suggestion [16] to minimize the sum of the squared absolute values of the relative errors [ ] (j ): This can be achieved by changing the weights to forming the diagonal weighting matrix which we apply to the system matrices before (9) [ ] The error sum does not converge monotonically, why the iteration is aborted as soon as [ ] ≥ [ −1] considering [ −1] as the best solution. Finally, after the calculation of all coefficients 0 , . . . , , 1 , . . . , the determined frequency response function (j ) can be transformed to a transfer function ( ) by the substitution = j .
It is important to note, that in case the given dataset (j ) was derived by evaluating a transfer function ( ), e. g. (1), the result will most likely not equal, but approximate the input so that ( ) ≈ ( ).

Stability Enforcement
The example system in (1) is stable. While all analytical rational approximations proposed in [4] naturally maintain this property, it is not guaranteed for Levy's Method. However, for our introductory example (2), only very high orders > 50 produce unstable poles. In the rare case these high orders are unavoidable, a zero with a positive real part nearby can be identified usually (otherwise must be modified). If so, one can flip the real parts of both, pole and zero, to the left half plane to stabilize the system while only slightly affecting the phase but not the magnitude. We call this technique pole-zero flipping and borrow it from [7], where it is used in the context of vector fitting, which we discuss next.
The results of Levy's method for different orders can be seen in Fig. 1. While for orders ≤ 4, like in Herzog et al. [12], the maximum amplitude and phase errors (0.7 dB, 4.2°) are not insignificant, they become negligible already for the reasonable order = 9. Although the pole-zero cancellation (PZC), which we will explain later in section V, does increase the errors, they are still very small (0.2 dB, 0.76°) and barely visualizable. But the algorithm effectively eliminates obviously dispensable pole-zero pairs and further We like to point out, that the direct fitting with e. g. = 21 leads to smaller errors than for the reduced system of initially = 27. However, the PZC results in a different weighting, focusing on the fractional-order component of the system, which is characterized by the steady alternation between poles and zeros known from ladder networks.

IV. VECTOR FITTING METHOD
In contrast to Levy's Method, Vector Fitting [7] uses a transfer function in its factorized form as the approximation, such that the poles and zeros are calculated directly. The ansatz is defined in the Laplace-domain in the rational form with the gain , the residues and the poles , the latter being either real or occurring as complex conjugate pairs. The maximum index equals the desired order of the denominator polynomial. The term + or even + +ℎ is optional, depending on whether = or = + 1, respectively, is desired.

Pole identification
The poles in (16) 1 Real poles are recommended for smooth functions. Complex conjugate pairs = − ∓ j with = 100 can avoid ill-conditioned system matrices for functions with distinct resonance peaks [7]. is valid, where ( ) is an unknown function in the form: The fitting parameters = [ 1 , . . . , In the original approach ∼ ∈ R is fixed to 1, however, it can be used as additional degree of freedom to improve the convergence ability of the poles, called relaxation [19]. By multiplying the original transfer function ( ) with ( ) inserting (17) and (18) and bringing all terms to one side, one gets the functional description For any given argument = j of the original function, one can form the sampled system description By concatenating = [ T 1 , . . . , T ] T analogously to (9) for all sample points and splitting both vectors into real and imaginary part, one finally obtains the overdetermined SLE which is extended by one line with the weighting vector .
It is calculated for every starting pole ∼ with with denoting the identity matrix and 1 a column vector with ones. In case the solution leads to complex conjugate pole pairs, we advice to balance them as laid out in [7].
The zeros = [ 1 , . . . , ] cannot be determined from the residues in a similar manner just yet, as they still correspond the starting poles ∼ , so the whole pole identification process has to be repeated with the new poles as starting poles. Besides, the following additional steps may be required and recommended before.

Iterative Fast Vector Fitting
The pole identification process can be repeated not just once, but multiple times, each time with the previously calculated new set of poles as starting poles, until it converges to an optimal solution [19]. It also avoids an ill-conditioning of the system matrix if the difference between the starting and solution poles is large or they are not real as we assumed. As abort criterion we chose: In each iteration step -except for the final residue identification -the residues are discarded, wasting computational effort while solving the SLE in (23). That is why Deschrijver et al. [20] propose an optional modification, they call Fast Vector Fitting, which omits the dispensable solutions 1 , . . . , already within the algorithm. Alternatively, Noda [16] suggests to determine the poles by a coefficient-based algorithm, similar to Levy's Method in section III, before he continues with the residue identification. He observed a more accurate determination of non-dominant poles in case the original frequency trace ( ) shows various peaks. However, for our smooth dataset we cannot identify any advantage.

Stability Enforcement
As with Levy's method, vector fitting may also produce unstable poles. As the zeros are not calculated yet, the suggested pole flipping, which inverts the sign of the real part of all unstable poles [7], may only have a marginal negative effect in this case [16]. The magnitude is not affected by the inversion and the change of the phase angle is mostly compensated by the yet to be determined zeros.

Residue Identification
After calculating the optimal set of poles , a last complete iteration of the pole identification algorithm, but with fixed poles, finally leads to an appropriate set of residues and therefore completes the rational function description in (16) in form of its partial fraction expansion. Unlike Levy's method, vector fitting does not generally utilize polynomial coefficients unless a transfer function description is required. Hence, it is less prone to coefficient-quantization errors and can be calculated with double-precision, which makes it considerably faster than Levy's method. Nevertheless, one may determine the transfer function (4) via its state-space description by with ∼ = diag( ), ∼ = 1 T , ∼ = and ∼ = 0. We note that before this step we switch over to highdigit VPA, because any polynomial representation is critical in respect to floating-point-precision. Otherwise, it may falsify subsequent calculations, like the determination of the zeros = [ 1 , . . . , ] or the discretization of the system. This difficulty especially arises for fractional-order (transcendental) systems, like in our case, as the occurring slopes and phase angles are no multiples of 20 dB/dec and 90°, respectively, which the poles and zeros naturally represent. In case a digital implementation is the final aim, one may also consider the Discrete-Time Vector Fitting as proposed by [21], where they directly compute the discrete poles and zeros in the Z-domain.
If VPA is not available, the quantization errors can be avoided by calculating the zeros numerically with e. g. Newton's method, where the expected alternation of poles and zeros can help to find appropriate starting points to find every zero. Then one may directly use poles and zeros to form a socalled second-order-sections system as explained in [13]. That way an error-prone high-order polynomial can be avoided, which is an unique advantage of the vector fitting method.
The fitting results of the Vector Fitting method in Fig. 2 are comparable to Levy's method, although slightly less accurate. However, for the low order = 4, the algorithm cannot identify a sufficiently well distributed pole-zero constellation, resulting in an unacceptable high maximum phase error of 31.3°. Still, for higher orders the regression errors become negligible. At the cost of tolerable errors, the pole-zero cancellation can effectively further reduce the system order.
V. POLE ZERO CANCELLATION WITH TRACKING ERROR COMPENSATION All numerical approximation methods discussed in this article as well as the analytical approximations we studied in [4] may produce pairs of poles and zeros, which are close to each other or even match. They increase the overall system order , without significantly increasing the accuracy of the model. Hence, it would be advantageous to remove these pairs and reduce the order without a significant loss in accuracy.
Analytical approaches for an order reduction, like variants of the Padé approximation [4] or Routh approximation [22] are not able to produce system descriptions with well distributed poles and zeros, which tend to cluster together for already reasonable low orders , as we will show later in Fig. 4. On the other hand, if one chooses the order too low, decisive poles may be missing. That is why the minimum practical order of these analytical methods is constrained, as the underlying mathematical laws do not provide the necessary degree of freedom. For that reason, we apply the known concept of pole-zero cancellation (PZC) not only to the unreduced and reduced analytical approximations in [4], but also the numerical approaches studied here, to further reduce the system order . To avoid potential deviations of the transfer function, we propose a tracking error compensation (TEC) to allow for usable results even for small orders, as illustrated in Fig. 3. In our implementation, we presume that all poles and zeros are either real or occur as complex conjugate pair (cc-pair). Furthermore, we require them to be stable.

Algorithm
Step 1: First we bring the initial function in the form to keep the gain ∼ constant during the process. From the poles and zeros we define the sets and and their union : As the algorithm solely relies on the absolute values | |, | | of every pole and zero, we treat cc-pairs as a single instance ≡ [ , ], which can only be canceled by another pair ≡ [ , ], though.
Step 2: Additionally, is defined as the set of all pole-zero pairs [ , ], or pairs of cc-pairs [ , ] ≡ [ , , , ] for that matter, where no other pole or zero lies within the bandwidth of these pairs 2 : Step 3: For all pole-zero pairs [ , ] ∈ we calculate the logarithmic distance between each other as = log 10 | | − log 10 | | .
We then form the first-order pole-zero transfer functions for real [ , ] or in case of cc-pairs [ , , , ] : Step 4: Tracking Error Compensation (TEC): Each (cc-) pole-zero pair [ , ] has one or most likely two adjacent neighbors + and − or pairs of cc-neighbors + ≡ [ , ] + and − ≡ [ , ] − . To identify the most appropriate to relocate, by maintaining the local error as small as possible, we define the neighbor functions or for cc-pairs, respectively, (37) Thereafter, between every pair of corrected neighbors the one ∼ is chosen, whose maximum magnitude error +|− over an is smaller, with ∼ + |− q ( ) = q ( ∼ + |− , ) as the corrected neighbor transfer function.
Step 5: Finally, from all pole-zeros pairs [ , ] (or pairs of cc-pairs) one selects the pair, with the minimum error for cancellation and correct the according neighbor with ∼ = ∼ + |− ( ). Although it would be theoretically possible to cancel multiple pairs at a time, we recommend to use the algorithm iteratively with a single cancellation at every step, as it influences adjacent pole-zero pairs.
The algorithm can be implemented in two manners. First, as a for-loop with a pre-determined target system order ∼ .
However, as it utilizes only existing poles and zeros, it cannot reach the accuracy of the previously discussed fitting algorithms. Therefore we recommend, secondly, to implement the pole-zero cancellation as a while-loop and iterate until no further pole-zero pairs [ , ] can be found, whose logarithmic distance is smaller than a certain tolerance. We chose tol = 10 −2 , so every pole-zero pair that is separated by less than 1 % of a frequency decade is canceled out. Furthermore, in case of complex values, one may analogously only permit pairs whose imaginary parts are very close log 10 ( ) − log 10 ( ) < tol = 10 −2 . In the previous sections, we applied the numerical approximations to the transcendental system in (2) directly derived from the diffusion equation and omitting perturbations of the magnetic circuit, like leakage, hysteresis and saturation. We did the same for various analytical approximations in [4]. In this section, we now compare all variants based on a measured frequency response corresponding to the transfer function of the coil in (1). Beforehand, we deduct the dead-time element for improved fitting, which is later reinstated. Fig. 4 shows the results of the comparison in terms of the magnitude and phase errors and the pole-zero constellations.
First, we like to note, that the analytical approximations are based on an analytical model, while the numerical approximations are a direct curve fit. Hence, it is to be expected, that the latter outperforms the former in terms of accuracy. However, in the later practical application within the actuator control, the difference is not decisive. With reasonable high all approximations are sufficiently accurate with maximum magnitude errors of ≤1 dB and phase errors 1.9. . . 4.2°.
The analytical approximations exclusively produce stable real poles and zeros, which is not the case for Levy's method and Vector Fitting, contrary to each methods introductory section. That suggests, that complex poles and zeros are necessary to accurately map the influence of the perturbations. Furthermore, one has to enforce the stability of some of the produced poles for their lack of a mathematical reference to the intrinsically stable solution of the diffusion equation, which describes the naturally stable system. Even though this does not pose a problem, it emphasizes the necessity to thoroughly inspect the fitting results for unphysical system behavior, e. g. poles which are physically not possible.
The rational form of a half-order integrator would be composed of evenly distributed poles and zeros, so for the present half-order low-pass one would expect a similar constellation, at least in parts. However, the poles and zeros of the analytical approximations SCFE (no order reduction) and PASR (equals SCFE with order reduction by Padé approximation) are badly distributed across the frequency range. They are concentrated around the main field time constant and frequency h = 1/(2 h ), which is a direct consequence of the underlying continued fraction expansions. Unlike the Padé approximation (cf. [4]), the proposed pole-zero cancellation can effectively reduce the order of the unreduced SCFE to = 9 leading to the expected pole-zero distribution, without worsening the accuracy. In spite of being computed analytically, the MAEIS relies on well distributed pre-defined sample-points. In theory,   it avoids the clustering of poles and zeros by design, similar to the numerical approaches, but the superposition of a halfand quarter-order system (cf. (3)) diminishes this advantage. Without these mathematical restrictions, the numerical approaches produce well distributed poles and zeros for appropriate orders ≈ 9. Only higher orders would cause clustering, which can be reversed by the pole-zero cancellation resulting in an optimized weighting of the fractional part of the function.

VII. CONCLUSION
This article adapts Levy's method and Vector Fitting as the most popular curve-fitting algorithms on the coil transfer function of a nonlaminated electromagnetic actuator. The system shows a fractional-order behavior, characterized by a magnitude slope between −15. . . −10 dB/dec and a steady alternation of poles and zeros, which generally requires higher system orders. This fact makes it especially prone to quantization errors, why all studied algorithms should be implemented with the consequent use of variable-precision arithmetic. That way, we were able to realize system orders of at least = 9 reducing the regression errors to a negligible level below 0.4 dB and 0.8°. On the other hand, compared to the previously studied analytical approximations with ≥ 26, the system order could be kept considerably smaller, allowing more efficient implementations.
In principle, the referenced analytical rational approximations, are able to perfectly match a transcendent function, in case a continued fraction expansion or Padé approximation is feasible. However, real empirical data for electromagnetic actuators is affected by nonlinear perturbations (e. g. saturation, hysteresis, fringing and leakage), which cannot be described exactly by mathematical laws like the diffusion equation. Numerical algorithms, on the other hand, are not bonded to these laws and hence allow for a significantly higher accuracy in such cases, while keeping the system order much lower. But they are data-dependent, which makes them barely adaptable to parameter changes and may require a fallback to analytical generated data. It becomes clear, that numerical and analytical approximations both have their distinct use cases and may be combined for optimal results.
In particular, the coefficient-based Levy's method tends to be the most accurate over the entire approximation bandwidth [ min , max ], especially for lower system orders < 9. In contrast, as Vector Fitting primarily relies on poles and zeros, it is less prone to quantization errors and may be calculated with 64 bit-floating-point arithmetic, which makes it computationally very efficient. For higher orders > 9 the fitting results of both methods are practically identical, why the latter is recommended as the more elaborated choice.
As an addition to all approximations, numerical as well as analytical, we propose a pole-zero-cancellation algorithm with tracking error compensation, which is particularly appropriate in case of fractional-order systems. To accurately map their atypical slopes and phase shifts, a higher number of poles and zeros is required. However, as a side product, the algorithms usually also generate dispensable pole-zero pairs in the non-fractional segments of the transfer function, which can only be avoided by a deliberate manual weighting process. Pole-zero-cancellation allows to remove these pole-zero pairs subsequently in an automized way, while the tracking error compensation maintains the accuracy of the fit.
It can be concluded that, if used wisely, numerical approximations are generally superior to their analytical counterparts, due to their higher degree of freedom. But the lack of mathematical and physical references, which complicates the adaption to parameter changes and requires an enforcement of the otherwise not guaranteed system stability, may prohibits their use in every application.