Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials

ARTICLE IN PRESS

JID: AMC

[m3Gsc;May 13, 2017;9:30]

Applied Mathematics and Computation 0 0 0 (2017) 1–17

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials Miroslav Novak∗, Jakub Eichler, Miloslav Kosek Technical University in Liberec, Studentska 2, 46117 Liberec 1, Czechia

a r t i c l e

i n f o

Article history: Available online xxx Keywords: 2d numerical derivative of the experimental data Ferromagnetics Preisach model First order reverse curve

a b s t r a c t The Preisach model can be used for detailed analysis of devices based on ferromagnetic materials, if its parameter, its weighting function, is well-known. Usually the weighting function is approximated by analytical formula. The second approach is to determine it directly from experimental data. Most widely used method to obtain the weighting function is the first order reversal curve method that is based on two partial derivatives of measured magnetization using a special excitation pattern beginning from deep material saturation. Since the derivative enhances the experimental error, a precision experiment is necessary. Furthermore, it is not easy to achieve the deep saturation with the required signal pattern. Therefore sophisticated data processing followed, in order to reduce experimental errors before performing the numerical derivative. The paper concerns measurements errors caused by insufficient saturation and also problems of negative values of the weighting function, partially due to the noise. Irrespective of measurement errors, the agreement between model and experiment is good and fully acceptable in technical praxis. © 2017 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

1. Introduction Since circuit elements using a ferromagnetic material are often used (as transformers first of all), a lot of work was devoted to make reliable models of these circuits. The most universal model is the Preisach one [1,2]. The key part of the model is the weighting function that characterizes magnetic properties of the material. It can be derived from measurements of the material. The basic procedure is presented in [1]. The two-dimensional (2D) weighting function is theoretically obtained by two partial derivatives of magnetization surface with respect to its variables [1]. However, due to experimental data inaccuracy, the derivatives exhibit large deviations [3,4]. That is the reason why this method is rarely used, for instance in [5,6]. No flat function was found for the soft magnetic materials [5]. The same is valid for grain oriented silicon steel [6]. This approach is termed as an experimental one. It follows from physical considerations that the Preisach weighting function has a strong maximum and its position is determined by the coercive magnetic field strength. Therefore the function can be approximated analytically by some of probability density functions. The method, called analytical, is used practically in the extended literature [5–9], for instance. ∗

Corresponding author. E-mail addresses: [email protected] (M. Novak), [email protected] (J. Eichler), [email protected] (M. Kosek).

http://dx.doi.org/10.1016/j.amc.2017.05.017 0 096-30 03/© 2017 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC 2

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 1. a) Elementary dipole, b) Elementary dipoles (hysterons) distribution in Preisach plane.

Usually the 2D weighting function is considered as the product of two 1D functions. Their parameters are determined from the major hysteresis loop by the least squares approach. As prototypes for the 1D function, authors consider: normal distribution, Cauchy distribution and lognormal distribution [8]. The best results from extended calculations are found for the product of the Cauchy–Cauchy distribution. The normal distribution is used also in [5,9]. Among others, the Lorentz function is considered in [7] and a more complicated function is applied in [6]. The comparison of several methods belonging to both the categories (experimental and analytical) is in [7], given for the material with almost rectangular hysteresis loop. To be complete, advanced methods are mentioned. The multiple order reverse curves (MORC) are considered in [10]. The use of the major hysteresis loop only is the subject of an extended theory followed by some applications in [11,12]. Since the analytical method is applied to the major hysteresis loop, worse agreement can be expected, if a low excitation is used. Therefore, the typical experimental method [1] should be preferred, since the weighting function is derived from all the excitations. In order to reduce experimental errors, a lot of hysteresis loops were measured. The neighboring firstorder reversal curves (FORC) [1,2] differ by about 0.4 A/m from one another, as for the excitation. The preliminary results are presented in this paper. 2. Theory Since the Preisach model is not frequently used, we will mention its principle at the beginning. The model is based on the assumption that the material consists of an infinite number of magnetic particles whose statistical distribution parameters determine permeability and therefore the shape of the hysteresis loop. The elementary magnetic particles are modeled in form of ideal dipoles (hysterons) that exhibit a rectangular hysteresis loop, see Fig. 1a). The hysteron has two parameters: switching magnetic field strengths Hu and Hd .1 The absolute value of the dipole momentum ms is not important; therefore it can be set to one. If the external field H becomes greater than level Hu (H > Hu ), the dipole magnetic momentum is switched to upper direction +ms . Analogically, if the external field H goes under the level Hd (H < Hd ), the momentum direction is down −ms . In the interval between the both levels (Hd < H < Hu ), the momentum direction depends on its previous state. If the dipole momentum was in the up direction (the external field decreases under Hu ), it remains polarized up. If the external field increases from the position under Hd , the dipole remains polarized down. In the Preisach model, the elementary dipoles are systematically arranged according to both levels Hu and Hd . Therefore the model is two dimensional, see Fig. 1b). Since the physical condition Hd ≤ Hu is valid, only the upper triangle of the whole matrix is occupied by hysterons. It is usually termed the Preisach triangle. The magnetization M(t) of the modeled sample is given by formula

1

Some authors use symbol α and β , or Hd and Hs .

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

ARTICLE IN PRESS

JID: AMC

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Preisach model, Mtot= −45 kA/m

Excitation waveform 100

100

50

50 H (A/m)

0

0

u

Field strength H [A/m]

3

−50

−50

−100

−100 0

50

100 Time t [ms]

150

−100

−50

0 H (A/m)

50

100

d

Fig. 2. Application of the Preisach model for the increasing external field.

M (t ) =

 Hu ≥Hd

ˆ (Hu , Hd )H (t )dHu dHd w(Hu , Hd )m

(1)

ˆ (Hu , Hd ) forms the elementary magnetic momenwhere H(t) is the input magnetic field strength. The hysteresis operator m ˆ (Hu , Hd ).H(t) means either tum for given field strength H(t) according to the explanation in Fig. 1a, for instance. Therefore, m +ms or –ms depending on the values of Hu , Hd and the direction of change of H(t): increase or decrease. The weighting function w(Hu , Hd ) that specifies the exact elementary magnetic momentum is determined by the material and it should be found. The integration is over all possible magnetic field strengths, from negative to positive saturation. The hysteron coordinates in the Preisach model are the following. The field Hd and its index increases in the row from left to right. However, the field Hu and its index increases in the column from the bottom upwards. Therefore, the direction of the row index in the Preisach model is opposite to the convention used in a mathematical matrix. The hysteron indices are in the sequence Hu , Hd . However, in the graphs, the field Hd is on the horizontal axis, while the Hu is on the vertical axis. According to this explanation, the hypotenuse in the Preisach model is termed the main diagonal and the normal to it is the secondary diagonal. Only half of it is important. Inspecting Fig. 1, we simply reveal that no hysteresis is on the main diagonal (Hu = Hd ). On the secondary diagonal, the loops exhibit central symmetry (Hu = −Hd ). The dynamics of the Preisach model can be explained, if we use the time dependent external magnetic field. In the following two figures, the field changes linearly, but its high and low limits decrease. We consider that the magnetic field was switched some time before and now its strength increases, see Fig. 2. The hysterons polarized upwards are marked by asterisks and the empty circle is used for downwards polarization. The instantaneous external field strength is given by the horizontal line in the Preisach triangle and moves from bottom to up. All hysterons under the instantaneous level are switched up. The hysterons above the level remain unchanged. Since the magnetization is the sum of the dipole momentums, it increases for the growing external field strength. The second case is for that time moment, when the magnetic field strength decreases, see Fig. 3. The instantaneous external field strength is given by the vertical line in the Preisach triangle and moves from the right hand side to the left. All hysterons on the right hand side of the instantaneous vertical level are switched down. The hysterons on the opposite side of the vertical level remain unchanged. Since the magnetization is the sum of dipole momentums, it decreases for the decreasing external field strength. Mathematically, the dynamics of the Preisach model can be given by the integral formula [2]. Since we use a digital approach, the formula can be written in the form

M (t ) =

N  i 



 



ˆ Hui , Hdj H (t )Hu Hd w Hui , Hdj m

(2)

i=1 j=1

where H(t) is the input external magnetic field strength as an excitation, M(t) is the output magnetization as a response. Indices of the individual hysterons are (i, j), N is the (maximum) number of hysterons in the row or the column, Hu and Hd are the increments of the field strength, Hu = Hd . Other used symbols are digital forms on analogue ones explained in detail under the (1). The straightforward method for the determination of the weighting function is the Systematic Minor Loops (SML) method that systematically measures the magnetization in the Preisach triangle. It can be shown simply from the digital model (2) application that the demagnetized state is not defined exactly. The better approach is to use the state of negative saturation as the initial state. The starting field strength is therefore negative. The harmonic field is added to the initial negative Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC 4

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 3. Application of the Preisach model for the decreasing external field.

Fig. 4. Simulated important partial hysteresis loops. The curve order develops from the bottom to the top.

field (see Fig. 4 and subsequent figures), and the time dependence of the magnetization is measured and partial hysteresis loops are formed. In order to explain the idea, in Fig. 4 there are several important partial hysteresis loops obtained by simulation using the analytical weighting function from the Cauchy distribution. The shape of loops strongly depends on the used parameters of the analytical function. In the upper left corner inset of Fig. 4, there is the corresponding excitation waveform. It does not differ significantly in the amplitude, but its effect to partial loops is dominant. The maximum excitations are also in the Preisach triangle in the lower right corner. The key excitations are near the triangle center. The increasing parts of the partial hysteresis loops are identical. It follows from the dynamic model in Fig. 2 and from the Preisach triangle in Fig. 4. The horizontal line goes thru identical hysteron positions excepting horizontal lines above the last excitation. The only difference is, therefore, in their end points. However, the decreasing parts of the partial hysteresis loops have key importance. Let us denote them by symbols Mi (Hdj ), j = 1,2, …, i, for given Hui . The curve is the basic part of the computation. It is termed the First Order Reverse Curve (FORC).2 The FORC for the ith row is sampled into i equidistant samples by the step of Hd and the samples are inserted into ith row of the Preisach triangle.

2 The terminology is not stable. Some authors use Everest function instead FORC. Its two partial derivatives are termed distribution function instead of weighting function.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

ARTICLE IN PRESS

JID: AMC

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

5

Fig. 5. Electrical scheme of the apparatus. Table 1 Toroid geometrical and electrical parameters. Parameter Symbol and unit

Outer diameter D (mm)

Inner diameter d (mm)

Height h (mm)

Primary turns N1

Secondary turns N2

Field step H (A/m)

Minimum field strength Hmin (A/m)

Large toroid Small toroid

220 85

160 80

30 10

100 500

300 300

0.786 4.1

−686 −3900

In the case of the analogous basic Eq. (1), the formula for the weighting function w(Hu , Hd ) can be derived [1]



w Hu , Hd





2 1 ∂ M Hu , Hd =− 2 ∂ Hu ∂ Hd



(3)

Since we use a discretized model, numeric partial derivatives must be used. It is well-known that the numeric derivative enhances the experimental noise; therefore many authors [5–9] preferred analytical approximation of the weighting function to theoretically exact calculation [7]. 3. Computations The Preisach model cannot be applied without the use of a computer, if we neglect the experiment that is also computer controlled, however. All the computations were made in MATLAB because of efficient matrix operations, complex number arithmetic, simple and efficient graphics, etc. The application of the Preisach model itself, using formula (2), is very simple. The efficient matrix processing in MATLAB makes it possible to write the Preisach model function that contains only several rows of code. Much more complicated is the determination of the weighting function. In order to improve data for the numeric derivative in (3), several steps should be made: data mining from careful experiment, data preprocessing in order to reduce noise and other errors in data, creation of FORC surface in the Preisach triangle, performing the numeric derivative of FORC surface in order to get the weighting function and verification of the model correctness. Below, the steps are dealt in detail. 3.1. Data mining In the experiment we used standard approach [13] in accordance with IEC 60,404–4. The FORC exceptions mentioned in the previous section were modified to improve the measurement accuracy. Each wave of FORC waveform was measured separately so that the complete dataset consists of 1800 measurements. The harmonic current source was used instead of the voltage source. The frequency was 1 Hz, in order to reduce eddy currents and to ensure harmonic current excitation. The partial loop measurement started from the level close to the negative saturation. Due to small mathematical interest, the experiment is not presented in details. The scheme of the apparatus is in Fig. 5. The current source is made by the simplest way from the programmable voltage source (PS) Kikusui PCR20 0 0 LA by adding resistances R1 = 45  and R2 = 15  that have negligible inductance. The capacitor C = 7 μF is used for filtering of the noise coming from the programmable source. A core of the toroid transformer Tr is made of the tested material. The exciting harmonic current is derived from the voltage Vc and the response, the secondary voltage, is measured by another voltmeter Vsec . The measurement is fully automated by the control program written in MATLAB that controls namely the programmable voltage and reads the data from digital voltmeters. Time dependence of the primary current and the secondary voltage was scanned automatically at a high speed (200 kSps) by the computer controlled ADC NI USB 6212. Using toroid parameters (geometry and number of turns), the field strength can be calculated from the primary current, while the flux density is obtained by the integration of the secondary voltage. Two toroids were used, large and small, with parameters given in Table 1. A large toroid cannot ensure harmonic excitation for the high fields and the field was not uniform in it. On the other hand, the measurement was more accurate for low excitation. In the small toroid, among others, the possibility of an efficient negative saturation was investigated. The step H of the excitation field is limited in its minimum that can be defined from the given computer controlled voltage source resolution and the toroid winding. It is known [14] that the higher field strength step lowers maximum of the Preisach weighting function and causes the quantization noise. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 6. Very low exciting field strength and its filtering by FFT. Also the start and finish points of the waveform are shown by crosses. The total field strength change is about 4 A/m.

Fig. 7. Smoothing of the secondary voltage using FFT filtering. Starting and finishing peaks are well visible.

3.2. Data preprocessing Main tasks of data preprocessing are concentrated mainly to low excitations and should ensure the following basic requirements, among others: • • • •

Reduction of noise and other disturbing effects. Finding a useful signal. Elimination of drift of the secondary voltage. Integration of the secondary voltage.

Two basic methods for noise and other disturbing effects reduction is the filtering by the Fourier transform and the running average filter. As for the Fourier transform, we used Fast Fourier Transform (FFT) and low pass filtering of the FFT spectrum. The method was very fast, however, in some cases, deviation in the starting and finishing part of the waveform appeared. As for running average we used the Savitzky–Golay function from the MATLAB library [15]. In comparison with FFT it is usually slower. Smoothing was applied to the exciting field, first of all. As an example, two periods of very low harmonic current were applied to negative bias, see Fig. 6. The total change of the field strength is about 4 A/m. The effect of filtering by FFT is shown by the white line in Fig. 6. Almost all noise is eliminated. Total measurement time is 10 s and only small part of it (20%) belongs to the excitation field. The ideally constant starting and finishing parts should be used for the reduction of drift in the measured secondary voltage and in the magnetic flux calculated by means of the numeric integration. Unfortunately, the part without excitation is not ideally constant, and small waves are visible. The source of parasitic waves is the programmable voltage source PS in Fig. 5. The waves are triangular with a frequency of 1.06 Hz and they come from the internal offset compensation of the PS. The similar smoothing procedure can be applied to the secondary voltage. The result for considerably higher excitation (in comparison with Fig. 6) is in Fig. 7. In this case, the noise is high; in reality its waveform is about 16 times higher.3 Also the peaks at the start and finish of active waveform are well visible.

3

Using the same scale does not allow seeing the effect of filtering.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

[m3Gsc;May 13, 2017;9:30] 7

Fig. 8. Drift of flux density waveform reduction.

A small and approximately square wave is found in the expected constant starting and finishing areas of the secondary voltage in Fig. 7. It was found that these errors have the same reason as in Fig. 6. Both of them have the same period and are shifted by 90°. This square wave arises naturally as the derivative of the parasitic triangle wave with frequency of 1.06 Hz generated by the used programmable voltage source due to the Faraday’s induction law at an inductive reactance in the primary part of the toroid. Practically, this effect is disturbing especially for low excitations, where the experimental errors are dominant. The second task was to find the start and finish points of the field strength waveform. The simplest way was to use sharp peaks both in the current and the voltage waveforms. They are due to the used voltage source. The peaks on the voltage waveform are well visible in Fig. 7. However, this simple approach does not work well for very low excitations and, furthermore, the method is device dependent. Therefore, we used more robust approach; finding points that are by a small value higher than the DC part. The points in the starting and finishing parts are shown in Fig. 6 as white crosses. The distance of 1.81 s. between them is lower than the expected total time of 2 s. The actual waveform start and finish points were found by computing the average of the experimental starting and finishing times. Then we get the waveform time center and by subtracting and adding one period of time (1 s.) to the center, we find the actual active waveform. The procedure is precise, since the time is measured with a very high accuracy. The found start and finish points were compared with the results of methods using the peaks of excitation or response voltage waveform, when possible. The agreement was very good. Other alternative methods can be used, for instance finding maximum or minimum on the excitation waveform. Due to the noise, their application seems to be less correct. The secondary voltage is numerically integrated in order to get the desired flux density. Furthermore the integration has an efficient smoothing effect on noisy voltage in Fig. 7, for instance. Unfortunately, the unwanted effect of any integration is nonzero drift that is due to the small constant in the integrand. The waveform in the time domain is shifted up or down as in Fig. 8. The corresponding partial hysteresis loop is therefore not closed. The basic general method for drift reduction is based on the theoretically constant initial and final parts of the measurement. They are approximated by linear regression, when the part with sinusoidal waveform is excluded. Then the approximating polynomial is subtracted from the data and the corrected waveform is numerically integrated. The result is given in Fig. 8 where the integration of origin and corrected voltage waveform are shown. A strict examination reveals that all three minimums do not lie exactly on a horizontal line. 3.3. Experimental errors During the data preprocessing, some significant experimental errors were revealed. The current source did not produce harmonic excitation in all cases. Comparison of excitation and response waveforms for the maximum excitation is shown in Fig. 9. The deviation of the field strength from its theoretical shape is evidently near the zero axis crossing. The deviation coincides with a sharp flux density change and can be explained by electromagnetic induction [13]. When the flux density increases (at a quarter of period), a voltage is generated in the primary winding that is directed against the current. As a result, the exciting field increases more slowly than expected. In the area near three quarters of period, the flux density decreases very fast. Also in this case the voltage oriented against the current is generated in the primary winding. The result is that the field strength decreases more slowly than expected. The complication is in some cases experimentally reduced by an air gap in the toroid [11], but other problems appear. For instance, the correct material parameters are not measured. The real partial hysteresis loops constitute another complication. Their central part is in Fig. 10. There are important differences from the simulated partial loops in Fig. 4. The increasing parts are not identical, but differ by several A/m. The most critical deviation is in the peaks of the loops, up to about a maximum of 30 A/m. The maximums of excitation and response are not identical. In the simulation in Fig. 4, this situation does not take place. The discrepancies can be explained by a non-uniform excitation field in the large toroid. In the small toroid the deviations are smaller, but the effect is disturbed by low experiment accuracy. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 9. The field strength and flux density at maximum excitation.

Fig. 10. Central part of partial hysteresis loops. Loops are numbered from bottom to top. Maximum excitation field is in the legend.

3.4. Saturation state in preisach model Another serious experimental problem is the response accuracy. In Fig. 11 there are neighboring curves of the flux density in the central part of the period for strong excitation. The curve order and the maximum of the sinusoid field strength are given in the legend. A real curve order is 1, 4, 2, and 3. The distance between the close curves in Fig. 11 is about 0.5 mT. For low excitation and especially for the small toroid, the curves similar to those in Fig. 11, usually cross one another. The points at a half period are collected for a large number of loops with strong excitation; they are given in Fig. 12. The curves are termed as cuts, since the points are for a lot of curves similar to those in Fig. 11 in the time of t = 0.5 T, where T is the signal period. This time is in the upper left inset of Fig. 12. Of course, the cuts may be taken for other times. In the ideal case, the curve should be monotonically increasing as the ones for the field strength.4 In fact, there are quasi-periodical changes and the total local change may be almost 5 mT. However the histogram of deviation in the lower right inset of Fig. 12 shows that the deviations are of random nature. We have examined a possibility of polynomial approximation and correction of individual responses to get a monotonic increase. It happens for a selected time, but for

4

They are not shown here, because of lack of place.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

9

Fig. 11. Neighboring flux density curves for maximum excitation. Curves are numbered for increasing excitation maximums as stated in the legend.

Fig. 12. Cuts of flux density for the time of half period and high maximum excitation.

other time there is none, or a very small improving. The correlation between different cuts exists, but not for a wide range of time. The complete Preisach model supposes the full saturation state as its limits. However, the condition of a deep saturation for some materials, including the investigated one – the grain oriented steel, cannot be achieved by standard experimental equipment for those complicated waveforms as the FORC method needs. The practical fields may be hundred times lower than required. In some papers, there are measurements with insufficient saturation and its typical consequences are listed below, which are mistakenly attributed to properties of the material being measured – 2 kA/m [16], 0.5 kA/m [17]. Fig. 13a, shows the complete Preisach triangle containing the saturation limits Hsat . The large theoretical triangle contains the full information. The smaller triangle defines areas that can be controlled experimentally.5 We termed it as a working triangle. If the field increases from minimum Hmin to maximum Hmax , all the hysterons under the horizontal line Hmax are polarized up. The area is colored and marked by the plus sign in Fig. 13a. Now, if the field decreases to Hmin , all the hysterons in the area on the right hand side of the vertical line Hmin are polarized down. Only a part of the area with 5

The triangle area does not reflect the reality. In real conditions the small triangle will not be visible.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC 10

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 13. Theoretical and practical saturation states in the Preisach model and the application of real excitation.

negative sign is in Fig. 13a, since it is drawn for positive quasi-saturation. The area marked by the question mark cannot be achieved experimentally. It can be designated as a forbidden area. However, its contribution to the total magnetization is of the same importance as from other parts. In order to proceed further, the assumption must be made as for the initial state of the total Preisach triangle to be able to define the state of the forbidden area. We assume an ideal demagnetization; therefore the boundary between negative and positive dipole momentums in the big triangle has a form of finest stairs along the secondary diagonal. The stair height is given by the step of the digitized magnetic field strength H. The initial state of the negative saturation for H = Hmin , achievable practically, is then in Fig. 13b. The magnetization has the minimum value of M0 . Now we consider application of the harmonic excitation waveform starting from the experimentally attainable negative saturation Hmin . The state, when the excitation field achieves its local maximum Hu , is in Fig. 13c. Then the excitation decrease to the negative attainable saturation follows. Now, three cases are possible: (1) the decrease finishes at negative saturation level Hd = Hmin , (2) the attainable saturation level was not reached, Hd > Hmin and (3) the excitation field goes below the saturation level Hd < Hmin . All three possibilities are schematically sketched in Fig. 13d through 13f. The expected case of achieving just the level of attainable negative saturation Hd = Hmin , is in Fig. 1d. The magnetization has a value of M = M0 as in previous Fig. 13b. Fig. 13e shows the case when the operating negative saturation field Hmin was not reached. Let us suppose for simplicity that Hd = Hmin + H, where H is the “width” of the column. The last column in the operating triangle remains polarized positively and the total magnetization is higher, M = M0 +M, where M is the dipole momentum of the hatched part in the column in Fig. 13e. It starts in the main diagonal and finishes on the horizontal line for the local maximum of the excitation field. The case of the excitation field running across the level Hmin is in Fig. 13f. For simplicity, we suppose the field Hd = Hmin – H. Then the column just below the level Hmin becomes negatively polarized. The column is hatched in Fig. 13f. The total magnetization is now M = M0 − M, where M is the total momentum in the hatched column in Fig. 13f. As for the value of magnetization change in this case, there are two possibilities depending on the previous history. If the lower level was reached for the first time, the column was negatively polarized and the magnetization change M is the sum of the hysteron momentums in the whole column, from the main diagonal to the stair curve in the forbidden area. This case is shown in Fig. 12f. In the second case, the excitation should be periodic and in the second and further periods, it reaches only a value of Hmin . Then the column just below Hmin is being processed repetitively. In the second and higher loop runs it is polarized positively only from the main diagonal to the horizontal line for the maximum excitation. The magnetization change M is therefore lower than in the first period. This case is similar to the case of no-reaching the negative saturation level in Fig. 13e. The difference in the magnetization change in these two cases is only one hysteron momentum, as for absolute values. These theoretically described effects take place in reality. We summarize important experimental errors that may be responsible for them. We do not concentrate on material contribution, as Berkhausen noise, for instance. For the large toroid the column width in the Preisach triangle is 0.786 A/m. The noise extremes at the very low excitation as in Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

11

Fig. 14. Illustrative simulation of partial hysteresis loops respecting effect of a negative starting pulse.

Fig. 6 differ by about 1 A/m. The practical low saturation level is not defined correctly before a useful signal is coming. Also the local minimum inside the useful signal is not defined exactly. Both these conclusions follow from Fig. 6. These effects are present also at a higher excitation. Besides the noise in a whole excitation, the starting parts of the waveform differ from one another. Their mean values can differ by 2 A/m, which is more than the column width. Therefore the responses from the neighboring excitations can differ significantly. Effects described above are random cases. But also systematic impulses can exist. The used instruments are the source of systematic disturbing effects. In the time instant at the start of the current useful signal, a small negative peak appears. At the end of the useful signal there is a positive peak. The peak height is proportional to the useful signal amplitude; therefore the peak for the very low excitation in Fig. 6 is not visible. These peaks can be used conveniently for the exact determination of the beginning and end of the useful signal. On the other hand the starting level of the excitation field is lower than it was set up by its constant initial part. The result is the hysteresis loop distortion caused by these parasitic peaks. These effects partially contribute to disorder of the responses in Fig. 11 and to a no-monotonic increase of the flux density in Fig. 12. The extended systematic simulation can reveal how much this effect is important. The effect of parasitic impulses from experimental equipment on hysteresis loop was simulated numerically using the Preisach model. We suppose a negative peak of one or more column widths, which is followed by an ideal periodic signal. Also for this special case the explanation in Fig. 13 can be used with a little modification. The negative peak makes the initial field under the value of Hmin ; therefore the initial state is shown in Fig. 13f, one or more columns are magnetized negatively up to the stair boundary. Since the pulse is very short, the final state after one cycle of excitation is similar to that in Fig. 13d. The magnetization change is the sum of dipoles in columns on the left hand side of the vertical line Hd = Hmin in Fig. 13f. The illustrative simulation in Fig. 14 confirms the conclusion above. The major loop for the larger triangle is the boundary curve of real measurable loops. The loop for excitation with a small negative parasitic pulse, for instance, starts in point A. The position of this point close to the bottom part of the major loop is due to negative magnetization, one of the negative polarized columns is in Fig. 13f. In point B, the loop is almost identical with the major loop. Its maximum in point C is also almost on the major loop. It is explained for simulation of the partial loops in Fig. 4. After the cycle run, the finishing point D is above the starting point A and the starting loop is not closed. The magnetization difference is the sum of momentums of dipoles that are polarized negatively by the parasitic pulse. The next ideal periodic excitation leads to closed hysteresis loops. The starting part goes from point D and only a small part of the loop is shown. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 15. Corrected FORCs for rows in Preisach triangle. Curves are ordered from bottom to up.

The key points are also in the Preisach triangle in the upper left inset of Fig. 14. The present state differs from an ideal one by the negative dipole momentums in the marked part of the column under the minimum ideal excitation. An important conclusion is the fact that insufficient saturation during FORC measurements distorted the flux density function. The consequence is incorrect weighting function. 3.5. FORC surface and numeric derivative The descending part of the hysteresis loops (FORCs) resulted from preprocessing must be converted to points of the Preisach triangle. The total number of measurements was 1800, which is also the number of rows (and columns) in the Preisach triangle. The instruments define the distance between neighboring rows H = 0.786 A/m, see also Table 1. This value is given by voltage resolution of the programmable power supply and the sample geometry. The points for the Preisach triangle were chosen from FORCs at multiples of the modified step H. Since in each original time domain FORC there are 105 points, the rows are defined much more accurately. The lowest distance between the measured samples in rows is 56 points, practically almost two orders higher than in columns. The decrease of the distance between curves requires important changes in the experimental apparatus, the time of measurement increases rapidly as well as the data volume. In theory, the quality of a weighting function depends on the step between the experiment rows that are inversely proportional to the number of measurements. The rule is valid in general: the lower the step, the better the weighting function. In some experiments only 100 rows are used [16,17]. Since we use a high number of rows, we can expect a better weighting function. The FORC were smoothed before the numeric derivative. The corrected FORC in the selected rows are in Fig. 15. The smoothing has little effect in this case. The smoothing has much bigger effect for columns. Smoothed FORCs in the selected columns are in Fig. 16. The parameter in the legend is the maximum field strength Hdmax in a row. The Hdmax is near the zero value. Irrespective of smoothing, low noise in curves even for higher excitations is evident, in comparison with Fig. 15. A very simple method of differences was applied to the smoothed FORCs. The alternative method is to use the Savitzky– Golay function in MATLAB to perform the partial derivatives. Then no preliminary smoothing is necessary. We tested algorithm efficiency on the analytical functions with a known derivative and the relative deviation between the calculation and theory was in an order of 10−8 . Therefore the algorithm accuracy is perfect. The complicated shape of the weighting function is confirmed by Fig. 17, for instance. The maximum position defines coercive force of the material and lies on the secondary diagonal, which is marked by a white line. Selected cuts along rows near the maximum are in the inset. Because the obtained weighting function is complicated and has unexpected features, a big amount of work was devoted to reduce the errors. These possibilities were investigated: • The signal filtration for removing the negative weighting function values. • FORC rows resampling in accord with the maximal Hu value. • Ensuring the symmetry with respect to the secondary diagonal. On the basis of simple physical considerations, the weighting function should be nonnegative [1]. However, the negative values occur due to experimental noise, FORCe rows scattering (see Figs. 11 and 12) and at a close proximity of the main diagonal due to rounded tips of the hysteresis loops. Removing the negative areas leads to worse agreement with experiment since they do not contain only noise and errors but also useful information. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

[m3Gsc;May 13, 2017;9:30] 13

Fig. 16. Selected FORCs in columns in Preisach triangle. Curves are ordered from bottom to top.

Fig. 17. Weighting function near its maximum.

The idea was that the filtration uses energy of negative peaks of the noise and transfers this energy to the surrounding points. Therefore, it is not possible to use median filters that lose this information. The averaging filter or more sophisticated methods from [4,18] are appropriate. The non-sufficient saturation in conjunction with the experimental noise during FORC curves measurement leads to the row curve scattering as mentioned above, Figs. 11 and 12. The result is a horizontal ridge in the obtained Preisach weighting function, see main part of Fig. 17, at values Hd = 13 A/m situated left from the main peak. It is possible to suppress this effect by the FORC lines re-sampling in accord with the maximal obtained Hu value. The complete removal of this problem is indeed not possible. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC 14

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 18. Weighting function after symmetrization.

Another trouble lies in the lateral weighting function peak, see Fig. 17 in position Hu = Hd = −18 A/m. The values at the main diagonal Hu = Hd of the Preisach plane are strongly dependent on the method of finding maximum field strength intensity Hmax in each measurement and on the proper placement of the decreasing branch data to the plane. On the grounds of experimental noise and rounding tips of the measured hysteresis loop, it is not an unambiguous task. The FORC hysteresis loops in proximity of the lower knee of the hysteresis loop are the most sensitive for the loop tips rounding, see Fig. 10. It leads to the fact that the lateral peak is originated from experimental errors. Fortunately, most of the above mentioned problems are manifested below the secondary diagonal of the Preisach plane. In addition, we expect a symmetrical weighting function in case of the grain oriented steel. Logical solution is to use only the data above the secondary diagonal and copy them over this secondary diagonal. The resulting weighting function is in Fig. 18. 4. Results The results can be summarized into two groups: application of the Preisach model and finding the weighting function from experiment. The MATLAB scripts for Preisach model applications are relatively simple. All the well-known hysteresis effects can be demonstrated by the use of constant or analytical weighting functions. Also several animations of the Preisach model dynamics were made6 and in Figs. 2 and 3 there are frames from one of these animations. The use of simulation in the investigation of the saturation states in part 3.4. was very useful. The main task was to determine the weighting functions from experimental data.7 Typical results will be presented in this chapter. It was already stressed in the part of the chapter dealing with computations that there are several possibilities how to improve (or damage) the weighting function and, therefore, the agreement with experiment. At present time, this issue is the subject of computational experiments. Only a few results are presented here. Computation experiments showed that these four corrections have a dominant effect: smoothing before the derivative, resampling the data according to Hmax , method of finding Hmax and Mmax in signal and finally the mirroring of the data over the secondary diagonal. 6 7

The scripts for animation can be got from authors under request. There are also several scripts using the analytical weighting functions, but it is not a scope of the paper.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

[m3Gsc;May 13, 2017;9:30] 15

Fig. 19. Hysteresis loop of large toroid.

Fig. 19 presents the simulated hysteresis loop of the large toroid using symmetrized weighting function. The agreement in the used scale is good. The details in the knee and the saturation areas, shown as insets of Fig. 19, also show a good agreement. The only exception was the negative saturation, where the experimental loop is not closed. The model loop is not closed too, which confirms the model quality. The hysteresis loop for the small toroid is in Fig. 20. Note that there is about 5 time higher excitation; nevertheless the saturation is not present. The weighting function is symmetrized again. The agreement is very good also in this case. The shape is similar to the large toroid including the unclosed experimental loop. No data preprocessing is necessary for a good agreement. But in details, the results from such a model are not smoothed. They exhibit ripples or noise of different values or amplitudes, which can be explained just by lack of data preprocessing. Furthermore, the step in excitation is about 5 times higher. Only hysteresis loops for maximum excitation are presented in Figs. 19 and 20. The agreement with experiment is good in all areas and therefore the model is fully acceptable for praxis in the case of the grain oriented steel. Also the agreement for lower excitations is good. 5. Discussion The application of the Preisach model, outlined above, can be made by mathematical tools only that contain a lot of algorithms. Fortunately, thanks to MATLAB, all the functions are either built-in functions or they have a form of free scripts. We found that they work perfectly. Therefore, there is no need to prepare one’s own algorithms. Usually we used at least two methods in the solution of each problem. The solution offering a better output was selected. Application of partial derivatives onto FORCs is the most general and the most precise method for determination of the weighting function. The function is determined from FORCs values in all points of the Preisach triangle. Therefore its application must be in agreement with experiment for any type of excitation. The assumption for it is an accurate experiment, which is a very difficult task. Although this condition was not valid in our work, the agreement with experiment was good or acceptable for excitations with positive maximums. For lower excitations there are not reliable values of the weighting function and the agreement is unsatisfactory. For very low excitations, the experimental loops are far from reality in some cases and therefore unacceptable. The frequent parallel method of modeling uses analytical approximation that is realized from only one loop for the maximum excitation. The parameters of the analytical function are obtained by comparison of modeled and experimental curves and then the calculation of the sum of squares of deviations is made. The goal is to find optimum parameters that ensure the minimum sum of squares, i.e. least squares. There are very sophisticated methods of finding the optimum. In Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC 16

ARTICLE IN PRESS

[m3Gsc;May 13, 2017;9:30]

M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

Fig. 20. Hysteresis loop of small toroid.

MATLAB, the fminsearch function is an efficient tool for this purpose. The advantage of the method is a simple experiment that has high accuracy due to a maximum excitation. The disadvantage is that the result can be used only for excitations that are close to maximum. For low excitations there is almost no information. Some difficulty constitutes physical interpretation of the weighting function. A simple explanation is that the weighting function represents probability of an individual hysterons that have given parameters, field strengths Hd and Hu . In the case of a continuous model, it is the probability density. In this case the weighting function must be nonnegative [18]. The condition of symmetry is only strictly valid on the secondary diagonal. It is a practical knowledge resulted from a lot of applications of the model. From a strict mathematical point of view, the theorem is not valid. For instance, two maximums symmetric with respect to the secondary diagonal result in a shifted positive and negative part of the hysteresis loop. Nevertheless symmetry application usually improves the agreement with experiment. The main reason for a difference between model and experiment lies in the individual FORC curve positions. Referring to Fig. 11, we can see that they are not in a predicted order and usually they also crosses. However, Fig. 12 shows that the curve shift may be systematic. Unfortunately, the correlation revealed that the similarity exists but only in a small time interval. Therefore, the difference cannot be reduced by mathematical procedures. A very important task is to find the source of the FORCs disorder. Probably it may be caused by fact that saturation cannot be reached and experimental noise is present, since the distance between the neighboring curves is of the order of 1 mT. Another question is the uniformity of the magnetic field in the toroid. The Amper’s law computations say that the nonuniformity of the magnetic field in the large toroid cannot be neglected. The hysteresis loops of the small and the large toroid differ one from another. Unfortunately, the experimental accuracy for the small toroid is lower in the presented data. 6. Conclusion The paper shows that the determination of the weighting function from experiment leads to practically usable results using the standard present-day experiment. On the other hand the experimental errors are significant and therefore the special data processing is necessary. The data processing starts with the reduction of experimental errors in the time domain. Because of errors in the application of the numeric derivative there are several choices how to improve the weighting function. Therefore, the present time scripts do not allow a routine use. The symmetrical weighing function enables us to use only data from the most valid area above the secondary diagonal. The positive side effect reduces the number of the measured waveforms. Therefore, the FORC loops with Hu > 0 are only needed. The computation experiments based on the Preisach model verified the effect of experimental errors on the resulting weighting function. Especially, impossibility to achieve perfect saturation was investigated. Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017

JID: AMC

ARTICLE IN PRESS M. Novak et al. / Applied Mathematics and Computation 000 (2017) 1–17

[m3Gsc;May 13, 2017;9:30] 17

In order to increase the speed of computation, all measured data should be in one file instead of 1800 files. It requires a large processor memory. Parallel computing also significantly increases computational speed and allows more computational experiments. Optimum parameters can be obtained from a lot of computational experiments. These parameters together with the improved experiment can be the base for a routine system applicable in laboratory or industry praxis. Acknowledgment The work was supported by the Student Grant Competition of Technical University of Liberec. References [1] G. Bertotti, I. Mayergoyz, The Science of Hysteresis, Vol. 1, 2 and 3, 1st ed., Elsevier, 2006 ISBN 978-0–2-369431-7, Chapter 3. [2] I. Mayergoyz, Mathematical Models of Hysteresis and Their Applications, 1st ed., Elsevier, 2003 ISBN 0-12-480873-5. [3] L. Stoleriu, A. Stancu, Using experimental FORC distribution as input for a Preisach-type model, in: INTERMAG 2006 - IEEE International Magnetics Conference, San Diego, CA, 2006 pp. 342–342, doi:10.1109/INTMAG.2006.376066. [4] R.J. Harrison, J.M. Feinberg, FORCinel: an improved algorithm for calculating first-order reversal curve distributions using locally weighted regression smoothing, Geochem. Geophys. Geosyst. 9 (5) (2008) 1–11, doi:10.1029/2008GC001987. [5] G. Consolo, G. Finocchio, M. Carpentieri, E. Cardelli, B. Azzerboni, About identification of scalar Preisach functions of soft magnetic materials, IEEE Trans. Magn. 42 (4) (2006) 923–926. [6] J. Füzi, Analytical approximation of Preisach distribution functions, IEEE Trans. Magn. 39 (3) (2003) 1357–1360. [7] L.L. Rouve, Th. Waeckerle, A. Kedous-Lebous, Application of Preisach model to grain oriented steels: comparison of different characteristizations for the Preisach function p(α , β ), IEEE Trans. Magn. 31 (6) (1995) 3557–3559. [8] P. Pruksanubal, A. Binner, K.H. Gonschorek, Determination of distribution functions and parameters for the Preisach hysteresis model, in: 17th International of Zurich Symposium on Electromagnetic Compatibility, 2006, pp. 258–261. [9] M. Cirrincione, R. Miceli, G.R. Galluzzo, M. Trapanese, Preisach function identification by neural networks, IEEE Trans. Magn. 38 (5) (2002) 2421–2423. [10] A. Stancu, Identification procedure for Preisach-type models bassed on FORC diagrams, J. Photoelectron. Adv. Mater. 8 (5) (2006) 1656–1659. [11] M. Ruderman, T. Bertram, Identification of soft magnetic B-H characteristics using discrete dynamic Preisach model and single measured hysteresis loop, IEEE Trans. Magn. 48 (4) (2012) 1281–1284. [12] P. Andrei, A. Stancu, Identification method analyses for the scalar generalized moving Preisach model using major hysteresis loops, IEEE Trans. Magn. 36 (4) (20 0 0) 1982–1989. [13] F. Fiorillo, Characterization and measurement of magnetic materials, Elsevier Series in Electromagnetism, Elsevier, Amsterdam, 2004 ISBN: 978-0-12-257251-7. [14] X. Zhao, D. Heslop, A.P. Roberts, A protocol for variable-resolution first-order reversal curve measurements, Geochem. Geophys. Geosyst. 16 (5) (2015) 1364–1377, doi:10.10 02/2014GC0 05680. [15] S.J. Orfanidis, Introduction to Signal Processing, Prentice-Hall, 1995 ISBN 0-13-209172-0, Chapter 8. [16] Chang Seop Koh, Young Hwan Eum, Sun-ki Hong, Pan-seok Shin, Comparison of the identification methods of preisach model for a grain-oriented electrical steel sheet, in: 12th Biennial IEEE Conference on Electromagnetic Field Computation, Miami, FL, 2006 pp. 99–99, doi:10.1109/CEFC-06.2006. 1632891. [17] T. Yamaguchi, F. Ueda, K. Takeda, Automatic measurement of preisach distributions in soft magnetic materials, IEEE Trans. J. Magn. Jpn. 2 (5) (1987) 412–413, doi:10.1109/TJMJ.1987.4549469. [18] M. de Wulf, L. Vandevelde, J. Maes, L. Dupre, J. Melkebeek, Computation of the Preisach distribution function based on a measured Everett map, IEEE Trans. Magn. 36 (5) (20 0 0) 3141–3143, doi:10.1109/20.908713.

Please cite this article as: M. Novak et al., Difficulty in identification of Preisach hysteresis model weighting function using first order reversal curves method in soft magnetic materials, Applied Mathematics and Computation (2017), http://dx.doi.org/10.1016/j.amc.2017.05.017