Open access peer-reviewed chapter - ONLINE FIRST

Preisach Hysteresis Model. Some Applications in Electrical Engineering

Written By

Alfredo Bermúdez, Dolores Gómez and Pablo Venegas

Reviewed: 21 July 2021 Published: 03 November 2021

DOI: 10.5772/intechopen.99590

Modern Permanent Magnets - Fundamental and Application IntechOpen
Modern Permanent Magnets - Fundamental and Application Edited by Dipti Ranjan Sahu

From the Edited Volume

Modern Permanent Magnets - Fundamental and Application [Working Title]

Prof. Dipti Ranjan Sahu

Chapter metrics overview

260 Chapter Downloads

View Full Metrics


In this chapter we recall the well-known hysteresis Preisach model, widely employed in the area of magnetism. Some applications of this model in electrical engineering are also described, with a specific focus on the estimation of electromagnetic losses in electrical machines, the simulation of magnetization-demagnetization processes arising in magnetic particle inspection, and the mathematical modeling of batteries for electric vehicles.


  • magnetic hysteresis
  • Preisach model
  • electromagnetic losses
  • magnetic particle inspection
  • batteries for electric vehicles

1. Introduction

Hysteresis is a very complex nonlinear behavior affecting many physical phenomena. Systems affected by this behavior are characterized by the fact that the way they evolve in response to a stimulus depends not only on the cause of the stimulus, but also on the preceding states of the system. Thus, the same instantaneous values of the input can give different outputs depending on the history of the input applied, which gives rise to a relationship that is not only nonlinear but also multivalued, making it very difficult to model and control. This memory-based property is found in various areas of science and engineering such as mechanics, biology, economics and multiphase flow in porous media or magnetism, among others. It was precisely in the latter that the term was initially coined since most ferromagnetic materials exhibit hysteresis [1]. This means that when the material is subjected to the application of a magnetic field, the magnetic induction reached at each point of the material depends not only on the intensity of the applied field at a given instant of time, but also on its previous magnetic history. Despite the difficulties involved in its study, having a good hysteresis model is essential in a multitude of applications in the field of electrical engineering; an example of this is the estimation of energy losses in electrical machines.

The variety of systems with hysteresis is very broad and there is a large amount of bibliography devoted to hysteresis models in different communities of physicists, engineers, and mathematicians, where great efforts have been devoted to their development and analysis. As a result, nowadays there are several hysteresis models, ranging from simple to complex, each of them valid in some specific situation; see, e.g., the monographs and reviews [2, 3, 4, 5].

The mathematical approach tries to deal with this phenomenon under a common mathematical framework, but this is not always possible. We highlight the monograph by Krasnosel’skii˘ and Pokrovskii˘ [2], who have conceptually introduced the notion of hysteresis operator and carried out a systematic analysis of its properties. More recently, research on hysteresis models as well as their coupling with partial differential equations has been progressing; see, among others, the works by Visintin [5, 6] and Brokate and Sprekels [7]; from the physical point of view, Mayergoyz [8] and Bertotti [9] are the classical references.

Hysteresis models can be classified into two main classes according to whether or not the system response depends on the input velocity. In static or rate-independent models, output values are not affected by the velocity of the input but only by the values in the input range and by the order in which these values have been attained. Consequently, the model cannot reflect the frequency or waveform dependence of the field. This is a characteristic of the classical hysteresis phenomenon and in this sense some authors [6] consider hysteresis as a rate-independent memory only. On the opposite side are the so-called dynamic or rate-dependent models, which account for the effect of the velocity of changes in the applied input.

The so-called classical Preisach operator [10] is a rate-independent model originally designed to model hysteresis of ferromagnetic materials which is based on physical assumptions derived from the concept of magnetic domains [8, 9, 11, 12]. It is suitable for modeling scalar hysteresis and for this purpose, the magnetic field strength H is used as the input variable while the magnetization M (or the magnetic flux density strength, B) acts as the output variable. Its main advantage is the ability to describe not only the major hysteresis loop, but also inner loops and other complex characteristics of the magnetization processes. There are several extensions of this model that consider the rate of change of the stimulus to be able to consider phenomena where this factor influences the magnetization. These are generically referred to as dynamic Preisach models (see, for instance [13, 14, 15, 16, 17]).

Currently, the Preisach-type operators [18] are used not only in the area of magnetism, but also to describe hysteresis phenomena in other fields as diverse as fluid flow in porous media [19], elastoplasticity [20], solid phase transitions [7], shape memory alloys [21], or biology [22], to name a few.

The mathematical and numerical treatment of PDEs with hysteresis operators is a challenging issue because, despite the importance of the topic, still few results are available. Due to the strong interdisciplinary character of hysteresis phenomena, the interest for their study is continuously increasing in a great variety of applications. During the last few years, the authors of this chapter have been working on the mathematical modeling, numerical analysis and computation of PDE hysteresis-related problems motivated by real applications, with special emphasis on the hysteresis of ferromagnetic materials. This work is a survey of previous co-authored studies [17, 23, 24, 25], where we have intended to provide original steps into the mathematical and numerical treatment of parabolic problems with hysteresis. In all of them, the Preisach model was considered as hysteresis operator. The problems were addressed in different aspects: mathematical analysis, numerical methods, until convergence of approximate solutions and computational results to confirm the theory and to compare with experimental data. The rigorous mathematical framework or the numerical analysis are not included here, and we will be concerned with the results obtained for different industrial problems.

The chapter is organized as follows. In Section 2, we start by recalling the basic properties of the hysteresis operators. The well-known Preisach hysteresis model, both in the rate-independent and rate-dependent cases, is introduced in Section 3. Finally, in Section 4 we show some examples of the applications of the Preisach model in electrical engineering. These examples include the computation of hysteresis losses, the numerical simulation of magnetization and demagnetization processes in ferromagnetic workpieces and the simulation of batteries for electric vehicles.


2. Hysteresis operators: Basic properties

In this section, we summarize the basic knowledge of scalar hysteresis operators based on the description given in [6].

Most systems experiencing hysteresis phenomena display hysteresis loops. To illustrate this in a simple setting, let us consider a system like the one depicted in Figure 1 (left), whose state is characterized by two time-dependent scalar variables, u and w. Let us also suppose that u is the independent variable, and thus the evolution of w depends on u. Such diagrams implicitly show the memory effect inherent to the hysteresis: at any time t>0, wt depends on the evolution of u in 0t, rather than only on ut. In many cases, not all the information in 0t is used to compute the value of w but only the one which has not been wiped-out from the memory (see Section 3.1). In other words, the state of the system at present time depends on the previous history of its state. Consequently, the same instantaneous values of the input can give different outputs and when describing the relation between w and u in the uw plane a multi-branch nonlinearity appears with branch to branch transitions (see Figure 1, left).

Figure 1.

Hysteresis loop (left) and rate independence example (right).

According to Ref. [6], two main characteristics of hysteresis phenomena can be distinguished: causality and rate independence. To clarify these concepts, we consider the uw relation introduced above. Causality means that, at any time t, the value of wt only depends on the previous evolution of u, namely, wt depends on u0t. On the other hand, rate independence is translated by the condition that, at any instant t, wt depends just on the range of function u:0tR and on the order in which the values of u have been attained. In other words, w is independent of the velocity of u. An example is given in Figure 1 (right) where it is shown that three different functions u:0tR lead to the same wu curve.

We notice that, even in the most typical hysteresis phenomena, like plasticity, ferroelectricity or ferromagnetism, the memory effect is not completely rate-independent, since viscous-type effects are coupled with hysteresis. However, in the applications presented in the following sections, this effect is neglected, so that the rate-independent component prevails.

To provide a functional setting for hysteresis operators, we first notice that, at any instant t, wt will depend not only on the previous evolution of u (i.e., on u0t), but also on the “initial state” of the system. Due to the memory dependence of hysteresis processes, additional information is needed to make up for the lack of history when the process begins. This initial information must represent the “history” of the function u before t=0. Hence, not only the standard initial value u0w0 must be provided. In general, a variable ξ containing all the information about the “initial state” is considered. This can be expressed, for instance, as follows:


In the case of partial differential equations, it is necessary to define an operator F acting between suitable function spaces involving the space variable. Let Ω̂ be an open subset of RNN1; given a hysteresis operator F˜, we introduce, for any u:Ω̂×0TR and any ξ:Ω̂Y, with Y a suitable space, the corresponding space dependent operator F as follows (see [23] for a rigorous mathematical setting)


We notice that operator F˜ is local in space, i.e., the output Fuξxt depends on ux0t, but not on uy0t for yx.

The general setting that we have discussed so far, and the majority of the hysteresis models proposed in the literature, are scalar. Thus, they can be applied only to model unidirectional inputs. However, in many applications, the hysteresis is characterized by a vector input ut and vector output Fut; thus, vector hysteresis is encountered. This is the case, for instance, of electric devices such as actuators, transformers, or rotating machines in which the direction of the magnetic field is apriori unknown. The properties of vector hysteresis are often very different from the properties of scalar hysteresis, and the derivation of a general model of vector hysteresis remains an open question. Useful references on vector hysteresis models and their features are given in [12, 26, 27, 28, 29, 30].


3. The Preisach hysteresis model

The Preisach model [31] is the most common and probably the most important model to represent magnetic hysteresis phenomenona in the literature. It was originally proposed by the physicist F. Preisach [10] in 1935 in the context of ferromagnetism and later the formalism was more broadly generalized to describe hysteretic behaviors in different fields [6, 8]. Nowadays, it is recognized as a fundamental tool for describing hysteretic systems with complex behaviors.

To describe the hysteresis, the classical Preisach model assumes that each particle of the material has an associated elementary hysteresis operator, called relay operator, which characterizes its state. Thus, the system can be modeled as a sum of elementary relays, whose calculation can be performed in parallel, weighted by a distribution function μ, called Preisach density function, determining, in a certain sense, the local influence of each relay on the global model. Then, the relationship between the input ut and the output wt is given by means of the integral:


where μ is the weight function, ξ contains the information about the “initial state” of every point of the domain and γρ1ρ2 is the relay. The function μ is a non-negative function with compact support in the Preisach half-plane with ρ1ρ2 that identifies the system. In practice, it can be analytically approximated by, e.g., Lorentzian or Gaussian distributions (see [9]) or, as it will be explained in Section 3.3, estimated from experimental measurements using the so-called Everett function (see [8]). Depending on the properties of these relays and the relationships between them, two types of Preisach models are defined: the classical (or static) Preisach model based on the rate-independent relay, and the dynamic Preisach model introduced in [14] based on the dynamic relay.

3.1 The static or rate-independent scalar Preisach operator

In this model, the elementary relay operator is represented by a rectangular loop in the input–output plane with “up” and “down” switching values (see Figure 2 for an example). Each elemental relay, here denoted by hρ, is associated to a point ρρ1ρ2R2, with ρ1<ρ2, and, for any ρ, it has two states: “up” (hρ=1) and “down” (hρ=1) and two switching thresholds: ρ2 is the switch-up threshold, and ρ1 is the switch-down one.

Figure 2.

Scalar relay.

Formally, for any continuous function uC0T and initial condition ξ11, hρuξ is a real function defined in the time interval 0T such that,


Then, for any t0T, we consider the set Xut{τ(0,t]:uτ=ρ1orρ2}. This set keeps account of the previous time instants in which u presents the thresholds ρ1 or ρ2. Next, we define


If u0=0, then the following initial condition can be considered


When working with ferromagnetic materials, this initial configuration results in zero magnetic induction; thus, the material is often said to be “demagnetized” or in a “virginal” state. We remark that, in this situation, hρ can only take values ±1 and it changes instantaneously from its last value depending on the previous evolution of the system; more precisely, when u reaches the threshold ρ2 from below, it “switches-up” to value 1, and when it attains ρ1 from above, it “switches-down” to 1. Therefore, hρ is not a local in time mapping: hρuξρt not only depends on ut but on its past history. The classical Preisach operatorFS is then defined as


where we recall that variable ξρ contains all the information of the initial state at each point of the system and p>0 denotes the density function.

The integral in (7) is calculated in the so-called Preisach triangle (see Figure 3, right):

Figure 3.

An arbitrary input ut (left) and its corresponding map on the Preisach triangle (right).


being ρ0ρ0 the thresholds setting the minimum and maximum values of ut. By using this triangle, it is possible to provide a geometric interpretation of the Preisach operator. The following is a summary of the most important aspects. For an extended description, the reader can consult Refs. [8, 23]. For a given uC0T and any time t, the triangle is split into two sets (one of them possibly empty): Su+t and Sut containing the points ρ1ρ2 for which the associated relays hρu are positive or negative, respectively, i.e.,


(see Figure 3). The interface Lut between the two sets is a staircase line with vertices having coordinates ρ1ρ2 respectively coinciding with the local minimum and maximum values of u at previous time instants. At time t, Lut intersects the line ρ1=ρ2 at utut. Lut moves up as ut increases and from right to left as ut decreases (see Figure 4).

Figure 4.

Staircase line Lut moving right to left (left) and moving up (right).

Then, the Preisach operator (7) can be equivalently expressed as


which can be readily calculated once Su+t is obtained. Since the second integral in (10) is constant in time, the computation of the output wt depends on the effective computation of the first integral.

From (10), is clear that wt depends on the shape of the interface Lut and the latter is determined by the extremal values of ut, at previous instants of time. Notice that not all extremal input values are needed. In fact, considering the dependence on Lut, the Preisach operator exhibits a wiping-out property. This means that every time the input reaches a local maximum ut, Lut erases (“wipes out”) the previous vertices of the staircase having a ρ2 value lower than the value ut. Similarly, when an input reaches a local minimum ut, the memory curve wipes out all previous vertices with a ρ1 greater than the value of ut. Therefore, only dominant extreme values contribute to the model, while all the other local extreme of the input are eliminated. Figure 3 illustrates this in a particular case.

From the description above, we conclude that three basic steps characterize the model application:

  • Identification, from experimental data, of the Preisach density p or a suitable auxiliary function.

  • At each time step the set Su+t has to be characterized from the local maxima and minima of the input u. This allows us to update the memory state of the system.

  • The value of wt is obtained by computing the integral of p in the domain Su+t.

The first and third steps, namely, the identification of the Preisach function and the output calculation are difficult issues. To obtain an efficient procedure for the computation of wt, we use, as usual, the so-called Everett function identified with the First-Order Reversal Curves (FORC).

3.2 The dynamic or rate-dependent Preisach operator

In [14, 32], Bertotti introduces a rate-dependent generalization of the classical Preisach model, aiming to take into account the rate of change of the input u. This new operator, termed as dynamic Preisach operator, overcomes some limitations of the classical model, in particular, the fact that the frequency of the input was not reflected in the shape of the hysteresis diagram.

Contrary to the classical relay operator, where only the two states ±1 are possible, the dynamic relay can attain all the intermediate values in the interval 11, switching at a finite velocity which is assumed to be proportional to the difference between ut and the threshold values ρ1 and ρ2 (see Figure 5, right panel). The proportionality factor is a parameter k depending on the material.

Figure 5.

Input function ut=150sin2πft (left) for a fixed frequency f = 20 Hz. A dynamic relay and its slope are presented for switching values ρ1ρ2=50,100 and initial state ξ=1. The center panel shows the relay with slope k=50 whereas the right panel shows the corresponding diagram uηρ/t. Black arrows represent increasing values of u while blue arrows represent decreasing values.

In a formal manner, and inspired by [14], for a fixed ρ=ρ1ρ2R2, ρ1<ρ2, we define the dynamic relay operator ηρ such that, for any u and ξ11, ηρuξ:0T11 is the unique function y such that yt11 be the solution of the (nonlinear) Cauchy problem:


In the expressions above, the standard notations x+=maxx0andx=maxx0, have been used; thus x=x+x. We remark that in this dynamic model the initial state ξ can attain all the values in the interval 11. Moreover, definition (11) is consistent with that given in [2] because cases third and fourth in (11) cannot occur. Nevertheless, if one wants to perform a proper mathematical analysis, all cases must be considered. This mathematical analysis is out of the scope of the chapter; the interested readers are referred to reference [17].

Next, we consider two examples aiming to give an idea of the dynamic curve utyt. For the first one, illustrated in Figure 5, a sinusoidal input ut=150sin2πft is considered for frequency f=20 Hz and initial condition ξ=1. This figure also shows the curve utyt parameterized with respect to the time variable when the relay ηρ is characterized by switching values ρ1ρ2=50,100 and slope k=50. The right panel shows the slope ηρ/t, which represents Eq. (11). From the diagram it can be seen that when ηρ reaches value 1 or 1, the derivative is again equal to zero, as in the interval between switching values ρ1ρ2. The second example, illustrated in Figure 6, shows the dynamic relay when different frequencies of the input and slopes k are considered. From the center panel, it can be seen that the dynamic relay is rate-dependent. On the other and, the right panel shows curve utyt for slopes k=102 (dash-dotted line), k=104 (dashed line) and k=108 (solid line). Notice that the solid line is an approximation to the discontinuous static relay of the classical Preisach model (see Figure 2).

Figure 6.

On the left panel, input function ut=150sin2πft+75. The center graph shows the uy curve corresponding to the dynamic relays with k=50 and initial state ξ=1 for frequencies f=50, 500, 5000 Hz (dashed, dash-dotted and solid line, respectively). On the right panel, the uy curve is depicted for f= 5000 Hz and k=102,104,108 (dash-dotted, dashed and solid line), respectively.

From previous considerations, the dynamic Preisach operator FD can be defined as


Notice that, if ηρ is replaced by hρ the classical, rate-independent, Preisach model is obtained. Figure 7 shows the classical and dynamic relay configurations with respect to the input u. The Preisach triangle is characterized by ρ0=300, a constant k=50 and the demagnetized state as initial condition.

Figure 7.

Input function ut (left) defined in [0, 0.0045]. The isolines represent the corresponding dynamic relay values ηρ with k=50 and the classical relay hρ, at t=0.003 (center) and t=0.0045 (right).

3.3 Computation of the Preisach operator based on Everett function

In [8], Mayergoyz developed an approach for the computation of the Preisach model that does not require the Preisach density function p but the Everett function, which describes the effect of p on the hysteresis operator. The Everett function is obtained from the First-Order Reversal curves by a procedure described below. A FORC diagram is generated from a class of minor hysteresis loops referred to as first-order reversal curves. A FORC branch Fρ2 is associated to the threshold ρ2. The input u rises up from the “reset” state (every relay is in the “down” state, i.e., Sut=T and u=ρ0). The output value in the ρ2 point is called wρ2 (inversion point). Then, u is brought back to ρ0. The branch Fρ2 is drawn by taking the output value wρ1,ρ2 for any value u=ρ1 as is shown in Figure 8 (left). The branch ends for u=ρ0, when Sut=T. The Everett function is defined in [8] as.

which is half of the output variation along the FORC branch starting in ρ2. The Everett function and the Preisach function are related by the following integral:

Figure 8.

Left – First order reversal curves (solid line) and major loop (dashed line). Right – Triangle Tρ1ρ2.


where the integration domain Tρ1ρ2 is the triangle highlighted in Figure 8 (right). It should be noted that the integral on the triangle Tρ1ρ2 is a function of its upper-left corner and that Eρ1ρ2=0 if ρ1=ρ2 (Tρ1ρ2 degenerates in a single point). The introduction of the Everett integral simplifies the computation of the first integral on the right hand side of (10). First of all, we subdivide Su+t into k trapezoids Qkt so that Su+=kQk. Moreover, each trapezoid can be represented as the set difference of two triangles Tmk1Mk and TmkMk:


Taking into account that the integral on a generic triangle Tρ1ρ2 can be expressed by (15), the Preisach integral (10) becomes


where Eρ0ρ0=Tpρ. Form the previous procedure we obtain the following results. First, the output can be computed by using a simple linear combination of the E values in the memory state points represented by the vertices on Su+. Second, the value of the Preisach density p is not required, since the identification of E in the domain T is sufficient to apply the model. The identification of E is simple and it is defined by a repeatable and reliable procedure based on experimental data (FORC branch).


4. Some applications of the hysteresis modeling in electrical engineering

In this section we present some applications of hysteresis modeling in electrical engineering. The first two are motivated by the hysteresis observed when modeling the BH curve of a ferromagnetic material. The third one has to do with the hysteresis present in the State of Charge (SoC)-Open Circuit Voltage (OCV) mapping of Li-ion batteries.

4.1 Hysteresis power losses computation

For the performance analysis of electrical machines, an important factor to be considered is the power losses that occur in the ferromagnetic materials that make up the core of the machine. These losses, usually known as iron losses, can generally be divided into eddy current, hysteresis and excess losses. Eddy current losses are resistive losses due to the currents induced in the magnetic material by the time varying magnetic induction; its magnitude strongly depend on the size of the continuous conductive regions. That is why the core of the machine is usually assembled from thin steel sheets, insulated from each other by a non-conductive coating on the surface.

The hysteresis and excess losses are essentially related to the microscopic properties of ferromagnetic materials. These consist of small magnetic domains, which tend to align along the external (time – varying) applied magnetic field. As a consequence of the molecular friction produced in this continuous movement, heat losses, generally referred to as magnetic hysteresis losses, are produced. Given a material point x, the density of hysteresis losses in the time interval t1t2 can be computed from the integral:


Nevertheless, these losses are especially difficult to compute, since for ferromagnetic materials, the magnetic induction in each point depends on the intensity of the present magnetic field to which it is exposed, and also on previous exposures to magnetic field intensity of each volume element. This causes differences in the magnetization curve under increasing and decreasing fields; therefore, hysteresis loops arise as the one illustrated in Figure 9 (left).

Figure 9.

Toroidal laminated media (left) and its meridian section (center); hysteresis loop measured and approached (right).

In the literature there are several simplified analytical expressions that are used to approximate the different components of the losses, such as those proposed by Bertotti [14], which are among the most widely used. However, the assumptions under which it is possible to apply these formulas are not met in most practical situations. In this context, numerical simulation is a viable option to overcome these limitations.

As an example, we consider an application consisting of the computation of hysteresis and eddy current losses in a laminated medium with toroidal geometry, as sketched in Figure 9. It consists of N circular sheets of rectangular section of thickness d surrounded by a coil. We denote by R1 and R2 the internal and external radius of the core, respectively. We will assume that the coil is infinitely thin so that it can be modeled as a surface current density (A/m).

As shown in [24], in this situation, it is possible to reduce the computational domain to the meridian section of one single sheet for which it is necessary to derive appropriate boundary conditions. More precisely, the axisymmetric transient eddy current problem to be solved reads:

FindHrztandBrztsuch that


where Ω=R1R2×0d, Γ its boundary, and σrzt,hrzt,ξrz and B0rz are given functions. We have employed the usual notation: H is the magnetic field, B the magnetic induction and σ the electrical conductivity. Moreover, it can be shown that (see [24] for details):


where It is the current intensity flowing through the coil at time t and ne denotes the number of winding turns. In order to model the B-H relation, the classical Preisach hysteresis operator was employed. The Everett function was approximated from experimental data. Figure 9 (right) shows a comparison between the B-H hysteresis curves measured experimentally and those computed from the Everett function thus assessing the validity of the approximation.

The energy losses can be evaluated by computing the magnetic field and magnetic induction solutions to (19)(22) with appropriate numerical techniques, such as Finite-Difference Time-Domain (FDTD) [33, 34] or Finite-Element method (FEM) (see, for instance [35, 36, 37, 38, 39, 40, 41, 42] where scalar and vector hysteresis models are considered, respectively). In particular, we have applied the latter and a fixed-point iteration (see [24], Section 3) and also [43, 44] for different fixed point iteratios applied to hysteresis models).

In order to compute the losses, we consider the energy balance in the axisymmetric setting, which reads [24]:


In this expression, the terms LE, LH and LT represent the eddy current, the hysteresis and the total losses, respectively. From the numerical point of view, these losses have been computed (in J/m3) as


The measurements reported in Table 1 were performed on an Epstein frame considering a material sheet of width 30 mm and thickness 0.5 mm. The sheet is subjected to sinusoidal flux excitation with peak induction levels Bm equal to 0.5, 0.9 and 1.4 T and frequencies f equal to 25 and 150 Hz. For each of these peak levels and frequencies, the physical measurements were the total electromagnetic losses per cycle and unit volume and the magnetic field on the boundary of the sheet. The following geometrical data were considered to simulate the experimental setting with our axisymmetric model: R1=100 m, R2=100.3 m, d=0.0005 m, σ=4064777 (Ohm/m)−1. Table 1 summarizes the results obtained for different frequencies and magnetic induction peak levels. Let us recall that equality LE+LH=LT holds at the continuous level and thus, the difference among the losses LEh+LHh and LTh (columns 6 and 7 in Table 1) is due to numerical approximation of the axisymmetric model.

f(Hz)Bm(T)LEhLHhLEh+LHhLThTotal (exp.)aError 1(%)bError 2(%)c

Table 1.

Total losses. Numerical versus experimental results (J/m3).

Total losses computed experimentally.

Relative error between experimental losses and LEh+LHh.

Relative error between experimental losses and LTh.

This table shows that the computed losses are close to the losses measured in the experiments. The largest differences are obtained for Bm=0.5 T and are probably due to the fact that our Everett function is unable to generate accurate B-H cycles for “small” values of B (see Figure 9, right), making the numerical approximation of the hysteresis losses more inaccurate. These discrepancies can be explained by the congruency properties of the classical Preisach model which states that B-H cycles between two fixed extreme magnetic fields H are independent of the induction level B. To avoid these discrepancies, a Moving Preisach hysteresis model can be applied. In this model the congruency property of minor loops is relaxed and it is expected to reproduce lower symmetric hysteresis loops (see [12]).

4.2 Hysteresis modeling in magnetic inspection particle processes

The technique known as Magnetic Particle Inspection (usually MPI for its acronym in English) is a non-destructive method for detecting flaws located on or near the surface of ferromagnetic parts. It exploits the fact that when a ferromagnetic sample is exposed to the influence of a magnetic field, the induced magnetic flux density accumulates inside the material. Then, if there is a crack, the magnetic field will be distorted, causing local magnetic leakage around the defect. Therefore, if fluorescent magnetic particles are sprayed on the magnetized sheet, they will concentrate on the cracks and produce deposits that are easy to identify under ultraviolet light. The purpose of MPI is to find faults in any direction. Because the breakings are easier to detect when orientated perpendicular to the specimen’s magnetic field, it is common to apply the magnetization in two orthogonal orientations, such as circular and longitudinal orientation (see Figure 10, left). Most times, after inspection, the workpieces must be demagnetized since residual magnetism could interfere with later processing.

Figure 10.

Field lines in circular (left) and longitudinal (center) magnetization. Computational domain used in tests (right).

In [25], the authors introduced two models for the numerical simulation of the entire magnetization and demagnetization procedures entailed in an MPI test. In particular, the remanent flux density in specimens with cylindrical symmetry is computed. The method is based on scalar models that include magnetic hysteresis and is carried out in three steps: long-term magnetization, circular magnetization and final demagnetization. To achieve circular magnetization (and demagnetization), the workpiece is clamped between two electrical contacts between which an electric current is circulated (see Figure 10). For modeling purposes, all fields are assumed to be θ-independent. It is also assumed that the current density has the form JxtJzrztez, that is, that it traverses the piece along its axial direction. Consequently, the magnetic field stands HxtHθrzteθ. As proven in [25], the problem can be spatially reduced to the meridional section Ω of the workpiece, and reads

where Γ1, Γ2 denote the boundaries clamped to the electrical contacts, L is the workpiece length in the axial direction, Sz any cross-section transversal to this direction and RSz its radius. As in previous sections, FS is the classical Preisach hysteresis operator and ξ the initial magnetization state.

For the longitudinal magnetization and demagnetization processes, the piece is placed inside a conducting coil carrying an alternating current in the azimuthal direction (see Figure 10). To avoid the use of a vector hysteresis law, an infinitely thin conducting surface ΓS of radius RS carrying a surface current density JSxtJSrteθ is employed to approximate the coil. Let Ωc be the part to be inspected, Rc representing its radius; moreover, let us denote Ω0 the air surrounding it which is assumed to be artificially bounded in the r-direction, R being its outer radius. Then problem reduces to (see [25])


where Aθ denotes the azimuthal of component magnetic vector potential. In (34)r=RS represents the jump across the surface r=RS. Notice that the source is the surface current density JS which is equal to the jump of H×n at r=RS.

It is worth noting that when solving (32)(35), the effective computation of the inverse of the hysteresis operator is avoided by using an iterative algorithm of a fixed point type, based on the properties of the Yosida regularization for maximal monotone operators (see [45]). Moreover, the fields H and B have only one non-null component in the eθ direction in the circular case, and in the ez direction in the longitudinal one. This allows using a scalar constitutive law for the nonlinear material, in both the circular and the longitudinal formulations.

For the sake of brevity, only the numerical results corresponding to the circular case are included. The ferromagnetic piece is characterized by the initial magnetization curve and the major hysteresis loop, as shown in Figure 11, which was obtained by using an artificial weight function p. The magnetization stage is performed by using a harmonic surface current density on the cylindrical surface Γs. When the source is turned off, the workpiece retains a residual field or remanence. In order to remove this remanent field, a surface current density on Γs is considered with amplitude equal to the one used for magnetizing the piece and continuously reversing while gradually decreasing to zero (see Figure 11). The numerical results give important information about the source needed to successfully demagnetize the piece. In particular, the efficiency of the demagnetization process strongly depends on the source and the number of cycles. In particular, the greater the number of source cycles, the better the demagnetization.

Figure 11.

Left – Major hysteresis loop and (blue) and initial magnetization curve (red). Center – Demagnetizing source. Right – Circular demagnetization. Remanent flux density at z=0.2 m vs. radius. Number of cycles comparison.

4.3 Hysteresis model applied to battery modeling

It is well-known that hysteresis has a significant impact on the ability to monitor battery performance since it affects the voltage during charging and discharging cycles for different battery chemistries. Therefore, it must be considered for diagnosis algorithms that use the so-called open circuit voltage (OCV) as a measure of the state of charge (SoC) of the battery. SoC indicates the residual charge of the battery with respect to its maximum nominal and it is usually expressed as a percentage of the battery capacity. Since the SoC provides a measure of the available energy of the battery, as well as of its instantaneous power capacity, an accurate estimate is important to monitor battery performance. Unfortunately, it cannot be directly measured but must be estimated from other battery quantities (see, for instance [46, 47, 48, 49, 50], to mention some recent works). However, the non-linear electrochemical reactions involved in the battery system, the effects of temperature, aging and, specially, hysteresis effects make this task especially difficult.

For a battery cell, hysteresis means that the battery cell reaches different OCV values at the same SoC value and temperature between charge and discharge, depending on the previous charge–discharge history. In particular, from experiments it is shown that, after discharging, the OCV always relaxes below the OCV after charging for the same SoC (see Figure 12, left). As a consequence, the SoC-OCV relation is a bundle of curves enclosed in a major loop. The loop consists of two SoC-OCV curves obtained by fully charging and discharging the battery, first to minimum voltage and then to maximum voltage while recording the battery voltage and accumulated ampere hours, for a complete battery cycle. Minor loops lying inside the major SoC-OCV loop can be obtained by conducting similar experiments but with partial cycles (see Figure 12, right). A simple approach is to compute the average of the major loop and to relate the battery’s OCV with its SoC through this average single-valued curve (see Figure 12, left), but usually this approximation suffers from significant errors in the real SoC estimation.

Figure 12.

Measurements along the major loop and the approximation given by the Preisach model (left). Experimental inner loops and approximation of the Preisach model for different SoC values (right).

Based on the work [51], in this section we consider a Preisach-type hysteresis model to approximate the SoC-OCV hysteresis loop (see [52, 53, 54] for different approaches to handle hysteresis in batteries). The battery dynamics is modeled by using an equivalent circuit model (ECM) as extensively used in the literature [55]. In particular, we consider a general j-th order circuit (see Figure 13 (left) for the case J = 2). It consists of the so-called equivalent series resistance and a number of J resistor-capacitor RjCj elements in series. The latter allow us the modeling of the diffusion voltages and are a good approximation of the so-called Warburg impedance. The notations employed in the sequel are the usual ones: Q is the maximum charge (or capacity) of the battery, it is the current intensity, ηt the Coulombic efficiency, Vt the voltage, zt the SoC, Rj the jth polarization resistance, Cj the jth polarization capacitance, and vRjt is the voltage between the ends of the the j-th resistor Rj. Thus, the problem to solve reads:

Figure 13.

Sketch of the 2RC model and parameters values.

Given functionsOCV̂,itandηt, constant parametersQ,R0,RjandCj, and initial conditionsz0andvRj,0,j=1,,J, find functionszt,VtandvRjt, satisfying,


where functionOCV̂gives the OCV from the state of chargez.

From a macroscopic standpoint, the SoC-OCV hysteresis relation can be considered rate-independent; this implies that the OCV depends on the SoC history and not on the velocity (battery current rate) of SoC. Our aim is now to apply a hysteresis model to a lithium-iron-phosphate (LFP) cell. It is known that the classical Preisach operator cannot be directly applied to model the SoC-OCV hysteresis loop. Thus, the hysteretic behavior is modeled with a modified Preisach operator in which the SoC is assumed to be the independent variable. In particular we consider


where Sz+ is defined in (9) and OCVmin is the minimum OCV value. As usual, the computation of the hysteretic term in (41) is done by means of the Everett function instead of the Preisach hysteresis function p. This function is identified from experimental data considering only the charging SoC-OCV curves of the battery, which is less time consuming than other approaches proposed in the literature. The SoC-OCV curves computed with the Preisach model (41) are compared with experimental data. The approximations are reported in Figure 12, where the experimental data and values simulated with the Preisach model are plotted. From this figure, it can be seen that the model approximates the measurements of both the major loop and the inner loops.

Next, we consider a realistic vehicle driving profile to quantify the accuracy of the ECM with hysteresis. We have solved (36)(39) and compared the computed voltage (40) with the experimental values. The values of OCV are obtained with the average open circuit voltage curve (see Figure 12, left) and with the Preisach hysteresis model (cf. (41)). Figure 14 shows the results obtained when a second-order circuit model (2RC) is considered. More precisely, this figure shows, on the left, the measured voltage and on the right the error between the experimental voltage and the voltage computed with the average curve (in red) and the Preisach model (in blue), respectively. The result has been excellent when the Preisach model is considered, having a relative error of 0.4%.

Figure 14.

Percentage error in the voltage approximation (left). Measured and computed voltage with ECM and Preisach model (right).


5. Conclusions

In this chapter we summarize some applications of the hysteresis Preisach model in electrical engineering. First, we have introduced the classical or rate-independent Preisach hysteresis model and the rate-dependent extension of this model, the so-called dynamic Preisach model. Moreover, the identification of this operator based on the Everett function is explained. Finally, we provide some examples of application of the Preisach model in electrical engineering which highlight the importance of having a good hysteresis model in those processes that are characterized by this complex behavior.



Work partially founded by FONDECYT-Chile project 1211030, by Xunta de Galicia project GI-1563 ED431C 2021/15 and by FEDER, Ministerio de Economía, Industria y Competitividad-AEI under grant project MTM2017-86459-R.


  1. 1. J. A. Ewing. X. experimental researches in magnetism. Philos. Trans. R. Soc. Lond., 176:523–640, 1885
  2. 2. M.A. Krasnosel’skii and A.V. Pokrovskii. Systems with Hysteresis. Springer, Berlin, 1989 (Russian edition, Nauka, Moscow (1983))
  3. 3. J. W. Macki, P. Nistri, and P. Zecca. Mathematical Models for Hysteresis. SIAM Review, 35(1):94–123, 1993
  4. 4. V. Hassani, T. Tjahjowidodo, and N. D. Thanh. A survey on hysteresis modeling, identification and control. Mech. Syst. Signal Proc., 49(1-2, SI):209–233, 2014
  5. 5. A. Visintin. Chapter 1 – Mathematical Models of Hysteresis. In Giorgio Bertotti and I. D. Mayergoyz, editors, The Science of Hysteresis, pages 1–123. Academic Press, Oxford, 2006
  6. 6. A. Visintin. Differential models of hysteresis, volume 111 of Applied Mathematical Sciences. Springer-Verlag, Berlin, 1994
  7. 7. M. Brokate and J. Sprekels. Hysteresis and Phase Transitions. Springer, Berlin, 1996
  8. 8. I. D. Mayergoyz. Mathematical models of hysteresis. Springer-Verlag, New York, 1991
  9. 9. G. Bertotti. Hysteresis in Magnetism. Academic Press, 1998
  10. 10. F. Preisach. Über die magnetische nachwirkung. Zeitschrift für Physik, 94:277–302, 1935
  11. 11. D. A. Philips, L. Dupré, J. Cnops, and J. A. A. Melkebeek. The application of the Preisach model in magnetodynamics: theoretical and practical aspects. J. Magn. Magn. Mater., 133(1):540–543, 1994
  12. 12. E. Della Torre. Magnetic Hysteresis. IEEE Press, New York, 1999
  13. 13. I. D. Mayergoyz. Dynamic Preisach Model of Hysteresis. IEEE Trans. Magn., 24(6):2925–2927, 1988
  14. 14. G. Bertotti. Dynamic generalization of the scalar Preisach model of hysteresis. IEEE Trans. Magn., 28(5):2599–2601, 1992
  15. 15. Y. Yu, Z. Xiao, N. G. Naganathan, and R. V. Dukkipati. Dynamic Preisach modelling of hysteresis for the piezoceramic actuator system. Mech. Mach. Theory, 37(1):75–89, 2002
  16. 16. L. Dupré, G. Bertotti, V. Basso, F. Fiorillo, and J.A.A Melkebeek. Generalisation of the dynamic Preisach model toward grain oriented FeSi alloys. Physica B, 275(1):202–206, 2000
  17. 17. A. Bermúdez, D. Gómez, and P. Venegas. Mathematical analysis and numerical solution of models with dynamic Preisach hysteresis. J. Comput. Appl. Math., 367:112452, 2020
  18. 18. O. Bottauscio, M. Chiampi, D. Chiarabaglio, and M. Repetto. Preisach-type hysteresis models in magnetic field computation. Physica B, 275(1-3):34–39, 2000
  19. 19. B. Schweizer. Hysteresis in porous media: Modelling and analysis. Interface Free Bound., 19(3):417–447, 2017
  20. 20. V. A. Lubarda, D. Šumarac, and D. Krajcinovic. Preisach model and hysteretic behaviour of ductile materials. European J. Mech. A Solids, 12(4):445–470, 1993
  21. 21. A. Rao and A.R. Srinivasa. A two species thermodynamic Preisach model for the torsional response of shape memory alloy wires and sprin under superelastic conditions. Int. J. Solids Struct., 50(6):887–898, 2013
  22. 22. P. Krejčí, J. P. O’Kane, A. Pokrovskii, and D. Rachinskii. Properties of solutions to a class of differential models incorporating Preisach hysteresis operator. Phys. D, 241(22):2010–2028, 2012
  23. 23. A. Bermúdez, D. Gómez, R. Rodríguez, and P. Venegas. Mathematical analysis and numerical solution of axisymmetric eddy-current problems with Preisach hysteresis model. Rend. Sem. Mat. Universitá e Politecnico Torino, special volume, 72(1-2):73–117, 2014
  24. 24. A. Bermúdez, L. Dupré, D. Gómez, and P. Venegas. Electromagnetic computations with Preisach hysteresis model. Finite Elem. Anal. Des., 126:65–74, 2017
  25. 25. A. Bermúdez, D. Gómez, M. Piñeiro, P. Salgado, and P. Venegas. Numerical simulation of magnetization and demagnetization processes. IEEE Trans. Magn., 53(12):1–6, 2017
  26. 26. M. Kuczmann. Vector Preisach hysteresis modeling: Measurement, identification and application. Physica B, 406(8):1403–1409, 2011
  27. 27. M. Kuczmann. Dynamic extension of vector Preisach model. Physica B, 549:47–52, 2018. 11th International Symposium on Hysteresis Modeling and Micromagnetics (HMM 2017)
  28. 28. I. D. Mayergoyz. Mathematical Models of Hysteresis and Their Applications. Elsevier, 2003
  29. 29. E. Cardelli and A. Faba. Modelling of vector hysteresis at macromagnetic scale: Open questions and challenges. Physica B, 486:130–137, 2016
  30. 30. E. Cardelli, E. Della Torre, and A. Faba. A General Vector Hysteresis Operator: Extension to the 3-D Case. IEEE Trans. Magn., 46(12):3990–4000, 2010
  31. 31. K. Preis. A contribution to eddy current calculations in plane and axisymmetric multiconductor systems. IEEE Trans. Magn., 19:2397–2400, 1983
  32. 32. G. Bertotti and M. Pasquale. Physical interpretation of induction and frequency dependence of power losses in soft magnetic materials. IEEE Trans. Magn., 28(5):2787–2789, 1992
  33. 33. F. Bertoncini, F. Beux, E. Cardelli, S. Di Fraia, and B. Tellini. Stable FDITD formulation for electromagnetic field diffusion in soft magnetic materials. IEEE Trans. Magn., 39(3):1681–1684, 2003
  34. 34. S. Quondam Antonio, A. Faba, H. P. Rimal, and E. Cardelli. On the analysis of the dynamic energy losses in NGO electrical steels under non-sinusoidal polarization waveforms. IEEE Trans. Magn., 56(4):1–15, 2020
  35. 35. R. Van Keer, L.R. Dupré, and J.A.A. Melkebeek. On a numerical method for 2D magnetic field computations in a lamination with enforced total flux. J Comput. Appl. Math., 72(1):179–191, 1996
  36. 36. N. Takahashi, S.-H. Miyabara, and K. Fujiwara. Problems in practical finite element analysis using Preisach hysteresis model. IEEE Trans. Magn., 35(3):1243–1246, 1999
  37. 37. E. Fallah and V. Badeli. A new approach for modeling of hysteresis in 2-D time-transient analysis of eddy current using FEM. IEEE Trans. Magn., 53(7):1–14, 2017
  38. 38. Y. Li, L. Zhu, and J. Zhu. Core loss calculation based on finite-element method with Jiles–Atherton dynamic hysteresis model. IEEE Trans. Magn., 54(3):1–5, 2018
  39. 39. L. R. Dupré, R. Van Keer, and J. A. A. Melkebeek. Complementary 2-D finite element procedures for the magnetic field analysis using a vector hysteresis model. Int. J. Numer. Meth. Eng., 42(6):1005–1023, 1998
  40. 40. J. B. Padilha, P. Kuo-Peng, N. Sadowski, and N. J. Batistela. Vector Hysteresis Model Associated to FEM in a Hysteresis Motor Modeling. IEEE Trans. Magn., 53(6), 2017
  41. 41. M. Tousignant, F. Sirois, G. Meunier, and C. Guérin. Incorporation of a vector Preisach–Mayergoyz hysteresis model in 3-D finite element analysis. IEEE Trans. Magn., 55(6):1–4, 2019
  42. 42. C. Guérin, K. Jacques, R. V. Sabariego, P. Dular, C. Geuzaine, and J. Gyselinck. Using a Jiles-Atherton vector hysteresis model for isotropic magnetic materials with the finite element method, Newton-Raphson method, and relaxation procedure. Int. J. Numer. Model El., 30(5):e2189, 2017
  43. 43. E. Dlala and A. Arkkio. Analysis of the convergence of the fixed-point method used for solving nonlinear rotational magnetic field problems. IEEE Trans. Magn., 44(4):473–478, 2008
  44. 44. P. Zhou, D. Lin, C. Lu, M. Rosu, and D. M. Ionel. An adaptive fixed-point iteration algorithm for finite-element analysis with magnetic hysteresis materials. IEEE Trans. Magn., 53(10):1–5, 2017
  45. 45. A. Bermúdez and C. Moreno. Duality methods for solving variational inequalities. Comput. Math. Appl., 7:43–58, 1981
  46. 46. M. A. Hannan, M. S.H. Lipu, A. Hussain, and A. Mohamed. A review of lithium-ion battery state of charge estimation and management system in electric vehicle applications: Challenges and recommendations. Renew. Sust. Energ. Rev., 78:834–854, 2017
  47. 47. J. P. Rivera-Barrera, N. Muñoz Galeano, and H. O. Sarmiento-Maldonado. SOC estimation for lithium-ion batteries: Review and future challenges. Electronics, 6(4), 2017
  48. 48. R. Zhang, B. Xia, B. Li, L. Cao, Y. Lai, W. Zheng, H. Wang, and W. Wang. State of the art of lithium-ion battery SOC estimation for electrical vehicles. Energies, 11(7), 2018
  49. 49. M. U. Ali, A. Zafar, S. H. Nengroo, S. Hussain, M. Junaid Alvi, and H.J. Kim. Towards a smarter battery management system for electric vehicle applications: A critical review of lithium-ion battery state of charge estimation. Energies, 12(3), 2019
  50. 50. Y. Wang, J. Tian, Z. Sun, L. Wang, R. Xu, M. Li, and Z. Chen. A comprehensive review of battery modeling and state estimation approaches for advanced battery management systems. Renew. Sust. Energ. Rev., 131:110015, 2020
  51. 51. P. Venegas, D. Gómez, M. Arrinda, M. Oyarbide, H. Macicior, and A. Bermúdez. Kalman filter and classical Preisach hysteresis model applied to the state of charge battery estimation. Submitted
  52. 52. V. Franzitta, A. Viola, and M. Trapanese. Description of Hysteresis in Lithium Battery by Classical Preisach Model. Adv. Mat. Res., 622-623:1099–1103, 2013
  53. 53. F. Baronti, N. Femia, R. Saletti, C. Visone, and W. Zamboni. Hysteresis modeling in Li-ion batteries. IEEE Trans. Magn., 50(11):1–4, 2014
  54. 54. L. Zhu, Z. Sun, H. Dai, and X. Wei. A novel modeling methodology of open circuit voltage hysteresis for LiFePO4 batteries based on an adaptive discrete Preisach model. Appl. Energ., 155:91–109, 2015
  55. 55. G. L. Plett. Battery Management Systems: Volume 1, Battery Modeling. Artech House, Boston, 2015

Written By

Alfredo Bermúdez, Dolores Gómez and Pablo Venegas

Reviewed: 21 July 2021 Published: 03 November 2021