Determination of a Hysteresis Model Parameters with the Use of Different Evolutionary Methods for an Innovative Hysteresis Model

Abstract

For precise modeling of electromagnetic devices, we have to model material hysteresis. A Genetic Algorithm, Differential Evolution with three different strategies, teaching–learning-based optimization and Artificial Bee Colony, were used for testing seven different modified mathematical expressions, and the best combination of mathematical expression and solving method was used for hysteresis modeling. The parameters of the hysteresis model were determined based on the measured major hysteresis loop and first-order reversal curves. The model offers a simple determination of the magnetization procedure in the areas between measured curves, with the only correction of two parameters based on only two known points in the magnetization process. It was tested on two very different magnetic materials, and results show good agreement between the measured and calculated curves. The calculated curves between the measured curves have correct shapes. The main difference between our model and other models is that, in our model, each measured curve, major and reversal, is described with different parameters. The magnetization process between measured curves is described according to the nearest measured curve, and this ensures the best fit for each measured curve. In other models, there is mainly only one curve, a major hysteresis or magnetization curve, used for the determination of the parameters, and all other curves are then dependent on this curve. Results confirm that the evolutionary optimization method offers a reliable procedure for precise determination of the parameters.

1. Introduction

Precise modeling of electromagnetic devices helps us to develop energy-efficient devices. Modeling of modern electromagnetic devices is done mainly with the use of numerical models, such as the finite element model. To get a precise model that allows accurate modeling of electromagnetic devices, the description of the materials’ properties should be as good as possible. Often nonlinear magnetic material is described using a magnetization curve. To make a better material description and consider hysteresis losses, material hysteresis also has to be modeled. In the past, different hysteresis models were presented by different authors.

To understand the meaning of the hysteresis model better, the magnetization of the motor cross-section using hysteresis is presented schematically in Figure 1. The value of the excitation current (Figure 1a) changes from the value presented with point 1 to values in points 2, 3, … to point 8. Magnetization changes according to the curves presented between points marked at the hysteresis loop (Figure 1b) and appropriate magnetic flux density B at a cross-section of the motor for magnetization at each presented point (Figure 1c).

The magnetization process starts with demagnetized material at point 1. Magnetization from zero to the maximum current, point 2, is made using the original magnetization curve. From point 2 to point 4, magnetization is made according to the left part of the major hysteresis loop. From point 4 to point 5, magnetization is made according to the right part of the major hysteresis loop. For smaller current at points 6, 7 and 8, the magnetization is made according to the first-order reversal magnetization curves. The current can have any value, so the magnetization can be anywhere inside the hysteresis loop. Some of the first-order curves inside the hysteresis loop are measured, but most of them should be obtained using a hysteresis model.

Mathematical models can be classified into two main categories: physical and phenomenological. The former (physical model) is created from the physics of a phenomenon, while the latter merely describes the end result with mathematical functions. Therefore, the latter is usually less complex, but it loses the connection to the physics of the material. Among the physical models is the Stoner–Wohlfarth model [1], which is often used to describe events in ferromagnetic thin layers. The material is treated as a single magnetic domain. Jiles-Atherton [2,3,4], on the other hand, also considers inter-domain interactions. Their main advantage is the connection with the physical parameters. More numerically oriented is the Preisach model of hysteresis [5,6,7]. It enables the description of different forms of hysteresis loops that are constructed from the so-called hysterons. They can switch between −1 and 1, depending on whether the input data exceed the start-up or shut-off threshold. Other less known hysteresis models can be found, such as the Mayergoyz Vector Hysteresis model [8], or Vector Magneto–Hysteretic H&S model [9], and many modifications and corrections of existing models. The phenomenological, numerically oriented approach is presented in the paper.

The history of the magnetization process must also be known in the case of magnetic field calculation considering hysteresis. For the Finite Element Method (FEM) calculation [10,11,12,13,14,15,16,17,18,19,20], the magnetization process could be different for almost every element. Due to this fact, it is favorable that the description of the magnetization history is as simple as possible. It is expected that the model is suitable for a wide range of different materials. It is a challenging problem to determine the parameters of the hysteresis models.

In our work, a phenomenological scalar hysteresis model is made, which offers a simple description of the magnetization process. Different mathematical expressions are tested to get the most appropriate.

This test is possible only if a reliable tool to determine the parameters is available. Different evolutionary optimization methods were used for the determination of the parameters, with the aim to determine the best for such a problem. The improper determination of model parameters due to the improper selection of the optimization method can lead to an incorrect evaluation of model quality. The used methods were genetic algorithm (GA) [21,22,23,24,25,26], differential evolution (DE) with three different strategies, which are DE/rand/1/exp, DE/rand/2/exp, DE/best/1/bin [27,28,29,30,31,32,33,34,35,36], teaching–learning-based optimization (TLBO) [37,38,39,40,41,42,43,44,45] and artificial bee colony (ABC) [46,47,48,49,50,51,52,53,54,55]. If an inappropriate solving method is used, the adequacy of the mathematical expression will not be determined correctly.

Two measured hystereses of materials 9S20 and Alnico with first-order reversal curves were used for the tests, and both are presented in Figure 2. The materials 9S20 and Alnico are magnetic, with very different hysteresis curves.

For the material 9S20, the parts of major hysteresis are marked with the left/right part of the measured major hysteresis loop (SML/SMR), decreasing first-order reversal curves are marked as SD1–SD4, and increasing first-order reversal curves are marked with SI1–SI4. For the material Alnico, the parts of major hysteresis are marked with the left/right part of the Alnico measured major hysteresis loop (AML/ AMR), decreasing curves are marked as AD1–AD3, and increasing curves are marked with AI1–AI3. The measured major hysteresis curve and first-order reversal curves were used to determine the parameters of the model. After that, only two points, the point at which the excitation is changed (for example, C–F in Figure 2) and the point at the ends of the hysteresis (for example, A and B in Figure 2), were used to determine the magnetization process in the hysteresis area between the measured curves.

No comparison between different mathematical expressions and their use on cases of very different hysteresis have been found in the literature, to know which of them is the most appropriate for the hysteresis modeling. Our contributions in this work are:

• Modification of the mathematical expressions from the literature, such as expressions from [56], by adding additional terms, with the aim to get a more suitable expression.
• Determination of a mathematical expression, selected between seven modified expressions, and determination of an appropriate method for the calculation of the parameters, selected between six methods (four methods and three different strategies of DE), to get the combination of expression and solving method which is suitable for a wide range of different materials. The appropriate method for the determination of the parameters is suggested, and tested by using it for the determination of parameters of a major hysteresis loop and more reversal curves inside the hysteresis loop for two very different materials.
• An innovative description of the magnetization process inside the hysteresis areas between the measured curves is presented, which offers a fast determination of the magnetization process and, with that, is suitable for use in the scope of the FEM calculation. It is impossible to make so many measurements that the calculation of curves in the area between the measurements would not be necessary.

Other authors have tried to make models appropriate for FEM calculation. In [7] M. Kuczmann developed an inverse vector Preisach hysteresis model, and it was inserted in the finite element procedure. Y. Zhang et al. in [57] presented a novel hysteresis core loss model for magnetic laminations, which is used in the scope of the finite element model. Like us, L. Petreswcu et al. in [56] also used Sigmoid’s functions in hysteresis phenomena modeling. The main difference between our model and other models is that, in our model, each measured curve, major and reversal, is described with different parameters, obtained using the optimization method, and the magnetization process between measured curves is described according to the nearest measured curve. In other models there is mainly only one curve, a major hysteresis or magnetization curve, used for the determination of the parameters, and all other curves are then dependent on this curve.

The paper consists of six Sections. In Section 2, different mathematical expressions are tested using standard optimization methods and evolutionary methods, and the best combination is presented of expression and method for determination of the parameters. The use of the model for the material 9S20 is presented in Section 3, and for the material Alnico in Section 4. The limits of the parameters are presented in Section 5 and the use of the model is presented in Section 6. In the last, Section 7, conclusions are given considering the proposed model.

2. Tested Mathematical Expressions and Solving Methods

Many engineering problems were solved using optimization methods. Sometimes the number of function evaluations is severely limited by time or cost. In [58] Jones et al. presented a response surface methodology to address this challenge. These surfaces can then be used for visualization, tradeoff analysis, and optimization. We were dealing with a problem presented in [59] by Jesenik et al., where the Newton–Raphson method was used to solve it. It is combined with approximation planes for a two-dimensional problem, obtained with the time consuming Finite Element Model. In the presented work, we are dealing with a curve-fitting problem, and function evaluations are not costly. Therefore, approaches for decreasing the number of fitness evaluations are not used, such as the approach with response surface methodology.

2.1. Tested Mathematical Expressions

Hysteresis occurs in various fields of engineering work [60,61,62]. In the past, different mathematical expressions were used for hysteresis modeling. The magnetization curve, and also hysteresis, can be modeled in the easiest way by using a Sigmoid function. In the literature, different expressions of the Sigmoid function can be found [56], and it plays a very important role in hysteresis modeling. A Sigmoid function is a mathematical function, which usually has an “S” shape. It is real-valued and differentiable. It is used in many fields, such as in the field of Neural Networks, to describe the processes in the human brain, in mathematical psychology, etc. We selected seven mathematical expressions of Sigmoid function from [56]. The Langevin function is related to the Boltzmann law of magnetic moment distributions. Other functions are mathematical expressions in which parameters have less physical meaning. These are used to adjust mathematical expressions to hysteresis curves.

Seven mathematical expressions from [56] are modified appropriately by adding additional terms into the expressions, and tested on the basis of two very different hystereses of a material named 9S20, with coercive magnetic field strength 4096 A/m, and a material named Alnico, with coercive magnetic field strength 47,289 A/m, and both materials are presented in Figure 2. Measurements were made in our laboratory (the Laboratory of Applied Electromagnetics at the Faculty of Electrical Engineering and Computer Science, University of Maribor). Based on our experience in working with magnetic materials, we added the following parts in each used expression:

-

After each H, a parameter was added in the form (H + P_a), which offers the possibility to move a curve with respect to the H axis.

-

Parameter + P_b was added, which offers the possibility to move a curve with respect to the B axis.

-

Part P_c∙μ₀∙(H + P_a) was added for better adjustment of saturation.

The used modified expressions are:

Langevin:

$$B(H)=P_1\left(\coth \left(\frac{H+P_3}{P_2}\right)-\frac{P_2}{H+P_3}\right)+P_4{\mu }_0\left(H+P_3\right)+P_5.$$

(1)

Grompertz:

$$B(H)=2P_1{e}^{-\exp \left(-P_2\left(H+P_3\right)\right)}-P_1+P_4{\mu }_0\left(H+P_3\right)+P_5.$$

(2)

Hyptan (Hyperbolic tangents):

$$B(H)=P_1\tanh \left(P_2\left(H+P_3\right)\right)+P_4{\mu }_0\left(H+P_3\right)+P_5.$$

(3)

Algebraic:

$$B(H)=\frac{P_1\left(H+P_3\right)}{\sqrt{P_4{\left(P_2\left(H+P_3\right)\right)}^2}}+P_5{\mu }_0\left(H+P_3\right)+P_6.$$

(4)

Logistic:

$$B(H)=\frac{2P_1}{P_4+e^{-P_2\left(H+P_3\right)}}-P_1+P_5{\mu }_0\left(H+P_3\right)+P_6.$$

(5)

Sigmoid:

$$B(H)=\frac{P_1{e}^{P_2\left(H+P_3\right)}}{P_4+e^{P_2\left(H+P_3\right)}}+P_5{\mu }_0\left(H+P_3\right)+P_6.$$

(6)

Elliot:

$$B(H)=\frac{P_1\left(H+P_3\right)}{P_4+P_2\left|H+P_3\right|}+P_5{\mu }_0\left(H+P_3\right)+P_6.$$

(7)

2.2. Calculation of the Parameters Using Standard Optimization Methods

Standard optimization methods can be used for the determination of the parameters. We tested the Levenberg-Marquardt algorithm and nonlinear regression to determine their usefulness for solving the described problem.For the test of the standard optimization methods, we used Elliot’s expression (7) and SML of the material 9S20, presented in Figure 2.

The Objective Function (OF) is written in (8).

$$\mathrm{OF}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(B_{\mathrm{calculated}\_i}-B_{\mathrm{measured}\_i}\right)}^2},$$

(8)

where n is the number of measured points, B_measured are the measured points, and B_calculated are the calculated points.

An initial guess is needed to start standard optimization methods. The Levenberg-Marquardt (LM) algorithm is often used for solving curve-fitting problems, which is also evident from the literature [63,64]. A set of nonlinear equations was solved in the least-square sense. The minimization of a sum of squares of equation residuals was made using the LM algorithm.

Unfortunately, the local minima caused that different initial values led to different results. The lack of change of the parameters was the stopping condition of the algorithm (the maximum change must be less than 10⁻¹²). Figure 3 shows the tests for the SML curve, and the values of the parameters which were used for both the initial and final curves are presented in

Table 1. The OF and parameters’ values calculated with LM and used for both the initial and final analytical curves presented in Figure 3.

	OF	Parameters
(a) Initial curve	8.02·10−1	P1 = 7; P2 = 20; P3 = 0; P4 = 2000; P5 = 15; P6 = 0
(a) Final curve	6.34·10−4	P1 = 6.89; P2 = 4.47; P3 = 428.8; P4 = 1999.6; P5 = 21.4; P6 = 1.60·10−3
(b) Initial curve	3.00·10−1	P1 = 19; P2 = 12; P3 = −400; P4 = 6000; P5 = 21; P6 = 0
(b) Final curve	1.55·10−1	P1 = 18.13; P2 = 12.28; P3 = −256.6; P4 = 6000; P5 = 24.09; P6 = 1.36·10−1
(c) Initial curve	4.61·10−1	P1 = 10; P2 = 12; P3 = −400; P4 = 6000; P5 = 21; P6 = 0
(c) Final curve	1.22·10−1	P1 = 11.58; P2 = 9.07; P3 = −190.8; P4 = 6000; P5 = 65.41; P6 = 1.36·10−1
(d) Initial curve	5.40·10−1	P1 = 20; P2 = 12; P3 = −400; P4 = 2000; P5 = 21; P6 = 0
(d) Final curve	6.34·10−4	P1 = 7.60; P2 = 4.93; P3 = 428.9; P4 = 2204.5; P5 = 21.24; P6 = 1.60·10−3

.

The results presented in Figure 3b,c are not correct (OF is 0.1548 and 0.1215). The OF for cases (a) and (d) is 6.342·10⁻⁴. Based on the results presented as (b) and (c) in Table 1 and Figure 3, it can be seen that, for two initial curves, the obtained results are in two different local minima (the final curve is not similar to the measured values), and they are not the global minimum. The results presented in (a) and (d) are global minima.

LM is very fast. A very good result can be obtained by using appropriate initial values. Due to the fact that LM is not able to bypass the local minima, it is not suitable to be used as a stand-alone method. Also, NonLinear regression using the Matlab function “nlinfit” and using analytical derivatives were not suitable for the presented problem.

2.3. Calculation of the Parameters Using Evolutionary Methods

Evolutionary optimization methods [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55] are population-based search algorithms, where the space of all possible solutions is searched through many solutions with various suitability (fitness). Better solutions have a higher probability to survive into the next generation, where they might be altered by mutation and/or crossover. Evolutionary optimization methods have a remarkable ability to balance exploration and exploitation [25] and, as such, avoid local optima. They are often used for technical problems, and also to determine the parameters of different models [65] describing real material properties or real devices. In this work, we used and compared GA [21,22,23,24,25,26], DE [27,28,29,30,31,32,33,34,35,36], TLBO [37,38,39,40,41,42,43,44,45] and ABC [46,47,48,49,50,51,52,53,54,55]. In GA 0.5 was used for the fraction of the kept population, and 0.2 was used for the mutation rate. For all three used DE strategies, which were DE/rand/1/exp (Algorithm 1), DE/rand/2/exp and DE/best/1/bin, 0.8 was used for the crossover probability, and 0.6 was used for the amplification of the differential variation. Different strategies are used, because it is known that the appropriateness of the strategy depends on the problem being addressed. The duplication elimination phase was omitted in the TLBO used in our work. Therefore, the number of Fitness Evaluations (FEs) which were consumed for the calculation was determined statically as FEs = 2 * population size * iterations. The limit value used in ABC was set at 100, and was used as a control parameter for the bee population. A scout bee might not be used for every iteration [49], therefore the number of FE’s used in ABC cannot be determined statically.

Algorithm 1. Pseudo-code for DE/rand/1/exp.

1.  Setting control parameters (maximum number of FEs, population size, F, Cr)

2.  Randomly generates the initial population of vectors in the D-dimensional search space

3.  repeat

4.  for each vector i in the population do

5.    select three mutually exclusive random vectors ${\overrightarrow{x}}_1\ne {\overrightarrow{x}}_2\ne {\overrightarrow{x}}_3$

6.    generate a donor vector according to (9)

7. $${\overrightarrow{v}}_i={\overrightarrow{x}}_1+F\cdot \left({\overrightarrow{x}}_{-1}{\overrightarrow{x}}_3\right)$$ (9)

8.    $L=0$

9.    do

10.    $L=L+1$

11.    while (rand(0,1) <= Cr AND L <= D)

12.    N = a random integer in the range of [1,D]

13.    generate a trial vector ${\overrightarrow{u}}_i$ using exponential crossover by (10)

14. $$u_{ij}=\{\begin{array}{ll}{v}_{ij}\mathrm{for}\mathrm{each}\mathrm{dimension}j=& \left(N\mathrm{mod}D\right),\left(N+1\mathrm{mod}D\right),\\ & \dots ,\left(N+L-1\mathrm{mod}D\right)\\ x_{ij}\mathrm{for}\mathrm{all}\mathrm{other}\mathrm{dimensions}& j\in \left[1,D\right]\end{array}$$ (10)

15.    evaluate the candidate vector ${\overrightarrow{u}}_i$ using (8)

16.    replace ${\overrightarrow{x}}_i$ with ${\overrightarrow{u}}_i$, if fitness of ${\overrightarrow{u}}_i$ is better than fitness of ${\overrightarrow{x}}_i$

17.   end-for

18.  until a termination condition is met

For all tests we used 83,160 FEs, which were appropriate in [66] for nine parameters, and in this work we are dealing with five (D = 5) or six parameters (D = 6). The population size is ten times the number of parameters, fifty for five parameters and sixty for six parameters. The OF which was used for the analytical curves’ quality estimation was defined previously in (8).

To compare the ability of the presented expressions (1)–(7) to express hysteresis analytically, they were tested at the left part of the 9S20 measured major hysteresis loop (SML) and AML, which are presented in Figure 2. The number of independent runs was 30 for each expression for GA, three DE strategies, (DE/rand/1/exp, DE/rand/2/exp and DE/best/1/bin), TLBO and ABC. The best (B), worst (W), mean (M) OF and standard deviation (SD) are written in

Table 2. Best (B), worst (W), mean (M) and standard deviation (SD) for 30 independent runs of different methods using expressions (1) to (7) for SML.

Expression				Method
		GA	DE/Rand/	DE/Rand/	DE/Best/	TLBO	ABC
			1/Exp	2/Exp	1/Bin
	B	1.11·10−3	9.75·10−4	9.75·10−4	9.75·10−4	9.75·10−4	9.79·10−4
LA 1	W	3.64·10−1	9.75·10−4	9.75·10−4	4.76·10−1	9.75·10−4	4.55·10−3
	M	7.41·10−2	9.75·10−4	9.75·10−4	6.49·10−2	9.75·10−4	1.76·10−3
	SD	1.25·10−1	3.25·10−19	3.25·10−19	1.61·10−1	3.25·10−19	8.61·10−4
	B	7.55·10−2	3.88·10−3	3.88·10−3	3.88·10−3	2.96·10−2	2.97·10−2
GR 1	W	1.23	5.99·10−3	2.79·10−2	4.76·10−1	2.99·10−2	1.57·10−1
	M	1.63·10−1	3.95·10−3	7.42·10−3	1.99·10−1	2.96·10−2	3.73·10−2
	SD	2.76·10−1	3.78·10−4	7.34·10−3	2.17·10−1	7.73·10−5	2.29·10−2
	B	4.21·10−3	1.77·10−3	1.77·10−3	1.77·10−3	1.77·10−3	1.77·10−3
HY 1	W	1.37	4.73·10−1	1.42·10−2	4.76·10−1	4.76·10−1	3.30·10−2
	M	5.67·10−1	5.19·10−2	1.13·10−2	3.68·10−1	1.60·10−1	5.47·10−3
	SD	5.93·10−1	1.34·10−1	5.27·10−3	1.96·10−1	2.08·10−1	6.98·10−3
	B	1.34·10−3	1.13·10−3	1.14·10−3	1.14·10−3	1.14·10−3	1.14·10−3
AL 1	W	1.70·10−1	1.1368·10−3	1.14·10−3	1.70·10−1	3.61	4.30·10−3
	M	5.41·10−2	1.13·10−3	1.14·10−3	1.24·10−2	3.45·10−1	1.82·10−3
	SD	7.44·10−2	8.67·10−19	4.90·10−12	4.20·10−2	1.11	7.75·10−4
	B	1.23·10−1	1.42·10−2	6.01·10−2	2.96·10−2	3.05·10−2	2.26·10−2
LO 1	W	2.01·10−1	6.01·10−2	6.01·10−2	4.76·10−1	1.70·10−1	1.50·10−1
	M	1.64·10−1	5.85·10−2	6.01·10−2	3.03·10−1	1.41·10−1	7.25·10−2
	SD	1.51·10−2	8.22·10−3	2.78·10−17	2.13·10−1	4.70·10−2	3.62·10−2
	B	1.97·10−1	4.83·10−2	7.49·10−2	1.78·10−2	2.21·10−3	1.77·10−3
SI 1	W	2.02·10−1	6.13·10−2	4.77·10−1	3.42·10−2	1.47	2.03·10−3
	M	2.00·10−1	5.49·10−2	2.11·10−1	2.35·10−2	2.36·10−1	1.83·10−3
	SD	1.04·10−3	5.32·10−3	1.88·10−1	7.58·10−3	3.06·10−1	5.99·10−5
	B	2.23·10−2	6.34·10−4	6.34·10−4	6.34·10−4	6.34·10−4	6.42·10−4
EL 1	W	4.95·10−1	6.34·10−4	6.34·10−4	1.83·10−1	6.34·10−4	3.89·10−3
	M	3.71·10−1	6.34·10−4	6.34·10−4	6.32·10−2	6.34·10−4	1.50·10−3
	SD	1.98·10−1	4.34·10−19	1.80·10−13	7.68·10−2	3.62·10−11	8.57·10−4

¹ LA-Langevin, GR-Grompertz, HY-Hyptan, AL-Algebraic, LO-Logistic, SI-Sigmoid, EL-Elliot.

for SML, and in

Table 3. B, W, M and SD for 30 independent runs of different methods using expressions (1) to (7) for left part of the Alnico measured major hysteresis loop (AML).

Expression				Method
		GA	DE/Rand/	DE/Rand/	DE/Best/	TLBO	ABC
			1/Exp	2/Exp	1/Bin
	B	1.82·10−2	4.36·10−4	4.36·10−4	4.36·10−4	4.36·10−4	4.36·10−4
LA 1	W	6.61·10−1	4.36·10−4	4.36·10−4	5.61·10−1	3.55·10−1	4.37·10−4
	M	2.86·10−1	4.36·10−4	4.36·10−4	2.04·10−1	2.41·10−2	4.36·10−4
	SD	1.66·10−1	2.71·10−19	2.71·10−19	2.55·10−1	8.84·10−2	2.63·10−7
	B	4.20·10−2	1.83·10−3	1.83·10−3	4.07·10−2	4.06·10−2	4.17·10−2
GR 1	W	6.56·10−1	5.61·10−1	3.56·10−1	5.61·10−1	5.61·10−1	4.86·10−2
	M	3.46·10−1	1.11·10−1	2.01·10−1	4.12·10−1	1.92·10−1	4.33·10−2
	SD	1.61·10−1	9.59·10−2	1.75·10−1	1.34·10−1	1.78·10−1	1.57·10−3
	B	5.72·10−2	1.93·10−3	1.93·10−3	1.93·10−3	1.93·10−3	1.93·10−3
HY 1	W	5.65·10−1	1.16·10−1	1.93·10−3	5.61·10−1	5.61·10−1	4.93·10−2
	M	3.78·10−1	5.74·10−3	1.93·10−3	2.56·10−1	9.47·10−2	4.54·10−3
	SD	9.97·10−2	2.05·10−2	8.67·10−19	2.05·10−1	1.33·10−1	8.56·10−3
	B	1.33·10−2	8.82·10−4	8.82·10−4	8.82·10−4	8.82·10−4	8.82·10−4
AL 1	W	1.48·10−2	8.82·10−4	8.82·10−4	5.61·10−1	8.82·10−4	9.07·10−4
	M	1.37·10−2	8.82·10−4	8.82·10−4	1.13·10−1	8.82·10−4	8.88·10−4
	SD	3.82·10−4	4.34·10−19	8.22·10−12	2.24·10−1	3.16·10−11	5.69·10−6
	B	4.17·10−2	2.04·10−1	2.04·10−1	4.06·10−2	4.23·10−2	1.94·10−3
LO 1	W	5.26·10−1	4.56·10−1	2.04·10−1	5.61·10−1	5.61·10−1	2.17·10−3
	M	3.25·10−1	2.12·10−1	2.04·10−1	3.17·10−1	2.16·10−1	2.02·10−3
	SD	1.39·10−1	4.53·10−2	8.33·10−17	1.85·10−1	1.45·10−1	7.30·10−5
	B	3.17·10−1	6.713·10−3	1.98·10−2	1.93·10−3	1.96·10−3	1.15·10−2
SI 1	W	5.43·10−1	2.68·10−1	6.24·10−2	1.96·10−3	1.80·10−2	1.74·10−1
	M	5.16·10−1	3.98·10−2	3.41·10−2	1.93·10−3	2.83·10−3	4.12·10−2
	SD	4.88·10−2	6.94·10−2	1.24·10−2	7.32·10−6	2.84·10−3	3.17·10−2
	B	5.86·10−2	4.85·10−4	4.85·10−4	4.85·10−4	4.85·10−4	4.88·10−4
EL 1	W	2.16·10−1	4.85·10−4	4.85·10−4	5.61·10−1	7.78·10−2	1.24·10−3
	M	1.57·10−1	4.85·10−4	4.85·10−4	6.78·10−2	3.90·10−3	6.41·10−4
	SD	3.65·10−2	1.63·10−19	1.63·10−19	1.42·10−1	1.45·10−2	1.53·10−4

¹ LA-Langevin, GR-Grompertz, HY-Hyptan, AL-Algebraic, LO-Logistic, SI-Sigmoid, EL-Elliot.

for AML. In Table 2 and Table 3 the lowest best OF and the lowest mean OF are marked with bold, and the second-lowest best OF and the second-lowest mean OF are marked with italics.

For each combination of expression and method from Table 2, appropriate parameters obtained using results with the best OF are shown in

Table 4. Parameters of expressions for the results (best OF) from Table 2 for SML.

Expression				Method
		GA	DE/Rand/	DE/Rand/	DE/Best/	TLBO	ABC
			1/Exp	2/Exp	1/Bin
	P1	1.31	1.26	1.26	1.26	1.26	1.26
LA 1	P2	192.7	180.3	180.3	180.3	180.3	179.0
	P3	424.8	436.8	436.8	436.8	436.8	433.9
	P4	49.46	61.21	61.21	61.21	61.21	61.42
	P5	3.61·10−3	−3.55·10−3	−3.55·10−3	−3.55·10−3	−3.55·10−3	−1.54·10−3
	P1	1.05	1.03	1.03	1.03	7.09·10−1	7.14·10−1
GR 1	P2	13.82	2.92·10−3	2.92·10−3	2.92·10−3	13.35	5.97
	P3	5.34·103	1.70	1.70	1.70	5.31·103	2.37·103
	P4	49.49	99.19	99.19	99.19	186.4	186.3
	P5	3.94·10−1	6.15·10−2	6.15·10−2	6.15·10−2	1.33	6.48·10−2
	P1	1.25	1.05	1.05	1.05	1.05	1.05
HY 1	P2	1.40·10−3	1.97·10−3	1.97·10−3	1.97·10−3	1.97·10−3	1.97·10−3
	P3	428.6	441.6	441.6	441.6	441.6	441.4
	P4	36.70	93.84	93.84	93.84	93.84	95.05
	P5	1.81·10−2	−7.25·10−3	−7.25·10−3	−7.25·10−3	−7.25·10−3	−7.18·10−3
	P1	18.56	22.16	22.78	22.81	22.65	22.71
AL 1	P2	16.04	−19.42	−19.96	19.99	19.84	20.00
	P3	415.6	437.9	437.9	437.9	437.9	437.2
	P4	6922·104	9633·104	1018·105	1021·105	1006·104	1004·104
	P5	70.91	−4.45·10−3	−4.45·10−3	−4.45·10−3	75.49	77.31
	P6	5.40·10−3	−4.45·10−3	−4.45·10−3	−4.45·10−3	−4.45·10−3	−3.46·10−3
	P1	6.01	2.75	2.52	2.98	2.31	3.00
LO 1	P2	−2.48·10−1	7.35·10−3	−20.00	5.46	−7.71	9.34·10−3
	P3	−689.5	−258.5	−433.1	−408.7	−405.6	−229.5
	P4	−41.91	3.83	−5.00	4.21	−3.31	4.85
	P5	351.1	2.05	3.00	2.27	189.9	223.1
	P6	5.97	2.05	3.00	2.27	3.00	2.41
	P1	8.03·10−1	2.05	2.63	1.70	2.03	2.12
SI 1	P2	1.74·10−1	2.75·10−2	3.08·10−2	1.34·10−2	4.42·10−3	3.90·10−3
	P3	−76.53	631.0	215.8	423.7	727.9	609.8
	P4	2.31	7.26·10−2	4.55·10−2	3.83	3.86	1.92
	P5	199.9	−1.18	−1.21	−8.24·10−1	97.75	92.15
	P6	7.45·10−3	−2.61·10−1	−4.05·10−1	1.06	−1.05	−1.08
	P1	8.23·10−1	19.63	19.63	19.63	19.63	13.87
EL 1	P2	1.05	12.73	12.73	12.73	12.73	9.13
	P3	385.3	429.1	429.1	429.1	429.1	430.3
	P4	15.08	5705	5705	5705	5705	3953
	P5	165.8	21.13	21.13	21.13	21.13	24.36
	P6	2.80·10−2	1.59·10−3	1.59·10−3	1.59·10−3	1.59·10−3	−3.79·10−4

¹ LA-Langevin, GR-Grompertz, HY-Hyptan, AL-Algebraic, LO-Logistic, SI-Sigmoid, EL-Elliot.

. For each presented combination in Table 3, appropriate parameters obtained using results with the best OF are shown in

Table 5. Parameters of expressions for the results (best OF) from Table 3 for AML.

Expression				Method
		GA	DE/Rand/	DE/Rand/	DE/Best/	TLBO	ABC
			1/Exp	2/Exp	1/Bin
	P1	1.42	1.55	1.55	1.55	1.55	1.55
LA 1	P2	1138	1501	1501	1501	1501	1501
	P3	47537	48339	48339	48339	48339	48339
	P4	2.75	2.65	2.65	2.6	2.65	2.65
	P5	1.63·10−3	−1.76·10−1	−1.76·10−1	−1.76·10−1	−1.76·10−1	−1.76·10−1
	P1	1.19	1.43	1.43	1.24	1.24	1.22
GR 1	P2	1.10	2.90·10−4	2.90·10−4	1.22	1.24	1.20
	P3	53233	14.39	14.39	59112	60000	58138
	P4	6.01	3.29	3.29	5.20	5.25	5.74
	P5	5.54·10−1	4.80·10−2	4.80·10−2	5.18·10−1	5.26·10−1	5.47·10−1
	P1	1.22	1.42	1.42	1.42	1.42	1.42
HY 1	P2	1.12·10−2	2.05·10−4	4.05·10−4	2.05·10−4	2.05·10−4	2.05·10−4
	P3	47990	48282	48282	48282	48282	48282
	P4	4.07	3.42	3.42	3.42	3.42	3.42
	P5	2.26·10−3	−1.68·10−1	−1.68·10−1	−1.68·10−1	−1.68·10−1	−1.68·10−1
	P1	24.30	27.43	28.48	28.50	27.48	21.83
AL 1	P2	−16.18	−18.69	19.41	19.43	−18.73	14.88
	P3	47439	48311	48311	48311	48311	48313
	P4	6192·10 6	7155·10 6	7716·10 6	7728·10 6	7183·10 6	4541·10 6
	P5	1.51	−1.72·10−1	−1.72·10−1	−1.72·10−1	−1.72·10−1	3.11
	P6	5.22·10−5	−1.72·10−1	−1.72·10−1	−1.72·10−1	−1.72·10−1	−1.72·10−1
	P1	5.43	2.77	2.77	5.48·10−1	1.95	4.58
LO 1	P2	1.23·10−1	−20	−20	19.99	−20	4.16·10−4
	P3	−48611	−49233	−49233	−48500	−48485	−45473
	P4	4.30	−5.00	−5.00	4.43·10−3	−1.58	3.23
	P5	4.89	3.00	3.00	−8.81·10−1	5.26	3.49
	P6	3.99	3.00	3.00	−8.81·10−1	3.00	3.00
	P1	9.36	3.01	3.00	2.85	2.86	2.76
SI 1	P2	−6.82·10−5	2.49·10−4	1.65·10−4	4.10·10−4	3.99·10−4	2.50·10−4
	P3	999.8	225.7	400.0	368.2	400.0	625.6
	P4	−10.00	6.20·10−6	5.61·10−4	0.0	0.0	6.80·10−6
	P5	10.26	−1.51	−1.41	−1.39	3.34	4.20
	P6	1.09	−1.52	−1.17	−1.39	−1.49	−1.30
	P1	13.38	20.83	20.83	20.83	20.83	19.82
EL 1	P2	10.59	12.64	12.64	12.64	12.64	11.99
	P3	47423	48378	48378	48378	48378	48382
	P4	11.38	35799	35799	35799	35799	34509
	P5	3.53	2.11	2.11	2.11	2.11	2.08
	P6	2.63·10−3	−1.82·10−1	−1.82·10−1	−1.82·10−1	−1.82·10−1	−1.82·10−1

¹ LA-Langevin, GR-Grompertz, HY-Hyptan, AL-Algebraic, LO-Logistic, SI-Sigmoid, EL-Elliot.

.

The aim of the results presented in Table 2, Table 3, Table 4 and Table 5 is also to find out if different algorithms lead to different parameters. For Langevin, Grompertz, Hyptan and Elliot, it can be seen that, with methods that reached the same best OF, the calculated parameters are also the same. If the best OF is only a little different, the parameters are also just slightly different. In the case of algebraic, the parameters are different, although the best OF is the same. In the case of logistic and Sigmoid, the obtained best OF is different for almost every algorithm, and, consequently, the calculated parameters are also different. Although the same best OF is obtained using different methods, mean OF can be very different, which means that some methods occasionally give poor results on this problem, which is not desirable.

A comparison of different expressions is presented graphically in Figure 4 for SML and in Figure 5 for AML. Parameters used were obtained from the calculation with the best OF.

From Table 2 it can be seen that the best results were obtained using the Elliot expression and methods DE/rand/1/exp, DE/rand/2/exp and TLBO. The second-best result was obtained using Langevin and methods DE/rand/1/exp, DE/rand/2/exp and TLBO. On the other hand, based on Table 3, it can be observed that the best results were obtained using Langevin and methods DE/rand/1/exp and DE/rand/2/exp. The second-best result was obtained using Elliot and methods DE/rand/1/exp and DE/rand/2/exp. We can conclude that the best results were obtained using the Elliot and Langevin expressions and DE/rand/1/exp and DE/rand/2/exp solving methods. To select the best, additional tests were made based on all curves of the material 9S20 (SML, SMR, SD1–SD4, SI1–SI4) and Alnico (AML, AMR, AD1–AD3, AI1–AI3), presented in Figure 2.

In

Table 6. B, W, M and SD for 30 independent runs for all 9S20 curves using DE/rand/1/exp and DE/rand/2/exp for the Elliot and Langevin expressions.

Curve			Method
		DE/Rand/1/Exp	DE/Rand/2/Exp	DE/Rand/1/Exp	DE/Rand/2/Exp
		Elliot	Elliot	Langevin	Langevin
	B	6.34·10−4	6.34·10−4	9.75·10−4	9.75·10−4
SML	W	6.34·10−4	6.34·10−4	9.75·10−4	9.75·10−4
	M	6.34·10−4	6.34·10−4	9.75·10−4	9.75·10−4
	SD	4.34·10−19	1.80·10−13	3.25·10−19	3.25·10−19
	B	7.19·10−4	7.19·10−4	1.11·10−3	1.11·10−3
SMR	W	7.19·10−4	7.19·10−4	1.11·10−3	1.11·10−3
	M	7.19·10−4	7.19·10−4	1.11·10−3	1.11·10−3
	SD	3.25·10−19	3.25·10−19	4.34·10−19	4.34·10−19
	B	6.51·10−5	6.51·10−5	1.41·10−4	1.47·10−4
SD1	W	7.59·10−4	8.08·10−3	6.50·10−3	6.60·10−3
	M	4.10·10−4	2.41·10−3	5.39·10−4	7.21·10−4
	SD	3.23·10−4	2.18·10−3	1.47·10−3	1.49·10−3
	B	3.25·10−4	3.25·10−4	6.12·10−4	6.12·10−4
SD2	W	8.19·10−3	1.17·10−2	6.12·10−4	6.27·10−4
	M	2.69·10−3	7.27·10−3	6.12·10−4	6.12·10−4
	SD	3.60·10−3	3.16·10−3	2.17·10−19	2.71·10−6
	B	3.80·10−4	3.80·10−4	6.02·10−4	6.02·10−4
SD3	W	3.80·10−4	3.97·10−2	6.02·10−4	6.02·10−4
	M	3.80·10−4	5.55·10−3	6.02·10−4	6.02·10−4
	SD	1.63·10−19	1.32·10−2	3.25·10−19	3.25·10−19
	B	7.44·10−4	7.44·10−4	5.67·10−4	5.67·10−4
SD4	W	7.44·10−4	7.44·10−4	5.67·10−4	5.67·10−4
	M	7.44·10−4	7.44·10−4	5.67·10−4	5.67·10−4
	SD	2.17·10−19	2.49·10−13	2.17·10−19	2.17·10−19
	B	6.51·10−5	6.51·10−5	1.47·10−4	1.47·10−4
SI1	W	7.09·10−4	7.12·10−3	1.13·10−2	1.04·10−2
	M	3.66·10−4	3.29·10−3	8.93·10−4	1.14·10−3
	SD	3.21·10−4	2.45·10−3	2.79·10−3	2.70·10−3
	B	3.25·10−4	3.25·10−4	6.12·10−4	6.12·10−4
SI2	W	8.19·10−3	9.00·10−3	6.12·10−4	6.12·10−4
	M	3.73·10−3	6.79·10−3	6.12·10−4	6.12·10−4
	SD	3.90·10−3	3.13·10−3	2.17·10−19	2.17·10−19
	B	3.80·10−4	3.80·10−4	6.02·10−4	6.02·10−4
SI3	W	3.81·10−2	9.96·10−2	2.21·10−1	6.02·10−4
	M	4.14·10−3	4.24·10−3	7.95·10−3	6.02·10−4
	SD	1.13·10−2	1.16·10−2	3.96·10−2	3.25·10−19
	B	4.54·10−4	4.54·10−4	4.22·10−4	4.22·10−4
SI4	W	4.54·10−4	8.13·10−2	4.22·10−4	4.22·10−4
	M	4.54·10−4	3.15·10−3	4.22·10−4	4.22·10−4
	SD	1.63·10−19	1.45·10−2	5.42·10−20	5.42·10−20

, the results of 30 independent runs are presented for all 9S20 curves from Figure 2 using DE/rand/1/exp and DE/rand/2/exp for the Elliot and Langevin expressions. In

Table 7. B, W, M and SD for 30 independent runs for all Alnico curves using DE/rand/1/exp and DE/rand/2/exp for the Elliot and Langevin expressions.

Curve			Method
		DE/rand/1/exp	DE/Rand/2/Exp	DE/Rand/1/Exp	DE/Rand/2/Exp
		Elliot	Elliot	Langevin	Langevin
	B	4.85·10−4	4.85·10−4	4.36·10−4	4.36·10−4
AML	W	4.85·10−4	4.85·10−4	4.36·10−4	4.36·10−4
	M	4.85·10−4	4.85·10−4	4.36·10−4	4.36·10−4
	SD	1.63·10−19	1.63·10−19	2.71·10−19	2.71·10−19
	B	3.24·10−4	3.34·10−4	3.95·10−4	3.95·10−4
AMR	W	3.24·10−4	3.24·10−4	3.95·10−4	3.95·10−4
	M	3.24·10−4	3.24·10−4	3.95·10−4	3.95·10−4
	SD	0	0	5.42·10−20	5.42·10−20
	B	7.36·10−4	7.36·10−4	7.28·10−4	7.28·10−4
AD1	W	7.36·10−4	7.36·10−4	7.28·10−4	7.28·10−4
	M	7.36·10−4	7.36·10−4	7.28·10−4	7.28·10−4
	SD	1.08·10−19	1.08·10−19	2.17·10−19	2.17·10−19
	B	7.22·10−4	7.22·10−4	2.61·10−4	2.61·10−4
AD2	W	7.22·10−4	7.22·10−4	2.61·10−4	2.61·10−4
	M	7.22·10−4	7.22·10−4	2.61·10−4	2.61·10−4
	SD	3.25·10−19	3.25·10−19	1.08·10−19	1.08·10−19
	B	3.14·10−4	3.14·10−4	2.01·10−4	2.01·10−4
AD3	W	3.14·10−4	3.14·10−4	2.50·10−1	2.01·10−4
	M	3.14·10−4	3.14·10−4	8.52·10−3	2.01·10−4
	SD	1.63·10−19	1.63·10−19	4.48·10−2	1.08·10−19
	B	1.56·10−4	1.56·10−4	8.48·10−4	8.48·10−4
AI1	W	1.56·10−4	1.56·10−4	8.48·10−4	8.48·10−4
	M	1.56·10−4	1.56·10−4	8.48·10−4	8.48·10−4
	SD	5.42·10−20	2.71·10−20	5.42·10−19	5.42·10−19
	B	1.75·10−4	1.75·10−4	2.12·10−4	2.12·10−4
AI2	W	1.75·10−4	1.75·10−4	2.12·10−4	2.12·10−-4
	M	1.75·10−4	1.75·10−4	2.12·10−4	2.12·10−4
	SD	2.71·10−20	2.71·10−20	5.42·10−20	5.42·10−20
	B	2.78·10−4	2.78·10−4	2.83·10−4	2.83·10−4
AI3	W	2.78·10−4	2.78·10−4	2.83·10−4	2.83·10−4
	M	2.78·10−4	2.78·10−4	2.83·10−4	2.83·10−4
	SD	1.08·10−19	1.08·10−19	1.08·10−19	1.08·10−19

, the results of 30 independent runs are presented for all Alnico curves from Figure 2 using DE/rand/1/exp and DE/rand/2/exp for the Elliot and Langevin expressions.

The numbers of best results achieved for the combination of expression and method presented in Table 6 and Table 7 are shown in

Table 8. Numbers of lowest B and M achieved using DE/rand/1/exp and DE/rand/2/exp for the Elliot and Langevin expressions.

Expression and Method	Number of the Lowest B	Number of the Lowest M
Elliot using DE/rand/1/exp	12	9
Elliot using DE/rand/2/exp	12	6
Langevin using DE/rand/1/exp	6	7
Langevin using DE/rand/2/exp	6	9

.

Based on Table 8, it can be seen that the best results were obtained using a combination of the Elliot expression and DE/rand/1/exp for the determination of the parameters. Based on the presented results, we decided to use the Elliot expression, written in (7), to model the hysteresis curve, and to use DE/rand/1/exp for determination of the parameters. The meaning of parameters P₁, P₂ and P₃ are presented in Figure 6, and the meaning of parameters P₄, P₅ and P₆ is presented in Figure 7.

Parameter P₃ moves the curve according to the H axis, P₆ according to the B axis, and parameters P₁, P₂, P₄ and P₅ influence the shape and size of the curve.

In the continuation, the use of the Elliot expression is presented for materials 9S20 and Alnico.

2.4. Discussion of Calculated Values with Respect to Measurement Error

Measurements were made on the device for the characterization of semi and hard magnetic materials. A Magnet-Physic yoke MIR-5 was used for the measurements. A sample of massive material was made with a circular cross-section (diameter of the sample was 10 mm) was used and it was 80 mm long. A tangential Hall probe HS-TGB5-104005 with basic accuracy 0.25% (to 3T) was used for measurement of H. The Hall probe was connected with a Magnet-Physic TESLAMETER with basic accuracy 0.3% (without probe). So that the accuracy of the H measurement, considering the accuracy of the TESLAMETER and probe, was 0.55%. In the measurements we also have to consider that, between the probe which was placed fixed near the sample and the sample, a small air gap appears. The magnetic field in the cross-section of the tested material was not completely homogeneous. The real accuracy for the measurement of H was around 5%.

B was measured by induction in the measuring coil with 500 turns, coiled around the center of the sample and connected with a Fluxmeter. A Magnet-Physic Electronic Fluxmeter EF5 was used for the measurement of B. The basic accuracy of the Fluxmeter was 0.25%. The frequency of the measurement was 10 mHz. Although the measuring coil was coiled tightly around the sample, there was still some air between the coil and the sample. For the measurement of B, it should be considered that the magnetic field in the cross-section of the tested material was not completely homogeneous. Because of that, the real accuracy for the measurement of B was around 5%.

The measured and calculated values of B in point A (presented in Figure 5) for 9S20 material and point K (presented in Figure 6) for Alnico, are presented in

Table 9. Measured and calculated values of B at points A and K.

Material: 9S20; Considered Point A (Figure 4); Measured H = 4100 A/m; Measured B = 1.50 T (Measurement Accuracy Range from 1.43 T to 1.58 T Considering 5% Accuracy of the Measurement)
Expression	LA1	GR1	HY1	AL1	LO1	SI1	EL1
Calculated B (T)	1.55	1.60	1.58	1.56	1.58	1.58	1.52
Material: Alnico; Considered Point K (Figure 5); Measured H = 80400 A/m; Measured B = 1.76 T (Measurement Accuracy Range from 1.67 T to 1.85 T Considering 5% Accuracy of Measurement)
Expression	LA1	GR1	HY1	AL1	LO1	SI1	EL1
Calculated B (T)	1.79	1.80	1.81	1.80	1.82	1.81	1.77

¹ LA-Langevin, GR-Grompertz, HY-Hyptan, AL-Algebraic, LO-Logistic, SI-Sigmoid, EL-Elliot.

.

From Table 9 it can be seen that three calculated values at point A were inside of the measurement accuracy range, three vales were at the border, and one value was outside of the measurement accuracy range. It can also be seen that all calculated values at point K were inside of the measurement accuracy range. Points A and K are points where the deviation between the measured value and the calculation is higher, but they are important for Hysteresis modeling, as will be presented in Section 3.2 and Section 4.2. Taking into account the measurement accuracy and the results presented in Table 9, we can conclude that it is reasonable to search for expressions and methods to match the measured and calculated curves best (some calculated values were at the border or outside of the measurement accuracy range).

3. Modeling of the Material 9S20

For better understanding of the hysteresis modeling procedure, modeling is explained briefly in steps, and, after that, each step is explained in detail in the subsections. The steps are presented in Figure 8.

In step 1, presented in Figure 8a, parameters P₁ to P₆ are determined for the measured curve SML, based on the measured curve using Elliot’s expression and the DE/rand/1/exp method.

In step 2, presented in Figure 8b, parameters P₁ and P₆ are adjusted according to the appropriate points (in the case of SML, points A and B).

In step 3, presented in Figure 8c, all measured curves are replaced with analytical expressions repeating steps 1 and 2 for each measured curve.

In step 4, presented in Figure 8d, any required curve with respect to magnetization changes, as shown in Figure 1, can be calculated on the basis of two points, for example, points E1 and B.

3.1. Calculation of the Parameters Using DE/Rand/1/exp

Major hysteresis loop left and right parts (SML and SMR), first-order decreasing reversal curves (SD1–SD4) and first-order increasing reversal curves (SI1–SI4) were measured in the case of the material 9S20. The first step in the modeling procedure is to determine the Elliot expression parameters using DE/rand/1/exp. 30 independent runs were used for each measured curve. The parameters obtained with the best OF are presented in

Table 10. Calculated parameters of the 9S20 major hysteresis loop.

Parameter	Elliot		Langevin
	SML	SMR	SML	SMR
P1	19.63	22.70	1.26	1.27
P2	12.73	14.46	180.3	188.1
P3	429.1	−436.0	436.8	−443.7
P4	5704	6856	61.21	59.61
P5	21.13	17.52	−3.55·10−3	9.53·10−4
P6	1.59·10−3	−4.06·10−3

,

Table 11. Calculated parameters of the 9S20 decreasing first-order reversal curves.

Parameter	Elliot				Langevin
	SD1	SD2	SD3	SD4	SD1	SD2	SD3	SD4
P1	15.23	13.29	17.86	23.82	4.91·10−1	5.66·10−1	7.96·10−1	1.09
P2	17.55	18.44	19.06	18.27	335.0	128.5	104.9	139.9
P3	369.0	322.8	339.6	408.6	391.4	338.3	346.4	412.0
P4	22721	6437	4716	6094	81.83	115.0	116.9	92.03
P5	38.90	90.53	92.49	54.64	−6.83·10−1	−4.51·10−1	−2.23·10−1	−5.08·10−3
P6	−6.73·10−1	−4.36·10−1	−2.15·10−1	−5.04·10−3

and

Table 12. Calculated parameters of the 9S20 increasing first-order reversal curves.

Parameter	Elliot				Langevin
	SI1	SI2	SI3	SI4	SI1	SI2	SI3	SI4
P1	13.27	12.43	17.66	19.26	4.91·10−1	5.66·10−1	7.96·10−1	9.83·10−1
P2	15.28	17.25	18.86	16.67	335.0	128.5	104.9	121.5
P3	−369.0	−322.8	−339.9	−377.9	−391.4	−338.3	−346.4	−382.1
P4	19784	6021	4665	4705	81.83	115.0	116.9	103.0
P5	38.90	90.53	92.49	72.88	6.83·10−1	4.51·10−1	2.23·10−1	1.02·10−1
P6	6.73·10−1	4.36·10−1	2.15·10−1	9.90·10−2

. The values of B, W, M and SD were presented previously in Table 6. In Table 10, Table 11 and Table 12 parameters obtained with the best OF of the Langevin expression are also shown, in order to enable the analysis of the results.

3.2. Adjusted Parameters P₁ and P₆ for 9S20

From Figure 4 it can be seen that the calculated curve SML does not flow through the hysteresis final points A and B. To adjust the calculated curve to points A and B only parameters P₁ and P₆ should be adjusted, using the system of two equations written in (11) and (12).

$$B_{\mathrm{A}}=\frac{P_1\left(H_{\mathrm{A}}+P_3\right)}{P_4+P_2\left|H_{\mathrm{A}}+P_3\right|}+P_5{\mu }_0\left(H_{\mathrm{A}}+P_3\right)+P_6$$

(11)

$$B_{\mathrm{B}}=\frac{P_1\left(H_{\mathrm{B}}+P_3\right)}{P_4+P_2\left|H_{\mathrm{B}}+P_3\right|}+P_5{\mu }_0\left(H_{\mathrm{B}}+P_3\right)+P_6.$$

(12)

The calculated SMR first-order decreasing and increasing reversal curves also do not flow through the final points, and parameters P₁ and P₆ were adjusted as a result. The values of adjusted parameters are presented in

Table 13. P₁ and P₆ for SML and the right part of the measured major hysteresis loop (SMR).

Points and Parameters	Curve
	SML	SMR
Points	A, B	A, B
P1	19.67	22.77
P6	−2.66·10−2	2.47·10−2

,

Table 14. P₁ and P₆ for the 9S20 decreasing first-order reversal curves.

Points and Parameters	Curve
	SD1	SD2	SD3	SD4
Points	B, C	B, D	B, E	B, F
P1	14.72	12.45	17.04	23.22
P6	−6.96·10−1	−4.55·10−1	−2.27·10−1	−8.24·10−2

and

Table 15. P₁ and P₆ for the 9S20 increasing first-order reversal curves.

Points and Parameters	Curve
	SI1	SI2	SI3	SI4
Points	A, G	A, H	A, I	A, J
P1	12.79	11.61	16.83	18.65
P6	6.96·10−1	4.53·10−1	2.25·10−1	1.19·10−1

, where the points used for adjustment are also presented (points are shown in Figure 9 and Figure 10).

The parameters of each of the curves (decreasing and increasing), written in Table 11, Table 12, Table 14 and Table 15, were obtained using a measured curve, and, with that, they were adjusted to the form of the curve. In this way, no special adjustment is required with respect to increasing and decreasing curves. Comparing similar increasing and decreasing curves (for example SI2 and SD2, shown in Figure 2), it can be seen that parameters P₁, P₂, P₄ and P₅, which influence the shape and size, are similar. Parameter P₃, which moves the curve according to the H axis, is different, and also parameter P₆, which moves the curve according to the B axis, is different.

The calculated and measured major hysteresis loop and first-order decreasing reversal curves are presented in Figure 9, and the first-order increasing reversal curves are presented in Figure 10.

3.3. Calculation of Curves in the Area between the Measured Curves for 9S20

The main purpose of the presented model is to calculate the magnetization in the area between the measured curves simply. To calculate magnetization, only points where the excitation starts to decrease after increasing (for example point C marked in Figure 9 and point B marked in Figure 9) are needed to calculate the decreasing magnetization. The parameters of the nearest curve are used, obtained based on the measured curve, and parameters P₁ and P₆ are adjusted using a system of two equations, such as written in (11) and (12), to get the appropriate magnetization. Calculated curves are presented in Figure 11 and marked with SE1–SE5. For example, curve SE1 is calculated based on the nearest curve obtained based on measurements, which is SD4. Point E1, where the excitation starts to decrease, and point B, are used to calculate parameters P₁ and P₆, and parameters P₂ to P₅ are the same as used for the curve SD4.

Parameters of the calculated curves SE1–SE5, together with used points and the nearest curve obtained based on the measured curve, which is the basis for the curve calculation, are presented in

Table 16. Parameters of the calculated curves SE1–SE5.

Basis curve, Points and Parameters	Curve
	SE1	SE2	SE3	SE4	SE5
Basis	SD4	SD4	SD3	SD3	SD2
Points	E1, B	E2, B	E3, B	E4, B	E5, B
P1	22.54	21.77	18.24	15.52	13.53
P2	18.27	18.27	19.06	19.06	18.44
P3	408.6	408.6	339.6	339.6	322.8
P4	6094	6094	4716	4716	6437
P5	54.64	54.64	92.49	92.49	90.53
P6	−1.17·10−1	−1.55·10−1	−1.67·10−1	−3.01·10−1	−4.01·10−1

.

The same procedure can be used for increasing magnetization. The values of B and H in the points presented in Figure 9, Figure 10 and Figure 11 are determined in

Table 17. B and H in the points presented in Figure 9, Figure 10 and Figure 11.

Point	H(A/m)	B(T)	Point	H(A/m)	B(T)
A	4104	1.50	I	−790	−6.39·10−1
B	−4082	−1.50	J	−1297	−9.92·10−1
C	290	−3.81·10−1	E1	1170	1.01
D	466	1.03·10−1	E2	1030	9.11·10−1
E	795	6.39·10−1	E3	880	7.70·10−1
F	1849	1.18	E4	590	4.50·10−1
G	−290	3.81·10−1	E5	490	2.05·10−1
H	−466	−1.03·10−1

.

The calculated and measured major hysteresis loop, decreasing and increasing curves are presented in Figure 12.

4. Modeling of the Material Alnico

The calculation steps of the Alnico hysteresis model are the same as presented briefly in Section 3.

4.1. Calculation of the Parameters Using DE/Rand/1/exp

A major hysteresis loop’s left and right parts (AML and AMR), decreasing curves (AD1–AD3) and increasing curves (AI1–AI3) were measured in the case of the material Alnico (presented in Figure 2). The Elliot expression parameters were determined using DE/rand/1/exp. 30 independent runs were used for each measured curve. The parameters obtained with the best OF are presented in

Table 18. Calculated parameters of the Alnico major hysteresis loop.

Parameter	Elliot		Langevin
	AML	AMR	AML	AMR
P1	20.83	30.24	1.55	1.55
P2	12.64	18.46	1500	1469
P3	48,378	−47,725	48,339	−47,729
P4	35,799	50,282	2.65	2.68
P5	2.11	2.16	−1.76·10−1	1.82·10−1
P6	−1.82·10−1	1.86·10−1

,

Table 19. Calculated parameters of the Alnico decreasing curves.

Parameter	Elliot			Langevin
	AD1	AD2	AD3	AD1	AD2	AD3
P1	14.96	16.95	23.37	1.13	1.30	1.40
P2	11.61	11.75	15.44	2000	2000	1570
P3	46,672	46,162	46,395	46,664	46,064	46,368
P4	53,560	49,604	47,693	4.38	3.66	3.35
P5	3.21	2.62	2.54	−5.43·10−1	−4.04·10−1	−2.93
P6	−5.46·10−1	−4.08·10−1	−2.95·10−1

and

Table 20. Calculated parameters of the Alnico increasing curves.

Parameter	Elliot			Langevin
	AI1	AI2	AI3	AI1	AI2	AI3
P1	8.23	13.96	16.07	7.78·10−1	1.10	1.37
P2	8.12	11.56	11.00	2000	1982	1477
P3	−44,120	−45,114	−46,485	−44,362	−45,089	−46,491
P4	54,719	45,731	30,288	5.50	4.19	3.43
P5	3.52	3.38	2.75	7.78·10−1	5.47·10−1	3.10·10−1
P6	7.73·10−1	5.49·10−1	3.11·10−1

. The values of B, W, M and SD were presented previously in Table 7. In Table 18, Table 19 and Table 20 parameters obtained with the best OF of the Langevin expression are also shown in order to enable an analysis of the results.

4.2. Adjusted Parameters P₁ and P₆ for Alnico

To adjust curves to the start and end points, parameters P₁ and P₆ were corrected using the system of two equations, such as in (11) and (12). The values of adjusted parameters are presented in

Table 21. P₁ and P₆ for AML and AMR.

Points and Parameters	Curve
	AML	AMR
Points	K,L	K,L
P1	20.57	29.90
P6	−1.74·10−1	1.76·10−1

,

Table 22. P₁ and P₆ for the Alnico decreasing curves.

Points and Parameters	Curve
	AD1	AD2	AD3
Points	L, M	L, N	L, O
P1	14.57	16.68	23.06
P6	−5.13·10−1	−3.77·10−1	−2.77·10−1

and

Table 23. P₁ and P₆ for the Alnico increasing curves.

Points and Parameters	Curve
	AI1	AI2	AI3
Points	K, P	K, R	K, S
P1	8.12	13.75	15.95
P6	7.56·10−1	5.41·10−1	3.01·10−1

, where the points used for adjustment are also presented.

4.3. Calculation of Curves in the Area between the Measured Curves for Alnico

For Alnico, the curves in the areas between the calculated curves, based on measured curves, were determined using the point where the calculation started to decrease after increasing, and the point in the lower left part of the hysteresis. Calculated curves AC1–AC5 are presented in Figure 13. For example, curve AC1 was calculated based on the nearest curve obtained based on the measurement, which is AD3. Point C1, where the excitation started to decrease, and point L, were used to calculate parameters P₁ and P₆, using the system of two equations, such as in (11) and (12). Parameters P₂ to P₅ are the same as used for the curve AD3. The new magnetization curve was calculated simply, based on only two points.

The parameters of the calculated curves AC1–AC5, together with used points and the nearest curve obtained based on the measured curve, which was the basis for the curve calculation, are presented in

Table 24. Parameters of the calculated curves AC1–AC5.

Basis Curve Points and Parameters	Curve
	AC1	AC2	AC3	AC4	AC5
Basis	AD3	AD3	AD2	AD1	AD1
Points	C1, L	C2, L	C3, L	C4, L	C5, L
P1	23.33	22.31	15.79	13.62	12.34
P2	15.44	15.44	11.75	11.61	11.61
P3	46,395	46,395	46,162	46,672	46,672
P4	47,693	47,693	49,604	53,560	53,560
P5	2.54	2.54	2.62	3.21	3.21
P6	−2.60·10−1	−3.21·10−1	−4.45·10−1	−5.86·10−1	−6.83·10−1

.

The same procedure can be used for increasing magnetization. The values of B and H in the points presented in Figure 12 and Figure 13 are determined in

Table 25. B and H in the points presented in Figure 12 and Figure 13.

Point	H(A/m)	B(T)	Point	H(A/m)	B(T)
K	80,411	1.76	S	−54,567	−1.46
L	−81,255	−1.76	C1	60,600	1.55
M	46,922	1.06	C2	53,600	1.40
N	50,178	1.30	C3	47,800	1.15
O	55,967	1.50	C4	45,000	9.01·10−1
P	−46,900	−5.78·10−1	C5	43,600	6.92·10−1
R	−48,867	−1.00

.

5. Limits of the Parameters Used by Evolutionary Methods

When an evolutionary method is used, appropriate limits of the parameters should be set. Our test cases, 9S20 and Alnico, are materials with very different coercive magnetic field strength. Because of that, the limits of the parameters must contain a wide range of different values. We set limits in a way that they were appropriate for a wide range of different materials. The limits used are presented in

Table 26. Parameters’ limits used in the case of evolutionary calculations.

Parameter	Lower Limit	Upper Limit
P1	0	50
P2	−20	20
P3	−60,000	60,000
P4	−5	60,000
P5	0	200
P6	−3	3

.

Although the limits were set wide, evolutionary methods were able to calculate parameters.

6. The Use of the Hysteresis Model

The hysteresis model is intended to be used for the transient magnetic field calculation. The algorithm presenting the use of the hysteresis model is shown in Figure 14. The curves of Alnico, presented in Figure 2, are used in the algorithm for better understanding. The basis for the model is the mathematical description of each measured curve.

Different colors are used in Figure 14 to make the algorithm more transparent. The parts of the algorithm when a new magnetization curve must be calculated are marked with red. The parts of the algorithm when the magnetization curve remains unchanged due to the previous step or an original magnetization curve is used are marked with green. The calculation depends on the excitation change, and on the calculated values of H and B from the previous step, and these data must be memorized (marked with blue color).

7. Conclusions

Precise magnetic material hysteresis modeling is important for accurate dimensioning of electromagnetic devices using FEM. The aim of the paper was to find a mathematical expression that describes a wide range of different hystereses authentically. Seven different mathematical expressions were modified appropriately with added parts, and tested on two same test examples. The results presented in Table 6, Table 7 and Table 8 show that the combination of the Elliot expression and DE/rand/1/exp method for the determination of the parameters is the best between the tested expressions and methods. Both the highest number of the lowest best OF and the highest number of the lowest mean OF were obtained using the Elliot expression combined with DE/rand/1/exp.

The idea of the model is that both major and first-order reversal measured curves could be modeled using the Elliot expression. For each measured curve, parameters were determined using DE/rand/1/exp, and, with that, the best match was achieved between measured and calculated curves. After that, a simple procedure is presented for the calculation of the curves in the area between the measured curves. Therefore, DE/rand/1/exp offers a reliable procedure for the precise determination of the parameters.

A large number of measured first-order reversal curves allow the model to function well. To simulate hysteresis, the closest curve should be selected, along with its set of six parameters. Therefore, it is impossible to simulate a curve if a similar experimental curve was not fitted previously. The model is, therefore, a fit of experimental data, and, consequently, it has no predictive ability, which is the limitation of the model.

Considering the measurement accuracy (5% for measurement of H and 5% for measurement of B) and results presented in Table 9, it can be concluded that it is reasonable to search for expressions and methods to match best the measured and calculated curves.

It can also be concluded that different algorithms lead to the same results. For almost all expressions in the case of the same calculated OF, calculated parameters are also the same. Although the same best OF is obtained using different methods, mean OF can be very different, which means that some methods occasionally give poor results on this problem, which is not desirable. Because of this, we have identified the method that is best for the presented problem.

The presented model shows the successful use of the evolutionary method for the determination of a higher number of parameters and simple determination of further magnetization processes. Findings about the used evolutionary methods, which are presented in the paper, will be used in the future to get better solutions of the optimization problems in Electromagnetics. The integration of the presented model into the FEM software will also be one of the future studies.