*Article* **Accounting for Magnetic Saturation E**ff**ects in Complex Multi-harmonic Model of Induction Machine**

### **Tomasz Garbiec \* and Mariusz Jagiela**

The Faculty of Electrical Engineering, Automatic Control and Informatics, Opole University of Technology, 45–758 Opole, Poland; m.jagiela@po.edu.pl

**\*** Correspondence: t.garbiec@po.edu.pl

Received: 4 August 2020; Accepted: 7 September 2020; Published: 8 September 2020

**Abstract:** Computations of quasi-dynamic electromagnetic field of induction machines using the complex magnetic vector potential require the use of the so-called effective magnetization curves, i.e., such in which the magnetic permeability is proportional to the amplitudes of magnetic flux density *B* or magnetic field strength *H*, not their instantaneous values. There are several definitions of that parameter mentioned in the literature provided for the case when *B* or *H* are monoharmonic. In this paper, seven different methods of determining the effective magnetization curves are compared in relation to the use of a field-circuit multi-harmonic model of an induction machine. The accuracy of each method was assessed by computing the performance characteristics of a solid-rotor induction machine. One new definition of the effective permeability was also introduced, being a function of multiple variables dependent on amplitudes of all the harmonics considered. The analyses demonstrated that the best practical approach, even for the multi-harmonic case, is to express the effective magnetic permeability as the ratio of the amplitudes of the fundamental harmonics of the magnetic flux density and the magnetic field strength, and assuming the sinusoidal variation of the latter.

**Keywords:** magnetic permeability; effective parameters; induction motor; finite element method

### **1. Introduction**

Numerical models based on the use of the finite element method have become an indispensable tool in the design of electrical machines. As is well known, the most accurate approach for overall consideration of physical phenomena in an analyzed converter is the use of time-domain models. Unfortunately, this approach usually involves significant computational cost, sometimes making its practical application impossible. For that reason, it is justified to explore methods which allow, particularly in the analysis of steady states, to achieve similar results, but in a much shorter time. While analyzing induction machines, one of such methods is the use of the multi-harmonic field-circuit model [1–9]. The results presented in [1,2,9] prove that this type of the field-circuit model can be particularly useful when analyzing steady states of the special constructions of induction machines (single-phase machines, high-speed solid-rotor machines). As demonstrated by Garbiec et al., with an appropriate modification of the coupling scheme between individual sub-models, the multi-harmonic field-circuit model makes it possible to achieve comparable calculation results with those obtained using the time-domain model, at a calculation cost comparable to the monoharmonic model. The multi-harmonic model in this case is understood as the De-Gersem type model formulated using the magnetic vector potential, where permeance slot harmonics of the magnetic field distribution in the air gap are extracted via Fourier transform of an air-gap magnetic field. These harmonics can affect the power losses in the rotor and the electromagnetic torque developed by the machine. The magnetic field in the stator in these models is considered as a function of the fundamental time harmonic only.

The application of the above-mentioned approach essentially involves the problem of modeling the nonlinearity of the magnetic materials. Taking this into account in the monoharmonic field models based on the use of the complex magnetic vector potential method is not a new task and it has been described to date in several works, including [10–16]. They propose several different methods of defining the so-called effective magnetic permeability based, among others, on the development of nonlinear waveforms in the Fourier series [14,16], averaging the magnetic reluctivity [10,15,16] or equivalence of energy stored in ferromagnetic components [10–13]. However, analyzing the conclusions drawn in the mentioned works, it is not possible to clearly state which method of defining the effective magnetic permeability will be most appropriate when using the multi-harmonic model with a strong coupling, as proposed by Garbiec et al. Some of them assume the sinusoidal variation of the magnetic flux density in calculating the effective magnetic permeability (which facilitates the calculations using the magnetic vector potential) [12–15], while others assume that it should be determined assuming sinusoidal variation of the magnetic field strength, or that it is insignificant for the accuracy of the calculations [10,11,14]. The latter approach was adopted in the authors' previous work, and satisfactory results were achieved [9]. It was based on the results of the earlier work of the authors, where the effective magnetic permeability defined with the assumption of the sinusoidal variation of the magnetic field strength provided correct results, modeling a solid-rotor induction machine with the use of the multi-harmonic field-circuit model with a weak coupling [17]. To the best of authors' knowledge, the research providing an analysis of the accuracy of definitions of the effective magnetic permeability on the results of computations using multi-harmonic field-circuit models was not reported.

This paper is a continuation of works by Garbiec et al. that aim at an in-depth verification of the De-Gersem type multi-harmonic complex model of induction machines from the point of view of their application in designing and optimization of real machines. The authors focus on the analysis of the impact of formulating the effective magnetic permeability. For that purpose, the basic operating characteristics of an exemplary machine were calculated for seven different methods of formulating the effective magnetic permeability. The first six methods use existing definitions derived in literature. The available formulas use averaging of magnetic reluctivity in terms of mean or RMS value. Also, the Fourier series truncated to the fundamental component and assuming sinusoidal variation of the magnetic flux density or magnetic field strength is used [14–16]. The seventh method considered in this paper is a new approach proposed by the authors that relies upon formulating the so-called multidimensional effective magnetic permeability. The results of calculations were compared with those obtained with the time-domain model and with the results of measurements performed on the physical model under one of the previous works where the efficient computational method for fast determination of the performance characteristics of solid-rotor induction machines was presented [17].

### **2. Mathematical Model**

### *2.1. E*ff*ective Magnetic Permeability*

Definition of the effective magnetic permeability, as used in the previous works of the authors [9,14,16,17], is the reference point for the considerations presented herein:

$$
\mu\_{\rm I} \{ H\_{pk} \} = \frac{2}{\pi H\_{pk}} \int\_0^\pi \left( \mu\_{\rm DC} \{ H\_{pk} \sin \alpha \} H\_{pk} \sin \alpha \right) \sin \alpha d\alpha,\tag{1}
$$

where *Hpk*—amplitude of the sinusoidal magnetic field strength *H*, and μ*DC*—DC magnetic permeability. The above definition corresponds to the ratio of the fundamental harmonic amplitude of the magnetic flux density waveform (obtained on the assumption of the sinusoidal variation of the magnetic field strength) and *Hpk*. This approach requires a successive determination of the magnetic field

*Energies* **2020**, *13*, 4670

strength in a given iteration of the algorithm based on the knowledge of its value determined in the previous iteration.

Assuming the sinusoidal variation of the magnetic field strength, two further definitions can be introduced [16]:

$$
\mu\_{\rm II} \left( H\_{pk} \right) = \frac{\pi}{\int\_0^{\pi} \nu\_{\rm DC} \left( H\_{pk} \sin \alpha \right) d\alpha} \,\tag{2}
$$

$$
\mu\_{\rm III} \left( H\_{pk} \right) = \frac{1}{\sqrt{\frac{1}{\pi} \int\_0^{\pi} \nu\_{\rm DC}^2 \left( H\_{pk} \sin \alpha \right) d\alpha}},
\tag{3}
$$

where ν*DC* = <sup>1</sup> <sup>μ</sup>*DC* . The expressions (2) and (3) correspond to the inverse of the mean value and the inverse of the RMS value of the magnetic reluctivity obtained for the half of the period.

Similarly, the definitions of the effective magnetic permeability can be derived, assuming the sinusoidal variation of the magnetic flux density *B* with the amplitude *Bpk* [15,16]:

$$
\mu\_{\rm IV} \left( B\_{pk} \right) = \frac{\pi B\_{pk}}{2 \int\_0^\pi \left( \nu\_{\rm DC} \left( B\_{pk} \sin \alpha \right) B\_{pk} \sin \alpha \right) \sin \alpha d\alpha},\tag{4}
$$

$$
\mu\_{\rm V} \left( B\_{pk} \right) = \frac{\pi}{\int\_0^{\pi} \nu\_{\rm DC} \left( B\_{pk} \sin \alpha \right) d\alpha},
\tag{5}
$$

$$
\mu\_{\rm VI} \left( B\_{pk} \right) = \frac{1}{\sqrt{\frac{1}{\pi} \int\_0^{\pi} \nu\_{\rm DC}^2 \left( B\_{pk} \sin a \right) d\alpha}}. \tag{6}
$$

Each of those definitions assumes that the magnetic saturation coefficient of the ferromagnetic material is a function of the fundamental harmonic of the magnetic field only, expressed by *B* or *H.* The definitions (1), (3), (5) and (6) are taken from the literature, whereas (2) and (4) are the equivalents of (5) and (1), respectively.

### *2.2. Analyzed Machine*

In this work, the effect of how the effective permeability is taken into account was examined on the example of a high-speed induction machine with a uniform solid rotor analyzed previously by Garbiec et al. The basic data of the considered motor is presented in Figure 1 and Table 1.

**Figure 1.** Solid-rotor induction machine taken into consideration: (**a**) dimensions of the stator package (cross section area), (**b**) dimensions of the rotor.


**Table 1.** Basic specifications of tested induction motor.

A field-circuit multi-harmonic model with a strong coupling was developed for the analyzed machine. Based on the results of the analysis presented in by Garbiec et al., it was decided to reduce the number of the rotor models to three, i.e., those related to the fundamental harmonic of the magnetic field in the air gap and two most influential higher slot harmonics with ordinal numbers—11 and 13. Neglecting the rotor-end effects the system of equations describing the model takes the form [9]:

$$
\begin{bmatrix}
\mathbf{M}\_{11} \begin{pmatrix} \mu\_{eff} \end{pmatrix} & \mathbf{M}\_{12} & \mathbf{M}\_{13} \\
\mathbf{M}\_{21} & \mathbf{M}\_{22} & 0 \\
\mathbf{M}\_{31} & 0 & 0
\end{bmatrix}
\begin{array}{c}
\mathbf{P} \\
\bar{\mathbf{I}} \\
\bar{\lambda}\_{-} \\
\bar{\lambda}\_{-}
\end{array}
= 
\begin{bmatrix}
0 \\
\mathbf{E} \\
0 \\
\bar{0}
\end{bmatrix} \tag{7}
$$

where: **M**11—matrix describing the magnetic and electrical properties of materials, **M**<sup>12</sup> = <sup>−</sup>**M***<sup>T</sup>* <sup>21</sup>/(jω*lz*)—matrix describing the distribution and connection method of the stator winding, **M**<sup>13</sup> = **M***<sup>T</sup>* <sup>31</sup>—matrices describing coupling between the rotor models and the stator model, **M**22—stator winding impedance matrix, <sup>ϕ</sup>**\_** —vector of the nodal values of the complex magnetic vector potential, **I \_** —vector of the amplitudes of the loop currents in the stator winding, λ **\_** —vector of complex circulations of the magnetic field strength vector, **E \_** —vector of the complex voltage amplitudes in the loops in the stator winding circuit.

For a comparison, in this paper, calculations were also carried out using a developed complementary time-domain model which uses the DC magnetization curves. The equations for such a model after discretization via Galerkin procedure and the implicit Euler method take the form:

$$
\begin{bmatrix}
\mathbf{S}(\boldsymbol{\mu}\_{\rm DC}) + \mathbf{G}\boldsymbol{\Lambda}^{-1} & -\mathbf{D}^{\mathsf{T}}\mathbf{K}^{\mathsf{T}} \\
 l\_{\boldsymbol{z}}\mathbf{K}\mathbf{D}\boldsymbol{\Lambda}^{-1} & \mathbf{K}\Big{\Big{}}\mathbf{R} + \mathbf{L}\boldsymbol{\Lambda}^{-1}
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{\Psi} \\
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{\Psi} \\
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{G}\boldsymbol{\Lambda}^{-1} & \mathbf{0} \\
l\_{\boldsymbol{z}}\mathbf{K}\mathbf{D}\boldsymbol{\Lambda}^{-1} & \mathbf{K}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{-1}\mathbf{K}^{\mathsf{T}}
\end{bmatrix}^{\mathsf{T}-1}
\begin{bmatrix}
\boldsymbol{\Psi} \\
\end{bmatrix}^{\mathsf{T}-1} + \begin{bmatrix}
\boldsymbol{\Psi} \\
\end{bmatrix}^{\mathsf{u}-1},\tag{8}
$$

where: **S**—reluctivity matrix, **G**—conductivity matrix, **D**—matrix describing the winding, **K**—matrix describing the winding connection method, **R**—winding resistance matrix, **L**—winding leakage inductivity matrix, Δ*t*—integration step, ϕ—vector of nodal values of the vector magnetic potential, **i**—vector of instantaneous values of the stator currents, **e**—vector of the instantaneous supply voltages values. The rotational movement was modeled using the moving band technique. The computing cost for the multi-harmonic model (when choosing an appropriate coupling scheme) is of the order of a per cent compared to that using the time-domain model [9]. This is the main reason for developing the model described by the Equation (7), however, in this paper the authors focus on the accuracy rather than on the computing cost. The nonlinearity in both above formulated computational problems is solved iteratively by updating, respectively μ*eff* and μ*DC* from iteration to iteration unless the relative change of solution is below a prescribed value.

Based on the DC magnetization curves for the stator and rotor materials, the effective magnetization curves were determined according to (1)–(6) (Figure 2). For saturated material, each effective magnetization curve is noticeably characterized by higher magnetic flux density values compared to the characteristic for the direct current. At a low level of magnetic saturation, the effective curves ale very close to the DC curve. Furthermore, a relatively small difference is found for different curves determined assuming the sinusoidal variation of the magnetic flux density for the rotor material. Moreover, the curves calculated with the assumption of a sinusoidal variation of the magnetic flux density are closer to the DC curve than the case with the assumption of a sinusoidal variation of the magnetic field strength. Additionally, curves (4)–(6) are closer to each other than (1)–(3). This is due to the completely different shape of the magnetic flux density waveforms (magnetic field strength) in the case when the magnetic field strength (magnetic flux density) is sinusoidal.

**Figure 2.** Effective magnetization curves determined from the expressions (1)–(6): (**a**) stator, (**b**) rotor.

### **3. Calculation Results**

The magnetization curves shown in Figure 2 were used to determine three basic operational characteristics of the machine (electromagnetic torque, RMS phase current, power losses in the rotor) vs. rotational speed at two different supply voltages (phase voltages), i.e., 50 V (rated) and 75 V (150% of the rated value) without changing the frequency of power supply. Note that the effective magnetic permeability for the rotor models is determined in each iteration based on the knowledge of the magnetic field strength distribution (expressions (1)–(3)) or magnetic flux density ((4)–(6)) only for the rotor model associated with the fundamental harmonic of the magnetic field distribution. The characteristics calculated were compared with the results of calculations using time-domain model, as described in Section 2. It should be kept in mind that the considered multi-harmonic model, despite the associated mathematical burden, is a simplified approach. This simplification can be observed in e.g., the obtained steady-state distributions of the magnetic flux over the machine cross-section shown in (Figure 3). Based on the results of that comparison, two groups of characteristics were created

(Figure 4) showing the total relative percentage calculation error for the three determined quantities (electromagnetic torque, RMS phase current, power loss in the rotor) in relation to the calculation results obtained from the model (8):

$$
\varepsilon = \frac{\varepsilon\_{T+\mathcal{E}I+\mathcal{E}P}}{3} \times 100\% \tag{9}
$$

where: ε*<sup>T</sup>* = *Tt*−*Tp Tt* , <sup>ε</sup>*<sup>I</sup>* <sup>=</sup> *It*−*Ip It* , <sup>ε</sup>*<sup>P</sup>* <sup>=</sup> *Pt*−*Pp Pt* , *Tt*,*It*, *Pt*—respectively, electromagnetic torque, stator RMS phase current, power loss in the rotor calculated with the model formulated in the time domain, *Tp*, *Ip*, *Pp*—respectively, electromagnetic torque, stator RMS phase current, and power loss in the rotor calculated using the multi-harmonic model. Table 2 shows the averaged values of that indicator against the entire speed range considered.

**Figure 3.** Comparison of sample distributions of magnetic flux over cross section of the motor at rotor slip equal to 10%: (**a**) time-domain model, (**b**) multi-harmonic model with all harmonics considered (the magnetic field in the rotor area is a sum of the complex magnetic vector potential distributions obtained from solution of (7) for the individual rotor submodels).

**Figure 4.** Variation of relative error (9) vs. rotational speed using various methods of calculating the effective magnetic permeability for the RMS phase supply equal to: 50 V (**a**), and 75 V (**b**).


**Table 2.** Averages percent value of the relative error for *N* = 11 considered points of the curves.

Out of the definitions of the effective magnetic permeability considered, noticeably the best results were achieved using μ<sup>I</sup> (originally used in calculations). Slightly worse results (except for the locked rotor state) can be obtained with the definitions assuming the sinusoidal variation of the magnetic flux density. On the other hand, an interesting conclusion is that the method how the nonlinearities are taken into account has practically no impact on the results of calculations for near-synchronous rotational speeds. It is caused by a much lower saturation of the rotor surface, crucial from the point of view of the machine properties, and consequently use by the time-harmonic model (7) similar values of effective permeability (see Figure 2). An increased error in the range of nominal operation point (21–27 krpm) is an adverse trend observed. However, note that the waveforms of the global error indicator shown in Figure 4 provide information on the discrepancy in the results of calculations of power losses in the rotor (ε*P*). The results of those calculations are distinguished by the greatest discrepancy (for μ*<sup>I</sup>* and *N* = 11 considered points of the curve *<sup>N</sup> <sup>i</sup>*=<sup>1</sup> ε*Pi <sup>N</sup>* = 15.42% at 50 V and 19.32% at 75 V), thus overestimate the global error indicator. In the rotational speed range from 21,000 rpm to 27,000 rpm, the average electromagnetic torque error and the RMS current error is 7.56% and 1.43% at 50 V, 10.93% and 1.67% at 75 V, respectively.

Calculating the effective magnetic permeability based on the distribution of the magnetic field strength or magnetic flux density, determined only for the model associated with the fundamental field harmonic in the machine air gap, can raise doubts. To take into account the impact of the higher harmonics and to assess the error related to this assumption, a new definition of the multidimensional effective magnetic permeability was proposed:

$$\mu\_{VII} \left( H\_{1pk'} H\_{11pk'} H\_{13pk} \right) = \frac{\sqrt{B\_{1pk}^2 + B\_{11pk}^2 + B\_{13pk}^2}}{\sqrt{H\_{1pk}^2 + H\_{11pk}^2 + H\_{13pk}^2}} \tag{10}$$

where *H*1*pk*, *H*11*pk*, *H*13*pk* are the magnetic field strength amplitudes related to the fundamental, eleventh and thirteenth harmonics of the field and

$$\begin{aligned} B\_{\eta pk} &= \frac{2}{\pi} \int\_0^\pi \mu\_{D\zeta} \left( H\_{1pk} \sin a + H\_{11pk} \sin 11a + H\_{13pk} \sin 13a \right) \\ \left( H\_{1pk} \sin a + H\_{11pk} \sin 11a + \dots + H\_{13pk} \sin 13a \right) \sin nadaf &= 1,11,13 \end{aligned} \tag{11}$$

For zero values of the higher harmonics, the definition (10) is equivalent to (1). This approach only slightly complicates the model implementation (in each algorithm iteration the distribution of the magnetic field strength in all considered rotor models needs to be obtained). An advantage of this expression, however, is the fact that effective permeability can be expressed before starting the calculations, without knowing the individual amplitudes. For the expression (10), it is possible since *Bnpk* does not depend on the phase angles of individual harmonics. In general, this occurs when the frequencies of the higher harmonics are the multiples of the fundamental harmonic. The value of the multidimensional effective permeability is determined only once as a multidimensional table for the assumed variability ranges of individual harmonics of the magnetic field strength. Figure 5 shows the determined distribution of this quantity assuming that *H*13*pk* = 0 A/m. The multidimensional permeability determined was used again for calculations and compared with the results obtained for (1). The result of that comparison is shown in Table 3. As can be seen, the application of the multidimensional effective permeability brings only slight improvement to the accuracy of calculation results. This is due to relatively low magnetic field strength values in the areas of the models associated with the slot harmonics, as illustrated in Figure 6.

**Figure 5.** Relative multidimensional effective permeability calculated for the DC magnetization curve of the rotor, assuming zero value of the magnetic field strength of the thirteenth harmonic.

**Table 3.** Averages percent value of the relative error for *N* = 11 considered points of the curves.

**Figure 6.** Comparison of the magnetic field strength values in the individual elements creating FE rotor mesh for the rotor models associated with the considered harmonics for a 20% slip: (**a**) at 50 V, (**b**) at 75 V.

### **4. Physical Validation**

The analysis demonstrated that, from a practical point of view, it is sufficient to assume the first definition of the effective magnetic permeability, both due to the accuracy of calculations, and ease of implementation. It should be mentioned that, when adopting the sinusoidal variation of the magnetic field strength in determining the effective magnetic permeability, the fixed point method is applied in the multi-harmonic model to solve the nonlinearity problem, which is easy to implement, and it is distinguished by good convergence quality in case of time-harmonic magnetic fields. In the latter case, it is recommended to use the Newton–Raphson method [16].

μ μ

The model developed, together with the first definition of the effective magnetic permeability, was applied to calculate the RMS current of a real machine (supplied via a quasi-square wave voltage inverter, fundamental frequency of supplying voltage equal to 500 Hz) with three different solid rotors tested by Garbiec and to compare the results with the calculations performed using a model with a weak coupling the same definition of the effective magnetic permeability applied (see Figure 7).

**Figure 7.** Comparison of the ideas of the different multi-harmonic model: (**a**) model with the weak coupling, (**b**) model with the strong coupling.

As it can be seen in the figure the idea of using the model with the strong coupling is to simultaneously solve the problem composed of a few submodels (the stator model and rotor models associated with the considered harmonics) coupled by appropriate boundary conditions imposed on the interface boundary in the air gap. This formulation results in one system of equations and determination of the magnetic vector potential distribution for all models via single solution of (7). The interface conditions in the air-gap are formulated with the use of matrix operators of the Fourier transform and the inverse Fourier transform, as shown in Garbiec at. al, enforce the propagation of only the selected harmonics of the air-gap field distribution between the stator model and the considered rotor models. All rotor models act via the same constraint on the stator model. In the considered model with the weak coupling, the monoharmonic field-circuit model is used. The power loss due to slot harmonics is estimated by means of extraction of these harmonics from the calculated distribution of the magnetic vector potential. In this way it is possible to solve the field problem limited to only the area of rotor. In this approach the power factor, torque and stator current can be appropriately

corrected, although there is no backward relationship, i.e., the air-gap field remains unaffected by these effects. To this end, it was also necessary to take into account the rotor-end effect [18].

The physical test-stand and result of the comparison of computed and measured data for the machine with different rotors is shown in Figure 8. With the application of the multi-harmonic model with the strong coupling and with the modified sub-model coupling scheme, it was possible to achieve practically identical calculation results, yet in a shorter time by approx. 30%. The errors found result mainly from supplying the physical model via a quasi-square wave voltage inverter (in numerical models, the pure sine-wave voltage supply with an RMS value equal to the RMS value of the measured voltage is assumed).

(**a**)

**Figure 8.** *Cont.*

**Figure 8.** Comparison of calculation results with the multi-harmonic model with the weak coupling to those of calculations with the multi-harmonic model with the strong coupling, with the measurement results [17]. (**a**) Physical test-stand, (**b**) RMS current vs. rotor speed for the machine equipped with an uniform rotor, (**c**) RMS current vs. rotor speed for the machine equipped with a slitted rotor with uniform end regions, (**d**) RMS current vs. rotor for a machine equipped with a slitted rotor with slits through the whole length.

### **5. Conclusions**

The comparison of various methods for calculating the effective magnetic permeability for a multi-harmonic model of a solid-rotor induction machine carried out in the study proved that, from a practical point of view, the best results are obtained using the basic definition, described by Equation (1). Determining the effective magnetic permeability as the ratio of the magnitudes of the fundamental harmonics of the magnetic flux density and the magnetic field strength calculated assuming the sinusoidal variation of the latter allows for an easy implementation of the algorithm which is important from practical point of view. Furthermore, it was demonstrated that, in the studied case, the most influential magnetic saturation effects come from the fundamental harmonic of the magnetic field which means that the slot harmonics of the air-gap magnetic flux density do not have to be computed in that process. For the same definition of the effective magnetic permeability and an appropriate modification of the coupling scheme, it is possible to obtain practically the same results as in the case of the multi-harmonic model with the weak coupling, but yet in a shorter time. Another step to take will be to conduct detailed research using the developed model to analyze the magnetic saturation effects in the squirrel cage induction machines with both three- and single-phase windings.

**Author Contributions:** Conceptualization, T.G. and M.J.; methodology, T.G. and M.J.; software, T.G.; investigation, T.G.; resources, T.G.; data curation, T.G.; writing—original draft preparation, T.G.; writing—review and editing, M.J.; visualization, T.G.; supervision, M.J.; project administration, T.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Application of Surrogate Optimization Routine with Clustering Technique for Optimal Design of an Induction Motor**

**Aswin Balasubramanian \*, Floran Martin \*, Md Masum Billah, Osaruyi Osemwinyen and Anouar Belahcen \***

Department of Electrical Engineering and Automation, Aalto University, 02150 Espoo, Finland; md.billah@aalto.fi (M.M.B.); osaruyi.osemwinyen@aalto.fi (O.O.)

**\*** Correspondence: aswin.balasubramanian@aalto.fi (A.B.); floran.martin@aalto.fi (F.M.); anouar.belahcen@aalto.fi (A.B.)

**Abstract:** This paper proposes a new surrogate optimization routine for optimal design of a direct on line (DOL) squirrel cage induction motor. The geometry of the motor is optimized to maximize its electromagnetic efficiency while respecting the constraints, such as output power and power factor. The routine uses the methodologies of Latin-hypercube sampling, a clustering technique and a Box–Behnken design for improving the accuracy of the surrogate model while efficiently utilizing the computational resources. The global search-based particle swarm optimization (PSO) algorithm is used for optimizing the surrogate model and the pattern search algorithm is used for fine-tuning the surrogate optimal solution. The proposed surrogate optimization routine achieved an optimal design with an electromagnetic efficiency of 93.90%, for a 7.5 kW motor. To benchmark the performance of the surrogate optimization routine, a comparative analysis was carried out with a direct optimization routine that uses a finite element method (FEM)-based machine model as a cost function.

**Keywords:** induction motors; surrogate optimization; Box–Behnken design; Latin-hypercube sampling; clustering; particle swarm optimization; pattern search

### **1. Introduction**

Electrical machines have a wide range of use cases, from household utilities to industrial applications, which consume a huge share of all the generated electrical energy [1]. To reduce the global greenhouse gas emissions, it is important to design electrical machines with high energy efficiency. The characteristics of the electrical machines are usually analyzed with a finite element method (FEM)-based electromagnetic simulation for better accuracy. The output characteristics of the electrical machine are highly sensitive to the design variables and the global search optimization algorithms, such as particle swarm optimization (PSO) or genetic algorithm (GA), require many model evaluations to reach the desired optimal solution [2,3]. This causes the optimization process with a FEM-based machine model as a cost function to be computationally expensive [4]. To utilize the time and computational resources efficiently, surrogate optimization techniques are used to optimize the electrical machines, which requires only a few FEM simulations for evaluation.

Response surface methodology (RSM) is a technique used to develop a polynomial function for a complex FEM-based multi-physics model, which defines the relationship between design variables and the output response of an electrical machine [5–9]. This polynomial model can be used with an optimization algorithm to search for the optimal solution. The Box–Behnken design, one of the popular response surface approaches, is used in conjunction with the FEM-based machine model to generate second-order polynomial functions for objective and constraints of the electrical machine [10–14]. The range of the boundaries of the design variables affect the accuracy of the polynomial function and in turn the optimal solution of the electrical machine.

In this article, a novel surrogate optimization routine is proposed for optimizing a three-phase direct on line (DOL) squirrel cage induction motor. The rotor bars of the

**Citation:** Balasubramanian, A.; Martin, F.; Billah, M.M.; Osemwinyen, O.; Belahcen, A. Application of Surrogate Optimization Routine with Clustering Technique for Optimal Design of an Induction Motor. *Energies* **2021**, *14*, 5042. https:// doi.org/10.3390/en14165042

Academic Editor: Mario Marchesoni

Received: 8 July 2021 Accepted: 13 August 2021 Published: 17 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

induction machine are not skewed for the sake of the demonstration of the surrogate optimization routine. The aim of the new routine is to discretize the problem domain into a number of subdomains for improving the accuracy of the polynomial models for a better search of the optimal solution, while efficiently using the computational resources with a smaller number of FEM simulations. To utilize the full capacity of the computational resource, the routine is programmed in such a way as to handle 15 FEM simulations in parallel. The methodologies used in the optimization process are Latin-hypercube sampling for design of the experiments, a clustering algorithm for dividing the problem domain and a Box–Behnken design as a response surface approach. The particle swarm optimization (PSO) algorithm is used for optimizing the polynomial functions of the response surfaces, while a pattern search algorithm is used for fine-tuning the surrogate optimal solution from the PSO. For the purpose of visualization, the proposed surrogate optimization routine is demonstrated with a simplification of the design problem, leaving three design variables. Validation of the results from the proposed surrogate optimization routine for a multivariate design problem is performed by comparing it to the direct optimization routine, which uses FEM simulation as a cost function.

### **2. Optimization Problem**

The electrical machine analyzed in the optimization problem is a three-phase squirrel cage induction motor for a direct online industrial application. The electrical steel core material used in the motor is M400-50A and the rotor cage is made of aluminum. The goal of the optimization problem is to maximize the electromagnetic efficiency, *η*, satisfying the constraints of the output power, *P*out, and power factor, *PF*, for a given volume of the machine. The outer diameter of the stator, *D*se, and axial length of the machine, *l*, are fixed so that the volume of the machine remains constant throughout the optimization process. The rotor end-ring overhang length, *l*oh, is kept constant so that the cross-section area of the end ring depends only on the height of the rotor slot. The specifications and fixed parameters of the induction motor for the optimization problem are shown in Table 1. The objective and constraints of the optimization problem are specified in Table 2. The optimization variables and their ranges for the induction motor are shown in Figure 1 and Table 3. The analysis of the machine design is done with timestepping simulation of a 2D finite element solver software, FCSMEK, developed by the research group of electromechanics at Aalto university [15]. The simulation of the timestepping analysis computes the electromagnetic characteristics of electrical machines by solving the circuit and field equations with the Crank–Nicholson timestepping method. The time is discretized at short time intervals and the magnetic field, currents, and potentials of the windings are solved at successive instants of time. The rotation of the rotor is accomplished by changing the finite element mesh in the air gap. The non-linear system of equations obtained at each timestep is solved using the Newton–Raphson method. The core losses are evaluated using the modified Jordan loss equation with a two-component loss model, namely with eddy current loss and hysteresis loss. The excess losses are included in the dynamic eddy current loss computation [15]. For simplicity, the mechanical losses and other manufacturing losses are not considered for comparing the results of the optimal solutions. Hence, the electromagnetic efficiency *η* is computed as shown in Equation (1).

$$\eta = \frac{P\_{\text{in}} - P\_{\text{elec}}}{P\_{\text{in}}} \times 100\% \tag{1}$$

where *P*in is the input power of the induction motor and *P*elec is the electromagnetic loss of the induction machine, which is comprised of iron losses and copper losses of the stator and rotor. The steady state temperature of the stator and rotor are considered as 80 ◦C and 100 ◦C, respectively, for the simulation.

**Figure 1.** Optimization variables of the induction motor.



**Table 2.** Objective and constraints of the optimization problem.



**Table 3.** Optimization variables and their range.

### **3. Response Surface Methodology Optimization with Box–Behnken Design**

Response surface methodology (RSM) is a set of mathematical and statistical techniques used to draw a relationship between control variables (inputs) and output response. This relationship can be approximated into a polynomial model, which can be useful in predicting the response of the control variables, hypothesis testing and finding the optimal condition of the variable settings [16]. In practice, the response surface methodology can be applied to simulate experimental results or for constructing a surrogate function of a computationally expensive multi-physics model. Optimizing the geometric variables of the induction motor directly with the finite element model is computationally expensive. In this article, for boundary-constrained input variables the RSM is used for approximating the electrical quantities of an induction machine into a second-order polynomial model. The second-order polynomial model that describes the functional relationship of the RSM between the control variables and the output response is as shown in Equation (2) [17].

$$y = \beta\_0 + \sum\_{i=1}^{k} \beta\_i \mathbf{x}\_i + \sum\_{i=1}^{k} \beta\_{ii} \mathbf{x}\_i^2 + \sum\_{i=1}^{k-1} \sum\_{j=i+1}^{k} \beta\_{ij} \mathbf{x}\_i \mathbf{x}\_j \tag{2}$$

where *y* is the output response, *xi* and *xj* are the input control variables, *β*0, *βi*, *βii*, *andβij* are the coefficients of the input control variable terms, and *k* is the number of control variables. The coefficients are estimated as shown in Equation (3).

$$\mathcal{J} = \left[\mathbf{X}^{\mathrm{T}}\mathbf{X}\right]^{-1}\mathbf{X}^{\mathrm{T}}\mathbf{Y} \tag{3}$$

where *X* is the matrix of input control variables sampled at multiple points and *Y* is the corresponding output response vector. The surrogate function shown in Equation (2) works well for interpolation of design variables to predict the output response, but prediction of output response by extrapolation of design variables can be inaccurate. One of the commonly used designs for determining the response surface, as shown in Equation (2), is the Box–Behnken design codeveloped by Box and Behnken in 1960 [18]. If the control variable space is defined as a cube, then the sample points are taken at the geometric center of the cube and at the middle points of the edges of the cube. The Box–Behnken design sample points represented for a three-dimensional control variable space are shown in Figure 2. Positioning sample points in this way preserves a uniform variance within the definition of the hyper-cube [18]. The number of sample points, *N*, for a given number of control variables is shown in Equation (4) [19].

$$N = 2k(k-1) + \mathbb{C}\_0\tag{4}$$

where *k* is the number of design variables and *C*<sup>0</sup> is the number of center points.

**Figure 2.** Box–Behnken design represented for a 3-dimensional variable space.

The Box–Behnken design was applied to the problem of an induction machine as described in Section 2 for performing the surrogate optimization. The process flow of optimization with the Box–Behnken design is presented in Figure 3. For easier representation, the optimization variables presented in Table 3 are assumed as *x*1, *x*2, *x*3..., *x*<sup>8</sup> in their respective order. Based on the optimization variable boundaries from Table 3, the Box–Behnken design sample points were created for the variables as an eight-dimensional hypercube. These samples were simulated with FCSMEK finite element software for calculating their corresponding response characteristics, such as efficiency, *η*, output power, *P*out, and power factor, *PF*. The relationship between the optimization variables and the output response was established as a second-order polynomial function as shown in Equation (2). The coefficients of the polynomial terms were calculated from the predetermined output responses by FEM simulations sampled at Box–Behnken sample points as shown in Equation (3). The polynomial response functions for the surrogate optimization problem are presented in Appendix A. These surrogate functions were used as the cost function of the PSO algorithm. A population of 1000 particles of the PSO was initialized with the Latin-hypercube sampling method.

**Figure 3.** Flow chart: optimization with Box–Behnken design.

The objective and constraints of the problem are as shown in Table 2. The optimal motor design from Box–Behnken design is validated with FEM simulation as shown in Table 4. It is observed that a difference in the result of the objective electromagnetic efficiency, *η*, between the surrogate optima and its FEM validation is considerably high. Moreover, the optimal solution does not respect the constraints of the output power, *P*out, and the power factor, *PF*, coupled with a high margin of error. The accuracy of the surrogate functions is impacted by the application of the Box–Behnken design to a large design variable space. Hence, a new optimization routine is proposed in this article (Section 6) for improving the accuracy of the surrogate functions resulting in an improved optimal solution.


**Table 4.** FEM validation of the optimal solution from Box–Behnken design

### **4. Latin-Hypercube Sampling**

Latin-hypercube sampling is a statistical method for selecting near-random samples from the input variable space, proposed by McKay, Beckman, and Conover [20]. It uses a sampling scheme of stratification to improve the distribution of samples in the input variable space. To select *n* Latin-hypercube samples for a sampling function with *xi* = (*x*1, *x*2) as input variables, the range of each of the *xi* is stratified into *n* equiprobable intervals. One observation is selected at random from each of the *n* intervals. These observations corresponding to *x*<sup>1</sup> and *x*<sup>2</sup> are matched at random to form *n* Latin-hypercube samples. A set of five samples generated with the Latin-hypercube sampling method for the input variables of *xi* = (*x*1, *x*2) is shown in Figure 4. The Latin-hypercube sampling method is widely used in the design of experiments for various applications of computer modeling [21–23]. In this article, the Latin-hypercube sampling method is used in the proposed optimization routine for selecting samples that satisfy a set of criteria and for initializing the first swarm of the PSO algorithm.

**Figure 4.** Latin-hypercube sampling.

### **5. Clustering**

Clustering is a method involved in partitioning a given dataset into different groups or clusters. The data that are mapped to a particular cluster tend to have similar characteristics and follow a similar pattern [24]. Clustering helps in classifying and analyzing large datasets, which can be applied in fields of machine learning, data science, pattern recognition, image processing, and bioinformatics [25–27]. The k-means is one of the oldest computational techniques used in solving clustering problems, based on the algorithm proposed by Lloyd [28]. If an integer *k* is chosen for partitioning a dataset into *k* clusters and *n* is the number of data points of the dataset, the goal of Lloyd's (k-means) algorithm is to find *k* centroids so as to minimize the potential function, *γ*. The potential function *γ* is a measure of the total squared Euclidean distance between each data point and its closest centroid, as shown in Equation (5).

$$\gamma = \sum\_{j=1}^{k} \sum\_{i=1}^{m} \left\| \mathbf{x}\_{i} + \mathbf{c}\_{j} \right\|^{2} \tag{5}$$

where *m* is the number of data points in the *j*th cluster, *xi* and *cj* are the data points and centroid of the *j*th cluster. Since Lloyd's (k-means) algorithm involves selecting the initial *k* centroids uniformly at random from the dataset, it suffers from inconsistency and accuracy issues. Arthur and Vassilvitskii proposed a randomized seeding technique for selecting the initial centroids of the k-means and combined it with the original k-means algorithm to call it the k-means++ algorithm with improved speed and accuracy [29]. A set of randomly generated samples clustered with the k-means++ algorithm is shown in Figure 5. In the proposed surrogate optimization routine, the k-means++ algorithm is used in partitioning the samples created from design variables into different clusters. In a multi-variable clustering problem, the variables can have varying scales of magnitude and incomparable units. Thus, it is required to apply the feature scaling technique to normalize the data for standardization. Z-score transformation is one of the successful standardization techniques utilized before applying the k-means clustering method to a dataset [30]. Equation (6) is used to estimate the Z-score,

$$\mathbf{x}' = \frac{\mathbf{x} - \overline{\mathbf{x}}}{\sigma\_{\mathbf{x}}} \tag{6}$$

where *x* is the variable vector that needs to be standardized, *x* and *σ*<sup>x</sup> are the mean and standard deviation of the vector *x*, and *x* is the transformed variable vector that implies the z-score.

**Figure 5.** A set of randomly generated samples clustered with the k-means++ algorithm.

### **6. The Proposed Surrogate Optimization Routine**

To overcome the issue in the accuracy of the surrogate model as shown in Section 3, a new surrogate optimization routine is proposed for the induction machine problem. The concepts of Latin-hypercube sampling, clustering, and Box–Behnken design are used in the proposed optimization routine, and algorithms such as PSO and pattern search are used for optimizing the control variables of the problem statement. The proposed optimization routine is divided into two parts, namely a surrogate optimization part and a pattern search part, which are presented as flow charts in Figures 6 and 7, respectively. The implementation of the surrogate optimization routine was carried out with MATLAB programming.

**Figure 6.** Flow chart: proposed optimization routine—part 1 (surrogate optimization).

**Figure 7.** Flow chart: proposed optimization routine—part 2 (pattern search).

The surrogate optimization part begins with the design of experiments using the Latin-hypercube sampling method. The design variables shown in Table 3 were sampled for 500 Latin-hypercube samples within the boundaries of the variable design space. These samples were simulated with FCSMEK finite element software for calculating their output responses, including electromagnetic efficiency, *η*, output power, *P*out, and power factor, *PF*. A set of criteria (satisfying the threshold values of the output responses) was devised to select a set of samples from the design of experiment. The design variables of the selected samples were standardized as shown in Equation (6), and clustered with the k-means ++ algorithm into different groups or clusters. The optimality of the number of clusters was evaluated by the MATLAB built-in function *evalclusters* using the gap statistics criterion [31,32]. The domain of the optimization problem was divided into *n* clusters and within those clusters, new boundaries of the design variables were established. The surrogate optimization process with Box–Behnken design as explained in Section 3 was applied to each cluster. At the end of the surrogate optimization, *n* surrogate optimal

solutions were obtained. The accuracy of the surrogate optimal solutions was validated with timestepping FEM simulation.

The surrogate optimal solutions were used for initializing the pattern search optimization algorithm. The goal of the pattern search algorithm is to search for an optimal solution in the vicinity of the surrogate optimal solution, which improves the objective of the optimization problem while respecting the constraints. The FEM simulation was used as a cost function for the pattern search algorithm. The final improved optimal solution was obtained at the end of the pattern search optimization process.

### **7. Results**

### *7.1. Visualization of a 3-Variable Optimization*

A three-variable optimization problem for an induction machine is demonstrated to visualize the flow of the proposed surrogate optimization routine as shown in Figure 8. The stator inner diameter, *D*s, stator yoke width, *h*ys, and slip, *s*, from Table 3 were chosen as the design variables and the remaining variables were fixed to constant values. The objective and constraints of the optimization problem are presented in Table 2. A set of 500 samples of the optimization variables were generated using the Latin-hypercube sampling method as shown in Figure 8a. These samples were simulated with FCSMEK finite element software to calculate the respective output responses. Selection criteria based on the output response of the samples as presented in Table 5 were applied to the Latin-hypercube samples to pick the samples of interest, as shown in Figure 8b. These samples were clustered into different groups as shown in Figure 8c. The Box–Behnken domain of each of the groups was as shown in Figure 8d. The surrogate optimal solution was computed and validated with FEM simulation for each of the clusters. The best solution of the surrogate optimization from one of the clusters is presented in Table 6. It is seen that the difference in the computation of the surrogate optimal solution has been reduced considerably when compared with the results from Table 4, but the constraint of the output power *P*out is not respected. The pattern search algorithm with FEM simulation as the cost function was applied to the surrogate optimal solution to search for a better solution in its vicinity. The objective and design variables of the optimal solution from the surrogate optimization part and pattern search part are compared in Tables 7 and 8. It can be seen that the pattern search algorithm found a marginally better solution in the neighborhood of the surrogate optimal solution, while respecting the constraints of the optimization problem. The electromagnetic efficiency, power factor and electromagnetic losses are compared for various load points in Tables 9 and 10.


**Table 5.** Sample selection criteria from the Latin-hypercube sampling method.

**Table 6.** Validation of surrogate optimal solution of the proposed optimization routine (part 1—surrogate optimization) with FEM simulation 3-variable optimization problem.



**Table 7.** Surrogate optimal solution compared with the improved optimal solution (part 1—surrogate optimization vs. part 2—pattern search) of the proposed optimization routine 3-variable optimization problem.

**Figure 8.** Visualization of the process involved in the proposed surrogate optimization routine.

**Table 8.** Comparison of design variables of the optimal solution (part 1—surrogate optimization vs. part 2—pattern search) from the proposed surrogate optimization routine 3-variable optimization problem.



**Table 9.** Comparison of electromagnetic efficiency and power factor for various load points of the optimal design 3-variable optimization problem.

**Table 10.** Comparison of losses for various load points of the optimal design 3-variable optimization problem.


### *7.2. Multivariate Optimization*

The proposed surrogate optimization routine was applied to a multivariate optimization problem as specified in Section 2. The selection criteria of the samples based on its output response for the clustering process were as shown in Table 5. The output response of the surrogate optimal solution from one of the clusters was validated with FEM simulation as presented in Table 11. It was found that difference in the output responses has been reduced considerably when compared with the solution presented in Table 4. The decrease in difference of the output response leads to the surrogate optimal solution respecting the set of constraints of the optimization problem. The pattern search algorithm improves the objective of the surrogate optimal solution by searching in the vicinity of the surrogate design. The objective and design variables at the end of both the surrogate optimal part and the pattern search part are compared in Tables 12 and 13. It can be noted from Table 13 that to improve the electromagnetic efficiency, *η*, of the surrogate optimal solution, the values of the design variables, such as air gap width, *δ*, stator tooth width, *b*ds, and stator yoke width, *h*ys, have changed by a small margin. The electromagnetic efficiency, power factor, and electromagnetic losses are compared for various load points in Tables 14 and 15. The flux density distribution of the optimal solution (quadrant of the optimal induction machine) at the end of pattern search algorithm is shown in Figure 9.

**Table 11.** Validation of surrogate optimal solution of the proposed optimization routine (part 1—surrogate optimization) with FEM simulationdesign–multivariate optimization problem.


**Table 12.** Surrogate optimal solution compared with the improved optimal solution (part 1—surrogate optimization vs. part 2—pattern search) of the proposed optimization routine–multivariate optimization problem.


**Table 13.** Comparison of design variables of the optimal solution (part 1—surrogate optimization vs. part 2—pattern search) from the the proposed surrogate optimization routine–multivariate optimization problem.


**Table 14.** Comparison of electromagnetic efficiency and power factor for various load points of the optimal design–multivariate optimization problem.


**Table 15.** Comparison of losses for various load points of the optimal design 3-variable optimization problem.


**Figure 9.** Flux density distribution of the optimal solution from the proposed surrogate optimization routine–multivariate optimization problem.

### *7.3. Comparison: Proposed Surrogate Optimization Routine vs. Direct Optimization Routine*

In this section, the result from the proposed surrogate optimization routine is compared with the direct optimization routine, which uses FEM simulation as the cost function. The flow chart of the direct optimization routine is presented in Figure 10. The computational cost of the direct optimization routine is high since it uses FEM simulation as the cost function. The PSO algorithm and pattern search algorithm were used in the direct optimization routine with the same configuration as that of the proposed surrogate optimization routine. Due to the high computational cost of the FEM simulation and limitation in the computational capabilities of the research workstation, the size of the population was fixed to 30 particles for the PSO algorithm in the direct optimization routine. The optimal solution from the PSO algorithm was used to initialize the pattern search algorithm, which searches for a better solution in the vicinity. The objective and design variables of the optimal solution from the proposed surrogate optimization routine are compared with the optimal solution from the direct optimization routine in Tables 16 and 17. The electromagnetic efficiency, *η*, of the optimal solution from the proposed surrogate optimization routine reached closer to that of the direct optimization routine. The marginal difference in the design variables slip, *s*, rotor slot lower width, *b*5r, and rotor yoke width, *h*yr, between the routines impacts the electromagnetic efficiency, *η*, of the optimal solutions. The electromagnetic efficiency, power factor, and electromagnetic losses are compared for various load points in Tables 18 and 19. The advantage of using the proposed surrogate optimization routine is that it requires far fewer FEM simulations than the direct optimization routine, while maintaining an accurate evaluation of the optimal design. The optimization routines were performed in the computer with dual processors of Intel Xeon Silver 4114 CPU at a clock-rate of 2.2 GHz, which can handle parallel computations of 15 FEM simulations. The comparison of the number of FEM simulations in the proposed surrogate optimization routine and direct optimization routine is presented in Table 20.

**Figure 10.** Flow chart: direct optimization routine with FEM as the cost function.


**Table 16.** Optimal design from the proposed surrogate optimization routine compared with the direct optimization routine.

**Table 17.** Design variables of the optimal solution from the the proposed surrogate optimization routine compared with direct optimization routine.


**Table 18.** Comparison of electromagnetic efficiency and power factor for various load points of the optimal design proposed surrogate optimization routine vs. direct optimization routine.


**Table 19.** Comparison of losses for various load points of the optimal design-proposed surrogate optimization routine vs. direct optimization routine.


**Table 20.** Comparison of the number of FEM simulations in proposed surrogate optimization routine and direct optimization routine.


The reliability of the proposed surrogate optimization routine for the induction machine problem was assessed with 20 continuous runs. The electromagnetic efficiency, *η*, output power, *P*out, and power factor, *PF*, of the optimal solution from each run were analyzed to provide the probability distribution as presented in Figure 11. It can be seen that all of the optimal solutions from the proposed surrogate optimization routine respect the constraints specified in the optimization problem and that the range of the objective, electromagnetic efficiency, *η*, varies between 93.75% and 93.95%.

**Figure 11.** Probability distribution of the objective and constraints from 20 runs of the proposed surrogate optimization routine.

### **8. Conclusions**

This paper proposes a novel, efficient, and reliable surrogate optimization routine that can be applied to multiple design problems. The proposed clustering technique used in the routine enables improving the accuracy of the surrogate model while exploring promising subsets of the design variable range. The surrogate optimization routine was

applied to design an optimal three-phase induction motor, maximizing its efficiency for a given volume. The surrogate functions of the electromagnetic efficiency, output power, and power factor were constructed as a function of eight design variables and these functions acted as the objective and constraints of the optimization problem. A precision of 0.01 mm was considered for the optimization process, which is possible only with laser cutting of the electrical sheets at the prototyping level. A three-variable optimization problem was performed to demonstrate the discretization of the optimization problem into a few subdomains with the clustering algorithm for searching for the optimal solution. The results of the proposed surrogate optimization routine applied to the multivariate optimization problem show an improved optimal solution when compared with optimization with a simple Box–Behnken design. This proves the improvement of the accuracy of the surrogate functions by the application of the proposed surrogate optimization routine. To benchmark the proposed surrogate optimization routine, a direct optimization routine was applied to the induction motor problem, which uses FEM simulation as a cost function. Upon comparing the results of both routines, the optimal solution from the proposed surrogate optimization routine was shown to reach closer to that from the direct optimization routine. Additionally, the proposed surrogate optimization routine used 1364 FEM simulations compared with 75,208 FEM simulations of the direct optimization routine, thus greatly improving the computational efficiency. Future work on the proposed surrogate optimization routine will focus on performance evaluations on different types of machines and its application for multi-objective optimization problems with several constraints.

**Author Contributions:** Conceptualization, A.B. (Aswin Balasubramanian); methodology, A.B. (Aswin Balasubramanian), F.M. and A.B. (Anouar Belahcen); software, A.B. (Aswin Balasubramanian), F.M. and A.B. (Anouar Belahcen); validation, A.B. (Aswin Balasubramanian) and F.M.; writing—original draft preparation, A.B. (Aswin Balasubramanian); writing—review and editing, A.B. (Aswin Balasubramanian), F.M., O.O., M.M.B. and A.B. (Anouar Belahcen); visualization, A.B. (Aswin Balasubramanian) and M.M.B.; supervision, A.B. (Anouar Belahcen); project administration, A.B. (Anouar Belahcen); funding acquisition, A.B. (Anouar Belahcen). All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the Academy of Finland consortium grant 330747.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We would like to acknowledge our colleague Alireza Nematsaberi, Department of Electrical Engineering and Automation, Aalto University, for his contribution in brainstorming the idea of this article, and Devi Geetha Nair, R&D Engineer, ABB Oy, Finland for proofreading the article.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A**

**Table A1.** Quadratic closed form equation of Box–Behnken design—constant and first order terms.


**Table A2.** Quadratic closed form equation of Box–Behnken design—second-order squared terms.


**Table A3.** Quadratic closed form equation of Box–Behnken design—second-order product terms.



**Table A3.** *Cont.*

### **References**

