Data-driven techniques to estimate parameters in a rate-dependent ferromagnetic hysteresis model

Physica B 407 (2012) 1394–1398

Contents lists available at ScienceDirect

Physica B journal homepage: www.elsevier.com/locate/physb

Data-driven techniques to estimate parameters in a rate-dependent ferromagnetic hysteresis model Zhengzheng Hu a, Ralph C. Smith a,, Jon M. Ernstberger b a b

Department of Mathematics and Center for Research in Scientific Computation, North Carolina State University, Raleigh, NC 27695, USA Department of Mathematics, 601 Broad Street, LaGrange College, LaGrange, GA 30240, USA

a r t i c l e i n f o

a b s t r a c t

Available online 8 July 2011

The quantification of rate-dependent ferromagnetic hysteresis is important in a range of applications including high speed milling using Terfenol-D actuators. There exist a variety of frameworks for characterizing rate-dependent hysteresis including the magnetic model in Ref. [2], the homogenized energy framework, Preisach formulations that accommodate after-effects, and Prandtl–Ishlinskii models. A critical issue when using any of these models to characterize physical devices concerns the efficient estimation of model parameters through least squares data fits. A crux of this issue is the determination of initial parameter estimates based on easily measured attributes of the data. In this paper, we present data-driven techniques to efficiently and robustly estimate parameters in the homogenized energy model. This framework was chosen due to its physical basis and its applicability to ferroelectric, ferromagnetic and ferroelastic materials. & 2011 Elsevier B.V. All rights reserved.

Keywords: Hysteresis Ferromagnetic materials Data-driven estimation techniques

1. Introduction The quantification of rate-dependent effects, including accommodation or reptation, magnetic after-effects, and eddy current effects, is crucial when employing magnetic transducers in high drive regimes. There exist a number of modeling frameworks that incorporate these effects, for materials acting in hysteretic regimes, including the magnetic models of Ref. [2], homogenized energy models (HEM) [9,10], extended Preisach models [5,7], and rate-dependent Prandtl–Ishlinskii models [1]. A critical issue when employing any of these frameworks concerns the robust determination of model parameters based on least squares fits to data. In this paper, we present a data-driven technique to obtain initial parameter estimates for the homogenized energy model for ferromagnetic hysteresis. We employ this framework due to its energy basis at mesoscales, its ease of implementation and relative accuracy for a range of operating conditions, the degree to which model parameters can be correlated with physical properties of measured data, and its ubiquity for ferromagnetic, ferroelectric (e.g., PZT), and ferroelastic (e.g., SMA) materials. We focus on rate-dependent effects and creep behavior associated with magnetic after-effects [4] since, for the transducer designs and data under consideration, eddy current losses have been minimized. The extension of these models and data-driven parameter estimation techniques to incorporate eddy currents is under current investigation.

 Corresponding author. Tel.: þ 1 919 515 7552; fax: þ1 919 515 1636.

E-mail addresses: [email protected] (Z. Hu), [email protected] (R.C. Smith), [email protected] (J.M. Ernstberger). 0921-4526/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2011.06.084

A short description of the homogenized energy model is provided in Section 2. In Section 3, we detail the data-driven techniques to efficiently and robustly estimate parameters in the homogenized energy model. In Section 4, we compare the model fits to experimental data for both steel and nickel rods.

2. Homogenized energy model for magnetic hysteresis As detailed in Refs. [9,10], the homogenized energy model for ferromagnetic materials is constructed in two steps. In the first, a Gibbs energy is balanced with a relative thermal energy to construct a local average magnetization relation. The effects of material and field nonhomogeneities are subsequently incorporated by assuming that local coercive and interaction fields are manifestations of underlying densities rather than constant parameters. Stochastic homogenization in this manner yields a macroscopic model that is accurate for a variety of operating regimes and is efficient to construct and implement. For 1801 moment switching, we employ the Gibbs energy GðH,MÞ ¼ cðMÞHM,

ð1Þ

where H and M respectively denote the field and magnetization, and the Helmholtz energy is 8 1 > > ZðM þ MR Þ2 , M r MI > > 2 > > > <1 2 M Z MI cðMÞ ¼ 2ZðMMR Þ , > >  2  > > 1 M > > > : 2ZðMI MR Þ M I MR , jMjo MI :

Z. Hu et al. / Physica B 407 (2012) 1394–1398

Here Z ¼ dH=dM after switching and the remanent magnetization MR are two of the parameters to be identified. As detailed in Refs. [9,10], thermal relaxation is accommodated by balancing the Gibbs energy and relative thermal energy kT=V using the Boltzmann relation:

mðGÞ ¼ CeGV =kT ,

midpoint or Gaussian rules, can be used to approximate two integrals in Eq. (8). Following Ref. [6], we employ expansions

nI ðHI Þ ¼ c2

MeGðH,MÞV=kT dM

ð4Þ

yields the relation H þ HI

Z

MR

ð5Þ

ð6Þ

involving the likelihoods of moments switching from negative to positive ðp þ  Þ and conversely ðp þ Þ. As detailed in Ref. [3], we employed the likelihood relations

g1

g

1 , ð7Þ erfcxðH Þ pffiffiffiffi R 1 2 2 where erfcxðxÞ ¼ ex ð2= pÞ x et dt (the scaled complementary error function), H þ ffi ¼ g2 ðHc þp ðHffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ HI ÞÞ, ffi H ¼ g2 ðHc ðH þ HI ÞÞ, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g1 ¼ ð1=tÞ ð2V Z=pkT Þ and g2 ¼ ðV=2kT ZÞ. Here t is the relaxation time, or the reciprocal of the frequency at which moments attempt to switch. To incorporate the effects of polycrystallinity, material nonhomogeneities, and variable interaction fields, we subsequently assume that local coercive and interaction fields are manifestations of underlying densities rather than constant coefficients. This yields the macroscopic magnetization relation: Z 1Z 1 MðHðtÞÞ ¼ MðHðtÞ þHI ; Hc Þ  nc ðHc ÞnI ðHI Þ dHI dHc , ð8Þ

pþ ¼

erfcxðH þ Þ

0

,

p þ ¼

bk ck ðHI Þ,

ð10Þ

whose corresponding normal distribution has mean ln H cm and standard deviation sck ¼ 2k sc . As detailed in Section 3, one of the H cm (e.g., m¼2) can be obtained directly from the coercive field of the data. Once H c2 is determined, H c1 and H c3 can be determined by, for example, letting H c1 ¼ 1:5H c2 and H c3 ¼ 0:5H c2 . The interaction field basis functions are taken to be normal distributions: ! H2 1 pffiffiffiffiffiffi exp  I2 , ck ðHI Þ ¼ ð12Þ 2sIk sIk 2p

for the local average magnetization. Here x þ denotes the fraction of positively oriented moments which evolve via the differential equation x_ þ ¼ ðp þ  þ p þ Þx þ þ p þ ,

k1 X

where the basis functions fk,m ðHc Þ and ck ðHI Þ are lognormal and normal distributions and the coefficients ak,m and bk are determined through a least squares fit to data. The constants P P c1 ¼ ð kk1¼,Mka0 ,m ¼ 1 ak,m Þ1 and c2 ¼ ð kk1¼ k0 bk Þ1 ensure integration to unity. Examples of the two densities are shown in Fig. 1. The coercive field basis functions are taken to be lognormal distributions: ! 1 ðln Hc ln H cm Þ2 pffiffiffiffiffiffi exp  fk,m ðHc Þ ¼ , ð11Þ 2s2ck sck 2pHc

1

M ¼ 2MR x þ þ

ð9Þ

k ¼ k0

and MI

ak,m fk,m ðHc Þ,

k ¼ k0 ,m ¼ 1

ð2Þ

MI

Z

kX 1 ,Ma

nc ðHc Þ ¼ c1

where k is the Boltzmann constant and T is the temperature in degree Kelvin. The reference volume V reflects the definition of magnetization as magnetic moments per unit volume and yields a relative thermal energy based on the Gibbs energy density. Approximation of the relations Z 1 /M þ S ¼ c MeGðH,MÞV=kT dM ð3Þ

/M S ¼ c

1395

with mean 0 and standard deviation sIk ¼ 2k sI for sI where sI can be obtained from the data.

3. Data-driven parameter estimation 3.1. Optimization method The parameters to be estimated are q ¼ fZ,MR ,H c2 , sI , sc , g1 , g2 , ak,m , bk g:

1

ð13Þ

To ensure the accuracy of nested, biased minor loops, we take k ¼ 2,1,0,1, and m ¼ 1,2,3. For simpler operating regimes, fewer values can be used.

where nc ðHc Þ and nI ðHI Þ are densities associated with the coercive and interaction fields. Various quadrature techniques, including

–7 Interactive Field Density 4 x 10

–6 Coercive Field Density 1.4 x 10

3.5

1.2

3

1

2.5

0.8

2 0.6

1.5

0.4

1

0.2

0.5

0

0

1

2 MV/m

3

4

0 –4

–2

0 MV/m

2

4

Fig. 1. (a) Sampled lognormal density nc with 80 points. (b) Sampled normal density nI with 80 points.

Z. Hu et al. / Physica B 407 (2012) 1394–1398

For a data set with Nd measurements, the least squares functional to be minimized is JðqÞ ¼

Nd X

2

½MðHðti Þ; qÞMi  ,

ð14Þ

 For a fixed field H0, we note that Eq. (6) has the analytic solution: x þ ðt,H0 Þ ¼

i¼1

3.2. Data-driven estimation The optimization problem (14) is difficult due to the fact that the functional exhibits multiple local minima and is shallow near the global minimum (if it exits). Hence it is critical to develop data-driven techniques to provide reasonable initial parameter estimates. Major loop data. Major loop data can be used to obtain initial estimates for MR, H c2 , Z, sI , sc .

 The remanent magnetization MR and the coercive field H c can 2

  



be obtained directly from the major hysteresis loop as shown in Fig. 2(a). The parameter Z ¼ dH=dM is obtained using the maximal polarization after switch at point a and MR and is approximated by Ha =ðMðHa ÞMR Þ. The parameter sI is determined from the point b in Fig. 2(a). This is because switching begins when nI ðHI Þ intersects nc ðHc Þ with 2sI specifying the 94.5% confidence level. After determining sI , we obtain an estimate for parameter sc by using Hf, the magnetization where switching completes for increasing H. As shown in Fig. 2(b), for increasing input magnetic field, nI ðHI Þ quits intersecting nc ðHc Þ at a 94.5% confidence level, which yields sc ¼ 1=2lnðHf 2sI Þ=H c2 . Initial values for the coefficients ak,m and bk are simply taken to be 1.

Creep data. Creep data collected at two fixed fields can be used to determine initial estimates for the parameters g1 and g2 .

eðtt0 Þ½p þ ðH0 Þ þ p þ  ðH0 Þ :

ð15Þ

To obtain an algorithm to estimate initial values of g1 and g2 , we first assume that the decay behavior of x þ can be adequately approximated by neglecting interaction fields and consider the behavior about the measured coercive field Hc or Hc . It then follows that p þ  0 if H o 0 and p þ   0 if H 40. Hence Eq. (15) reduces to ( 1 þ½x þ ðt0 ,H0 Þ1eðtt0 Þp þ ðH0 Þ , H Z 0 x þ ðt,H0 Þ ¼ H o 0: x þ ðt0 ,H0 Þeðtt0 Þp þ  ðH0 Þ , We secondly assume that M and x þ have the same creep rate; that is x þ ðt3 ,H0 Þx þ ðt2 ,H0 Þ d3 M3 M2 ¼ ¼ : x þ ðt2 ,H0 Þx þ ðt1 ,H0 Þ d2 M2 M1

ð16Þ

 These two assumptions yield the relations g1 erfcx½g2 ðH0 Hc Þ

¼

  1 d2 ln , Dt d3

H o 0,

x 10–3 1

0.9

Magnetization (MA/m)

where Mi are magnetization measurements corresponding to field inputs Hðti Þ. We note that an important property of the model is that analytic Jacobians can be computed which significantly facilitates optimization.

  p þ ðH0 Þ p þ ðH0 Þ þ x þ ðt0 ,H0 Þ p þ ðH0 Þ þ p þ  ðH0 Þ p þ ðH0 Þ þp þ  ðH0 Þ

0.3 0.5

Density

1396

–0.3 νc

νI

 Let t0 designate the time at which the field is fixed and let M1, M2 and M3 be measured magnetization values at equally distanced subsequent times t1, t1 þ Dt, and t1 þ2Dt, and let d2 ¼ M2 M1 , d3 ¼ M3 M2 .

–0.9 –15

–10

–5

0

5

10

15

20

0 25

Magnetic Field (kA/m) 0.9

M

a

f

M

b

H –Hc

Hc

Hc

Ha

Decreasing H

Hf

H

Increasing H

Magnetization (MA /m)

MR

0.3 0.3

0.25

–0.3 0.2 25

30

35

Data

2 σI

H

c

~ H

c

Hc exp(2σc) 2σI Fig. 2. Experimental data marked by points used to obtain initial estimates for Z, MR, H c2 , sI , sc .

Model

–0.9 15

20

25

30

35

Time (s) Fig. 3. Model fit to the steel creep data in the (a) phase and (b) time domain. The vertical lines in (b) indicate the time where the magnetic field is held constant.

Z. Hu et al. / Physica B 407 (2012) 1394–1398

  g1 1 d2 ln , ¼ Dt d3 erfcx½g2 ðH0 E c Þ

H 40:

1397

ð17Þ

x 10–3 2

0.5

 For the fixed field H0 near H c , this yields g1 ¼ ð1=DtÞ

lnðd21 =d31 Þerfcx½g2 ðH1 H c2 Þ lnðd22 =d32 Þerfcx½g2 ðH2 H c2  ¼ 0

ð18Þ

for g2 . Note that Eq. (18) can be easily solved for g2 by plotting the solution as a function of g2 and approximating the root or by using a symbolic routine; e.g., the symbolic MATLAB command solve.m.

1.5

νc

–0.5

4. Experimental results

Magnetization (MA/m)

0.9

0 5 Magnetic Field (kA /m)

10

0.5

0

0

0.3

Data Model

–0.5 50

100

150

200

Time (s)

–0.3

–0.9 –15

Fig. 5. Model fit to the nickel biased minor loop data in the (a) phase and (b) time domain.

–10

–5 0 5 Magnetic Field (kA /m)

10

15

0.9

Magnetization (MA/m)

–5

νI

0.5

Magnetization (MA /m)

Here we illustrate the data-driven algorithms for experimental data from: (1) a cylindrical steel rod measuring 2.0 in. by 0.25 in., yielding the magnetization profiles shown in Figs. 3(a), 4(a), and (2) a rod comprised of Nickel 200 [8,10], yielding the magnetization profiles shown in Fig. 5(a).

1

0

Density



Magnetization (MA/m)

2

lnððM2 M1 Þ=ðM3 M2 ÞÞ. If two fixed fields available (e.g., H1 ,H2 o0Þ, with the corresponding magnetization differences d21 ,d31 , and d22 ,d32 , elimination of g1 yields the relation

0.3

–0.3

Data Model –0.9 15

20

25 Time (s)

30

35

Fig. 4. Model fit to the biased minor loop data in the (a) phase and (b) time domain.

The initial and estimated parameter values (excluding a’s and b’s) for the steel data are listed in Table 1. Note that the steel creep data shown in Fig. 3 was collected at a single fixed field, so we can only obtain initial value for g1 . The initial value for g2 is chosen so that the initial fit to data is reasonable. Also, the identified g1 value is significantly larger than the initial value suggested by the datadriven techniques detailed in Section 3.2. This is due to the small creep rate of steel, as shown in the inset of Fig. 3(b). The convergence of the optimization algorithm to obtain this larger value illustrates its robustness. In Fig. 3, the model constructed through a least squares fit to the creep data with extra weight to errors occurring immediately after the field was held constant is compared with magnetization data in the (a) phase domain, together with the two densities, and in (b) time domain with a detailed view of the creep region. It takes less than a minute to finish the optimization, and the maximum error between the data and the model is 0.078 MA/m, or 9% of the saturation magnetization. The estimated parameters are then used as initial values to fit the biased minor loop data, and the results are shown in Fig. 4. The estimated densities are similar to the ones shown in Fig. 3(a) and are not shown. It takes about 1.5 min to optimize and the maximum error between the data and the model is 0.043 MA/m, or 5% of the saturation magnetization. To further demonstrate the capability of the model, the same techniques are applied for a second data set collected from a nickel rod. The initial and estimated parameter values (excluding

1398

Z. Hu et al. / Physica B 407 (2012) 1394–1398

Table 1 Initial and identified parameter values obtained through a least squares fit to the creep and biased minor loop steel data. Parameters

Z

MR (MA/m)

H c2 (kA/m)

sc

sI (kA/m)

g1

g2

Initial Identified with creep data Identified with minor loop data

3.8368  104 1.0540  105 1.1474  105

0.4438 0.7091 0.7416

1.554 NA NA

0.9309 NA NA

2.5 NA NA

9.0775  10  1 1.1053  103 5.534  102

1.0  10  2 1.0403  10  2 1.1599  10  2

Table 2 Initial and identified parameter values obtained through a least squares fit to the biased minor loop nickel data. Parameters

Z

MR (MA/m)

H c2 (kA/m)

sc

sI (kA/m)

g1

g2

Initial Identified

1.2854  105 1.1947  105

0.3528 0.3653

2.592 NA

0.4197 NA

1.0 NA

10 1.4320

1.0  10  2 8.7036  10  2

a’s and b’s) are listed in the second row of Table 2. The resulting magnetization data and model fit are shown in Fig. 5. We observe that the model accurately characterizes both the major loop and biased minor loops. 5. Conclusions We note that the homogenized energy model, together with the data-driven initialization techniques, accurately and efficiently characterizes rate-dependent behavior of both biased minor loop operation and the magnetization response to a fixed input field following dynamic operation. The former is important for accurate material/device characterization whereas the latter is crucial for characterizing devices required to maintain a fixed magnetization or position for timescales comparable to the relaxation time of the material. Because the after-effects or creep rates in steel are relatively small, algorithms to provide initial estimates for relaxation parameters are relatively insensitive compared to analogous algorithms for piezoelectric compounds. However, this insensitivity also provides robustness for optimization algorithms which provide accurate fits while converging quickly. Finally, we note that the extension of these techniques to include eddy current effects is under investigation.

Acknowledgments This research of ZH was supported through the NSF Grant DMS-0636590, EMSW21-RTG program while the research of RCS was supported in part by the Air Force Office of Scientific Research through the Grant AFOSR-FA9550-08-1-0348. References [1] [2] [3] [4] [5] [6]

[7] [8]

[9] [10]

M. Al Janaideh, C.-Y. Su, S. Rakheja, Smart Mater. Struct. 17 (2008) 1. G. Bertotti, Hysteresis in Magnetism, Academic Press, San Diego, CA, 1998. T.R. Braun, R.C. Smith, J. Intell. Mater. Syst. Struct. 19 (11) (2008) 1295. S. Chikazumi, Physics of Ferromagnetism, second ed., Carendon Press, Oxford, 1997. E. Della Torre, Magnetic Hysteresis, IEEE Press, New York, 1999. J.M. Ernstberger, High Speed Parameter Estimation for a Homogenized Energy Model, Ph.D Dissertation, North Carolina State University, Raleigh, NC, 2008. I.D. Mayergoyz, Mathematical Models of Hysteresis, Springer-Verlag, New York, 1991. A.P. Mortensen, Characterization, Modeling and Dynamic Implementation of Terfenol-D Particulate Composites, M.S. Thesis, The Ohio State University, Columbus, OH, 2005. R.C. Smith, Smart Material Systems—Model Development, SIAM, Philadelphia, 2005. R.C. Smith, M.J. Dapino, T.R. Braun, A.P. Mortensen, IEEE Trans. Magn. 42 (7) (2006) 1747.