Open access peer-reviewed chapter - ONLINE FIRST

Preisach Hysteresis Model. Some Applications in Electrical Engineering

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

Submitted: February 10th 2021Reviewed: July 21st 2021Published: November 3rd 2021

DOI: 10.5772/intechopen.99590

Downloaded: 28


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 staticor rate-independentmodels, 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 dynamicor rate-dependentmodels, 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, uand w. Let us also suppose that uis the independent variable, and thus the evolution of wdepends on u. Such diagrams implicitly show the memory effectinherent to the hysteresis: at any time t>0, wtdepends on the evolution of uin 0t, rather than only on ut. In many cases, not all the information in 0tis used to compute the value of wbut only the one which has not been wiped-outfrom the memory (see Section 3.1). In other words, the state of the system at present time depends on the previous historyof its state. Consequently, the same instantaneous values of the input can give different outputs and when describing the relation between wand uin the uwplane 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: causalityand rate independence. To clarify these concepts, we consider the uwrelation introduced above. Causalitymeans that, at any time t, the value of wtonly depends on the previous evolution of u, namely, wtdepends on u0t. On the other hand, rate independenceis translated by the condition that, at any instant t, wtdepends just on the range of function u:0tRand on the order in which the values of uhave been attained. In other words, wis independent of the velocityof u. An example is given in Figure 1 (right) where it is shown that three different functions u:0tRlead to the same wucurve.

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, wtwill 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 ubefore t=0. Hence, not only the standard initial value u0w0must 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 Facting 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:Ω̂×0TRand any ξ:Ω̂Y, with Ya suitable space, the corresponding space dependent operator Fas follows (see [23] for a rigorous mathematical setting)


We notice that operator F˜is local in space, i.e., the output Fuξxtdepends on ux0t, but not on uy0tfor 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 utand 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 utand the output wtis 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ρ2is the relay. The function μis a non-negative function with compact support in the Preisach half-plane with ρ1ρ2that 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 modelbased on the rate-independent relay, and the dynamic Preisach modelintroduced 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: ρ2is the switch-up threshold, and ρ1is the switch-down one.

Figure 2.

Scalar relay.

Formally, for any continuous function uC0Tand initial condition ξ11, hρuξis a real function defined in the time interval 0Tsuch 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 upresents the thresholds ρ1or ρ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 ±1and it changes instantaneously from its last value depending on the previous evolution of the system; more precisely, when ureaches the threshold ρ2from below, it “switches-up” to value 1, and when it attains ρ1from above, it “switches-down” to 1. Therefore, hρis not a local in time mapping: hρuξρtnot only depends on utbut on its past history. The classical Preisach operatorFSis then defined as


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

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

Figure 3.

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


being ρ0ρ0the 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 uC0Tand any time t, the triangle is split into two sets (one of them possibly empty): Su+tand Sutcontaining the points ρ1ρ2for which the associated relays hρuare positive or negative, respectively, i.e.,


(see Figure 3). The interface Lutbetween the two sets is a staircase line with vertices having coordinates ρ1ρ2respectively coinciding with the local minimum and maximum values of uat previous time instants. At time t, Lutintersects the line ρ1=ρ2at utut. Lutmoves up as utincreases and from right to left as utdecreases (see Figure 4).

Figure 4.

Staircase lineLutmoving 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+tis obtained. Since the second integral in (10) is constant in time, the computation of the output wtdepends on the effective computation of the first integral.

From (10), is clear that wtdepends on the shape of the interface Lutand 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, Luterases (“wipes out”) the previous vertices of the staircase having a ρ2value lower than the value ut. Similarly, when an input reaches a local minimum ut, the memory curve wipes out all previous vertices with a ρ1greater 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 por a suitable auxiliary function.

  • At each time step the set Su+thas 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 wtis obtained by computing the integral of pin 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 dynamicPreisach 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 ±1are 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 utand the threshold values ρ1and ρ2(see Figure 5, right panel). The proportionality factor is a parameter kdepending on the material.

Figure 5.

Input functionut=150sin2πft(left) for a fixed frequency f = 20 Hz. A dynamic relay and its slope are presented for switching valuesρ1ρ2=50,100and initial stateξ=1. The center panel shows the relay with slopek=50whereas the right panel shows the corresponding diagramuηρ/t. Black arrows represent increasing values ofuwhile 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 uand ξ11, ηρuξ:0T11is the unique function ysuch that yt11be 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πftis considered for frequency f=20Hz and initial condition ξ=1. This figure also shows the curve utytparameterized with respect to the time variable when the relay ηρis characterized by switching values ρ1ρ2=50,100and 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 kare 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 utytfor 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 functionut=150sin2πft+75. The center graph shows theuycurve corresponding to the dynamic relays withk=50and initial stateξ=1for frequencies f=50, 500, 5000 Hz (dashed, dash-dotted and solid line, respectively). On the right panel, theuycurve is depicted for f= 5000 Hz andk=102,104,108(dash-dotted, dashed and solid line), respectively.

From previous considerations, the dynamic Preisach operator FDcan 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=50and the demagnetized state as initial condition.

Figure 7.

Input functionut(left) defined in [0, 0.0045]. The isolines represent the corresponding dynamic relay valuesηρwithk=50and the classical relayhρ, att=0.003(center) andt=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 pbut the Everett function, which describes the effect of pon 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ρ2is associated to the threshold ρ2. The input urises up from the “reset” state (every relay is in the “down” state, i.e., Sut=Tand u=ρ0). The output value in the ρ2point is called wρ2(inversion point). Then, uis brought back to ρ0. The branch Fρ2is drawn by taking the output value wρ1,ρ2for any value u=ρ1as 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 – TriangleTρ1ρ2.


where the integration domain Tρ1ρ2is the triangle highlighted in Figure 8 (right). It should be noted that the integral on the triangle Tρ1ρ2is a function of its upper-left corner and that Eρ1ρ2=0if ρ1=ρ2(Tρ1ρ2degenerates 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+tinto ktrapezoids Qktso that Su+=kQk. Moreover, each trapezoid can be represented as the set difference of two triangles Tmk1Mkand TmkMk:


Taking into account that the integral on a generic triangle Tρ1ρ2can 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 Evalues in the memory state points represented by the vertices on Su+. Second, the value of the Preisach density pis not required, since the identification of Ein the domain Tis sufficient to apply the model. The identification of Eis 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 t1t2can 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 Ncircular sheets of rectangular section of thickness dsurrounded by a coil. We denote by R1and R2the 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,ξrzand B0rzare given functions. We have employed the usual notation: His the magnetic field, Bthe magnetic induction and σthe electrical conductivity. Moreover, it can be shown that (see [24] for details):


where Itis the current intensity flowing through the coil at time tand nedenotes 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, LHand LTrepresent 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 Bmequal to 0.5, 0.9and 1.4T and frequencies fequal to 25and 150Hz. 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=100m, R2=100.3m, d=0.0005m, σ=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=LTholds at the continuous level and thus, the difference among the losses LEh+LHhand 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.5T 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 propertiesof 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, Γ2denote the boundaries clamped to the electrical contacts, Lis the workpiece length in the axial direction, Szany cross-section transversal to this direction and RSzits radius. As in previous sections, FSis 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 ΓSof radius RScarrying a surface current density JSxtJSrteθis employed to approximate the coil. Let Ωcbe the part to be inspected, Rcrepresenting its radius; moreover, let us denote Ω0the air surrounding it which is assumed to be artificially bounded in the r-direction, Rbeing its outer radius. Then problem reduces to (see [25])


where Aθdenotes the azimuthal of component magnetic vector potential. In (34)r=RSrepresents the jump across the surface r=RS. Notice that the source is the surface current density JSwhich is equal to the jump of H×nat 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 Hand Bhave only one non-null component in the eθdirection in the circular case, and in the ezdirection 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 Γsis 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 atz=0.2m 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 Jresistor-capacitor RjCjelements 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: Qis the maximum charge (or capacity) of the battery, itis the current intensity, ηtthe Coulombic efficiency, Vtthe voltage, ztthe SoC, Rjthe jth polarization resistance, Cjthe jth polarization capacitance, and vRjtis 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 OCVminis 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 2017/60 and by FEDER, Ministerio de Economía, Industria y Competitividad-AEI under grant project MTM2017-86459-R.


chapter PDF

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Alfredo Bermúdez, Dolores Gómez and Pablo Venegas (November 3rd 2021). Preisach Hysteresis Model. Some Applications in Electrical Engineering [Online First], IntechOpen, DOI: 10.5772/intechopen.99590. Available from:

chapter statistics

28total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us