Preisach Hysteresis Model – Some Applications in Electrical Engineering

Abstract

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.

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’ski$\breve{\mathrm{i}}$ and Pokrovski$\breve{\mathrm{i}}$ [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$, $w\left(t\right)$ depends on the evolution of $u$ in $\left[0,t\right]$, rather than only on $u\left(t\right)$. In many cases, not all the information in $\left[0,t\right]$ 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 $\left(u,w\right)$ 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 $\left(u,w\right)$ relation introduced above. Causality means that, at any time $t$, the value of $w\left(t\right)$ only depends on the previous evolution of $u$, namely, $w\left(t\right)$ depends on ${\left.u\right|}_{\left[0,t\right]}$. On the other hand, rate independence is translated by the condition that, at any instant $t$, $w\left(t\right)$ depends just on the range of function $u:\left[0,t\right]\to \mathbb{R}$ 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:\left[0,t\right]\to \mathbb{R}$ lead to the same $w-u$ 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$, $w\left(t\right)$ will depend not only on the previous evolution of $u$ (i.e., on ${\left.u\right|}_{\left[0,t\right]}$), 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 $\left(u\left(0\right),w\left(0\right)\right)$ must be provided. In general, a variable $\xi$ containing all the information about the “initial state” is considered. This can be expressed, for instance, as follows:

$$\left(u,\xi \right)\to w=\stackrel{\sim }{\mathcal{F}}\left(u,\xi \right).$$E1

In the case of partial differential equations, it is necessary to define an operator $\mathcal{F}$ acting between suitable function spaces involving the space variable. Let $\widehat{\mathrm{\Omega }}$ be an open subset of ${\mathbb{R}}^N\left(N\ge 1\right)$; given a hysteresis operator $\stackrel{\sim }{\mathcal{F}}$, we introduce, for any $u:\widehat{\mathrm{\Omega }}\times \left[0,T\right]\to \mathbb{R}$ and any $\xi :\widehat{\mathrm{\Omega }}\to Y$, with $Y$ a suitable space, the corresponding space dependent operator $\mathcal{F}$ as follows (see [23] for a rigorous mathematical setting)

$$\left[\mathcal{F}\left(u,\xi \right)\right]\left(x,t\right)≔\left[\stackrel{\sim }{\mathcal{F}}\left(u\left(x,\cdot \right),\xi \left(x\right)\right)\right]\left(t\right),\forall t\in \left[0,T\right],\mathrm{a}.\mathrm{e}.\text{in}\widehat{\mathrm{\Omega }}.$$E2

We notice that operator $\stackrel{\sim }{\mathcal{F}}$ is local in space, i.e., the output $\left[\mathcal{F}\left(u,\xi \right)\right]\left(x,t\right)$ depends on ${\left.u\left(x,\cdot \right)\right|}_{\left[0,t\right]}$, but not on ${\left.u\left(y,\cdot \right)\right|}_{\left[0,t\right]}$ for $y\ne x$.

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 $\mathit{u}\left(t\right)$ and vector output $\mathcal{F}\left(\mathit{u},t\right)$; 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 $\mu$, 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 $u\left(t\right)$ and the output $w\left(t\right)$ is given by means of the integral:

$$w\left(t\right)={\iint }_{{\rho }_2\ge {\rho }_1}{\gamma }_{{\rho }_1{\rho }_2}\left(u\left(t\right),\xi \right)\mu \left({\rho }_1,{\rho }_2\right)\mathrm{d}{\rho }_1\mathrm{d}{\rho }_2,$$E3

where $\mu$ is the weight function, $\xi$ contains the information about the “initial state” of every point of the domain and ${\gamma }_{{\rho }_1{\rho }_2}$ is the relay. The function $\mu$ is a non-negative function with compact support in the Preisach half-plane with ${\rho }_1\ge {\rho }_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_{\rho }$, is associated to a point $\rho ≔\left({\rho }_1,{\rho }_2\right)\in {\mathbb{R}}^2$, with ${\rho }_1<{\rho }_2$, and, for any $\rho$, it has two states: “up” ($h_{\rho }=1$) and “down” ($h_{\rho }=-1$) and two switching thresholds: ${\rho }_2$ is the switch-up threshold, and ${\rho }_1$ is the switch-down one.

Figure 2.

Scalar relay.

Formally, for any continuous function $u\in C\left(\left[0,T\right]\right)$ and initial condition $\xi \in \left\{1,-1\right\}$, $h_{\rho }\left(u,\xi \right)$ is a real function defined in the time interval $\left[0,T\right]$ such that,

$$h_{\rho }\left(u,\xi \right)\left(0\right)≔\left\{\begin{array}{ll}-1& \text{if}u\left(0\right)\le {\rho }_1,\\ \xi & \text{if}{\rho }_1<u\left(0\right)<{\rho }_2,\\ 1& \text{if}u\left(0\right)\ge {\rho }_2.\end{array}\right.$$E4

Then, for any $t\in \left(0,T\right]$, we consider the set $X_u\left(t\right)≔\{\tau \in (0,t]:u\left(\tau \right)={\rho }_1\mathit{\text{or}}{\rho }_2\}$. This set keeps account of the previous time instants in which $u$ presents the thresholds ${\rho }_1$ or ${\rho }_2$. Next, we define

$$h_{\rho }\left(u,\xi \right)\left(t\right)≔\left\{\begin{array}{llll}{h}_{\rho }\left(u,\xi \right)\left(0\right)& \text{if}{X}_u\left(t\right)=\varnothing ,& & \\ -1& \text{if}{X}_u\left(t\right)\ne \varnothing & \text{and}& u\left(max{X}_u\left(t\right)\right)={\rho }_1,\\ 1& \text{if}{X}_u\left(t\right)\ne \varnothing & \text{and}& u\left(max{X}_u\left(t\right)\right)={\rho }_2.\end{array}\right.$$E5

If $u\left(0\right)=0$, then the following initial condition can be considered

$$h_{\rho }\left(u,\xi \right)\left(0\right)≔\left\{\begin{array}{ll}-1& \text{if}{\rho }_1+{\rho }_2>0,\\ 1& \text{if}{\rho }_1+{\rho }_2<0.\end{array}\right.$$E6

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_{\rho }$ can only take values $\pm 1$ and it changes instantaneously from its last value depending on the previous evolution of the system; more precisely, when $u$ reaches the threshold ${\rho }_2$ from below, it “switches-up” to value $1$, and when it attains ${\rho }_1$ from above, it “switches-down” to $-1$. Therefore, $h_{\rho }$ is not a local in time mapping: $h_{\rho }\left(u,{\xi }_{\rho }\right)\left(t\right)$ not only depends on $u\left(t\right)$ but on its past history. The classical Preisach operator${\mathcal{F}}_S$ is then defined as

$$w\left(t\right)=\left[{\mathcal{F}}_S\left(u,\xi \right)\right]\left(t\right)={\int }_{\mathcal{T}}\left[h_{\rho }\left(u,\xi \left(\rho \right)\right)\right]\left(t\right)p\left(\rho \right)\mathit{d\rho },$$E7

where we recall that variable $\xi \left(\rho \right)$ 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 $u\left(t\right)$ (left) and its corresponding map on the Preisach triangle (right).

$$\mathcal{T}≔\left\{\rho =\left({\rho }_1,{\rho }_2\right)\in {\mathbb{R}}^2:-{\rho }_0\le {\rho }_1\le {\rho }_2\le {\rho }_0\right\},$$E8

being $\left(-{\rho }_0,{\rho }_0\right)$ the thresholds setting the minimum and maximum values of $u\left(t\right)$. 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 $u\in C\left(\left[0,T\right]\right)$ and any time $t$, the triangle is split into two sets (one of them possibly empty): $S_u^{+}\left(t\right)$ and $S_u^{-}\left(t\right)$ containing the points $\left({\rho }_1,{\rho }_2\right)$ for which the associated relays $h_{\rho }\left(u\right)$ are positive or negative, respectively, i.e.,

$$S_u^{-}\left(t\right)=\left\{\rho \in \mathcal{T}:\left[h_{\rho }\left(u\right)\right]\left(t\right)=-1\right\}\text{and}{S}_u^{+}\left(t\right)=\left\{\rho \in \mathcal{T}:\left[h_{\rho }\left(u\right)\right]\left(t\right)=1\right\}$$E9

(see Figure 3). The interface $L_u\left(t\right)$ between the two sets is a staircase line with vertices having coordinates $\left({\rho }_1,{\rho }_2\right)$ respectively coinciding with the local minimum and maximum values of $u$ at previous time instants. At time $t$, $L_u\left(t\right)$ intersects the line ${\rho }_1={\rho }_2$ at $\left(u\left(t\right),u\left(t\right)\right)$. $L_u\left(t\right)$ moves up as $u\left(t\right)$ increases and from right to left as $u\left(t\right)$ decreases (see Figure 4).

Figure 4.

Staircase line $L_u\left(t\right)$ moving right to left (left) and moving up (right).

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

$$w\left(t\right)={\int }_{S_u^{+}\left(t\right)}p\left(\rho \right)\mathit{d\rho }-{\int }_{S_u^{-}\left(t\right)}p\left(\rho \right)\mathit{d\rho }=2{\int }_{S_u^{+}\left(t\right)}p\left(\rho \right)\mathit{d\rho }-{\int }_{\mathcal{T}}p\left(\rho \right)\mathit{d\rho },$$E10

which can be readily calculated once $S_u^{+}\left(t\right)$ is obtained. Since the second integral in (10) is constant in time, the computation of the output $w\left(t\right)$ depends on the effective computation of the first integral.

From (10), is clear that $w\left(t\right)$ depends on the shape of the interface $L_u\left(t\right)$ and the latter is determined by the extremal values of $u\left(t\right)$, at previous instants of time. Notice that not all extremal input values are needed. In fact, considering the dependence on $L_u\left(t\right)$, the Preisach operator exhibits a wiping-out property. This means that every time the input reaches a local maximum $u\left(t\right)$, $L_u\left(t\right)$ erases (“wipes out”) the previous vertices of the staircase having a ${\rho }_2$ value lower than the value $u\left(t\right)$. Similarly, when an input reaches a local minimum $u\left(t\right)$, the memory curve wipes out all previous vertices with a ${\rho }_1$ greater than the value of $u\left(t\right)$. 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 $S_u^{+}\left(t\right)$ 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 $w\left(t\right)$ is obtained by computing the integral of $p$ in the domain $S_u^{+}\left(t\right)$.

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 $w\left(t\right)$, 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 $\pm 1$ are possible, the dynamic relay can attain all the intermediate values in the interval $\left[-1,1\right]$, switching at a finite velocity which is assumed to be proportional to the difference between $u\left(t\right)$ and the threshold values ${\rho }_1$ and ${\rho }_2$ (see Figure 5, right panel). The proportionality factor is a parameter $k$ depending on the material.

Figure 5.

Input function $u\left(t\right)=150sin\left(2\mathit{\pi ft}\right)$ (left) for a fixed frequency f = 20 Hz. A dynamic relay and its slope are presented for switching values $\left({\rho }_1,{\rho }_2\right)=\left(\mathrm{50,100}\right)$ and initial state $\xi =-1$. The center panel shows the relay with slope $k=50$ whereas the right panel shows the corresponding diagram $\left(u,\partial {\eta }_{\rho }/\partial t\right)$. Black arrows represent increasing values of $u$ while blue arrows represent decreasing values.

In a formal manner, and inspired by [14], for a fixed $\rho =\left({\rho }_1,{\rho }_2\right)\in {\mathbb{R}}^2$, ${\rho }_1<{\rho }_2$, we define the dynamic relay operator ${\eta }_{\rho }$ such that, for any $u$ and $\xi \in \left[-1,1\right]$, ${\eta }_{\rho }\left(u,\xi \right):\left[0,T\right]\to \left[-1,1\right]$ is the unique function $y$ such that $y\left(t\right)\in \left[-1,1\right]$ be the solution of the (nonlinear) Cauchy problem:

$$\frac{\mathit{dy}}{\mathit{dt}}\left(t\right)=F\left(t,y\left(t\right)\right)≔\left\{\begin{array}{ll}k{\left(u\left(t\right)-{\rho }_2\right)}^{+}& -k{\left(u\left(t\right)-{\rho }_1\right)}^{-}\text{if}-1<y\left(t\right)<1,\\ 0& \text{if}y\left(t\right)=-1\text{and}u\left(t\right)\le {\rho }_2,\\ k\left(u\left(t\right)-{\rho }_2\right)& \text{if}y\left(t\right)=-1\text{and}u\left(t\right)\ge {\rho }_2,\\ k\left(u\left(t\right)-{\rho }_1\right)& \text{if}y\left(t\right)=1\text{and}u\left(t\right)\le {\rho }_1,\\ 0& \text{if}y\left(t\right)=1\text{and}u\left(t\right)\ge {\rho }_1,\end{array}\right.$$E11

$$y\left(0\right)=\xi .$$E12

In the expressions above, the standard notations $x^{+}=max\left\{x,0\right\}\text{and}{x}^{-}=max\left\{-x,0\right\},$ have been used; thus $x=x^{+}-x^{-}$. We remark that in this dynamic model the initial state $\xi$ can attain all the values in the interval $\left[-1,1\right]$. 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 $\left(u\left(t\right),y\left(t\right)\right)$. For the first one, illustrated in Figure 5, a sinusoidal input $u\left(t\right)=150sin\left(2\mathit{\pi ft}\right)$ is considered for frequency $f=20$ Hz and initial condition $\xi =-1$. This figure also shows the curve $\left(u\left(t\right),y\left(t\right)\right)$ parameterized with respect to the time variable when the relay ${\eta }_{\rho }$ is characterized by switching values $\left({\rho }_1,{\rho }_2\right)=\left(\mathrm{50,100}\right)$ and slope $k=50$. The right panel shows the slope $\partial {\eta }_{\rho }/\partial t$, which represents Eq. (11). From the diagram it can be seen that when ${\eta }_{\rho }$ reaches value 1 or $-1$, the derivative is again equal to zero, as in the interval between switching values $\left({\rho }_1,{\rho }_2\right)$. 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 $\left(u\left(t\right),y\left(t\right)\right)$ for slopes $k={10}^2$ (dash-dotted line), $k={10}^4$ (dashed line) and $k={10}^8$ (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 $u\left(t\right)=150sin\left(2\mathit{\pi ft}\right)+75$. The center graph shows the $\left(u,y\right)$ curve corresponding to the dynamic relays with $k=50$ and initial state $\xi =-1$ for frequencies f=50, 500, 5000 Hz (dashed, dash-dotted and solid line, respectively). On the right panel, the $\left(u,y\right)$ curve is depicted for f= 5000 Hz and $k={10}^2,{10}^4,{10}^8$ (dash-dotted, dashed and solid line), respectively.

From previous considerations, the dynamic Preisach operator ${\mathcal{F}}_D$ can be defined as

$$w\left(t\right)=\left[{\mathcal{F}}_D\left(u,\xi \right)\right]\left(t\right)={\int }_{\mathcal{T}}\left[{\eta }_{\rho }\left(u,\xi \left(\rho \right)\right)\right]\left(t\right)p\left(\rho \right)\mathit{d\rho }.$$E13

Notice that, if ${\eta }_{\rho }$ is replaced by $h_{\rho }$ 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 ${\rho }_0=300$, a constant $k=50$ and the demagnetized state as initial condition.

Figure 7.

Input function $u\left(t\right)$ (left) defined in [0, 0.0045]. The isolines represent the corresponding dynamic relay values ${\eta }_{\rho }$ with $k=50$ and the classical relay $h_{\rho }$, 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_{{\rho }_2^{\prime }}$ is associated to the threshold ${\rho }_2^{\prime }$. The input $u$ rises up from the “reset” state (every relay is in the “down” state, i.e., $S_u^{-}\left(t\right)=\mathcal{T}$ and $u=-{\rho }_0$). The output value in the ${\rho }_2^{\prime }$ point is called $w_{{\rho }_2^{\prime }}$ (inversion point). Then, $u$ is brought back to $-{\rho }_0$. The branch $F_{{\rho }_{2^{\prime }}}$ is drawn by taking the output value $w_{{\rho }_1^{\prime },{\rho }_2^{\prime }}$ for any value $u={\rho }_{1^{\prime }}$ as is shown in Figure 8 (left). The branch ends for $u=-{\rho }_0$, when $S_u^{-}\left(t\right)=\mathcal{T}$. The Everett function is defined in [8] as.

Figure 8.

Left – First order reversal curves (solid line) and major loop (dashed line). Right – Triangle $\mathcal{T}\left({\rho }_1,{\rho }_2\right)$.

$$E\left({\rho }_1^{\prime },{\rho }_2^{\prime }\right)≔\frac{w_{{\rho }_2^{\prime }}-w_{{\rho }_2^{\prime },{\rho }_1^{\prime }}}{2},$$E14

which is half of the output variation along the FORC branch starting in ${\rho }_2^{\prime }$. The Everett function and the Preisach function are related by the following integral:

$$E\left({\rho }_1,{\rho }_2\right)={\int }_{\mathcal{T}\left({\rho }_1,{\rho }_2\right)}p\left(\rho \right)\mathit{d\rho }\forall \left({\rho }_1,{\rho }_2\right)\in \mathcal{T},$$E15

where the integration domain $\mathcal{T}\left({\rho }_1,{\rho }_2\right)$ is the triangle highlighted in Figure 8 (right). It should be noted that the integral on the triangle $\mathcal{T}\left({\rho }_1,{\rho }_2\right)$ is a function of its upper-left corner and that $E\left({\rho }_1,{\rho }_2\right)=0$ if ${\rho }_1={\rho }_2$ ($\mathcal{T}\left({\rho }_1,{\rho }_2\right)$ 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 $S_u^{+}\left(t\right)$ into $k$ trapezoids $Q_k\left(t\right)$ so that $S_u^{+}={\cup }_k{Q}_k$. Moreover, each trapezoid can be represented as the set difference of two triangles $\mathcal{T}\left(m_{k-1},M_k\right)$ and $\mathcal{T}\left(m_k,M_k\right)$:

$${\int }_{Q_k\left(t\right)}p\left(\rho \right)\mathit{d\rho }={\int }_{\mathcal{T}\left(m_{k-1},M_k\right)}p\left(\rho \right)\mathit{d\rho }-{\int }_{\mathcal{T}\left(m_k,M_k\right)}p\left(\rho \right)\mathit{d\rho }.$$E16

Taking into account that the integral on a generic triangle $\mathcal{T}\left({\rho }_1,{\rho }_2\right)$ can be expressed by (15), the Preisach integral (10) becomes

$$w\left(t\right)=2\sum _{k=1}\left(E\left(m_{k-1},M_k\right)-E(m_k,M_k)\right)-E\left(-{\rho }_0,{\rho }_0\right),$$E17

where $E\left(-{\rho }_0,{\rho }_0\right)={\int }_{\mathcal{T}}p\left(\rho \right)\mathit{d\rho }$. 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 $S_u^{+}$. Second, the value of the Preisach density $p$ is not required, since the identification of $E$ in the domain $\mathcal{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 $\left[t_1,t_2\right]$ can be computed from the integral:

$${\int }_{t_1}^{t_2}\frac{\partial \mathbf{B}}{\partial t}\left(x,t\right)\cdot \mathbf{H}\left(x,t\right)\mathit{dt}.$$E18

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 $R_1$ and $R_2$ 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:

Find$H\left(r,z,t\right)$and$B\left(r,z,t\right)$such that

$$\frac{\partial B}{\partial t}-\frac{\partial }{\partial r}\left(\frac{1}{\mathit{\sigma r}}\frac{\partial \left(\mathit{rH}\right)}{\partial r}\right)-\frac{\partial }{\partial z}\left(\frac{1}{\sigma }\frac{\partial H}{\partial z}\right)=0\text{in}\mathrm{\Omega }\times \left(0,T\right),$$E19

$$B\left(r,z,t\right)={\mu }_0\left(H\left(r,z,t\right)+\left[{\mathcal{F}}_S\right(H,\xi \left)\right](r,z,t)\right)\text{in}\mathrm{\Omega }\times \left(0,T\right),$$E20

$$H\left(r,z,t\right)=h\left(r,z,t\right)\mathrm{on}\mathrm{\Gamma }\times \left(0,T\right),$$E21

$$B\left(r,t,0\right)=B_0\left(r,z\right)\text{in}\mathrm{\Omega },$$E22

where $\mathrm{\Omega }=\left[R_1,R_2\right]\times \left[0,d\right]$, $\mathrm{\Gamma }$ its boundary, and $\sigma \left(r,z,t\right),h\left(r,z,t\right),\xi \left(r,z\right)$ and $B_0\left(r,z\right)$ are given functions. We have employed the usual notation: $H$ is the magnetic field, $B$ the magnetic induction and $\sigma$ the electrical conductivity. Moreover, it can be shown that (see [24] for details):

$$h\left(r,z,t\right)≔\frac{n_{e}I\left(t\right)}{2\mathit{\pi r}},$$E23

where $I\left(t\right)$ is the current intensity flowing through the coil at time $t$ and $n_e$ 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]:

$$\underset{{\mathcal{L}}_H}{\underbrace{{\int }_{\mathrm{\Omega }}\left({\int }_{t_1}^{t_2}\frac{\partial B}{\partial t}\mathit{Hdt}\right)r\mathit{\text{drdz}}}}+\underset{{\mathcal{L}}_E}{\underbrace{{\int }_{t_1}^{t_2}\left({\int }_{\mathrm{\Omega }}\frac{{\left|\mathbf{curl}\left(H{\mathbf{e}}_{\theta }\right)\right|}^2}{\sigma }r\mathit{\text{drdz}}\right)\mathit{dt}}}=\underset{{\mathcal{L}}_T}{\underbrace{{\int }_{t_1}^{t_2}\left({\left.\left(\mathit{rH}\right)\right|}_{\mathrm{\Gamma }}\frac{\partial }{\partial t}\left({\int }_{\mathrm{\Omega }}B\mathit{\text{drdz}}\right)\right)\mathit{dt}}}.$$E24

In this expression, the terms ${\mathcal{L}}_E$, ${\mathcal{L}}_H$ and ${\mathcal{L}}_T$ represent the eddy current, the hysteresis and the total losses, respectively. From the numerical point of view, these losses have been computed (in $\left(\mathrm{J}/{\mathrm{m}}^3\right)$) as

$${\mathcal{L}}_E^h≔\frac{1}{\mathit{\pi d}\left(R_2^2-R_1^2\right)}\left\{{\int }_0^T\left(2\pi {\int }_{\mathrm{\Omega }}\frac{{\left|curl\left(H_h{\mathbf{e}}_{\theta }\right)\right|}^2}{\sigma }\mathit{\text{rdrdz}}\right)\mathit{dt}\right\},$$E25

$${\mathcal{L}}_H^h≔\frac{1}{\mathit{\pi d}\left(R_2^2-R_1^2\right)}\left\{2\pi {\int }_{\mathrm{\Omega }}\left({\int }_0^T\frac{\partial B_h}{\partial t}{H}_h\mathit{dt}\right)\mathit{\text{rdrdz}}\right\},$$E26

$${\mathcal{L}}_T^h≔\frac{1}{\mathit{\pi d}\left(R_2^2-R_1^2\right)}\left\{{\left.2\pi {\int }_0^T{R}_2{H}_h\right|}_{\mathrm{\Gamma }}\left(\frac{d}{\mathit{dt}}{\int }_{\mathrm{\Omega }}{B}_h\mathit{\text{drdz}}\right)\mathit{dt}\right\}.$$E27

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 $B_m$ 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: $R_1=100$ m, $R_2=100.3$ m, $d=0.0005$ m, $\sigma =4064777$ (Ohm/m)⁻¹. Table 1 summarizes the results obtained for different frequencies and magnetic induction peak levels. Let us recall that equality ${\mathcal{L}}_E+{\mathcal{L}}_H={\mathcal{L}}_T$ holds at the continuous level and thus, the difference among the losses ${\mathcal{L}}_E^h+{\mathcal{L}}_H^h$ and ${\mathcal{L}}_T^h$ (columns 6 and 7 in Table 1) is due to numerical approximation of the axisymmetric model.

$f$(Hz)	$B_m$(T)	${\mathcal{L}}_E^h$	${\mathcal{L}}_H^h$	${\mathcal{L}}_E^h+{\mathcal{L}}_H^h$	${\mathcal{L}}_T^h$	Total (exp.)a	Error 1(%)b	Error 2(%)c
25	$0.5$	9.91	126.23	136.14	134.60	121.26	12.27	11.00
	$0.9$	32.77	269.50	302.27	300.32	300.25	0.67	0.02
	$1.4$	91.37	585.73	677.10	676.60	638.93	5.97	5.90
150	$0.5$	48.25	146.10	194.34	186.49	167.55	15.10	11.30
	$0.9$	178.53	282.28	460.81	4570.42	459.16	0.36	0.46
	$1.4$	506.40	588.19	1094.58	1095.70	1090.18	0.40	0.51

Table 1.

Total losses. Numerical versus experimental results (J/${\mathrm{m}}^3$).

^a

Total losses computed experimentally.

^b

Relative error between experimental losses and ${\mathcal{L}}_E^h+{\mathcal{L}}_H^h$.

^c

Relative error between experimental losses and ${\mathcal{L}}_T^h$.

This table shows that the computed losses are close to the losses measured in the experiments. The largest differences are obtained for $B_m=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 $\theta$-independent. It is also assumed that the current density has the form $\mathit{J}\left(\mathit{x},t\right)\equiv J_z\left(r,z,t\right){\mathit{e}}_z$, that is, that it traverses the piece along its axial direction. Consequently, the magnetic field stands $\mathit{H}\left(\mathit{x},t\right)\equiv H_{\theta }\left(r,z,t\right){\mathit{e}}_{\theta }$. As proven in [25], the problem can be spatially reduced to the meridional section $\mathrm{\Omega }$ of the workpiece, and reads

$$\frac{\partial {\mathcal{F}}_S\left(H_{\theta },\xi \right)}{\partial t}{\mathbf{e}}_{\theta }+\mathbf{curl}\left(\frac{1}{\sigma }\mathbf{curl}\left(H_{\theta }{\mathbf{e}}_{\theta }\right)\right)=\mathbf{0}\text{in}\mathrm{\Omega }\times \left[t_0,T\right],$$E28

$$H_{\theta }\left(0,z,t\right)=0\mathrm{on}\left(0,L\right)\times \left[t_0,T\right],$$E29

$$H_{\theta }\left(R_{\mathcal{S}}\left(z\right),z,t\right)=\frac{I\left(t\right)}{2\pi R_{\mathcal{S}}\left(z\right)}\mathrm{on}\left(0,L\right)\times \left[t_0,T\right],$$E30

$$\frac{\partial H_{\theta }}{\partial z}\left(r,z,t\right)=0\mathrm{on}\left({\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\right)\times \left[t_0,T\right],$$E31

where ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$ denote the boundaries clamped to the electrical contacts, $L$ is the workpiece length in the axial direction, ${\mathcal{S}}_z$ any cross-section transversal to this direction and $R_{\mathcal{S}}\left(z\right)$ its radius. As in previous sections, ${\mathcal{F}}_S$ is the classical Preisach hysteresis operator and $\xi$ 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 ${\mathrm{\Gamma }}_{\mathrm{S}}$ of radius $R_{\mathrm{S}}$ carrying a surface current density ${\mathit{J}}_{\mathrm{S}}\left(\mathit{x},t\right)\equiv J_{\mathrm{S}}\left(r,t\right){\mathit{e}}_{\theta }$ is employed to approximate the coil. Let ${\mathrm{\Omega }}_c$ be the part to be inspected, $R_c$ representing its radius; moreover, let us denote ${\mathrm{\Omega }}_0$ the air surrounding it which is assumed to be artificially bounded in the $r$-direction, $R_{\infty }$ being its outer radius. Then problem reduces to (see [25])

$$\sigma \frac{\partial A_{\theta }}{\partial t}-\frac{\partial }{\partial r}\left({\mathcal{F}}_S{}^{-1}\left(\frac{1}{r}\frac{\partial }{\partial r}\left({\mathit{rA}}_{\theta }\right),\xi \right)\right)=0\text{in}\left(0,R_c\right)\times \left[t_0,T\right],$$E32

$$\frac{\partial }{\partial r}\left(\frac{1}{{\mu }_0}\frac{1}{r}\frac{\partial }{\partial r}\left({\mathit{rA}}_{\theta }\right)\right)=0\text{in}\left[\left(R_c,R_{\mathrm{S}}\right)\cup \left(R_{\mathrm{S}},R_{\infty }\right)\right]\times \left[t_0,T\right],$$E33

$${\left[\frac{1}{{\mu }_0}\frac{1}{r}\frac{\partial }{\partial r}\left({\mathit{rA}}_{\theta }\right)\right]}_{r=R_{\mathrm{S}}}=J_{\mathrm{S}}\mathrm{on}\left\{R_{\mathrm{S}}\right\}\times \left[t_0,T\right],$$E34

$$A_{\theta }\left(0,t\right)=0\text{and}\frac{1}{r}\frac{\partial }{\partial r}\left({\mathit{rA}}_{\theta }\right)\left(R_{\infty },t\right)=0,\forall t\in \left[t_0,T\right],$$E35

where $A_{\theta }$ denotes the azimuthal of component magnetic vector potential. In (34)${\left[\cdot \right]}_{r=R_{\mathrm{S}}}$ represents the jump across the surface $r=R_{\mathrm{S}}$. Notice that the source is the surface current density ${\mathit{J}}_{\mathrm{S}}$ which is equal to the jump of $\mathit{H}\times \mathit{n}$ at $r=R_{\mathrm{S}}$.

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 $\mathit{H}$ and $\mathit{B}$ have only one non-null component in the ${\mathbf{e}}_{\theta }$ direction in the circular case, and in the ${\mathbf{e}}_z$ 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 ${\mathrm{\Gamma }}_{\mathrm{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 ${\mathrm{\Gamma }}_{\mathrm{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 $\left(R_j,C_j\right)$ 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, $i\left(t\right)$ is the current intensity, $\eta \left(t\right)$ the Coulombic efficiency, $V\left(t\right)$ the voltage, $z\left(t\right)$ the SoC, $R_j$ the $j-$th polarization resistance, $C_j$ the $j-$th polarization capacitance, and $v_{R_j}\left(t\right)$ is the voltage between the ends of the the $j$-th resistor $R_j$. Thus, the problem to solve reads:

Figure 13.

Sketch of the 2RC model and parameters values.

Given functions$\widehat{\mathit{OCV}}$,$i\left(t\right)$and$\eta \left(t\right)$, constant parameters$Q,R_0,R_j$and$C_j$, and initial conditions$z_0$and$v_{R_{j,0}}$,$j=1,\dots ,J$, find functions$z\left(t\right)$,$V\left(t\right)$and$v_{R_j}\left(t\right)$, satisfying,

$$\frac{\mathit{dz}}{\mathit{dt}}\left(t\right)=-\frac{1}{Q}\eta \left(t\right)i\left(t\right)$$E36

$$\frac{{\mathit{dv}}_{R_j}}{\mathit{dt}}\left(t\right)+\frac{1}{R_j{C}_j}{v}_{R_j}\left(t\right)=\frac{1}{C_j}i\left(t\right),j=1,\dots ,J$$E37

$$v_{R_j}\left(0\right)=v_{R_{j,0}}j=1,\dots ,J$$E38

$$z\left(0\right)=z_0$$E39

$$V\left(t\right)=\widehat{\mathit{OCV}}\left(z\left(t\right)\right)-\sum _{j=1}^J{v}_{R_j}\left(t\right)-R_{0}i\left(t\right)$$E40

where function$\widehat{\mathit{OCV}}$gives the OCV from the state of charge$z$.

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

$$\widehat{\mathit{OCV}}\left(z\left(t\right)\right)={\int }_{S_z^{+}\left(t\right)}p\left(\rho \right)\mathit{d\rho }+{\mathit{OCV}}_{\mathit{min}}$$E41

where $S_z^{+}$ is defined in (9) and ${\mathit{OCV}}_{\mathit{min}}$ 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.

Acknowledgments

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.