A Self-Sensing Method for Electromagnetic Actuators with Hysteresis Compensation

: Self-Sensing techniques are a commonly used approach for electromagnetic actuators since they allow the removal of position sensors. Thus, costs, space requirements, and system complexity of actuation systems can be reduced. A widely used parameter for self-sensing is the position-dependent incremental inductance. Nevertheless, this parameter is strongly affected by electromagnetic hysteresis, which reduces the performance of self-sensing. This work focuses on the design of a hysteresis-compensated self-sensing algorithm with low computational effort. In particular, the Integrator-Based Direct Inductance Measurement (IDIM) technique is used for the resource-efﬁcient estimation of the incremental inductance. Since the incremental inductance exhibits a hysteresis with butterﬂy characteristics, it ﬁrst needs to be transformed into a B-H curve-like hysteresis. Then, a modiﬁed Prandtl–Ishlinskii (MPI) approach is used for modeling this hysteretic behavior. By using a lumped magnetic circuit model, the hysteresis of the iron core can be separated from the air gap, thus allowing a hysteresis-compensated estimation of the position. Experimental studies performed on an industrial switching actuator show a signiﬁcant decrease in the estimation error when the hysteresis model is considered. The chosen MPI model has a low model order and therefore allows a computationally lightweight implementation. Therefore, it is proven that the presented approach increases the accuracy of self-sensing on electromagnetic actuators with remarkable hysteresis while offering low computational effort which is an important aspect for the implementation of the technique in cost-critical applications.


Introduction
In the field of actuation technology, electromagnetic actuators represent an established solution for several decades. In particular, actuators based on the reluctance principle like solenoid actuators are widely used for positioning tasks such as in hydraulic valves. Their main merits lay in their high forces and large strokes as well as simple construction, which allows the mass production of such actuators at low costs. By applying different geometries for plunger and stator, their force characteristic can be easily adapted to a switching or proportional behavior [1], allowing various application scenarios even under open loop operation. In the past years, the monitoring and control of such actuators raised interest because it allows to drive the actuators in a more precise, energy-efficient and robust manner compared to the open loop operation. Within the concept of Industry 4.0, Condition Monitoring and Predictive Maintenance become additionally interesting for such kind of actuators. An important measurement quantity for such control and monitoring tasks is the actuator position. Usually, this information is measured by mechanical sensors such as encoders, hall sensors or linear variable differential transformers (LVDTs). The disadvantage of such sensor usage is the increased system cost, system size and complexity. Nevertheless, the position information can be also obtained by self-sensing techniques that allow the implementation of control and monitoring techniques even in low-cost 2 of 19 applications [2][3][4][5][6]. Besides the avoidance of position sensors, a self-sensing algorithm in combination with position sensors increases the functional safety of the actuation system due to the obtained redundancy [7].
Works using the inductance information exploit the dependency of that parameter, which is present even at standstill with high sensitivity. The inductance can only be identified when a suitable persistently exciting voltage is applied. Some approaches inject directly a dedicated excitation signal into the actuator [5,7,10]. The demerit of this approach is the presence of acoustic noise and force/torque ripples caused by the injection. Due to this limitation, other works use the current ripple caused by a pulse-width-modulated (PWM) voltage. As most controlled electromagnetic actuators are driven by switching power electronics, a PWM voltage and the resulting current ripple are inherently present inside the actuator and allows a continuous persistent excitation. Many techniques evaluate the current using a numerical derivative [11][12][13][14][15][16][17], which usually leads to a low signal-to-noise ratio (SNR) because of noisy measurements. More sophisticated works make usage of an oversampled current measurement and apply linear regression methods [18][19][20]. While these approaches offer a high SNR and a precise identification, they require high-performance AD converters and computation units with high computational power. The last category of inductance-based techniques uses analog signal pre-processing for increasing the SNR of the measurements while reducing the computational power [21][22][23][24]. Those techniques provide a suitable tradeoff between estimation precision, SNR, and computational costs, therefore allowing the implementation of such techniques in low-cost applications.
The estimation using the incremental inductance has one demerit which lays in its position dependency which in some actuators is non-monotonic and therefore does not allow a unique inversion. This issue leads to ambiguities in the estimation. Unlike the inductance, the eddy current losses exhibit a strictly monotonous behavior over the position. Works evaluating the eddy currents [24][25][26][27] model them as a resistor in parallel to the main inductor, representing the leakage that is caused by eddy currents. All works concede that this parameter has a low sensitivity as well as a high measurement variance, thus they combine the eddy-current based estimation with an inductance-based estimation. In particular, a merging of information can be done by binary rules [25,26] or by means of neural networks [27].
The state-of-the-art works show that a robust estimation of the actuator position over the entire speed range including standstill requires the estimation through the incremental inductance. Being based on the relative permeability and, therefore, the B-H curve of the used ferromagnetic materials, the incremental inductance shows an electromagnetic hysteresis which influences the self-sensing algorithm significantly. Most of the works neglect this hysteresis, which still provides good estimation results when high accuracy is not required or when ferromagnetic materials with small hysteresis such as ferrites are used. Nevertheless, to increase the produced reluctance force, ferromagnetic materials with a high relative permeability µ r and high saturation flux density are commonly used in reluctance actuators [1]. These materials usually show a remarkable hysteresis, leading to high position estimation errors. This makes a hysteresis compensation necessary in self-sensing algorithms when an acceptable estimation accuracy and robustness is required. Hysteresis modeling is well known in the state-of-the-art, as it is mostly applied in controllers to compensate for the hysteresis phenomena [28]. For the case of reluctance actuators, the work [29] shows the modeling of such kinds of actuators by means of the Generalized Preisach model for hysteresis. Similarly, the work [30] applies a Preisach model for observing the magnetic flux in the core and the air gap. Self-sensing with hysteresis compensation is shown in [31] by applying a pseudo-inverse Preisach model on the eddy current resistor. Moreover, the work presented in [7] shows the usage of a modified Prandtl-Ishlinskii model for self-sensing based on the differential inductance that is obtained by a numerical derivative.
This work aims at developing a self-sensing algorithm with hysteresis compensation that provides a good estimation accuracy with high SNR as well as a low computational effort, making an implementation on cost-critical systems feasible. To obtain an estimate with high SNR and low computational power, the Integrator-Based Direct Inductance Measurement (IDIM) technique [23] is used since it allows the estimation of the incremental inductance through the usage of analog pre-processing. Since the incremental inductance shows a butterfly hysteresis, which cannot be modeled by classical hysteresis models, an algorithm will be shown which allows one to transform the butterfly hysteresis into a classical B-H curve characteristic. This is done using an integrator approach that is optimized for implementation. The obtained characteristic is modeled using a modified Prandtl-Ishlinskii (MPI) model that has been proven to provide a good accuracy with low computational power [32] as it is based on the linear superposition of simple elementary hysteresis operators. Finally, a lumped magnetic circuit model allows the separation of the hysteretic iron core reluctance from the non-hysteretic air gap reluctance, thus allowing to isolate and compensate the hysteretic behavior.
This work is divided as follows: after an overview of hysteresis phenomena in reluctance actuators and the definition of important quantities, the MPI model is briefly described in Section 2. Section 3 introduces shortly the IDIM technique for inductance estimation. Section 4 presents the integrator approach that is required to transform the butterfly hysteresis and explains the magnetic equivalent circuit that is needed for hysteresis compensation. An approach to position estimation is proposed. Section 5 shows the implementation and evaluation of the self-sensing approach on an industrial solenoid actuator and compares the performance of the estimator with and without hysteresis compensation. Finally, conclusions about accuracy and implementation effort are drawn before an outlook on further research aspects is provided.

Hysteresis Phenomena and Their Modeling
Electromagnetic actuators such as the reluctance actuator shown in Figure 1 consists of a coil, a back-iron, an air gap with length x, and a movable plunger made of ferromagnetic material.

Electromagnetic Hysteresis
Such soft-magnetic ferromagnetic materials exhibit a characteristic behavior, which is described by their so-called B-H curve. Such a B-H curve is illustrated in Figure 2 and shows the dependency of the flux density B on the magneto-motive force (MMF) denoted as H. The characteristic shows a dependency on the history of the previous working points, which is called electromagnetic hysteresis. For cycling currents between the minimum and maximum current, the major hysteresis loop is obtained. When operating at a particular working point, a small signal excitation causes the presence of a minor loop around that working point.  Figure 2. B-H curve of a ferromagnetic material with the definition of the incremental and differential permeability, adopted from [33].
By applying Maxwell's equations, it can be observed that the B-H curve resembles the characteristic of the flux ψ over the current I, also known as ψ-I curve [1]: with W being the number of windings and A denoting the cross-section of the actuator. Electromagnetic actuators (EMA) are typically driven with a pulse-width-modulated (PWM) voltage provided by switching power electronics. Such a PWM voltage switches at a high frequency and allows to set an almost constant average voltage due to the lowpass behavior of inductive loads like an EMA. Nevertheless, this switching component creates inherently a current ripple, that charges and discharges the inductor around its actual working point. Figure 3 shows an example of such a PWM voltage and its resulting current ripple in an inductive load. The current ripple causes a minor loop in the B-H curve around the actual working point on the major loop.  The permeability in this minor loop represents another quantity than the static permeability and is often described as incremental permeability [33]: with B ∆ being the flux density inside the minor loop. The incremental permeability, therefore, represents the average slope of the minor loop, as depicted in Figure 2. From the incremental permeability, the incremental inductance can be obtained: with ∆I being the current change in the minor loop and ψ ∆ being the flux inside the minor loop. Contrary to this, the differential permeability is described [33] as: and represents the tangent of the actual working point on the major loop. The differential inductance can be obtained as: Unlike the classical hysteresis that is present in a B-H curve, the characteristic of the incremental inductance shows a symmetric butterfly hysteresis [33], which is shown in Figure 4. I L Δ Figure 4. Incremental inductance over current, which shows a butterfly hysteresis.

Modified Prandtl-Ishlinskii Model for Hysteresis Modeling
Hysteresis modeling has been addressed by a large amount of scientific works in the past. The underlying model nature can be generally divided into two groups: physical models and phenomenological models. A short review of different models is provided here, while the work [34] provides a detailed overview. Physical models like the Jiles-Atherton model [35], Carpenter model [36], Globus model [37] and the models described by Stoner-Wohlfahrt [38] and Müller-Aschenbach-Seelecke [39] provide a good physical description of the underlying microscopic effects and working principles. On the other side, phenomenological models proposed by Preisach [40], Bouc-Wen [41] and Coleman-Hodgdon [42] as well as the Prandtl-Ishlinkskii model and its modified form [32,43] provide better modeling and identification results. In this work, the modified Prandtl-Ishlinkskii (MPI) model is chosen due to its good accuracy and its low computational effort. Moreover, it benefits in practical implementation since the identification can be obtained by solving a quadratic optimization problem and the model itself can be inverted by algebraic manipulation [43].
The modified Prantdl-Ishlinskii (MPI) model is based on elementary hysteresis operators, which allow modeling mathematically nonlinearities with memory. In particular, the play operator, which is depicted in Figure 5, is the basic hysteresis operator. The play operator H(x, y, r H ) is based on the function [43] H Energies 2021, 14, 6706 with a parameter named r H , which defines the width of the hysteresis loop. The operator can be applied onto a system with input x, output y and initial state y 0 [43]: To increase the model accuracy, a linear superposition of play operators can be used in the so-called threshold-discrete implementation [43] H with w H being a vector containing the weights, H r H being a vector containing a certain number of play operators with different threshold parameters r H and z H0 being the initial states of these operators.
x y Figure 5. Play operator characteristic.
The so-called superposition operator allows to model memory-free nonlinearities and is shown in Figure 6. It resembles a one-sided dead-zone function and can be described as [43]: with r S being the threshold of the dead-zone. The operator can be applied on an input x: Similar to the play operator, the model accuracy can be increased using a linear superposition of the operators [43]: with w S being the vector of weight and S r S being the vector of superposition operators with different thresholds r S . The modeling of hysteretic nonlinearities can be achieved by the function composition of both operators [43]: with Γ being the model output of the MPI model. This expression can be analytically inverted for model inversion [43]: where the index denotes the inverse weights.
x y r s = 0 r s > 0 r s < 0 Further information about the model consistency and the calculation of the weights as well as the inverse weights can be found in [32]. To train the MPI model, the hysteretic characteristic needs to be excited using a time-changing signal which follows the complete major loop. It is important to remark that the initial conditions need to be zero for identification, so in the case of the EMA, the actuator needs to be fully demagnetized [43]. This can be achieved by using a proper demagnetization signal, which is usually an alternating signal with an initial amplitude that is high enough to reach saturation. The amplitude is then reduced until zero is reached [33]. After the measurements are performed, the MPI model can be trained using a solver for quadratic problems [43].

IDIM Technique for Estimation of the Incremental Inductance
For position self-sensing, the estimation of the incremental inductance is required. In this work, the Integrator-Based Direct Inductance Measurement (IDIM) method is used and briefly described in this section. This technique estimates the incremental inductance by processing the current ripple caused by a PWM-driven electromagnetic actuator. This information is then used in Section 4 for the hysteresis-compensated position estimation of such actuators. For further reference, the IDIM technique is described and derived thoroughly in the works [23,27].
Electromagnetic actuators under PWM operation can be modeled by the electrical equivalent circuit shown in Figure 7. This lumped-parameter model consists of a series resistance R Σ , the position-dependent incremental inductance L ∆ and the parallel resistance R p representing losses such as eddy current or hysteresis losses [27]. The applied voltage is denoted as u(t) while the current in the actuator is represented by i s (t).
By applying the Kirchhoff rules: the electrical model can be mathematically expressed as [27]: where resistive losses are summarized into the total series resistance R Σ [27]: This resistance consists of the copper resistance R s as well as components which are induced during motion. Since a reluctance actuator does not contain permanent magnets, the induced voltage during motion only depends on the applied current i s and, therefore, can be modeled by a resistive component [27]. 0 8 of 20 By applying the Kirchhoff rules: the electrical model can be mathematically expressed as [27]: where resistive losses are summarized into the total series resistance R Σ [27]: This resistance consists of the copper resistance R s as well as components which are induced during motion. Since a reluctance actuator does not contain permanent magnets, the induced voltage during motion only depends on the applied current i s and, therefore, can be modeled by a resistive component [27].
The voltage applied on the phase terminals, shown in Figure 3, is pulse-width modulated (PWM) and can be mathematically described as with U DC being the bus voltage, α being the applied duty cycle and t PW M being the chosen PWM period. The driving of the actuator with such a voltage leads inherently to the presence of a current ripple, as shown in Figure 3. This current ripple contains the information about the lumped parameters R s , L ∆ and R p and can be used for the estimation of those [27]. In particular, the IDIM technique estimates those parameters by processing the sensed current by an analog circuitry, which is depicted in Figure 8. The voltage applied on the phase terminals, shown in Figure 3, is pulse-width modulated (PWM) and can be mathematically described as with U DC being the bus voltage, α being the applied duty cycle and t PW M being the chosen PWM period. The driving of the actuator with such a voltage leads inherently to the presence of a current ripple, as shown in Figure 3. This current ripple contains the information about the lumped parameters R s , L ∆ and R p and can be used for the estimation of those [27]. In particular, the IDIM technique estimates those parameters by processing the sensed current by analog circuitry, which is depicted in Figure 8.
where resistive losses are summarized into the total series resistance R Σ [27]: This resistance consists of the copper resistance R s as well as components which are induced during motion. Since a reluctance actuator does not contain permanent magnets, the induced voltage during motion only depends on the applied current i s and, therefore, can be modeled by a resistive component [27].  By applying the Kirchhoff rules: the electrical model can be mathematically expressed as [27]: where resistive losses are summarized into the total series resistance R Σ [27]: This resistance consists of the copper resistance R s as well as components which are induced during motion. Since a reluctance actuator does not contain permanent magnets, the induced voltage during motion only depends on the applied current i s and, therefore, can be modeled by a resistive component [27].
The voltage applied on the phase terminals, shown in Figure 3, is pulse-width modulated (PWM) and can be mathematically described as with U DC being the bus voltage, α being the applied duty cycle and t PW M being the chosen PWM period. The driving of the actuator with such a voltage leads inherently to the presence of a current ripple, as shown in Figure 3. This current ripple contains the information about the lumped parameters R s , L ∆ and R p and can be used for the estimation of those [27]. In particular, the IDIM technique estimates those parameters by processing the sensed current by an analog circuitry, which is depicted in Figure 8. The voltage applied on the phase terminals, shown in Figure 3, is pulse-width modulated (PWM) and can be mathematically described as with U DC being the bus voltage, α being the applied duty cycle and t PW M being the chosen PWM period. The driving of the actuator with such a voltage leads inherently to the presence of a current ripple, as shown in Figure 3. This current ripple contains the information about the lumped parameters R s , L ∆ and R p and can be used for the estimation of those [27]. In particular, the IDIM technique estimates those parameters by processing the sensed current by an analog circuitry, which is depicted in Figure 8.
From the sensed current signal, the fundamental component is subtracted by means of a sample-and-hold (S/H) stage and a subtraction stage. In a next step, the offset-eliminated signal is fed to a fast resettable integrator, which can be reset at every PWM instant. The Figure 8. Electronic circuitry implementing the IDIM technique, adopted from [27].
From the sensed current signal, the fundamental component is subtracted using a sample-and-hold (S/H) stage and a subtraction stage. In the next step, the offset-eliminated signal is fed to a fast resettable integrator, which can be reset at every PWM instant. The trigger signal r(t), which controls the S/H stage and the integrator, is generated by a timer and has the following specifications with the timings The parameter t r is a tunable waiting time, which should be long enough to ensure a proper sampling in the S/H stage as well as a full reset of the integrator. Finally, the output of the integrator Q(t) can be sampled by an analog-digital (AD) converter at the time instants t + s , t + e , t − s and t − e and can be digitally processed in a microcontroller. Figure 9 shows schematically the current ripple and the integral together with the timings. The output of the circuit can be written as where the generic time t x stands for t + s in the case the positive voltage pulse of the PWM is being observed and t − s in case the negative PWM pulse is being observed. Inserting the electrical model shown in Equation (18) as well as the PWM voltage from Equation (20) yields to Integral signal Q(t) Trigger signal r(t) Figure 9. Applied PWM voltage u(t), resulting current i s (t), trigger signal r(t) as well as output of the IDIM circuit Q(t), schematically shown for the defined trigger times, adopted from [27].
This equation can be solved analytically by integration. Nevertheless, during the time of integration, it is assumed that the position-dependent parameters are constant. This assumption holds since usually the mechanical time constant of electromagnetic actuators is considerably larger than the electrical time constant and the chosen PWM frequency. The initial condition of the integral can be set to zero since the reset mechanism of the integrator forces the integration to start at an initial value of zero. Furthermore, it can be assumed that R Σ R p , since electromagnetic actuators are designed to have low resistive losses represented by R Σ and low eddy current losses, which leads to a high value of the parallel leakage resistor R p . Under these assumptions, the expression of the integral can be solved and simplified to [27]: for the positive PWM voltage pulse and with the matrix A being the matrix of the measurements: For the estimation of the incremental inductance, the system of linear equations can be solved [27]:

Hysteresis Compensation Based on Magnetic Circuit Model
In the previous section, an incremental inductance estimation technique has been shown, which can be used for position estimation. Nevertheless, as already shown in Figure 4, the incremental inductance exhibits a hysteretic behavior, which leads to significant position estimation errors in case the hysteretic nature is not considered. Section 2 gives a short overview of existing hysteresis models and presents the MPI model as a suitable approach for hysteresis modeling. However, the hysteretic nature of the incremental inductance has a characteristic butterfly behavior, which cannot be modeled by standard hysteretic models since its characteristic does not follow the standard B-H curve hysteresis. Therefore, the butterfly characteristic first needs to be transformed into a B-H curve characteristic. Since the incremental inductance is based on the difference quotient of the flux inside a minor loop, an integration approach will be shown in the following, which allows us to estimate the flux inside a minor loop, which from now will be called incremental flux ψ ∆ .
Since the current inside the actuator exhibits a ripple and consequently changes during a PWM period, a mean current I is defined that is constant during one PWM period and represents a constant working point on the B-H curve, at which the minor loop occurs. Considering Equation (4), the incremental flux can be obtained through integration of the incremental inductance over the current: Continuously evaluating this equation allows us to track the incremental flux present inside a minor loop at a certain working point on the B-H curve. As it can be seen in Section 5, the incremental flux estimated by this technique has a B-H curve-like hysteretic characteristic which can be modeled by the MPI model. However, using the integration approach has several implementation issues. First of all, the initial condition of the integral Equation (30) has to be known exactly. Since knowing the initial point on the B-H curve under the consideration of remanence effects can be cumbersome, it is desirable to fully demagnetize the actuator using a demagnetization signal to have an initial point at the origin of the B-H curve. Since such an initialization procedure anyways needs to be implemented for the MPI model, no further implementation effort is needed. Secondly, Equation (30) is usually implemented on a time-discrete system such as a microcontroller, thus a numerical integration method needs to be applied. For a good trade-off between accuracy and implementation effort, the trapezoidal rule is used [44]: where j stands for the measurement sample. Finally, the practical implementation of such integrators suffers from continuous drifting since measurement and calculation errors sum up until the integral diverges. This is a major problem for this approach. Therefore, the integrator drift needs to be actively compensated using the knowledge of well-known points on the B-H respectively the ψ-I curve. Because of the hysteresis, this must be done at unique operating points of the actuator. For the maximum current, the actuator is fully in saturation and no hysteresis is present. Thus, when the actuator reaches the points [−I max , I max ], the integrator can be forced to the value of the incremental flux. For increasing the accuracy, the same can be applied to the origin, when the moving direction of the current is known. Applying the offset compensation allows us to limit the drift as every quarter cycle of the major loop the integrator output gets corrected.
To compensate for the effect of hysteresis during self-sensing, the hysteretic part of the incremental inductance needs to be separated from the non-hysteretic position-dependent part. This can be achieved by modeling a magnetic equivalent circuit, which has been proposed by many works such as [7,[29][30][31]. The magnetic circuit shown in Figure 10, consists of an MMF source as well as a magnetic reluctance for the iron core R mFe and the position-dependent air gap reluctance R mx . Only the iron reluctance exhibits a hysteretic characteristic and therefore, hysteretic and non-hysteretic components can be separated. In this model, leakage fluxes and parasitic air gaps are not considered.
initialization procedure anyways needs to be implemented for the MPI model, no further implementation effort is needed. Secondly, Equation (30) is usually implemented on a timediscrete system such as a microcontroller, thus a numerical integration method needs to be applied. For a good trade-off between accuracy and implementation effort, the trapezoidal rule is used [44]: where j stands for the measurement sample. Finally, the practical implementation of such integrators suffers from continuous drifting since measurement and calculation errors sum up until the integral diverges. This is a major problem for this approach. Therefore, the integrator drift needs to be actively compensated using the knowledge of well-known points on the B-H respectively the ψ-I curve. Because of the hysteresis, this must be done at unique operating points of the actuator. For the maximum current, the actuator is fully in saturation and no hysteresis is present. Thus, when the actuator reaches the points [−I max , I max ], the integrator can be forced to the value of the incremental flux. For increasing the accuracy, the same can be applied to the origin, when the moving direction of the current is known. Applying the offset compensation allows to limit the drift, since every quarter cycle of the major loop the integrator output gets corrected. In order to compensate the effect of hysteresis during self-sensing, the hysteretic part of the incremental inductance needs to be separated from the non-hysteretic positiondependent part. This can be achieved by modeling a magnetic equivalent circuit, which has been proposed by many works such as [7,[29][30][31]. The magnetic circuit shown in Figure 10, consists of a MMF source as well as a magnetic reluctance for the iron core R mFe and the position dependent air gap reluctance R mx . Only the iron reluctance exhibits a hysteretic characteristic and therefore, hysteretic and non-hysteretic components can be separated. In this model, leakage fluxes and parasitic air gaps are not considered. A reluctance represents a resistance in the magnetic path and can be calculated as [1]: where v is the magnetic voltage. Under the knowledge of the number of windings W, the magnetic voltage can be calculated from the applied current I. For sake of brevity and comprehension, the number of windings will be now assumed equal to one.
The total reluctance in the magnetic path can be summarized as A reluctance represents a resistance in the magnetic path and can be calculated as [1]: where v is the magnetic voltage. Under the knowledge of the number of windings W, the magnetic voltage can be calculated from the applied current I. For sake of brevity and comprehension, the number of windings will be now assumed equal to one. The total reluctance in the magnetic path can be summarized as Only the part of the air gap reluctance is position-dependent and non-hysteretic, thus it can be described by an invertible function f (x): When the air gap equals zero (x = 0), the air gap reluctance disappears and therefore, the total reluctance only consists of the iron core reluctance: Under this condition, the core reluctance can be easily identified since no position dependency is present. At zero position, the identification of the MPI model can be done, which yields an estimate of the incremental flux under zero air gap. By knowing the incremental flux under these circumstances as well as the current, the reluctance can be estimated according to Equation (32) to: While the part of the core reluctance gets estimated in an online way, the IDIM technique, together with the integrator approach shown in Equation (30), allows to measure the actual value of the total reluctance: By applying Equation (33), the air gap reluctance can be calculated as the difference between the two reluctances and finally the position x can be obtained by inverting the function f : Such a function f is found by a pre-identification of the air gap reluctance and modeled using polynomial or physical models. The function must be strictly monotonous to be inverted. Theoretically, the function only depends on the air gap x. In practice, modeling errors in the MPI model as well as the negligence of parasitic leakage fluxes leads to a current dependency of the air gap reluctance. To increase the accuracy of the estimation and consider the number of windings W, the current dependency should be considered during modeling. Thus, f is a function of position x and current I. Figure 11 summarizes schematically the obtained self-sensing algorithm. it can be described by an invertable function f (x): When the air gap equals zero (x = 0), the air gap reluctance disappears and therefore, the total reluctance only consists of the iron core reluctance: Under this condition, the core reluctance can be easily identified since no position dependency is present. At zero position, the identification of the MPI model can be done, which yields to an estimate of the incremental flux under zero air gap. By knowing the incremental flux under these circumstances as well as the current, the reluctance can be estimated according to Equation (32) to: While the part of the core reluctance gets estimated in an online way, the IDIM technique, together with the integrator approach shown in Equation (30), allows to measure the actual value of the total reluctance: By applying Equation (33), the air gap reluctance can be calculated as the difference between the two reluctances and finally the position x can be obtained by inverting the function f : Such a function f is found by a pre-identification of the air gap reluctance and modeled using polynomial or physical models. The function must be strictly monotonous in order to be inverted. Theoretically, the function only depends on the air gap x. In practice, modeling errors in the MPI model as well as the negligence of parasitic leakage fluxes leads to a current dependency of the air gap reluctance. In order to increase the accuracy of the estimation and consider the number of windings W, the current dependency should be considered during modeling. Thus, f is a function of position x and current I. Figure 11 summarizes schematically the obtained self-sensing algorithm.

Experimental Results
In the following, the technique derived in Section 4 is experimentally verified. First, the performance of the MPI model will be validated when the actuator is at zero position. Then, the position-dependent characteristic of the actuator will be shown, modeled, and used for self-sensing. The results are then compared to a high-precision encoder. Figure 12 shows the experimental test-bench. The actuator under test is coupled to a linear positioning table with a resolution of 200 nm that allows us to block the plunger up to forces of 50 N. Additionally, the table with its high precision encoder is used as a validation sensor. The actuator under test is a switching solenoid actuator from the Type GTC A 40 from Magnet-   Figure 13 shows the experimentally obtained incremental inductance when the air gap is zero (x = 0mm). It resembles the butterfly characteristic shown in Figure 4. Nevertheless, the measured curve is not strictly symmetric compared to the idealized behavior. This phenomenon is also observable for different materials [33] and can be explained by a preferred magnetization direction, which can be caused e.g., by fabrication processes. The obtained characteristics are then fed to the integrator shown in Equation (30) to calculate the incremental flux at a certain working point. Figure 14 shows the obtained incremental flux over the current. It can be seen that the integral of the butterfly-shaped incremental inductance resembles a B-H curve with a saturation area at higher currents. The MPI model explained in Section 2 is used to model this behavior and is identified using an excitation signal and quadratic optimization. As a good trade-off between accuracy and computational effort, the model is parameterized with 6 play-operators and 10 superposition-operators. A further increase in the number of operators does not con-siderably decrease the identification error. The online estimation using the MPI model is shown in Figure 14 along with the original measurement. The model achieves a good performance especially at low and high currents with a maximum relative error of 1.6%. In the following, only the results for positive currents are used since reluctance actuators operate in a unipolar driving mode.  Figure 15 shows the measured incremental inductance for different positions of the actuator. It can be seen that for zero air gap (x = 0 mm), the characteristic from Figure 13 is obtained. For increasing air gaps, the total value of the inductance decreases, and the width of the butterfly hysteresis decreases. For the maximum air gap (x = 8 mm), the ferromagnetic core is almost removed from the coil and therefore only little hysteretic behavior can be observed.  Figure 16 illustrates the calculated total magnetic reluctance R mtot over the current for different positions. The total reluctance is obtained by applying Equation (37) on the measurements shown in Figure 15. It can be seen that the total reluctance is strongly hysteretic and position-dependent. Its absolute value rises for increasing air gap sizes. For certain positions, the curves overlap due to the hysteresis, making a unique estimation of the position impossible. Based on Equation (38), the air gap reluctance can be analytically calculated. The result shown in Figure 17 is the hysteresis-compensated air gap reluctance, which is theoretically only positioned dependent. Compared to the total reluctance shown in Figure 16 Another important aspect is the presence of ambiguities in the obtained inductance characteristic, which are highlighted by an ellipse in Figure 17. Those ambiguities are caused by the behavior of the B-H curve of the material. In particular, there is an extremum in the incremental permeability, which leads to a non-monotonic behavior. Many works have reported this behavior [1,[24][25][26][27]31] and it can be either addressed by the avoidance of those working points or by the consideration of a different parameter to be identified, such as the eddy current resistor R p . Since this ambiguity problem is not within the scope of this paper, the working points showing this phenomenon are not considered in this work and are neglected.

Position Estimation Results Using Self-Sensing with Hysteresis Compensation
By avoiding the ambiguous working point in the ellipse, the remaining points are unique over the position range. Thus, the mapping is strictly monotonous and can be inverted. In particular, the characteristic is inverted and identified using a polynomial model of an order of 3 for the reluctance and an order of 3 for the current. The resulting function f −1 (R mx , I) is implemented for validation. Figure 18 shows the performance of the obtained self-sensing algorithm compared to the reference sensor in case the hysteresiscompensated air gap reluctance is used and in case the non-compensated total reluctance is used. Additionally, the actual current I of the working point and the relative error referred to as the maximum stroke of the actuator are presented. In the experiment, current loops from zero to maximum current are applied, while the ambiguous working points are left out. The current gets applied with a finer step size around the origin to prove the model's performance at the initial curve. It can be generally seen that using the hysteresiscompensated algorithm significantly increases the accuracy of the estimated position. For the uncompensated estimator, a maximum relative error of 49.23% is achieved while the maximum relative error for the compensated algorithm amounts to 9.02%. Especially for small positions, the accuracy is increased since the influence of the iron core at small air gaps is more significant. For larger air gaps, both errors decrease with one of the compensated algorithms being smaller. For small positions, the deviation of the compensated algorithm also shows a considerable error of 9% which is mainly caused by the ambiguities and low sensitivity of the incremental inductance at these working points. It is visible that the uncompensated algorithm shows repetitive peaks in the error diagram when a certain current on the major loop is reached. In the compensated algorithm, these peaks are less evident.

Conclusions
This work presents a self-sensing algorithm for electromagnetic actuators with hysteresis compensation. The self-sensing approach is based on the incremental inductance, a parameter that has a good sensitivity over the position but exhibits a strong butterfly hysteresis. To model this hysteretic butterfly characteristic, an integration approach is shown and practically implemented using a drift compensation. The obtained characteristic shows the original behavior of a standard B-H curve hysteresis, thus it can be modeled by standard hysteresis models. In particular, the modified Prandtl-Ishlinksii (MPI) model is chosen due to its accuracy as well as low implementation and computation effort. By using a magnetic equivalent circuit, the hysteretic behavior can be theoretically separated from the characteristic of the air gap, therefore allowing a position estimation without the influence of the hysteresis.
During the design of the method, special focus is laid on the computational effort. The Integrator-Based Direct Inductance Measurement (IDIM) technique that is used for the estimation of the incremental inductance performs only 7 measurements per PWM period and allows the calculation of the inductance value within a closed equation. Furthermore, the MPI model allows a computationally lightweight calculation since it is based on simple hysteresis operators which are linearly superposed. This allows the implementation of the technique on low-cost systems which are preferred in the case of solenoid actuators.
Experimental results on an industrial switching actuator validate the concept of the proposed algorithm. The discussed approach can reduce the relative error on the actuator under test from 49.23% in case the hysteresis is not compensated to 9.02% with hysteresis compensation. Nevertheless, a remaining hysteresis is still present in the estimation of the air gap reluctance, which stems from identification errors present inside the MPI model as well as model uncertainties in the magnetic equivalent circuit. Therefore, a further research goal is to increase the model accuracy, e. g. with the Generalized Prandtl-Ishlinskii model [45]. During the experiments, ambiguities are present which are not caused by the hysteresis, but by the non-monotonous behavior of the incremental permeability. To address the problem of finding a unique position estimate at those working points, another parameter like the eddy current resistor R p needs to be taken into account. While the work [27] offers a possible solution, it needs to be combined in the future with the hysteresis compensation presented in this work. Further investigations are required concerning the temperature-dependency of the incremental inductance which can affect the position estimation at varying temperature ranges significantly.

Patents
The IDIM technique has been submitted for patenting and references can be found under [22].