Analysis of Higher-Order Bézier Curves for Approximation of the Static Magnetic Properties of NO Electrical Steels

Abstract

Adequate mathematical description of magnetization curves is indispensable in engineering. The accuracy of the description has a significant impact on the design of electric machines and devices. The aim of this paper was to analyze the capability of Bézier curves systematically, to describe the nonlinear static magnetic properties of non-oriented electrical steels, and to compare this approach versus the established mathematical descriptions. First, analytic functions versus measurements were analyzed. The Bézier curves were then compared systematically with the most adequate analytic functions. Next, the most suitable orders of Bézier curves were determined for the approximation of nonlinear magnetic properties, where the influence of the range of the input measurement dataset on the approximation process was analyzed. Last, the extrapolation capabilities of the Bézier curves and analytic functions were evaluated. The general conclusion is that Bézier curves have adequate flexibility and significant potential for the approximation and extrapolation of nonlinear properties of non-oriented electrical steels.

1. Introduction

The magnetic properties of soft magnetic materials are highly nonlinear and reflect intricate underlying phenomena in hysteresis behavior. In engineering applications, these properties are often described using basic magnetization curves, e.g., the anhysteretic curve, the first magnetization curve, the major loop, and the minor loops. For example, the major loop represents the input data of several well-known hysteresis models [1,2,3]. Examples of using the anhysteretic curve as input data are the widespread Jiles–Atherton hysteresis model and its variations [4,5] and various finite element simulation tools. Hysteresis models are used to describe nonlinear magnetization and power loss in iron cores, whereas the nonlinear behavior is indispensable in the design of electromagnetic machines and devices [6,7]. For these reasons, it is important to find adequate methods or analytic functions to describe such nonlinear curves with high accuracy and with low computational power.

Researchers have applied various approaches to model all the discussed magnetization curves. On the one hand, researchers try to establish adequate links between the physical parameters of materials and various mathematical descriptions, which would result in adequate physics-based material models [4]. This is a very complex task, and essential for understanding the intricate phenomena within soft magnetic materials, and, consequently, for their design and optimization processes. The complexity arises from the fact that hysteresis is a macroscopic property that emerges from the collective behavior of microscopic entities, and this behavior can be influenced by a multitude of factors [8]. On the other hand, researchers also try to model nonlinear magnetization in a more practical way to predict the magnetization behavior, which is based on measured datasets. These models are phenomenological, i.e., they simply attempt to describe nonlinear relationships with the assumption that the modeled relationships extend past the measured values without a deep understanding of the underlying mechanisms [9,10,11,12]. These models are equally as important as the physical ones, and are essential in applied engineering, where the existing soft magnetic materials are applied in electromagnetic devices and machines. The most important requirements are accuracy and computational simplicity. A common and straightforward method is to use mathematical functions [2,10,13] or curves [14] that approximate the data obtained from measurements. This method has its difficulties, as it needs costly and specific instruments for measurement, and the data obtained may be limited in both amount (i.e., volume) and quality. To measure the magnetization curves, the measurement equipment can reach different levels of saturation, depending on its capability. The equipment also has different degrees of accuracy that affect the noise and errors in the obtained curves. The quality of the approximations depends on several factors: the size and noise level of the input dataset, the properties of the functions used, and the type of the fitting method. For example, an analysis based on the application of sigmoid functions is described in [15]. The authors investigated the relevant sigmoid functions and used them to construct the major hysteresis loop. Recently, the sigmoid functions from [15] were modified and fitted using different optimization algorithms to match measured data [10]. The modified analytic functions offered more degrees of freedom during the fitting process. In this paper, we examined different methods for estimating and extending the data obtained from measurements.

An approach used in many fields for the description of a wide range of different shapes and curves is (classical) Bézier curves. Bézier curves enable the construction of free-form curves and surfaces. The shape of a Bézier curve depends on the placement of the so-called control points. This feature makes them intuitive, easy to implement, and flexible. It is possible to construct very complex shapes only by placing control points in adequate positions. Bézier curves were first used for the shape optimization of automobile bodies [16]. Nowadays, these curves are used widely in computer graphics and animation. Examples of using them in computer graphics are font generation (Chinese [16] and Arabic [17]), image processing, and path design in various animations. However, they are not limited to the field of computer graphics. Due to their simplicity and continuous curvature, Bézier curves are being used in path planning for ground vehicles [18,19,20] and for unmanned aerial vehicles [21]. Their free-form ability made them popular in the field of shape design. An example of shape optimization with Bézier curves is the modeling of laminated plates of variable stiffness [22]. To expand their usability in shape optimization, Bézier curves are used to construct circular arcs in CAD applications (e.g., SolidWorks) [23] using mainly cubic Bézier curves [22]. The free-form ability to represent various shapes and geometries enhanced and is useful in the field of additive manufacturing, e.g., for the topology optimization of electrical machines [24,25]. Their ability to fit any dataset smoothly also makes them useful in the field of electrical engineering. The applications include the modeling of the flux linkage characteristics of a switched reluctance motor [26] and the modeling of PV module characteristics [27,28].

Classical Bézier curves are often used for simple smooth curve approximations. The number of control points defines the order of the corresponding Bézier curve. The order can be set to any arbitrary number for more flexibility, but it will increase the complexity due to an increasing number of control points [22]. Furthermore, Bézier curves constructed by a large number of control points are numerically unstable [19], and the calculation is slower since additional coordinates have to be calculated. For this reason, an alternative approach is to join low-order Bézier curves together in a smooth way to form a composite curve. This concept is called piecewise Bézier curves, and it builds a curve as a union of Bézier curves connected end to end [29]. Using this concept, even more complex shapes can be described by joining a larger number of lower orders of Bézier curves. A few rules must be followed to obtain a smooth curve when joining Bézier curves, e.g., the so-called parametric continuity rules [22,30]. Bézier curves are also used for spline interpolations, especially in the known B-spline group. Any piecewise Bézier curves with arbitrary order can be converted into a B-spline, and any B-spline can be converted into one or more Bezier curves [29]. The B-spline curve has one main disadvantage compared with piecewise Bézier curves. The B-spline concept requires more computation, and the mathematics used to describe them are more complicated [29].

Based on the discussed flexibility that Bézier curves offer, especially the smooth fitting of datasets, we performed an analysis of the Bézier curve fitting process on various fundamental magnetization curves of non-oriented (NO) electrical steels (i.e., the anhysteretic curve, the major loop, and the first magnetization curve). This work is an extension of the preliminary study [31], and the aim was to investigate the approximation capabilities of classical Bézier curves systematically. We limited this study to classical Bézier curves because of their straightforward approximation of datasets without the need to implement additional rules, e.g., the continuity rules [22,30]. To demonstrate their ability to approximate magnetization curves, we identified established analytic functions [10,15] in parallel and performed a systematic comparison based on the deviation from the measured curves. We analyzed the deviation from the measured data in different regions of individual curves, e.g., the high-permeability region and the region of high saturation. This analysis resulted in the identification of the most suitable orders of classical Bézier curves and the most adequate analytic functions for different input datasets. The high-saturation region is, in most cases, very challenging to evaluate experimentally; therefore, we further analyzed the extrapolation capabilities of the identified Bézier curves and analytic functions. The main challenge during the curve-fitting process was to find the position of the control points, which determine the shape and curvature of the curve [18]. The main goal was to find a precise and elegant way to represent the observed data mathematically beyond the range of observation and to examine how the size of the input data sets affected the results.

The article is split into six sections. Section 2 presents the considered analytic functions, as well as the description of the theoretical background and properties of the Bézier curves. Two methods are presented for calculating the Bézier curves. The materials and the methodology used are presented in Section 3. Section 4 contains the obtained results. The influence of the order of Bézier curves on the approximation of the measured data is analyzed. Additionally, the extrapolation capabilities of Bézier curves of selected orders and analytic functions are analyzed. In Section 5, a discussion is provided. Finally, in Section 6, the concluding remarks are presented with a description of future work.

2. Mathematical Descriptions of Nonlinear Magnetic Properties

Mathematical functions, which are adequate for the description of the nonlinear curves in hysteresis, must have specific properties. They must be continuous, differentiable, single valued, and monotone; have odd symmetry; approach saturation polarization $J_{\mathrm{s}}$, as magnetic field strength $H$ tends toward negative or positive infinity (e.g., in ideal cases they have asymptotes); and must have a knee point where the second derivative reaches its minimal value [4].

Natural candidates that have those properties are the so-called sigmoid functions. The aim of this work is to introduce and analyze a novel concept for the approximation of nonlinear magnetic properties that is based on curves instead of functions, where classical Bézier curves of different orders are considered.

2.1. Analytic Functions

Analytic functions that describe a curve shaped in the form of the letter “S” are called sigmoid functions. Besides their shape, sigmoid functions have further important properties. Sigmoid functions are real valued, differentiable, and have a pair of horizontal asymptotes [15]. Due to their intrinsic properties, sigmoid functions are suitable candidates to describe nonlinear magnetic properties. The described favorable basic properties are the reason for the frequent use of sigmoid functions in the description of the discussed magnetization curves [4]. The most adequate (and, for the discussed purpose, applied) sigmoid functions [10,15] are presented in

Table 1. Mathematical equations of selected analytic functions.

Analytic Function	Equation
Logistic	$\frac{2a}{d+{\mathrm{e}}^{-bx-c}}-a$	(1)
Hyperbolic tangent	$a·\tanh \left(b·\left(x+c\right)\right)$	(2)
Elliot	$\frac{a·\left(x+c\right)}{d+b·\left|x+c\right|}$	(3)
Gompertz	$2a·{\mathrm{e}}^{-{\mathrm{e}}^{-bx-c}}-a$	(4)
Langevin	$a·\left(\coth \left(\frac{x+c}{b}\right)-\frac{b}{x+c}\right)$	(5)
Inverse hyperbolic sine	$a·\mathrm{asinh}\left(b·\left(x+c\right)\right)$	(6)

with Functions (1)–(5). In addition, many functions exist that have a similar shape to sigmoid functions but do not have all their properties, e.g., do not have horizontal asymptotes. Such functions are still applicable to the discussed problem and are widely used in applied engineering. During the preliminary analysis [31], we determined that the inverse hyperbolic sine function is promising, especially in the saturation region, because of its logarithmic character for higher input values. Therefore, the inverse hyperbolic sine of Function (6) was included in the presented analysis.

Original analytic functions were modified with additional parameters to enable better agreement with the measured data during the approximation process. Functions (1)–(6) all have parameters:

• $a$—determines the value of the asymptote;
• $b$—controls the slope of the curve;
• $c$—enables the shifting of the function along the x-axis in an xy-coordinate system [15].

Functions (1) and (3) have a fourth parameter $d,$ which also influences the slope [15].

Functions (1)–(4) are analytic functions without any physical background. In contrast to this, the Langevin function (Function (5)) has a physical origin. The Langevin function is a special case of the quantum mechanical Brillouin function when the function is limited to infinity [32]. The function was originally derived for paramagnets. Because of its physical foundation, it is used in many models to describe the magnetization curves, e.g., in the Jiles–Atherton hysteresis model [5]. Further, in applied engineering, the so-called Frölich–Kennelly approximation is often used for the description of nonlinear properties, especially for extrapolation purposes [6,33]. The Frölich–Kennelly equation has an equivalent mathematical description in the Elliot function (Function (3)). Therefore, all results obtained for the Elliot function also apply for the Frölich–Kennelly description.

In addition, we analyzed the inverse hyperbolic sine function (Function (6)). Generally, the inverse hyperbolic sine function has all the properties needed for adequate approximation of nonlinear magnetic properties. The difference compared with sigmoid functions (Functions (1)–(5)) is that the inverse hyperbolic sine function does not have asymptotes. However, the function value of the inverse hyperbolic sine function increases almost linearly (with a small slope), as the input x values increase at very high values. For that reason, parameter $a$ has a scaling role for the inverse hyperbolic sine function.

Figure 1 presents a graphical comparison of the functions from Table 1. The parameters were set to $a=1$, $b=1$, $c=0$, and $d=1$ for the sigmoid functions. Such a choice ensured that Functions (1)–(5) approached horizontal asymptotes at −1 and 1, as the x-value tended to negative and positive infinities, respectively. For the sake of graphic comparison, the parameters of the inverse hyperbolic sine function were set to $a=1/\mathrm{asinh}\left(5\right)$, $b=1$, and $c=0$ to ensure that the first and last point of the function reached values of −1 and 1 in the presented interval $x\in [-5, 5]$ in Figure 1.

2.2. Bézier Curves

Bézier curves are parametric curves with Bernstein polynomials as their basis functions. Paul de Casteljau contributed significantly to their first development, but they are named after Pierre Bézier, who invented them in 1962 and made them widely popular for the shape design of automobile bodies while working for Renault [16,19]. Parametric curves are not defined as $y\left(x\right)$ dependencies, but they describe the shape of a curve based on a (local) input variable $t$ that is defined along the curve. Variable $t$ is 0 at the start of the curve and 1 at the end. For each value of $t$, two values are defined (i.e., $x\left(t\right)$ and $y\left(t\right)$), which, combined, define the coordinates ($x\left(t\right),y\left(t\right)$) of all the points on the curve. This means that the coordinates of a Bézier curve are obtained using two polynomials $x\left(t\right)=B_{\mathrm{x}}\left(t\right)$ and $y\left(t\right)=B_{\mathrm{y}}\left(t\right)$ [34]. Due to this property, parametric curves enable the description of highly complex curves. The free-form ability of Bézier curves, explained in Section 1, is demonstrated in Figure 2. It is possible to form a single-valued curve, as shown in Figure 2a; self-intersecting curves, as shown in Figure 2b,c; and complex shapes, as shown in Figure 2d.

Bézier curves are defined by their order $n$ and a corresponding set of ($n+1$) control points $\mathbf{P}=({\mathrm{P}}_0,{\mathrm{P}}_1,\dots ,{\mathrm{P}}_{\mathrm{n}})$. Each control point has its x and y coordinates ${\mathrm{P}}_{\mathrm{i}}\left(x_{\mathrm{i}},y_{\mathrm{i}}\right)$ ($i=0, 1,\dots ,n$) in the 2D Cartesian coordinate system. Considering order $n$ and corresponding vector of control points $\mathbf{P}$, the Bézier curve $B\left(t\right)=\left(B_{\mathrm{x}}\left(t\right),B_{\mathrm{y}}\left(t\right)\right)$ can be expressed in polynomial form by Equations (7) and (8)

$$B_{\mathrm{x}}\left(t\right)=x\left(t\right)=\sum _{i=0}^n{b}_{i,n}\left(t\right)x_i, 0\le t\le 1,$$

(7)

$$B_{\mathrm{y}}\left(t\right)=y\left(t\right)=\sum _{i=0}^n{b}_{i,n}\left(t\right)y_i, 0\le t\le 1,$$

(8)

where $b_{i,n}\left(t\right)$ represents the Bernstein basis polynomials and ($x_{\mathrm{i}}$, $y_{\mathrm{i}}$) are the coordinates of control point ${\mathrm{P}}_{\mathrm{i}}$. The Bernstein basis polynomials of order $n$ are defined by Equation (9) [18]:

$$b_{i,n}\left(t\right)=\left(\begin{array}{c}n\\ i\end{array}\right)t^i{\left(1-t\right)}^{n-i}, i=0,\dots ,n$$

(9)

where $t\in \left[\mathrm{0,1}\right]$, in which 0 is the starting point of the curve and 1 is the end point of the curve [18], and $\left(\begin{array}{c}n\\ i\end{array}\right)$ is a binomial coefficient defined by Equation (10):

$$\left(\begin{array}{c}n\\ i\end{array}\right)=\frac{n!}{i!\left(n-i\right)!}$$

(10)

Besides Bézier curves, Bernstein polynomials are also known as basis functions for splines [35]. The basis functions of the Bézier curve can be enhanced to be able to reproduce various shapes, such as those presented in [30]. De Casteljau’s recursive algorithm is used to evaluate Bézier curve polynomials because of its numerical stability [18].

Another feature of Bézier curves is their matrix representation, which allows for a straightforward calculation of the polynomial coefficients for functions $B_{\mathrm{x}}\left(t\right)$ and $B_{\mathrm{y}}\left(t\right)$. For a Bézier curve of order $n$, the matrix equation is defined by Equation (11) [36]:

$$B_{\mathrm{x}}\left(t\right)=\left[\begin{array}{ccc}{t}^{\mathrm{n}}& t^{\mathrm{n}-1}& \begin{array}{cc}\dots & \begin{array}{cc}t& 1\end{array}\end{array}\end{array}\right]{\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}\left[\begin{array}{c}{x}_0\\ x_1\\ \begin{array}{c}⋮\\ \begin{array}{c}{x}_{\mathrm{n}-1}\\ x_{\mathrm{n}}\end{array}\end{array}\end{array}\right]$$

(11)

The $B_{\mathrm{x}}\left(t\right)$ polynomial of the Bézier curve is obtained with Equation (11). To obtain the $B_{\mathrm{y}}\left(t\right)$ polynomial of the Bézier curve, the y coordinates of control points ${\mathrm{P}}_{\mathrm{i}}$ must be considered in Equation (11). The matrix ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ is an $\left(n+1\right)\times \left(n+1\right)$-sized coupling matrix between the variable $t$ and the control point’s coordinates $\mathbf{P}$. By evaluating the matrix ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ for Bézier curves of orders of up to 6 [36,37], the general form of matrix ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ for curves of order $n$ is determined by Equation (12) [38]:

$${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}=\left[\begin{array}{ccccc}\left(\begin{array}{c}n\\ 0\end{array}\right)\left(\begin{array}{c}n\\ n\end{array}\right)& \left(\begin{array}{c}n\\ 1\end{array}\right)\left(\begin{array}{c}n-1\\ n-1\end{array}\right)& \dots & \left(\begin{array}{c}n\\ n-1\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right)& \left(\begin{array}{c}n\\ n\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)\\ \left(\begin{array}{c}n\\ 0\end{array}\right)\left(\begin{array}{c}n\\ n-1\end{array}\right)& \left(\begin{array}{c}n\\ 1\end{array}\right)\left(\begin{array}{c}n-1\\ n-2\end{array}\right)& \dots & \left(\begin{array}{c}n\\ n-1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)& 0\\ \left(\begin{array}{c}n\\ 0\end{array}\right)\left(\begin{array}{c}n\\ n-2\end{array}\right)& \left(\begin{array}{c}n\\ 1\end{array}\right)\left(\begin{array}{c}n-1\\ n-3\end{array}\right)& & 0& 0\\ ⋮& ⋮& & ⋮& ⋮\\ \left(\begin{array}{c}n\\ 0\end{array}\right)\left(\begin{array}{c}n\\ 1\end{array}\right)& \left(\begin{array}{c}n\\ 1\end{array}\right)\left(\begin{array}{c}n-1\\ 0\end{array}\right)& & 0& 0\\ \left(\begin{array}{c}n\\ 0\end{array}\right)\left(\begin{array}{c}n\\ 0\end{array}\right)& 0& \dots & 0& 0\end{array}\right]$$

(12)

Each non-zero element in ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ is calculated as the product of two binomial coefficients by applying Equation (10). The elements within ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ have a positive or negative sign depending on the matrix order. A general rule is that the main counter diagonal elements of matrix ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ always have a positive sign, and then the sign alternates between negative and positive for all other diagonals above the main counter diagonal. The formation of matrix ${\mathbf{M}}_{\mathbf{B}\mathbf{é}\mathbf{z},\mathbf{n}}$ is presented in [36,37].

2.2.1. Properties of Bézier Curves

The basic properties of Bézier curves are:

• First, control point ${\mathrm{P}}_0$ defines the beginning, and the last control point ${\mathrm{P}}_{\mathrm{n}}$ defines the end of the curve.
• The tangent vector formed by the first two control points ${\mathrm{P}}_0$ and ${\mathrm{P}}_1$ $\left(\overline{{\mathrm{P}}_0{\mathrm{P}}_1}\right)$ defines the initial direction, and the tangent vector formed by the last two points $\overline{{\mathrm{P}}_{\mathrm{n}-1}{\mathrm{P}}_{\mathrm{n}}}$ defines the ending direction of the curve.
• The Bézier curve lies within the area that is formed by the control points [18,19,20].

Furthermore, an important property of Bézier curves is the straightforward calculation of its derivative in terms of the variable $t$. If the Bézier curve for x coordinates is $x\left(t\right)=B_{\mathrm{x}}\left(t\right)$ and for y coordinates is $y\left(t\right)=B_{\mathrm{y}}\left(t\right)$, then the first derivative $\mathrm{d}y/\mathrm{d}x$ is calculated by Equation (13):

$$\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}$$

(13)

and the second derivative ${\mathrm{d}}^{2}y/\mathrm{d}{x}^2$ by Equation (14):

$$\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{x}^2}=\frac{\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathrm{d}y/\mathrm{d}\mathrm{x}\right)}{\mathrm{d}x/\mathrm{d}t}$$

(14)

The equations for the first and second derivatives allow for the calculation of the curvature of Bézier curves [18,19,39]. Bézier curves offer the possibility to form symmetric curves. Symmetric Bézier curves are generated by placing pairwise control points into symmetrical positions [40].

2.2.2. Impact of the Control Point’s Placement on the Bézier Curve

In this subsection, we demonstrate the impact of the positioning of control points on the shape of a cubic (i.e., order of $n=3$) Bézier curve. As presented in Figure 2, Bézier curves enable the formation of various shapes—continuous single-valued curves, as shown in Figure 2a; self-intersecting curves, as presented in Figure 2b,c; and closed shapes, as shown in Figure 2d. Figure 3 shows the various possibilities regarding the change in shape of a cubic Bézier curve from Figure 2a only by changing the position of the control points ${\mathrm{P}}_1$ and ${\mathrm{P}}_2$ (the first ${\mathrm{P}}_0$ and last ${\mathrm{P}}_3$ control points have fixed positions).

We chose to fix the first and last points for the presented comparison to show the variety of curves that start and end in the same points. As explained in Section 2.2.1, the starting and ending directions of the curve were defined by the tangents $\overline{{\mathrm{P}}_0{\mathrm{P}}_1}$ and $\overline{{\mathrm{P}}_2{\mathrm{P}}_3}$, respectively. The cases shown in Figure 3a,b are examples of changing the slope and curvature by moving the control points but preserving the property analogous to a single-valued function. In the case presented in Figure 3c, when the control point ${\mathrm{P}}_1$ was below control point ${\mathrm{P}}_0$ and ${\mathrm{P}}_2$ was above ${\mathrm{P}}_3$, we obtained an overshoot over the normalized value of 1. That behavior is not adequate for the description of magnetization curves. Another nonphysical behavior was observed when the absolute values of the x coordinates of points ${\mathrm{P}}_1$ and ${\mathrm{P}}_2$ were higher than the x values of the last and first control points, respectively, as shown in Figure 3d. The obtained curve was not single valued in respect to x and thus not appropriate for the description of magnetization curves. Control points are decisive for the curvature and shape of the final curve [26].

3. Methodology

3.1. Assumptions and Limitations

The discussed analytical functions in Section 2.1 have limited ability to approximate various curves with high accuracy because they have only three or four adjustable parameters. In contrast to this, the Bézier curves can be based on many theoretically arbitrary control points; therefore, they offer higher flexibility for the approximation of measured data.

The focus of the presented analysis was on a systematic analysis of which functions and curves would result in the best fit based on different sets of measured data. Further, the extrapolation capabilities of the presented functions and curves were analyzed. During the analyses, we assumed that different sets of measured data for two NO steels were available. We limited the analyses to the application of the presented analytical functions and classical Bézier curves of different orders. One of the goals was to identify the most adequate orders for approximation of the discussed curves.

3.2. Measured Data

The analysis was performed systematically based on measured characteristic magnetization curves $J\left(H\right)$, where $J$ is the magnetic polarization and $H$ is the magnetic field strength. In the analysis, we considered three of the most basic curves that describe the nonlinear properties of soft magnetic materials:

(1)

First magnetization curves;

(2)

Anhysteretic curves;

(3)

Descending branches of the major loops.

We determined the discussed curves for two NO electrical steels with different thicknesses, i.e., 0.27 mm (NO27) and 0.35 mm (NO35), which are both commonly used in contemporary electrical machines. The measurement setup was based on a single sheet tester within a computer-aided measurement setup, in accordance with the International Standard IEC 60404-3. The measurements were performed at low excitation frequencies (i.e., under the so-called quasi-static conditions) up to $H=50 \mathrm{k}\mathrm{A}/\mathrm{m}$ at 5000 $J$-equidistant measurement points. Based on the measured curves, it was assumed that $H=H_{\mathrm{s}}=50 \mathrm{k}\mathrm{A}/\mathrm{m}$ was the saturation point of both materials. Consequently, the available measured dataset was ${\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}=\left[-H_{\mathrm{s}},H_{\mathrm{s}}\right]$, where $\left[-H_{\mathrm{s}},H_{\mathrm{s}}\right]$ denotes the interval of the discrete measured dataset, starting from $-H_{\mathrm{s}}$ and ending with $H_{\mathrm{s}}$.

The anhysteretic curve was calculated as the average value of polarization $J$ between the ascending and descending branches of the major loop.

3.3. Definition of the Data Subsets and Evaluation (Sub)Regions

All the discussed functions and curves were evaluated systematically on different input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ of the available measured data ${\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$ (i.e., ${\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}=\left[-H_{\mathrm{s}},H_{\mathrm{s}}\right]$). The input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ was, in all cases, a subinterval of the measured data ${\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$, i.e., ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}\subseteq {\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$. The evaluation subregions ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$ were defined within the discussed input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$. Special attention was paid to evaluation of the following characteristic subregions:

(1)

The high-permeability (HP) subregion (${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},2}=\left[-H_1,H_1\right]$);

(2)

The saturation (SAT) subregion (${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},3}=\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},{-H}_1\right]\cup \left[H_1,H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$);

(3)

The extrapolation (EXT) subregion (${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},6}=\left[{{-H}_{\mathrm{s}},-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]\cup \left[H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{s}}\right]$).

The HP and SAT subregions were important for the first part of the analysis, i.e., evaluation of the approximation capabilities, whereas the EXT subregion was crucial for the analysis of the extrapolation capabilities. These subregions are very important in applied engineering and are presented schematically in Figure 4.

The input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ used in the process of fitting analytical functions and Bézier curves is plotted in blue in Figure 4. The measured data between $\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$ represented the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ in individual cases for approximation, as presented in Figure 4a. It was slightly different for the case of extrapolation, where, additionally, the saturation points $H_{\mathrm{s}}$ were added to both sides of the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$, highlighted with the blue cross in Figure 4b. The characteristic subregions of the magnetization curves were identified in both cases. Considering those subregions, the measured data were divided accordingly into evaluation subsets, denoted with ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{1}}$ to ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{6}}$ and presented in Figure 4a,b. This approach and division enabled us to perform a systematic analysis via local evaluation of the goodness of fit of the approximated curves in specific subregions individually as well as a global evaluation on the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$.

It is important to note that $H_{\mathrm{s}}$ and $H_1$ were fixed at $H_{\mathrm{s}}=50 \mathrm{k}\mathrm{A}/\mathrm{m}$ and $H_1=1 \mathrm{k}\mathrm{A}/\mathrm{m}$, respectively. In contrast to this, $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$ was varied (i.e., 1, 2, 5, 10, 15, 20, 30, 40, and 50 kA/m) and determined various input datasets for approximation and extrapolation of the discussed functions and curves.

3.4. Curve Fitting

The parameters ($a$, $b$, $c$, and $d$) of the analytic function from Table 1 were fitted to the input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ of the measured magnetization curves using the Matlab function fit by applying the nonlinear least squares method.

For the process (i.e., the placement of control points) of fitting Bézier curves to the input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$, we applied the differential evolution (DE) algorithm [10,41]. We set the value of the mutation factor to 0.5 and the crossover factor to 0.7 and used the DE algorithm for 15,000 iterations for each curve. In each iteration, the x and y coordinates of the control points were calculated to minimize the deviation of the calculated Bézier and measured magnetization curves. Within the process, the first and last control points (which were assumed to be equal to the first and last points, respectively, of the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$) were fixed according to the adequate position of the measured curve. The goal was to determine the intermediate control points. To find the optimum values, we minimized the objective function $F$ in Equation (15):

$$F={\sum }_{i=1}^{N_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\left[J_{\mathrm{i}\mathrm{n},\mathrm{i}}\left(H_{\mathrm{i}\mathrm{n},\mathrm{i}}\right)-J_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c},\mathrm{i}}\left(H_{\mathrm{i}\mathrm{n},\mathrm{i}}\right)\right]}^2,$$

(15)

where $J_{\mathrm{i}\mathrm{n},\mathrm{i}}\left(H_{\mathrm{i}\mathrm{n},\mathrm{i}}\right)$ are the input measured values of the magnetic polarization at the corresponding magnetic field strengths $H_{\mathrm{i}\mathrm{n},\mathrm{i}}$, $J_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c},\mathrm{i}}\left(H_{\mathrm{i}\mathrm{n},\mathrm{i}}\right)$ is the magnetic polarization approximated by a Bézier curve at $H_{\mathrm{i}\mathrm{n},\mathrm{i}}$, and $N_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the number of elements of the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$.

In general, $2·\left(n-1\right)$ parameters had to be determined for an arbitrary order $n$ of a Bézier curve. Therefore, higher orders of Bézier curves naturally led to longer computational effort in the fitting procedure.

3.5. Evaluation of the Goodness of Fit

To evaluate the goodness of fit for the calculated curves, we applied the statistical measure normalized root mean square (NRMS) $\epsilon$, defined by Equation (16) [7].

$$\epsilon =\sqrt{\frac{1}{N_{\mathrm{m}\mathrm{a}\mathrm{x}}}\sum _{i=1}^{N_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\left(\frac{J_{\mathrm{i}\mathrm{n},\mathrm{i}}\left(H_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{i}}\right)-J_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c},\mathrm{i}}\left(H_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{i}}\right)}{{\Delta J}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}}\right)}^2}$$

(16)

In Equation (16), $J_{\mathrm{i}\mathrm{n},\mathrm{i}}\left(H_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{i}}\right)$ are the measured values of magnetic polarization at corresponding magnetic field strengths $H_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{i}}$ in the observed subregion for the evaluation of ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$. Furthermore, $J_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c},\mathrm{i}}\left(H_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{i}}\right)$ are the corresponding calculated values of magnetic polarization, and ${\Delta J}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}$ is the range of the evaluated subregion ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$. This range ${\Delta J}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}$ is defined as the difference between the maximum $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}$ and minimum $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{i}\mathrm{n}}$ values of such a subregion, i.e., ${\Delta J}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}=J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}-J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{i}\mathrm{n}}$.

During the analysis, we evaluated both the approximation and the extrapolation capabilities of the presented functions and curves. Individual curves were evaluated by Equation (16) in different subregions individually (i.e., HP, SAT, and EXT regions) by selecting the adequate evaluation region ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$.

3.5.1. Approximation Capabilities

We analyzed the goodness of fit and evaluated the approximation capabilities of the analytic functions and Bézier curves of different orders in specific subregions within the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$. The available measured data ${\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$ were divided into nine input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ between $\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, where $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$ equaled 1, 2, 5, 10, 15, 20, 30, 40, and 50 kA/m, respectively. An exemplary case of an input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}},$ where $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=30\mathrm{k}\mathrm{A}/\mathrm{m},$ is presented in Figure 4a with the blue curve. The first magnetization curves were fitted on the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}=\left[0,H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, whereas all the other curves were fitted on a symmetric subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}=\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right].$ In the first step, individual analytic functions were fitted, and Bézier curves of different orders were obtained for each individual input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$. The evaluation of the goodness of fit $\epsilon$ was carried out in three subregions of each subset, where the evaluation was performed in the following ranges ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$ (${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}\subseteq {\mathbf{H}}_{\mathbf{i}\mathbf{n}}$):

(a)

The input subset ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{1}}={\mathbf{H}}_{\mathrm{i}\mathrm{n}}=\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$;

(b)

The HP subregion (within the input subset), i.e., ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{2}}=\left[{-H}_1,H_1\right]$;

(c)

The SAT region (within the input subset), i.e., ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{3}}=\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},{-H}_1\right]\cup \left[H_1,H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$.

The evaluation subregions ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{1}}$ to ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{3}}$ are depicted graphically in Figure 4a. It is important to note that $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}$ and $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{i}\mathrm{n}}$ for case (c) were determined at $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$ and $H_1$, respectively.

The result of this analysis was the identification of the most suitable analytic functions and the most adequate orders of Bézier curves for the approximation of nonlinear magnetic properties for different input data subsets. These were included in the next analysis, i.e., the analysis of the extrapolation capabilities.

3.5.2. Extrapolation Capabilities

In the analysis of the extrapolation capabilities, the fitting was carried out analogously with the previous analysis. The difference was that the positive and negative saturation points $J_{\mathrm{s}}\left(H_{\mathrm{s}}\right)$ were added to all symmetric input subsets between $\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$. Therefore, the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ consisted of the negative saturation point $-H_{\mathrm{s}}$, the measured data between $\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, and the positive saturation point $H_{\mathrm{s}}$, i.e., ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}=-H_{\mathrm{s}}\cup \left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]\cup H_{\mathrm{s}}$. This input subset for a specific case, where $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=30 \mathrm{k}\mathrm{A}/\mathrm{m},$ is presented with the blue line and blue crosses in Figure 4b. The result of the fitting process was analytic functions and Bézier curves that ranged from $-H_{\mathrm{s}}$ to $H_{\mathrm{s}}$. They also enabled the analysis of the goodness of fit in the assumed EXT subregion. The evaluation of the goodness of fit $\epsilon$ was carried out again in three evaluation $\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$ ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$, where ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}\subseteq {\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$:

(d)

The full measured dataset up to saturation ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{4}}={\mathbf{H}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{s}}$;

(e)

The subregion which contained the measured data ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{5}}=\left[{-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$;

(f)

The EXT subregion (${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{6}}=\left[{{-H}_{\mathrm{s}},-H}_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]\cup \left[H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{s}}\right]$).

Analogous to the previous subsection, $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}$ and $J_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l},\mathrm{m}\mathrm{i}\mathrm{n}}$ for case (f) were determined at $H_{\mathrm{s}}$ and $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$, respectively.

4. Results

The available range (defined by $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$) of the (measured) input data ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}=\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$ has, in general, a high impact on the goodness of fit of the approximated functions and curves. The impact of the range on the approximated curves using both analytic functions and Bézier curves is presented in Figure 5. The curves in Figure 5a,b were approximated on the input subset limited by $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=0.5$ $\mathrm{k}\mathrm{A}/\mathrm{m}$ and in Figure 5c,d on the input subset limited by $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=15$ $\mathrm{k}\mathrm{A}/\mathrm{m}$. The curves approximated on a bigger subset ($H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=15$ $\mathrm{k}\mathrm{A}/\mathrm{m}$) deviated more in the HP region in comparison to those approximated on the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}},$ which coincided with the observed region ($H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=0.5$ $\mathrm{k}\mathrm{A}/\mathrm{m}$). The approximated curves in Figure 5 confirm that the size of the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ (i.e., the range of measured input data defined by $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$) had a high impact on the approximation accuracy.

4.1. Approximation with Analytic Functions

In general, all analytic functions produced the overall shape to describe the discussed nonlinear relationships, as demonstrated, e.g., in Figure 5a,c. For the sigmoid functions (Functions (1)–(5)) similar conclusions are presented in [10,15].

4.1.1. Anhysteretic Curves

To determine the most adequate analytic function for the approximation of the anhysteretic curve, we calculated the NRMS deviations $\epsilon$ on different regions within different input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ for all the discussed functions, as described in Section 3.5.1. The obtained results are presented in Figure 6.

First, we evaluated the global goodness of fit (i.e., within different symmetric input subsets ${{\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{1}}=\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ of the measured data) systematically, as described in case (a) of Section 3.5.1. Based on the presented results in Figure 6a,b, it was observed that the Langevin function had the best overall fit for the smaller input subsets up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=5 \mathrm{k}\mathrm{A}/\mathrm{m}$. Further, the Elliot function had the best goodness of fit for the measured anhysteretic curve of both materials when a symmetric input subset from $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=5 \mathrm{k}\mathrm{A}/\mathrm{m}$ up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=10 \mathrm{k}\mathrm{A}/\mathrm{m}$ was applied. These results agree with the findings in [10]. The goodness of fit with the Langevin function was slightly worse compared to the one evaluated with the Elliot function in the latter region. However, if the input subsets exceeded $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=10$ $\mathrm{k}\mathrm{A}/\mathrm{m}$, the global goodness of fit of the inverse hyperbolic sine function was the best.

To support the obtained results further, we analyzed the deviations in two specific subregions of the anhysteretic curve. The first local goodness of fit was defined between ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{2}}=\left[-H_1,H_1\right]$, where $H_1=1 \mathrm{k}\mathrm{A}/\mathrm{m}$ (the so-called HP region), as explained in case (b) of Section 3.5.1. The local goodness of fit for the HP region is shown in Figure 6c,d. The obtained results supported the premise that, as the input subset ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ increased into the saturation region, the deviations were more significant. This was expected due to the high nonlinearity of the input data, whereas individual subregions were not considered separately, and, consequently, could not be approximated with equal accuracy. In the HP region, the best results were again obtained using the Elliot and Langevin functions, whereas in the cases of larger input subsets above $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=20 \mathrm{k}\mathrm{A}/\mathrm{m},$ the inverse hyperbolic sine function again became the best option. If limited measured data not deep into saturation are available, the Elliot and Langevin functions are the best choice for approximation.

The second local subregion ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{3}}$ (SAT region) was defined from $1 \mathrm{k}\mathrm{A}/\mathrm{m}$ to the maximum value in the individual input subset $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$, i.e., when the materials were significantly saturated (i.e., case (c) in Section 3.5.1). The values of $\epsilon$ within this analysis are presented in Figure 6e,f. In this region the inverse hyperbolic sine function had superior results to other analytic functions, especially for input subsets that were measured deeper into saturation.

4.1.2. Major Loop Curves

Next, we analyzed the goodness of fit of the approximated analytic functions on the descending branch of the major loop. The calculated NRMS values $\epsilon$ are presented in Figure 7. The results obtained for the major loop are in line with those obtained for the anhysteretic curve. This could be attributed to the fact that the anhysteretic curve was calculated as the average value between both branches of the major loop, and, consequently, closely related to the major loop. It was concluded that the offset in $H$ of the major loop curves did not introduce a significant change in the results.

4.1.3. First Magnetization Curves

The results obtained were slightly different in the case of the first magnetization curves. The calculated NRMS errors $\epsilon$ for the case of the first magnetization curve are shown in Figure 8. In the input subset between 0 and $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$, i.e., ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}=\left[0,H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, the Elliot and Langevin functions gave best results in the case of sheet NO35, as shown in Figure 8a. In the case of NO27, presented in Figure 8b, the Elliot and Langevin functions gave the best results up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=30$ $\mathrm{k}\mathrm{A}/\mathrm{m}$. Beyond this subregion, the inverse hyperbolic sine function became the most adequate analytic function. In the HP region, the Langevin function had the lowest deviation from the measured data in the case of NO35, as presented in Figure 8c. Figure 8d shows the results for NO27, where the Langevin function gave better results for the intermediate fields (up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=20$ $\mathrm{k}\mathrm{A}/\mathrm{m}$), and, after that region, the Langevin and Elliot functions had similar goodness of fit. In the case of the SAT region (i.e., from 1 kA/m to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$), the results for both materials in Figure 8e,f were similar; overall, the inverse hyperbolic sine function yielded the best results. The Elliot and Langevin functions can be used for lower fields of up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=10$ $\mathrm{k}\mathrm{A}/\mathrm{m}$.

4.2. Approximation with Higher-Order Bézier Curves

Next, we analyzed the ability of Bézier curves of different orders systematically to describe the measured magnetization curves. We performed the analysis starting with Bézier curves of the order $n=3$ (order three is commonly used in many applications, e.g., in splines [42]) and increased the order up to $n=12$.

Based on the results that are presented in Figure 5b,d we concluded that lower-order Bézier curves (e.g., the third and fourth orders) had limited ability to approximate magnetization curves adequately. Hence, we excluded these orders from further analysis. We applied the analysis presented in cases (a)–(c) in Section 3.5.1 to the approximation with Bézier curves.

4.2.1. Anhysteretic Curves

The obtained values of $\epsilon$ for the anhysteretic curves are presented in Figure 9. This part of the analysis was performed with Bézier curves from the 5th to 12th orders. For the sake of clarity, we did not plot $\epsilon$ for the Bézier curves of all the discussed orders but only representative ones. From the $\epsilon$ obtained for the individual symmetrical input subsets ${{\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{1}}=\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ (as explained in Section 3.5.1 case (a)) in Figure 9a,b, we concluded that the most suitable were the eighth- and ninth-order Bézier curves. In general, it is reasonable to select a curve with the lowest number of control points. The obtained results showed that all Bézier curves starting from order six already had significantly improved goodness of fit compared to the analytic functions in all the analyzed scenarios.

For the HP region (case (b) in Section 3.5.1) in Figure 9c,d, the ninth-order Bézier curve offered the best approximation, which is in line with the results in Figure 5b,d. Along with the 9th order, the 7th-, 8th-, and 10th-order Bézier curves yielded a very good approximation of this region. In this case, we observed that, as ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ increased, the fit in the HP region became less accurate, which was analogous with the results for the analytic functions. The same orders of Bézier curves were estimated as the best for the remaining SAT region of the curve (case (c) in Section 3.5.1) presented in Figure 9e,f.

4.2.2. Major Loop Curves

The results for the descending branch of the major loop are similar to the results obtained for the anhysteretic curve. The results are presented in Figure 10. The ninth-order Bézier curve had the lowest deviation from the measured major loop in all cases. The results obtained with the ninth-order curve were followed closely by the eighth-order Bézier curve.

4.2.3. First Magnetization Curves

In the case of approximating the first magnetization curve, Bézier curves were superior to analytic functions, whereas the obtained NRMS values $\epsilon$ were significantly lower. The calculated NRMS deviation values $\epsilon$ for the Bézier curves are presented in Figure 11. For this specific curve, even the fourth-order Bézier curve demonstrated low deviation from the measured data. For input subsets up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=30 \mathrm{k}\mathrm{A}/\mathrm{m},$ shown in Figure 11a,b, the sixth- and seventh-order Bézier curves yielded the best goodness of fit, and, above 30 kA/m, the eighth- and ninth-order Bézier curves yielded the best results. In the HP region, presented in Figure 11c,d, the seventh-order Bézier curve had the best overall goodness of fit. Those results are followed closely by the results obtained with the eighth-order Bézier curve. In the SAT region, presented in Figure 11e,f, the seventh-order Bézier curve had the best goodness of fit. Additionally, in the region with magnetic field strength values of up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=30 \mathrm{k}\mathrm{A}/\mathrm{m},$ the sixth-order Bézier curves demonstrated low deviation from the measured data. In regions with high magnetic field strength values beyond 30 kA/m, the eighth- and ninth-order Bézier curves had the lowest NRMS deviation values $\epsilon$.

Based on the results in this subsection, we concluded that suitable orders of Bézier curves for approximation of the discussed nonlinear magnetic properties were the sixth to ninth.

4.3. Analysis of Extrapolation Capabilities

It is often not feasible to measure the hysteresis properties deep into saturation. Therefore, only measured magnetization curves in limited regions are available, as presented schematically in Figure 12a. The reasons for this are the highly nonlinear properties and the corresponding complexity and challenges in the measurement procedures. The measurement accuracy is decreased significantly when measurements are performed at values of high magnetic field strength $H$. Because of that, it is important to have a reliable method to determine the magnetic curves up to high saturation values, even if measured data are not available for the whole magnetization curve. This can be performed by extrapolation, which is the technique used most often for this purpose in applied engineering. Two of the most used functions for this purpose are the Langevin and the Elliot functions (the latter is often referred to as the Frölich–Kennelly extrapolation).

In Figure 12b, an example of limited available measured data is presented schematically. In such cases, a significant part of the magnetization curve is unknown and must be estimated with extrapolation. In the literature, several methods have been proposed for the extrapolation of magnetization curves [43]. Bézier curves could offer an additional reliable and good estimation of the extrapolated region of the magnetization curves if the saturation point is known. For that reason, we analyzed the capabilities of Bézier curves for the extrapolation of measured magnetic curve data and compared the results with the analytic functions. In this analysis, we included Bézier curves of the sixth to ninth orders, as well as the representative analytic functions, i.e., the Elliot, Langevin, and inverse hyperbolic sine functions. These were selected based on the results in Section 4.1 and Section 4.2.

The analysis was performed using input subsets ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$, and we assumed that the saturation point was known, i.e., we fixed the last point as the saturation point $J_{\mathrm{s}}\left(H_{\mathrm{s}}\right)$. Then, we shortened the input subset gradually. In this way, we could analyze the influence of the measured input subset’s size on the extrapolated part of the curve. The methodology is presented in Section 3.5.2 in cases (d)–(f). The measured and extrapolated curves of the NO27 electrical steel, where the input data were limited to the region between $\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, where $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=15 \mathrm{k}\mathrm{A}/\mathrm{m},$ are shown in Figure 13.

4.3.1. Anhysteretic Curves

First, we analyzed the goodness of fit for the full measured dataset of the anhysteretic curve, i.e., ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{4}}={\mathbf{H}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$ (case (d) in Section 3.5.2). The goodness of fit for these regions is presented in Figure 14a,b. The eighth- and ninth-order Bézier curves had the best goodness of fit in the whole data range, followed closely by the seventh-order and sixth-order Bézier curves. The NRMS deviation $\epsilon$ became constant when the input subset reached a certain range. This range was approximately $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=20 \mathrm{k}\mathrm{A}/\mathrm{m}$ for both discussed materials. As the ${\mathbf{H}}_{\mathbf{i}\mathbf{n}}$ increased (and the extrapolation region decreased), $\epsilon$ decreased, which supports the conclusion from the previous analysis.

In the next step, we evaluated the deviations in two separate subregions of the approximated curves. The first was the subregion of the input data ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{5}}=\left[-H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}},H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}\right]$, as described in Section 3.5.2 (e). The second was the EXT subregion, as presented in Section 3.5.2 (f). The goodness of fit for ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{5}}$ is presented in Figure 14c,d, and for the EXT subregion ${\mathbf{H}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l},\mathbf{6}}$ in Figure 14e,f. For both regions, the eighth- and ninth-order Bézier curves generated the lowest NRMS error $\epsilon$ and offered the best estimation of the measured curve. In both subregions, the inverse hyperbolic sine function had the lowest deviation from the measured data out of the chosen analytic functions and was comparable to the results obtained with the sixth-order Bézier curve.

4.3.2. Major Loop Curves

Figure 15 shows the calculated NRMS values $\epsilon$ obtained for the descending branch of the major loop. The conclusions are analogous to the analysis of the extrapolation of the anhysteretic curve. Overall, the eighth- and ninth-order Bézier curves offered the best goodness of fit, followed by the seventh-order Bézier curve. The sixth-order Bézier curve and the inverse hyperbolic sine function again had very similar results.

4.3.3. First Magnetization Curves

The results of the analysis of the first magnetization curve are presented in Figure 16. The Bézier curves again had superior results over the analytic functions in the case of the first magnetization curve. All orders from six to nine were suitable for the description of the first magnetization curve in all subsets and subregions in individual subsets. The results for the full region up to saturation are presented in Figure 16a,b, where the Bézier curves of orders seven to nine offered the best results. Slightly worse results were obtained using the sixth-order curve. In the region with the measured data (i.e., $\left[0,H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}} \right]$) in Figure 16c,d, the results were in line with the full region. The extrapolation region in Figure 16e,f was the only region where one of the analytical functions, namely, the inverse hyperbolic sine function, came close to the accuracy of the Bézier curves. Among the Bézier curves, in the low-field-strength input datasets of up to $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}=5 \mathrm{k}\mathrm{A}/\mathrm{m},$ the sixth-order Bézier curve had the lowest deviation, whereas in the more extended datasets, the seventh- to ninth-order Bézier curves again offered the best goodness of fit.

5. Discussion

Based on the results obtained within the performed analyses, the common takeaway is that both approaches, i.e., using adequate analytic functions or adequate Bézier curves, are suitable to describe or filter noisy measured data or extrapolate the nonlinear magnetic properties. Overall, Bézier curves of adequate orders (in most cases, between six and nine) offered significantly lower deviations from the measured curves compared with the analytic functions. The advantages and disadvantages of both approaches will be highlighted in this section.

Bézier curves are less sensitive to the size and range of the input dataset (limited by $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$) than analytic functions. This is especially highlighted in Figure 9 and Figure 10, where it is shown that the Bézier curves had almost constant values of the NRMS deviation $\epsilon$, whereas the Elliot function did not reach a constant value at different $H_{\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}$. Further, as the order of Bézier curves increased (at a fixed maximum number of iterations of DE), starting from order nine, the goodness of fit became worse, which can be observed in Figure 9, Figure 10 and Figure 11 (order 12 had a higher NRMS deviation $\epsilon$ than order 9). Therefore, not all higher orders are adequate for approximating or extrapolating the measured subsets.

A downside of analytic functions is that their final point does not reach the saturation value. This property can be observed in Figure 13. They are limited by their mathematical description. It is possible to force the functions to go through the saturation point exactly but with the tradeoff of increasing the deviations in other subregions of the curve significantly, to the extent that the description in them is useless. A special and limiting feature of Bézier curves is that their starting points and endpoints can be fixed exactly without losing the agreement with the measured data in any region (i.e., the first and last control points are interpolated, whereas the intermediate points are used to manipulate the shape and curvature of the curve). Furthermore, Bézier curves allow the slope at which the curve approaches saturation to be prescribed, as explained in Section 2.2.1. This can be both an advantage and a disadvantage. For example, a very accurate extrapolation is possible only when the saturation point of the material is known (e.g., can be calculated or estimated theoretically).

The straightforward calculation of derivatives in any point of the curve and, consequently, the curvature of the curve, allows for the precise determination of the shape of a curve under construction. Measured magnetic curves can have a lot of noise in the HP region as the maximum value of magnetic field density increases. Bézier curves have the potential to make the extraction of measured data easier. The user could set the theoretically known saturation value, the slope at which the saturation is reached, and the curvature in individual regions to obtain a symmetrical curve and, consequently, a symmetrical loop. In this way, the potential non-physical behavior of the obtained curves can be avoided (similar to the behavior of the curves in Figure 5b,d). In this work, the only prescribed limitation was not to exceed the value of magnetic polarization at saturation $J_{\mathrm{s}}\left(H_{\mathrm{s}}\right)$, i.e., the control points could have been placed anywhere below that point (and above the negative saturation $-J_{\mathrm{s}}\left(-H_{\mathrm{s}}\right)$). Additionally, pairwise symmetrical control points can be prescribed during the DE calculation to obtain symmetrical Bézier curves [40].

6. Conclusions

From the obtained results, we concluded that Bézier curves are promising candidates to adequately model the nonlinear properties of NO soft magnetic materials. Through the flexibility of curve order and adequate placement of control points, Bézier curves can approximate the shape of the magnetization curves with significantly lower deviation compared to analytic functions.

With the performed analysis, we concluded that Bézier curves from order six and higher are suitable for the approximation of the discussed magnetization curves. The approximation with those orders offers significantly improved fit based on the measured data compared with approximations with all the discussed analytic functions. In our specific case, the ninth order Bézier curve offered the best goodness of fit in all the analyzed subregions.

Bézier curves can approximate the high saturation region accurately (i.e., the extrapolation region). However, the theoretical point of saturation is required, in addition to the measurements in the HP subregion. To achieve high accuracy, it would be necessary to prescribe further boundary conditions at which the curve is generated. Examples are prescribing the slope for reaching the saturation region and adjusting the smoothness of the curvature along individual segments of the magnetization curve while calculating the control points. Our future work will focus on analyzing the conditions to generate symmetrical and accurate magnetization curves.