Comparison of a quasi Newton method using Broyden’s update formula and an adjoint method for determining local magnetic material properties of electrical steel sheets

The local magnetic properties of electrical steel sheets are critical for predicting the performance in applications such as transformers, motors and generators (Li et al., 2017). These properties include parameters like magnetic permeability, saturation magnetization and hysteresis characteristics, which vary locally across the material due to the influence of different manufacturing processes. In particular, cutting techniques are an essential part in the manufacturing chain to create the shape and geometry of the electrical steel sheets needed in different applications. As a result of these processes, changes in the microstructure and residual stresses at the cutting edges may occur, leading to a deterioration of the magnetic material properties of the electrical sheets. The extent of this deterioration at the cutting edge depends on the cutting process, e.g. punching, laser cutting or water-jet cutting, and the cutting parameters, e.g. laser intensity, blade sharpness or cutting speed (Schoppa et al., 2000; Hofmann et al., 2015; Sundaria et al., 2019).

Accurate and efficient determination of the local properties is a challenging task and remains a subject of active exploration within the scientific community. One widely used approach for assessing the impact of cutting edges involves the examination of the ratio between the cutting length and the overall bulk material. This is accomplished by dividing an electrical steel sheet into multiple narrower strips, such that the combined width of these strips, when placed side by side, mirrors the dimensions of the original sheet. Modifying the strip width yields distinct combinations of cutting length to bulk material ratio, which are subsequently evaluated using a single sheet tester (SST) or an Epstein frame (Sundaria et al., 2020; Bali et al., 2017).

A destructive and a nondestructive method for locally measuring the influence of cutting edges are presented in earlier studies (Nakata et al., 1992; Loisos and Moses, 2005; Lewis et al., 2018; Gmyrek, 2016). The destructive method includes the method of placing search coils near the cutting edge by drilling holes in sheets to measure the magnetic flux. The nondestructive method uses the needle probe method.

Our approach involves the integration of various scientific methodologies to comprehensively investigate the influence of cutting edges on the magnetic material properties. Therefore, a combination of measurement, numerical simulation and inverse modeling techniques are used. To gather measurement data, a sensor-actuator (SA) system is used to magnetize electrical steel sheets locally and measure the magnetic field above the sample. In addition to experimental measurements, numerical data are generated by solving the magneto-static problem employing the finite element (FE) method. Therefore, an appropriate model for the SA system and electrical steel sheets is used. The combination of measured and numerical data forms the basis for applying inverse schemes to determine the parameters of the defined material model, including the degradation of the properties due to cutting.

In this work, the ideas of Gschwentner et al. (2023) are taken up and extended by introducing the adjoint method for solving the inverse problem. In Section 2, a brief summary of the sensor-actuator system as well as the electrical sheet model is given. The inverse problem can be written as a minimization problem, containing the error term between the measured and simulated data as well as a regularization term. The minimization problem is then solved via a quasi Newton method (QNM) using Broyden’s update formula to approximate the Jacobian and via the adjoint method. A detailed description of the methods is presented in Section 3. In Section 4, both methods are tested for differently pronounced cutting edge effects. The goal of this work is to compare the QNM and the adjoint method regarding accuracy and efficiency.

The sensor-actuator system depicted in Figure 1 represents an assembly of stacked iron sheets, which is subjected to excitation from two coils. This system possesses the capability to magnetize electrical steel sheets locally while concurrently measuring the local magnetic flux density. The measurement process involves employing a sensor array equipped with S Hall and/or GMR sensors, which can accurately detect the x−, y− and z− components of the magnetic field above the electrical steel sheets. Two electrical steel sheets, denoted as Sample 1 and Sample 2, are placed in proximity along the cutting edge. This arrangement is necessary to ensure magnetization of the cutting edges with the sensor-actuator system and thus achieve a corresponding sensitivity for the inverse scheme (Gschwentner et al., 2023). Under the assumption that both samples originate from the same batch and identical cutting process parameters are maintained, it is valid to consider symmetrical and identical material behavior. To gather measurement data, Sample 1 and Sample 2 are measured at various positions along the x-direction with P being the number of positions. Notably, due to the assumed large variations in material behavior near the cutting edges, a higher density of measurement positions is concentrated in this region compared to the bulk material. Furthermore, to maintain consistency, the sheets are demagnetized between each measurement position to eliminate any residual magnetism, which is an essential assumption for subsequent numerical simulations. The resulting data set contains measurements for the three magnetic field components at each sensor and measurement position, denoted as $Bx,i,jmeas$⁠, $By,i,jmeas$ and $Bz,i,jmeas$ with i = 1, 2,...., S the sensor positions and j = 1, 2,....., P the measurement positions. Subsequently, based on this data set, the magnetic field density amplitude is computed using the Euclidean norm, represented as $||Bi,jmeas||2=(Bx,i,jmeas)2+(By,i,jmeas)2+(Bz,i,jmeas)2$⁠, leading to the final measurement data $Bmeas=(B1,1meas,B2,1meas,…,BS,Pmeas)T$⁠.

The degradation of material characteristics due to cutting becomes noticeable within a narrow span of millimeters near the cutting edges. The extent of this degradation depends on the particular cutting technique and the parameters used in the cutting process. To accurately model the significant material changes that occur in this region during the simulations, each electrical steel sheet is divided into M nonequidistant subdomains, denoted by Ω_m. The size of these subdomains is significantly smaller in the immediate vicinity of the cutting edges compared to the bulk material (Figure 2).

In the numerical simulation, the electrical steel sheets are described by a material model, which can encompass linear, nonlinear or hysteretic characteristics. The specific material model involves a variable number of parameters that must be determined to align the model’s behavior with the actual material properties. In this work, a linear material model v = v_r v₀, with v_r the relative reluctivity and v₀ the reluctivity of vacuum, is assumed. To take the influence of cutting edges into account, the chosen linear material model is assigned to each subdomain, allowing independent selection of model parameters for each subdomain. Consequently, the searched-for parameter vector can be expressed as p = [v_r,1, v_r,2,…, v_r,M]^T. The advantage of this approach is that no adaptation of the material model is necessary to take into account factors influencing the magnetic material behavior, e.g. residual stresses, microstructure, etc., as these are inherently included in the model parameter for each subdomain. Furthermore, this approach can be extended to applications that also lead to a change in magnetic material behavior, e.g. forming or heat-treatment of electrical steel sheets.

The inverse scheme calculates the searched-for parameter vector p based on the measured magnetic flux densities B^meas and the simulated magnetic flux densities B^sim. Therefore, a nonlinear least squares minimization problem has to be solved to find the optimal parameter p^opt, such that the error norm between B^meas and B^sim is minimized. Due to the inevitable measurement noise in the data, difficulties in solving the nonlinear least squares problem occur. More precisely, small perturbations in the measurement data have a pronounced negative effect on the computed parameters and cause the solution strategy to diverge. From a mathematical point of view, this can be stated as an ill-posed problem. To overcome this problem, a Tikhonov regularization (Tikhonov et al., 1995) is applied to ensure convergence. In doing so, the minimization problem reads as follows:

with $Fi(xj,u(p),p)=Bisim(xj,p)−Bimeas(xj)$⁠, P the number of measurement positions, S the number of sensors position, $Bisim(xj,p)$ the simulated magnetic flux density, $Bimeas(xj)$ the measured magnetic flux density, α the regularization parameter, p = [v_r,1,v_r,2,…,v_r,M]^T the relative magnetic reluctivities in each subdomain of the steel sheet, u the magnetic vector potential and J the electric current density. Finding the optimal parameter p^opt of the minimization problem, equation (1) is solved iteratively, using a QNM with Broyden’s update formula (see Section 3.1) and an adjoint method (see Section 3.2).

The first method is based on quasi Newton method with Broyden’s update formula (Nocedal and Wright, 2006). In doing so, equation (1) can be written as follows:

with I the identity matrix, q the search direction, p^ref the a priori information, λ the line search parameter (determined by Armijo rule) and $B$ the approximated Jacobian using Broyden’s update formula:

For the initialization of the Jacobian $Binit$⁠, the central difference method is used. To start at iteration step k = 0 with a good approximation, the finite difference method is applied to compute the Jacobian:

with $ΔBi,jsim=Bi,jsim(p+)−Bi,jsim(p−)$ and $Δpm=pm+−pm−$⁠, whereby $pm+=νr,m(1+ϵ)$ and $pm−=νr,m(1−ϵ)$⁠.

The adjoint method enables the direct computation of the gradients of our parameter vector p (Hinze et al., 2009). In doing so, we rewrite equation (1) as follows:

with the functional J being:

In a next step, we introduce the Lagrange function $L$⁠:

which is totally equal to our functional J, as we just added zero. In equation (9), λ denotes the vector of Lagrange multipliers. In a next step, we compute the gradient of the Lagrange function with respect to the parameters p resulting in the following:

Now, equation (10) can be rearranged as follows:

When we now set the terms in the parenthesis to zero, ∂u/∂p is no longer needed in the computation of the gradient of the functional J with respect to the parameters p. In doing so, we obtain the following equation for the vector of Lagrange multipliers:

which is also named adjoint equation. Finally, we arrive at the following:

The solution of the magnetostatic field is performed by the FE method, which discretizes the weak formulation via the Galerkin method. As our operator $A$ is self-adjoint, the weak form of the left-hand side of (12) reads as follows:

with the test function v. For the right-hand side of equation (12), we explore the Gateux – derivative and obtain for each position j of the sensor actuator system:

As each sensor has a finite volume, we evaluate equation (15) via an integral over each sensor volume Ω_i, and arrive at the weak form of the adjoint equation for each fixed sensor actuator position j:

which is also solved via the Galerkin FE formulation. Please note that the adjoint formulation results in the same system matrix as the forward simulation. Therefore, using a direct solver, the conjugation for the adjoint solution is highly efficient. To obtain the gradients of the parameter vector p, the term ∂A/∂p is needed [see equation (13)], which computes by the following:

where ∂v_i/∂v_j is one for i = j and zero else. Finally, the term ∂J/∂p in equation (13) calculates according to equation (8) by the following:

These results allow the evaluation of equation (13), which is the gradient q for adapting the material parameter p. The iterative procedure reads as follows:

with λ the line search parameter (determined by Armijo’s role).

To check the accuracy of the calculated gradient using the adjoint method, it is compared with the gradient obtained from using the finite difference method. For the sake of simplicity, only a variation of the magnetic reluctivity in the subdomain Ω₁ is assumed and the sensor-actuator system is positioned as shown in Figure 2. The relative error and the resulting gradients are shown in Table 1.

As the QNM and the adjoint method are iterative solution strategies, a stopping criterion has to be defined. In doing so, the following error norm is used:

Furthermore, the choice of regularization parameter is crucial for achieving an optimal solution during the iterative process. If the regularization parameter is set too high, the solution prioritizes the regularization term, while setting it too low can lead to divergence of the iterative process. According to the accuracy and resolution of the sensors, an a priori upper bound β for the error norm is available:

In equation (22), B^exact denotes the exact data without noise, and the discrepancy principle of Morozov (1968) is used. Therefore, starting from an initial regularization parameter α^init, the regularization parameter is reduced by each iteration step:

until the following condition is fulfilled:

For all computations, a = 0.5 and α^init = 1 has been chosen.

For the comparison of the two methods, electrical steel sheets with different cutting-edge influence are considered. It is assumed that the magnetic material properties decrease exponentially toward the cutting edge, using the empirical formula (Bali et al., 2014):

with v_r(x) the relative reluctivity, γ the degradation factor, δ the degradation skin depth and $νrglobal$ the relative reluctivity of the bulk material. Based on this formula, the initial reluctivity v^init(x) (initial configuration for the inverse scheme), reference reluctivity v^ref(x) (used for the Tikhonov regularization) and exact reluctivity v^exact(x) are computed for the case of small (see Figure 3) and for the case of large (see Figure 4) cutting-edge effect. The parameters for exact reluctivity v^exact(x) were chosen by the authors based on the investigations of Nakata et al. (1992), Bali et al. (2014), Holopainen et al. (2017) and Sundaria et al. (2020) and should reflect the real material behavior under the influence of cutting edges. The parameters are listed in Table 2.

In the simulation, each electrical steel sheet is divided into five subregions. The discrete reluctivities $νiinit$⁠, $νiref$ and $νiexact$ are calculated by integrating the corresponding continuous curves over the region length (see Table 3, where the reciprocal value is given, which corresponds to the relative permeability).

The reluctivities (⁠$νr,iinit$⁠, $νr,iref$ and $νr,iexact$⁠) for the subdomain i = 5 in both cases are equal. This is based on the assumption, that the material properties of the bulk material are known from SST or Epstein measurements. Thus, the parameter v_r,₅ is excluded from the optimization procedure and remains constant during the iterative procedure. The continuous and discrete material distribution for small and large cutting-edge effect are shown in Figures 3 and 4.

Measuring the magnetic flux density above the steel sheets, Hall sensors (in total seven sensors) are used, uniformly distributed along the line [(−3,0.4,0),(3,0.4,0)] in mm. In total, the electrical steel sheets are measured at six different positions, whereby the first measurement were taken such that sensor s₄ (see Figure 1) was directly above the cutting edge. The additional measurements were performed such that sensor s₄ was in the middle of each subregion.

The measurement data B^meas are generated artificially by forward simulations solving the magnetic field for the magneto-static case considering the exact material values. Furthermore, the generated data are overlaid by a Gaussian white noise $N(0,σ2)$⁠. For the sake of completeness, this means that the measurement noise is added to each measured flux density value for each sensor at each sensor-actuator position. To avoid an inverse crime, different meshes are used for the generation of the measurement data [Figure 5(a)] and for the inverse procedure [Figure 5(b)].

Based on the given data, the convergence behavior for the searched-for parameter vector p is investigated for the QNM and the adjoint method. In a first step, we demonstrated the convergence of both methods to the exact solution while analyzing measurement data without Gaussian white noise. Therefore, we consider the case with large cutting-edge effects and no a priori information. The outcomes of the QNM and adjoint method are presented in Figures 6, .

A few remarks pertaining to the results depicted in Figure 6 may be drawn. The convergence behavior for the QNM exhibits a smoother and faster progression when compared to the adjoint method. This can be clarified through the following observations. First, the QNM uses an initialization strategy using the finite difference method for the Jacobian, which results in a well-defined approximation of the Jacobian. Consequently, during the iterative process, significant adaptations through Broyden’s update formula are unnecessary. In addition, the QNM uses information from the second derivative through an approximation of the Hessian matrix. On the other hand, the adjoint method depends only on first-order derivative information, similar to an optimization method resembling the steepest descent approach. It is worth noting that a reduction in the initial regularization parameter, denoted as α^init, holds the potential to enhance the convergence behavior of the adjoint method. We therefore direct the reader’s attention to Figure 6(b). Until iteration step 5, no significant alterations in the parameters are recognizable, which suggests that the regularization term $||p−pref||22$ exerts a more pronounced influence in the initial iterations when compared to the error term $||Fi(xj,p)||22$ (see (1)).

To evaluate the performance of the proposed solution strategies under more realistic conditions, measurement data are perturbed with Gaussian white noise having a standard deviation of 10%. In the subsequent analysis, we compare the convergence behavior of the searched-for parameter vector p along with the error norm ε, as defined in equation (21). Furthermore, to evaluate the stability of the optimization procedure, we allow them to run for 25 iterations. For the sake of completeness, we denote the iteration step at which the stopping criterion would have been accomplished by a vertical black line in the subsequent plots.

In a first case, we make no a priori assumptions about the reluctivity and set p^ref in equation (1) to zero. The outcomes for both methods under conditions characterized by large cutting-edge effects and small cutting-edge effects are illustrated in Figures 7 and 8.

As explained previously, the adjoint method does not bring significant improvements during the first iterations, due to the initial regularization parameter. As the number of iterations increases, the stability of both methods is observed. In an overall evaluation, the adjoint method shows slightly reduced errors for both the searched-for parameter p and the error norm ε, while the QNM displays a relatively faster convergence with respect to the iteration steps at which the stopping criterion is satisfied.

The second investigation involves the incorporation of a priori information $νiref$⁠, detailed in Table 3. Similarly, both methods and cases are used as in the previous scenario. The results for the large and small cutting-edge effects are presented in Figures 9 and 10.

For the sake of completeness, the computation times and the total number of iterations (including all the iteration steps used for the line search) for the QNM and the adjoint method are listed in Table 4. Despite a similar number of total iterations, the simulation time of the adjoint method is about two times longer than that of the QNM. This is due to the currently inefficient implementation of the adjoint method in the finite element solver openCFS (Kaltenbacher, 2015), where we do not take advantage of the fact that the system matrix of the forward problem and the adjoint problem is the same. By optimizing the implementation, a significant improvement in simulation times for the adjoint method can be expected. It should also be mentioned here that the optimizers currently in use are implemented by the authors. Using established optimizers has the potential to improve convergence in general, thus reducing simulation time and the total number of iterations.

Incorporating a priori information produces notable differences. First, there is an improvement in the iteration step at which the stopping criterion is satisfied. Moreover, during the first iteration, there is a reduction in the deviation of the searched-for parameter vector p and the residual error ε compared to the results obtained without the inclusion of a priori information. However, it is noteworthy that taking into account a priori information does not seem to have a significant effect on the overall stability of the methods, nor does it significantly affect the residual error at higher iteration steps.

A crucial aspect that deserves discussion concerns the computation of the residual norm and the implications for the inverse procedure. Evaluating the previous results concerning the iteration steps at which the stopping criterion is satisfied and the convergence behavior of the searched-for parameter vector p, it becomes clear that this criterion is mainly satisfied when the error for larger subdomains, e.g. v₄ and v₃, is small, while the material error for the small subdomains, e.g. v₂ and v₁, has little influence. This behavior can be explained by the relationship between the change in material parameters and the effect on the magnetic field. A variation in the material parameter associated with larger subdomains leads to a pronounced change in the magnetic field, which is detected by more sensors than a variation of material parameter associated with small subdomains. Consequently, the procedure tends to optimize the material parameter of large subdomains due to the pronounced change in the residual norm.

In this work, the ideas of Gschwentner et al. (2023) based on a sensor-actuator model and the QNM using Brodyen’s update formula to locally determine the magnetic material behavior are taken up. In doing so, a dedicated adjoint method for solving the inverse problem is introduced and described in detail. Both methods are tested numerically for different cutting-edge effects, by generating the measurement data artificially by forward simulations. These generated data are overlaid by a Gaussian white noise with 10% standard deviation. Overall, both methods show a similar and fast convergence behavior for the investigated cases. Thereby, the importance of an a priori knowledge of the expected magnetic reluctivity values has been demonstrated. The investigation considering a priori information of the expected reluctivity values resulted in a strongly faster convergence. However, it has to be noted that the two proposed inverse schemes can also cope with the situation of no a priory knowledge. In future work, both methods will be extended, such that nonlinear and even hysteretic material models can be considered and tested with real world measurements.

The work is supported by the joint DFG/FWF Collaborative Research Centre CREATOR (CRC – TRR361/F90) at TU Darmstadt, TU Graz and JKU Linz.

Erratum: It has come to the attention of the publisher that the article Gschwentner, A., Kaltenbacher, M., Kaltenbacher, B. and Roppert, K. (2024), “Comparison of a quasi Newton method using Broyden’s update formula and an adjoint method for determining local magnetic material properties of electrical steel sheets”, COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. ahead-of-print No. ahead-of-print. https://doi.org/10.1108/COMPEL-11-2023-0566, contained several symbol errors in the equations and text. These errors were introduced during the production process and have now been corrected. The publisher sincerely apologises for this error and for any confusion caused.