A simple scheme for the inversion of a Preisach like hysteresis operator in saturation conditions

The modeling of rate-independent hysteresis is crucial in the description of many physical phenomena, e.g., nonlinear elasticity in solids, magnetic characterization of ferromagnets, and, more recently, coupled phenomena in smart materials for energy harvesting, sensing or actuation purposes.^1–4 The paradigm of hysteresis modeling is represented by the Preisach model, i.e., a weighted sum of ideal relays, able to describe rate-independent phenomena in any system fulfilling wiping out and congruency properties,⁵ where, usually, the magnetic field H is assumed as input, while magnetic Polarization J (measured in [T]) is the output field. Despite its complex structure, the Preisach model has been shown to be invertible.⁶ It is stressed that the trivial exchange of field variables defines a relationship with different properties, and this implies that the operator with $H=H(B)$ only provides an approximation of the $B=B(H)$ operator, known as pseudo-compensator.⁷ Therefore, the Preisach operator is able to provide the magnetization of the sample, once the H input history is given. This framework is suitable also when the output variable of interest is B rather then J, at least when the sample does not experience saturation conditions, i.e., negligible variations of J in response to large H− field changes. Unfortunately, also in such conditions, this picture is not suitable in several applications. For example, in numerical modeling of electric circuits or electrodynamics it is preferable to consider the flux density B as input, while magnetic field is assumed as output variable. This requires the inversion of the operator which is usually performed numerically, with all problems related to the rate-independent memory properties of the operator. To address the issue an algorithm able to update the Preisach state by using the output B field rather than H was proposed.⁸ A second way to address the issue, was to employ Prandtl-Ishilinskii (P-I) hysteresis operators.^1,9 They are simplified versions of Preisach operators and are characterized by the availability of the inverse in a closed form. This latter approach seemed effective in many applications and for this reason it is worth to be further investigated. In particular, the availability of simplified versions of the Preisach operator based on P-I operators enables the deployment of hysteresis operators, whatever the input variable is, in numerical electrodynamics, circuit theory and control engineering at least, as specified above, without saturation or in weak saturation conditions. Despite its complex structure, the Preisach model has been shown to be invertible and under simplified assumptions, the expression of the inverse operator can be given in a closed form.¹⁰ When the material experiences saturation conditions, J and B can sensibly differ each other, holding the condition

and the modeling of the B − H characteristic with a single and invertible operator is no more suitable.

In this work, we show that this problem can be addressed through a simple change of the P-I based operator, which, after simple manipulations, is formulated assuming the flux density B as input, as detailed in the sequel.

Let us preliminary recall the Prandtl-Ishilinskii hysteresis operators, which are weighted sums of a continuum of Stop or Play operators, the behavior of which is sketched in Figs. 1 and 2.^1,7 In particular, the classical P-I operator formulation is as follows:

with P_r[x] the Play operator with threshold r, and θ a density function, to be identified on the available measured data. Further, an alternative formulation could be given, as follows:

with S_r the Stop operator, and ξ a suitable density function. From θ and ξ, it is possible to define the following functions:

If ϕ = ψ⁻¹ then the operators π and σ are the inverse of each other, i.e., σ = π⁻¹. From this result,⁹ it is possible to define PI-based operators which admit the inverse in a closed form. A class of generalized P-I operators is defined as follows:¹⁰ 

which admits the inverse in a closed form. In the present analysis we choose $F(x)=x$⁠. The usual assumption exploited when soft Fe-based materials is concerned, is to state the model in an implicit form considering as input the field H_e = H + ν J, usually referred to as effective magnetic field^11,12 Finally, the model takes the form:

with ν is a parameter to be tuned. For Fe samples¹¹ suitable functions are:

where $G(x)$⁠, rather than in past applications¹¹ is not a saturating function but carries a linear term to take into account deep saturating conditions, as those experienced in usual transient behavior of power circuits. Of course, more accurate identification procedures for those functions could be adopted, but this is out the scope of the paper. By exploiting the existence of π and Ginverse, eqn. (7) can be easily arranged as follows:

where the new field J is, from the formal point of view, unknown. In order to bypass the problem, some discussion can be carried out, observing that magnetic polarization J and flux density B do coincide for low fields, while in saturation conditions they differ by the term μH. Further, for large H-field, magnetization approaches the saturation, say J_s, providing a constant value to the mean-field contribution in eqn. (10). Such connection allows a restatement of the previous model as follows:

being $P$ a saturating function. The model defined so far, in conclusion, represents the H − B constitutive model of a soft Fe-based magnetic material and could be easily implemented with low computational cost. The measured and modeled magnetic core characteristic is shown in Fig. 3. The global model’s behavior is quite good since the basic sample magnetic behavior is caught. However a visible but still acceptable discrepancy on the knee of the descending branches is observed. Such error can be reduced by improving the identification of model’s parameters, which is not in the scope of this study. The interesting issue is the handling easiness of the model which is able to provide a fair description of the material’s behavior, with minimal formulation and identification effort.

The model described so far can be easily coupled to a set of ODEs governing equations, where the simplest example is represented by the eqns. given below, modeling a nonlinear RL circuit where the iron core nonlinear inductor describes the primary of a transformer with the secondary in open circuit, as shown in Fig. 4. In the same circuit R models the global coils resistance. The whole circuit models an experimental test bench, made up of an Epstein Frame, a voltage generator and a precision resistor, as detailed in next section.

being Φ, the magnetic flux linked to the primary coil, while i is the primary current. Moreover, Φ(t) = N₁SB(t), with B the flux density, while $i(t)=lN1H(t)$⁠. The parameters N₁, S and l, represent the primary coil turns number, the core cross section and length of the magnetic circuit, respectively. Finally, the relationship between flux density and magnetic field, modeled by eqn. (11) takes the form of eqn. (13) for current and flux variables.

The availability of such tool for magnetic hysteresis modeling allows to investigate the influence of hysteresis (i.e., remanent magnetization and coercive field) on the global circuit’s dynamic behavior.

The proposed experimental setup is shown in Fig. 5 and consists of two fundamental blocks, i.e., the Power Block (PB) and Acquisition and Generation Block (AGB). The former is composed by a Brockhaus-Messtechnik Epstein frame, having the double scope to characterize the ferromagnetic material and to mimic a transformer with the secondary in open circuit conditions. The measured static major loop has been employed to tune the model’s parameters, while minor loops are used to check the effectiveness of the model in describing magnetization phenomena for lower fields, as shown in Fig. 3. The magnetic circuit consists of four strips with dimensions 32.5mm × 0.6mm of a soft ferromagnetic Fe-based alloy, for a total length of the magnetic path of 940mm. The voltage generator is realized by means of a commercial voltage-controlled power supply (Kepco BOP 20-20), which provides a sinusoidal voltage V₁(t) with tunable amplitude and frequency. The AGB consists of a modular National Instruments USB system (NI CompactDAQ 9174 with NI 9215 modules characterized by 16-bit resolution and maximum sampling frequency of 100 kS/s) that allows to backup and manage the experimental data in Matlab environment. Each input channel has its own Analog to Digital converter enabling the simultaneous recording of the supply voltage V₁(t), the supply current I₁(t) and the induced voltage at the secondary winding E₂(t). The AGB output signal V_ref(t) is used to set both the amplitude and the frequency of the primary voltage. Finally, a custom software, developed within the Matlab environment, manages the input variables and data acquisition. In order to check the effectiveness of the proposed model, the experimental bench illustrated so far has been simulated through the circuit model shown in Fig. 4, assuming f = 50Hz, initial phase α₀ = 0[rad], coil resistance R = 0.56Ω, residual flux density B = 0. Simulations were carried out with voltage source amplitudes, E₀ = 5.0 V and E₀ = 6.5 V. The model by its nature, takes trivially into account the remanent magnetization of the material which strongly affects the whole circuit’s dynamics. The results of simulations, compared to the measured current in the circuit are reported in Fig. 6. It is quite evident that the model is able to foresee the correct value of the current peak and the overall behavior of circuit’s variables. A similar performance is illustrated by Fig. 7. The measured loop appears larger than that described by the model. This is due to eddy currents taking place into the sample, which are not described by the (rate-independent) model. The plots also show that the employment in eqn. (14) of the model proposed so far, or using the numerical inversion of eqn. (1) leads to equivalent results. This result is better illustrated in Fig. 8, where the difference between currents provided by the model and by the numerical inversion of the hysteretic characteristic in eqn. (1), is provided in two different working conditions. This allows to conclude that the approach described so far could be deployed for simulation of power circuits to take, for example, into account the inrush current phenomenon due to the insertion of the transformer in the grid, with a good accuracy and low formulation effort and numerical weight. The overall performances could be improved, in particular by improving the description of the hysteresis descending branches, shown in Fig. 3, but this would be the object of a further study.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.