Multi-Objective Electromagnetic Design Optimization of a Power Transformer Using 3D Finite Element Analysis, Response Surface Methodology, and the Third Generation Non-Sorting Genetic Algorithm

Abstract

This paper presents a multi-objective design optimization of a power transformer to find the optimal geometry of its core and the low- and high-voltage windings, representing the minimum power losses and the minimum core and copper weights. The optimal design is important because it allows manufacturers to build more efficient and economical transformers. The approach employs a manufacturer’s design methodology, which is based on the usage of the laws of physics and leads to an analytical transformer model with the advantage of requiring a low amount of computing time. Afterward, the multi-objective design optimization is defined along with its constraints, and they are solved using the Non-Sorting Genetic Algorithm III (NSGA-III), which finds a set of optimal solutions. Once an optimal solution is selected from the Pareto front, it is necessary to fine-tune it with the 3D Finite Element Analysis (FEA). To avoid the large computing times needed to carry out the 3D Finite Element (FE) model simulations used in multi-objective design optimization, Response Surface Methodology (RSM) polynomial models are developed using 3D FE model transformer simulations. Finally, a second multi-objective design optimization is carried out using the developed RSM empirical models that represent the cost functions and is solved using the NSGA-III. The numerical results of the optimal core and windings geometries demonstrate the validity of the proposed design methodology based on the NSGA-III. The used global optimizer has the feature of solving optimization problems with many cost functions, but it has not been applied to the design of transformers. The results obtained in this paper demonstrate better performance and accuracy with respect to the commonly used NSGA-II.

1. Introduction

The power transformer is an important element in an electrical network because it makes possible the transmission and distribution of electrical energy. Even though a transformer is passive equipment, i.e., with no rotating parts, it may be exposed to malfunctions, and its failure can generate energy interruptions that directly impact consumers and cause economic losses to the electrical energy companies. A transformer’s electrical and magnetic performance is governed by Maxwell equations [1]; it can also be described in a simpler way by using Faraday’s induction equation and equivalent electrical and magnetic circuit theory [2]. By employing simpler equations than those of the electromagnetic theory, transformer manufacturers have developed their own proprietary design computer programs, which they have improved over the years based on their experience manufacturing many transformers. Alternatively, Maxwell equations can be numerically solved either in 2D or 3D, giving more accurate results than the manufacturer’s design method, but this approach has the disadvantage of having high computing times when used directly in an optimization process; the use of symmetries can significantly reduce the computing time [3,4]. The manufacturer’s design approach has the drawback of being less accurate, but, due to its simplicity, it has the advantage of being quickly solved and having fast designs. In the design process of a transformer, it is desirable to have the best possible design that can fulfill constraints, and this can be achieved by employing optimization algorithms. Initially, determinist optimization algorithms were developed and used in designing transformers. In the last two decades, applying global optimization algorithms to improve the design of power transformers has been a topic of interest to academic and industrial researchers. Recently, global optimization algorithms that can handle multiple objectives have also been used in engineering [5]. Amoiralis et al. proposed a global transformer optimization based on decision trees and adaptively trained neural networks, which were combined with the Finite Element Method (FEM) and Boundary Element, where the selection of winding materials was achieved [6]. Amoilaris et al. also used a parallel mixed-integer programming technique combined with the branch and bound algorithm for the global design optimization of several transformer ratings, where the cost of its active parts was chosen as the objective function to be minimized [7]. In 2009, Amoilaris et al. also published a complete bibliographic survey related to transformer design carried out in the previous forty years, where they reported the appearance of stochastic methods that included artificial intelligence techniques and genetic algorithms [8]. Arjona et al. developed a methodology that they applied to the optimal design of a power transformer, where they combined the conventional transformer design with genetic algorithms at the start of the design process and then used a 2D Finite Element (FE) transformer model with the Gauss-Newton method [9]. Tran et al. proposed an adapted output space-mapping technique with a bi-objective optimization based on a 3D FE model of a single-phase transformer; the Pareto front was constructed for the loss and total mass of the transformer; this technique reduces the computation time needed when using direct 3D FE simulations [10]. Cheema et al. published a practical approach for the global optimization of a three-phase core-type distribution transformer design where the total owning cost was considered as the cost function, which was minimized using four different algorithms: direct search, differential evolution, simulated annealing, and random search [11]. Hultman Ayala et al. proposed a multi-objective optimal design for a safety insulating transformer where the metaheuristic wind-driven optimization algorithm was employed and the minimization of two cost functions was used: the total mass, iron loss, and copper loss; no Finite Element Analysis (FEA) was applied [12]. Xia et al. proposed an optimal design applied to the winding transposition of power transformers, where an adaptive co-Kriging surrogate model was constructed with 2D and 3D FEA and the binary particle swarm optimization algorithm was used; the single objective of circulating current winding loss was minimized [13]. Mohammed et al. proposed the use of evolutionary algorithms and FEA for the multi-objective transformer design where total loss and total ownership cost were defined as the cost functions to be minimized; the Non-Sorting Genetic Algorithm (NSGA) II, the genetic algorithm, and the differential evolution algorithm were used [14]. Liu et al. applied the NSGA-II to the optimal design of a distribution transformer with an amorphous core; they considered the total owning cost as the objective function to be minimized along with the transformer constraints [15]. Orosz et al. proposed an optimal design of a power transformer, where NSGA-II was applied to minimize the total cost of ownership and the geometric programming method was used to calculate the optimal winding layout of the transformer; 2D FEA was included in the optimization process [16]. Recently, a thermal analysis based on the electromagnetic and fluid dynamic modeling of a power transformer was presented by Stebel et al. [3].

Metamodeling and global optimization have also been applied to the optimal design of low-power transformer applications. Cove et al. applied response surface modeling and FEA to design small planar transformer windings, where heuristic models were developed for predicting leakage and magnetizing inductances, capacitances, and winding resistances; optimization algorithms were not applied [17]. Li et al. presented the optimal design of a high-frequency coaxial transformer where artificial neural networks and FEA were used; the design was focused on determining the ring core size, the winding structure using Litz wires, and the working parameters [18]. Xiaowei et al. used the NSGA-II for the multi-objective optimal design of a high-frequency transformer. Two cost functions were defined: the first one is related to the product of the magnetic window area and cross-section area of the core, and the second function represents the losses of the core and windings [19]. Tan et al. conducted a multi-response optimization of the thermo-hydraulic performance of a transformer [20].

The above-published research results demonstrate that the transformer industry can benefit from using powerful evolutionary multi-objective algorithms to design transformers that can fulfill complex goals and constraints. The advent of optimization algorithms that can solve real optimization problems with many objectives can help electrical engineering designers [21,22]. Multi-objective problems with at least four objectives are known as many objectives, even though some authors consider problems with at least three objectives to be many-objective problems [23]. Numerical methods have become mature transformer designer tools like the FEM [1]. Furthermore, multi-physics problems can be numerically solved to obtain an optimal solution to complex problems. In addition, the role of statistics in process optimization by using the design of experiments and surrogate modeling is also important [24]. Nonetheless, few publications related to multi-objective optimization of power transformers have been reported, and none of them have applied the NSGA-III algorithm, whose application is being proposed in this paper.

This paper presents an optimal design to find the geometry of the core and the high- and low-voltage windings of a power transformer. The multi-objective electromagnetic design minimizes the transformer’s electrical power losses and the copper and core weights with two equality constraints related to the apparent output power and the leakage reactance of the transformer. The problem is solved using the many-objective optimization algorithm NSGA-III, an analytical transformer model, Response Surface Methodology (RSM) polynomial models, and 3D FEA. An analysis of the optimal solution, convergence performance using the running metric, computing time, and the accuracy obtained with the optimizers NSGA-II and NSGA-III is reported.

2. Multi-Objective Transformer Design Methodology

The proposed analytical approach is illustrated in Figure 1. It starts by using a manufacturer’s approach to design a transformer. The procedure begins by defining the technical specification of the device (

Table 1. Transformer design specification.

Power Rating (kVA)	High Voltage (V)	Low Voltage (V)	Type	Connection	Frequency (Hz)	Reactance (%)
1500	13,200	220	Oil filled	Delta-wye	60	5.5

). Afterward, a multi-objective optimization is formulated and solved using the NSGA-III algorithm, which can handle any number of goals defined by the designer. When this step ends, one of the best optimal solutions must be selected and considered in the second step, where FEA is employed to solve a more detailed transformer model [25]. In the next step, a RSM polynomial model is constructed; this model is obtained through the design of experiments, where the FEA provides information on the desired performance of the transformer. Finally, a multi-objective optimization is carried out using the RSM models, which can speed up the numerical solution of the multi-objective optimization. The last step may also include verifying the selected optimal solution against FE simulations of the power transformer.

3. Transformer Modeling

The transformer modeling consists of simplified analytical closed-form equations, electromagnetic formulation, and response surface methodology.

3.1. The Manufacturer’s Transformer Design Approach

The initial problem in designing a transformer consists of determining its appropriate dimensions and performance. A summary of the transformer’s performance can be described by Equations (1)–(5) [26].

$$V=4.44Nf{B}_c{A}_c,$$

(1)

$$S=3VI,$$

(2)

$$P_{Cu}=\raisebox{1ex}{$J^2{\rho }_{Cu}{m}_{Cu}$}\!\left/ \!\raisebox{-1ex}{${\rho }_{dCu}$}\right.,$$

(3)

$$P_{core}=V_c{\rho }_c{P}_c,$$

(4)

$$X=\frac{{\left(2\pi \right)}^2{\mu }_{0}f{\left(VI\right)}_{base}}{{\left(\raisebox{1ex}{$V$}\!\left/ \!\raisebox{-1ex}{$N$}\right.\right)}^2{h}_m}\left(\frac{r_{m1}{b}_1}{3}+\frac{r_{m2}{b}_2}{3}+r_{mg}g\right),$$

(5)

where V is the winding voltage, N represents the winding turns, f stands for frequency, B_c is the flux density, A_c is the core area, S denotes apparent power, and I is the phase current. P_Cu stands for cooper loss, J is the winding current density, ρ_Cu is the copper resistivity, m_Cu is the winding weight, and ρ_dCu denotes the density of cooper. P_core represents the core loss, ρ_c is the density of oriented steel, P_c indicates specific iron loss, and V_c is the core volume. X represents the transformer reactance, h_m is the winding mean height, r_m denotes mean radius, b is radial winding depth, r_mg is the airgap mean radius, and g is the airgap length between primary and secondary windings. Subscripts c, 1, and 2 indicate the core, low voltage, and high voltage windings, respectively.

The transformer schematic diagram with the main dimensions of its active elements is presented in Figure 2; this information is the output of the transformer analytical model used by manufacturers. Figure 2 represents the geometry considered in the design of the power transformer, where a set of distances and diameters related to the electrical clearances of the main transformer’s electrical insulation and the magnetic circuit, along with its windings, can be seen. A side view of the transformer is shown in Figure 2a, illustrating the height, radial thickness, and inner diameters of the low- and high-voltage windings. Additionally, the core dimensions are indicated, such as its diameter and the core leg height; the distance between the centers of the core legs can be obtained from the winding’s dimensions and their insulation clearances. The gaps related to the main insulation of the transformer, such as the distance between the core and tank, the electrical clearance between the HV winding and tank, the gap between the high voltage windings, the distances between the yokes and tank, and the windings and core, are also illustrated. The expressions employed for the calculation of the height and width of the tank are indicated in Figure 2a. The top schematic view of the transformer is presented in Figure 2b, where the dimensions of the HV and LV windings are seen and the electrical insulation distance between the windings and core is also indicated. The algebraic expression employed to calculate the tank depth is also shown; this distance depends on the core, windings, and electrical clearance dimensions. The description of the dimensions indicated in Figure 2 is given in

Table 2. Description of the main transformer dimensions of Figure 2.

Dim.	Description	Dim.	Description	Dim.	Description
d0	Core diameter	h	Core leg height	h3	Yoke bottom tank clearance
d1	LV winding diameter	g	HV–LV windings gap	h4	Yoke top tank clearance
d2	HV winding diameter	g0	LV–Core clearance	tw	Tank width
h1	LV winding height	g1	LV–Yoke clearance	th	Tank height
h2	HV winding height	g2	HV–Yoke clearance	td	Tank depth
b1	LV winding radial depth	g3	HV–HV clearance
b2	HV winding radial depth	g4	HV–Tank clearance

.

3.2. Electromagnetic Model

The 3D magnetodynamic formulation includes the variation of the magnetic scalar potential (Ω) and the electric vector potential in conducting regions. The time variation is assumed to be sinusoidal, and complex numbers are used to represent it. The electric vector potential T is related to the vector current density J and to the magnetic field intensity H as follows (6)–(9) [1]:

$$\nabla \times \mathit{T}=\mathit{J},$$

(6)

$$\left(\mathit{T}-\nabla \mathsf{\Omega }\right)=\mathit{H}.$$

(7)

By using the Faraday and Gauss laws, the induced current in conducting regions is given by (8):

$${\nabla }^2\mathit{T}=-\sigma {\mu }_0{\mu }_r\left(\mathit{T}-\nabla \mathsf{\Omega }\right),$$

(8)

where σ is the electric conductivity, μ₀ is the free space permeability, and μ_r is the relative permeability. In other regions, the scalar magnetic potential is used (9).

$${\nabla }^2\mathsf{\Omega }=0.$$

(9)

3.3. Response Surface Methodology

Response Surface Methodology uses statistics and mathematics to develop heuristical models that are used to analyze and improve physical processes [24]. It has been applied to the industrial world, such as chemical, automotive, electromagnetics, and electronics. RSM allows the determination of models, which can be of the type of polynomial functions, Kriging, or neural networks. RSM analyzes the effect of several inputs (also known as factors, independent variables, or design variables) on the response of a system. This methodology has the advantage that it can represent systems with multiple inputs and multiple outputs, and hence can be applied to complex processes such as an electrical power transformer. The resulting models of applying RSM are not based on physical principles; they are also known as metamodels; hence, the methodology can be applied to different areas of engineering. RSM is mainly applied as a sequential procedure, as shown in Figure 3a. The first step consists of defining the design variables that may impact the outputs of interest; this phase may be based on the expertise of the design engineer who knows the process response. In the work presented here, the design factors considered in the RSM are the core diameter, core leg height, low-voltage winding radial thickness, high-voltage winding radial thickness, and the gap radial length between the windings. The second RSM step is related to determining the values of the design variables where the optimum or near-optimum operating point lies; this procedure in RSM is commonly carried out using deterministic optimization using the RSM-derived first-order polynomial models. However, in this paper, the optimum operating point of the power transformer is obtained using a global optimizer and an analytical model of the transformer, which was derived from the manufacturer’s design experience. The third step of a response surface is to find a heuristic model that represents the response with higher accuracy in a small region near the optimum. In this stage, the design of experiments is used to collect FEA-derived data, which are then used to develop RSM-based polynomial models. There are several methods for carrying out the experiments for data collection, such as the Central Composite Design (rotatable, face-centered, and spherical), Box-Behnken, full factorial, and fractional factorial, among others. The CCD rotatable fractional factorial design of experiments is used in this paper for two main reasons: It allows for a reduced set of experiments, and it has five-level design variables (−α, −1, 0, +1, +α), where the +α and −α represent the maximum and minimum values of each input parameter (Figure 3b), and level 0 is the mean value of that input. The center point and the axial points located at distances ±α from the center point allow the estimation of curvature and quadratic terms. The number of design points in the experiment is calculated as 2^(k−f) + 2k + 1, where k is the number of factors and f is the fraction of the full factorial experiment [27]. Furthermore, the rotatable CCD has the property of having an equal predicted variance for the samples that are equidistant from the center of the design. Furthermore, the larger value of α allows for including a larger region of the operating point, and it is calculated as α = 2^(k−f)/4. An empirical second-order model is used because it is recommended when the response near the optimum usually exhibits a curvature. In the fourth RSM step, the resulting model is analyzed to determine the optimum values of the factors that represent the core and winding dimensions of the power transformer. In the final step, a confirmatory experiment must be carried out, in which the objective is to confirm the identified optimum operating point. In RSM, a second-order polynomial model is given by (10) [24].

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

(10)

where x’s represent the input coded variables, $\beta$’s are the unknown parameters, and ε stands for the prediction error. The second term in the Equation (10) indicates the main effects, the third term shows the quadratic components, and the fourth term shows the interactions for factor x_i.

4. Multi-Objective Optimization

Genetic algorithms are search and optimization tools that work differently compared with classical methods and are capable of handling multi-objective problems with a wide scope of applications [28]. Since their introduction by Holland in 1975, they have evolved, passing through different extensions: real-coded, multimodal, multi-objective, and many-objective, which led them to solve different optimization problems. Genetic algorithms search simultaneously for multiple solutions on discontinuous regions next to the approximate Pareto front. Deb and Jain proposed a many-objective evolutionary algorithm whose basic framework remains similar to the NSGA-II algorithm [21]. The main difference is the substitution of crowding distance for a selection based on well-distributed reference points. Reference points of the original NSGA-III are uniformly distributed on a hyperplane to guide solutions to converge.

The NSGA-III algorithm starts with an initial population of feasible solutions of size N and a set of reference directions on a unit simplex. One of the requirements of this algorithm is the initial supply of a set of reference directions. A reference direction is constructed by a vector starting from the origin and connected to each of the reference points that belong to a set of points initialized on an M-dimensional unit simplex, where M is the number of objectives. Das and Dennis’s method is one of the most commonly used methods to place reference points [29]. However, with Das and Dennis’s method for a large number of objectives, very few points of a reasonable-sized reference set lie in the interior of the simplex. Blank et al.’s Riesz s-Energy method was used here to create a well-spaced set of points on the unit simplex [30], whose concept is to find the z-matrix that minimizes the s-Energy function subject to every z⁽ⁱ⁾ vector lying on the unit simplex (11):

$${\sum }_{i=1}^M{z}_m^{\left(i\right)}=1,$$

(11)

Inverted Generational Distance (IGD) is a performance indicator used to evaluate the quality of an obtained solution set in comparison with a pre-specified reference point set [31,32]. One of the issues with the IGD performance indicator is that it requires a set of solutions from the true Pareto-optimal (PO) front, which makes it difficult to apply to a real-world problem where PO solutions are not available before an evolutionary multiobjective optimization algorithm is run. Blank and Deb updated to IGD metrics with no use of PO front [33]. They keep track of how the realized ideal point (z*) and realized nadir point (z^nad) move from one generation to another, considering the maximum absolute difference of each component. The first one is a vector with the minimum of all objectives, and the latter is a vector with the maximum of non-dominated solutions in the objective space. The non-dominated solutions are accumulated from the initial generation to the current generation (τ) and then z* and z^nad points are calculated from the accumulated set. Afterwards, the normalized i-th objective value of the j-th non-dominated point at the t-th generation $P_i^{\left(j\right)}\left(t\right)$ using $P\left(\tau \right)$ for (0 ≤ t ≤ τ) is calculated by:

$${\overline{P}}_i^{\tau ,\left(j\right)}\left(t\right)=\frac{P_i^{\left(j\right)}\left(t\right)-z_i^{*}\left(\tau \right)}{z_i^{nad}\left(\tau \right)-z_i^{*}\left(\tau \right)}$$

(12)

Blank and Deb [33] state that once the non-dominated sets at generations 0 ≤ t ≤ τ are normalized with fixed z* and z^nad points, a performance metric that requires a reference set P* and an evolving non-dominated set Q(t) can now use ${\overline{P}}^{\tau }\left(\tau \right)$ and ${\overline{P}}^{\tau }\left(t\right)$, respectively, to evaluate the performance of the algorithm for the mentioned generations. Then, the IGD metric for 0 ≤ t ≤ τ can be calculated using (13):

$$IGD\left({\overline{P}}^{\tau }\left(t\right),{\overline{P}}^{\tau }\left(\tau \right)\right) \mathrm{for} 0\le t\le \tau ,$$

(13)

where t represents the t-th generation, τ stands for the accumulated non-dominated solutions from the initial to the current generation, ${\overline{P}}^{\tau }\left(t\right)$ is an evolving non-dominated set, and ${\overline{P}}^{\tau }\left(\tau \right)$ represents the reference set. To reduce the computational complexity and IGD computations, Blank and Deb [33] calculate the average improvement of the IGD metric over generations $\left(t-1\right)$ to $t$ as follows (14):

$${\mathsf{\Phi }}_t=\mathrm{IGD}\left({\overline{P}}^t\left(t-1\right),{\overline{P}}^t\left(t\right)\right),$$

(14)

The demonstration of the proposed new design methodology is given in this section. First, a multi-objective optimization is mathematically formulated based on the design approach used by manufacturers, where a cost vector function, f(x), is solved with the NSGA-III algorithm as indicated by (15)–(21).

$$\min \mathit{f}\left(\mathit{x}\right)=\left\{f_1\left(\mathit{x}\right),f_2\left(\mathit{x}\right),f_3\left(\mathit{x}\right)\right\},$$

(15)

s.t.

$$f_1\left(\mathit{x}\right)={\rho }_{dCu}{V}_{Cu},$$

(16)

$$f_2\left(\mathit{x}\right)={\rho }_c{V}_c,$$

(17)

$$f_3\left(\mathit{x}\right)=P_{Cu}+P_{core},$$

(18)

$$S\left(\mathit{x}\right)=S_{rated},$$

(19)

$$X\left(\mathit{x}\right)=X_{ref},$$

(20)

$${\mathit{x}}_l\le \mathit{x}\le {\mathit{x}}_u,$$

(21)

where x is the design variable vector, f₁ is the copper weight of the windings, f₂ stands for core weight, f₃ represents the copper and core electrical power losses, S(x) is the apparent electrical power, S_rated is the desired transformer’s apparent power, X(x) is the transformer’s reactance, and X_ref represents the specified transformer’s reactance. The copper and core volumes are represented by V_Cu and V_c, respectively. The subscripts l and u mean the lower and upper bounds, respectively.

5. Numerical Results

The population size and number of generations for the NSGA-III optimization power transformer problem with the three objectives and two equality constraints are 200 individuals and 1000 generations.

5.1. Optimization Results with the Manufacturer’s Transformer Design Approach

The Pareto front with three objectives minimized with the NSGA-III, i.e., f₁(x) for copper weight, f₂(x) for core weight, and f₃(x) for copper and core power losses, are shown in Figure 4 in normalized values. The used vector of design variables is defined as x = [B, J_HV, J_LV, d0, h, g]. The blue circles indicate the Pareto front with all possible optimal solutions, whereas the cyan region indicates the plane of the reference points. As mentioned in the previous section, the performance of the applied optimization was evaluated by using the running performance metric proposed by Blank and Deb [33]. The convergence of the multi-objective solutions plotted in Figure 5 sets out two phases of the running performance metric for the transformer manufacturer’s design approach. The criteria employed to choose the best solution are based on selecting the optimum operating point with the lowest power losses with the minimum weight of copper and grain-oriented silicon steel, such that it represents the lowest cost considering a price of 15.15 USD/Kg for the copper and 6 USD/Kg for the core. The selected optimal point obtained with the NSGA-III is shown in

Table 3. Optimal solution obtained with the transformer manufacturer’s design approach.

Design Variables and Transformer Quantities	Description	Optimum NSGA-II	Optimum NSGA-III	FEA of Optimum NSGA-III	Units
B	Flux density	1.48	1.69	-	T
JHV	HV current density	3.36	3.20	-	A/mm2
JLV	LV current density	3.25	2.81	-	A/mm2
d0	Core diameter	255.83	239.08	-	mm
h	Core leg length	457.60	493.21	-	mm
g	HV–LV winding gap	12.00	14.82	-	mm
S	Apparent power	1499.72	1499.84	1501.09	kVA
X	Leakage reactance	5.54	5.5	5.86	%
Ploss	Power losses	14,751	13,740	13,248	W
WCu	Copper weight	515.88	548.40	548.3	Kg
Wcore	Core weight	1261.15	1113.19	1124.9	Kg

.

Figure 5 shows the reference points used with NSGA-III to ensure diversity in the obtained solutions. Blank et al.’s Riesz-Energy method was used to place points on a normalized hyper-plane that is equally inclined to all objective axes and has an intercept of one on each axis [30]. The blue circles correspond to the obtained optimal solutions.

The running metric is used to visualize the performance of the NSGA-III algorithm. Figure 5 shows the proceedings with an interval of 50 generations. To compute the blue line in the early phase, non-dominated sets for all past generations where 0 ≤ t ≤ τ with τ = 50 are normalized using Equation (12), then IGD values are computed with Equation (13). After 50 more generations, the IGD plot for τ = 100 is shown in orange. The IGD values in the first 50 generations are different from the ones of the plot in the blue line because all non-dominant sets are re-normalized with z* and z^nad points from the 100th generation. From the final phase, it can be seen how NSGA-III improves and the IGD value approaches zero.

A comparison of the performance of both NSGA-II and NSGA-III in terms of the running metric is presented in Figure 6. It is clearly seen how the NSGA-III shows better convergence performance with an optimization problem defined by three objectives and two equality constraints. It is seen that the NSGA-II IGD metric does not get better values as the number of generations increases. Meanwhile, the NSGA-III running metric converges to lower values after generation 200. This demonstrates that NSGA-III has better convergence performance than NSGA-II.

The optimum solutions obtained with the NSGA-II and NSGA-III, using the analytical model for the transformer design, are shown in Table 3. The optimum dimensions of the core and the best values of the current densities of the LV and HV transformer windings are given in Table 3. The winding areas can be calculated using the rated transformer current and these optimal current densities. The core leg height defines the height of both LV and HV windings because they are calculated considering the insulation gaps between the core and the windings. In Table 3, it is seen that the current density values have relatively high values, as is expected because of the power loss minimization problem, and at the same time, the optimizer minimizes the copper weight to be used in the transformer. Alternatively, it can be observed that the core flux density has a value of 1.69 T (NSGA-III), which means better utilization of the core by minimizing the core weight. It can also be seen that the core diameter and the core leg height have values that the optimizer finds to satisfy both the cost functions and equality constraints. It is important to note that the independent variables, core diameter and core leg height, directly impact the windings’ volume because their diameters and heights depend on these two independent variables. Regarding the fulfillment of the equality constraints, Table 3 shows values that are very close to the desired apparent power and leakage reactance. However, it should be taken into account that the analytical model of the transformer does not have high accuracy when compared with a 3D FE model, yet this solution is important because it defines the region to be analyzed using RSM (Figure 2a). Regarding the accuracy of the optimal solution (NSGA-III) using the analytical model shown in Table 3, a 3D FE simulation was carried out using the optimal transformer dimensions; the results are shown in the same table. It is seen that the transformer leakage reactance is higher, i.e., 5.86%, and the apparent power is higher by 1.09 kVA than the specified transformer power. Differences in the power losses, copper, and core weights were also found, and these are due to the higher accuracy of the FE model when compared with the analytical model.

A comparison of the solutions obtained with NSGA-II and NSGA-III is also shown in Table 3. It is seen that the power loss is lower with NSGA-III than with NSGA-II; this means that the designed transformer will have a higher efficiency. Alternatively, the copper weight is higher, and the core weight is lower. It is important to point out that these optimum points belong to the Pareto front. Regarding the fulfillment of the equality constraints, NSGA-III has more accuracy in the leakage reactance than NSGA-II; meanwhile, the apparent transformer power is almost the same for both optimizers. The power loss and the core weight are lower. The differences may be attributed to the difference in convergence as indicated by the running metric shown in Figure 6. The computing time for the NSGA-II is 5349 s, while for the NSGA-III it is 3200 s using a computer with a CPU i7-7700 at 3.6 GHz and 16 GB of RAM; hence, it was found that the NSGA-III is faster.

5.2. Multi-Objective Optimization Results Using RSM Polynomial Models

In the development of the empirical models using RSM, the optimum point obtained with NSGA-III, shown in Table 3, is employed to construct the 3D FE model of the power transformer, which is voltage fed; this means that the current in its windings is dependent on the electrical load connected to the LV winding. Therefore, the radial thickness of both LV and HV windings and the core leg height were selected as the design variables. The windings’ height is a function of the core leg height and the insulation gaps. Choosing the above three independent variables is equivalent to selecting the current density in the windings. The additional independent variables are the core diameter and the air gap between the HV and LV windings. The core flux density is dependent on the voltage source of the FE 3D transformer model. The used vector of design variables is defined as x = [d₀, g, b₁, b₂, h]. The values of the five independent variables (d₀, g, b₁, b₂, and h) and their coded values, based on the rotatable CCD experimental design, are shown in

Table 4. Design variables used in the experimental design.

Levels	Core Diameter (mm), d0	HV–LV Gap (mm), g	LV Radial Depth (mm), b1	HV Radial Depth (mm), b2	Core Leg Height (mm), h
−α (min)	215.17	13.34	32.35	28.32	443.89
−1 (low)	227.13	14.08	34.14	29.90	468.55
0 (center)	239.08	14.82	35.94	31.47	493.21
1 (high)	251.03	15.56	37.74	33.04	517.87
+α (max)	262.99	16.30	39.53	34.62	542.53

.

Table 5. Rotatable CCD experimental design matrix and FEA transformer responses for each run.

Run	d0 (mm)	g (mm)	b1 (mm)	b2 (mm)	h (mm)	Power Loss (W)	Copper Weight (Kg)	Core Weight (Kg)	Sout (kVA)	X (%)
1	239.08	14.82	35.94	31.47	493.21	13,248.4	548.3	1124.9	1501.09	5.86
2	215.17	14.82	35.94	31.47	493.21	13,421.7	507.8	874.6	1447.99	5.41
3	262.99	14.82	35.94	31.47	493.21	12,997.5	588.8	1415.7	1501.85	6.26
4	239.08	13.34	35.94	31.47	493.21	13,248.6	546.1	1121.3	1502.08	5.62
5	239.08	16.3	35.94	31.47	493.21	13,245.5	550.5	1128.5	1500.71	6.09
6	239.08	14.82	32.35	31.47	493.21	13,243.9	513.3	1116.1	1501.03	5.6
7	239.08	14.82	39.53	31.47	493.21	13,232.8	584.0	1133.6	1500.86	6.12
8	239.08	14.82	35.94	28.32	493.21	13,237.6	517.9	1117.2	1502.15	5.69
9	239.08	14.82	35.94	34.62	493.21	13,256.9	579.2	1132.5	1500.26	6.01
10	239.08	14.82	35.94	31.47	443.89	13,155.9	485.3	1079.9	1500.16	6.53
11	239.08	14.82	35.94	31.47	542.53	13,325.3	611.3	1169.9	1502.74	5.3
12	227.13	14.08	34.14	29.90	517.87	13,427.9	523.8	1006.1	1495.79	5.06
13	251.03	14.08	34.14	29.90	468.55	13,064.5	503.2	1229.2	1502.40	6.04
14	227.13	15.56	34.14	29.90	468.55	13,321.0	468.9	968.8	1494.00	5.86
15	251.03	15.56	34.14	29.90	517.87	13,135.8	566.7	1282.8	1502.96	5.68
16	227.13	14.08	37.74	29.90	468.55	13,321.3	499.0	973.4	1494.08	5.89
17	251.03	14.08	37.74	29.90	517.87	13,143.3	602.7	1288.5	1502.91	5.69
18	227.13	15.56	37.74	29.90	517.87	13,380.8	561.9	1017.3	1493.31	5.54
19	251.03	15.56	37.74	29.90	468.55	13,066.3	539.3	1242.8	1501.34	6.58
20	227.13	14.08	34.14	33.04	468.55	13,322.8	494.5	972.4	1492.01	5.78
21	251.03	14.08	34.14	33.04	517.87	13,136.8	597.5	1287.3	1503.22	5.59
22	227.13	15.56	34.14	33.04	517.87	13,429.2	557.4	1016.3	1495.18	5.44
23	251.03	15.56	34.14	33.04	468.55	13,060.7	534.6	1241.6	1501.53	6.47
24	227.13	14.08	37.74	33.04	517.87	13,424.5	591.3	1020.9	1492.19	5.46
25	251.03	14.08	37.74	33.04	468.55	13,067.2	567.1	1247.3	1501.34	6.5
26	227.13	15.56	37.74	33.04	468.55	13,336.8	529.2	983.6	1493.92	6.3
27	251.03	15.56	37.74	33.04	517.87	13,161.2	638.7	1300.9	1502.26	6.09

shows the experimental design using the fractional rotatable CCD for the five factors (d₀, g, b₁, b₂, and h) and the five model outputs (power loss, copper weight, core weight, apparent power, and leakage reactance). The combinations of the input variables, the model responses, and the runs, executed using 3D FE transformer models, are also shown in Table 5. It is necessary to point out that variations of the five design variables modify the dimensions of the transformer tank because the insulation clearances must be taken into account. The matrix information in Table 5 is used to determine the second-order polynomial models. The three output functions (power loss, copper weight, and core weight) represent the three objectives of the multi-objective optimization problem. In addition, there are two output functions (apparent power and leakage reactance). The total computing time to obtain the RSM models, including the 3D FE transformer simulations, was 2061 s. A workstation with two Xeon E5-2695 2.3 GHz processors, 28 cores, and 128 GB of RAM was used. The full set of RSM second-order polynomial functions is then employed to find the optimal value of the transformer’s core and windings dimensions using the optimizer NSGA-III (

Table 6. The optimal solution obtained using the RSM polynomial model and NSGA-III.

Design Variables and Transformer Quantities	Description	Optimum Values Using RSM Models	FEA Validation	Units
d0	Core diameter	232.49	-	mm
g	HV–LV windings gap	14.51	-	mm
b1	LV winding radial depth	33.10	-	mm
b2	HV winding radial depth	28.57	-	mm
h	Core leg height	478.52	-	mm
S	Apparent power	1499.92	1500.22	kVA
X	Reactance	5.52	5.51	%
Ploss	Power losses	13,287.5	13,280.47	W
WeightCu	Copper weight	466.1	465.9	Kg
Weightcore	Core weight	1025.5	1025.4	Kg

).

The 3D representation of the surface response of the RSM polynomial models as a function of the combination of some independent variables is shown in Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. The surface response of the transformer power losses is illustrated in Figure 7a, where it is seen that the power losses decrease as the core leg height (h) and core diameter increase. The increase of the core diameter (d₀) makes the inner diameters of the LV and HV transformer windings increase the copper volume, while an increase in the core diameter generates an increment in the core cross-section, which originates a decrease in the magnetic flux density and hence the core power loss. Figure 7b shows how the transformer power loss increases with an increment of the core leg height, which means a core volume increment and, therefore, higher core losses. Alternatively, the increase in the radial thickness (b₁) of the LV transformer winding will increase the conductor’s area, and consequently, the LV current density increases. A similar response can be seen in Figure 7c using the factors h and b2. Figure 7d shows a slight variation in the power losses when the HV–LV gap (g) varies. Dots indicate the FEA results in Figure 7a–d, where they are close to the response surfaces.

Figure 8 shows the objective function that represents the copper weight obtained by RSM. Figure 8a shows the response of the copper weight, which increases when the core leg height increments; therefore, the transformer windings’ heights (h₁ and h₂) will also increase, as shown in Figure 2a. The surface model also shows how the copper weight increases with the increment of the core diameter (d₀); this will generate an increment in the inner diameters of the LV and HV transformer windings. Figure 8b illustrates the response of the function that represents the copper weight, which increases with the increment of the core leg height (h), which represents an increase of the low voltage winding (h₁) and also of the high voltage winding (h₂) of the power transformer. As the LV winding thickness increases, the copper volume increases. A similar response on the model’s surface is seen when the high-voltage winding thickness varies (Figure 8c). Figure 8d shows the model’s output, represented by the copper weight, as a function of the core leg height as the HV–LV gap (g) varies; an increase in the copper weight is seen as both design variables increase. Finally, it is seen that surface responses are close to the 3D FE transformer model results (indicated by dots).

The cost function, representing the core weight, is shown in Figure 9. It is seen that an increase in the core leg height (h), along with an increment of the core diameter (d₀), generates a higher volume of the grain-oriented silicon steel and hence a higher core weight (Figure 9a). While the core weight response exhibits slight variation with changes in the independent variables b₁, b₂, and g, as can be seen in Figure 9b–d, respectively, this is an expected result due to the changes in the winding dimensions. Finally, it is observed that the RSM empirical model response closely follows the dots (3D FE transformer model) and with acceptable accuracy in Figure 9b–d.

Figure 10 illustrates the output of the RSM polynomial model that represents the apparent power (S), which is used in the equality constraint of the multi-objective optimization problem. Figure 10a shows that when the core diameter decreases, a decrease in the apparent power is obtained due to the reduction in the core area that leads to a lower magnetic flux, which impacts the induced voltage and, consequently, the output power. While a little effect is seen on the response (S) due to the changes in factor h with variations of radial depth of the LV and HV windings (b₁, b₂), as can be seen in Figure 10b,c. A similar response is seen in Figure 10d, where the variation of the apparent power is plotted as a function of core leg height and the air gap between HV and LV windings.

The second equality constraint uses the leakage reactance of the power transformer, and its 3D response in the associated RSM heuristic model is shown in Figure 11. It is seen that the reactance is inversely proportional to the core leg height and, consequently, to the height of the transformer windings (h₁, h₂); the core diameter (d₀) directly affects the mean diameter of the LV and HV transformer windings; the model response is shown in Figure 11a. The model’s response as a function of the variations on the factors h and b₁ is shown in Figure 11b, where an increase is seen due to increments in the radial depth of the low voltage winding (b₁). Figure 11c shows a similar surface response to Figure 11b when there are changes in the independent variables h and b₂. Finally, Figure 11d shows the RSM polynomial model of reactance under variations of the factors h and g. In all figures, there is a close agreement with respect to the 3D FE transformer model results (indicated by the dots).

Once the RSM polynomial models have been obtained for the three objective functions (illustrated in Figure 7, Figure 8 and Figure 9) and the RSM models for the apparent power and leakage reactance, the multi-objective optimization defined by Equations (15)–(21) is solved using the NSGA-III. Figure 12 shows the Pareto front (blue points) with the normalized cost functions. The criteria employed to select the best solution consist of choosing the operating point with the lowest power losses (f₃(x)) and those weights of copper (f₂(x)) and core (f₂(x)) that represent the minimum cost. The optimal solution selected is indicated in Table 6. Figure 13 shows the running metric of the IGD that shows the convergence of the multi-objective problem.

The running IGD metric shown in Figure 13 oscillates, but in the last few generations, the performance suddenly gets better. It could be helpful to compare the true IGD values calculated with the true Pareto optimal points. However, this is one of the real-world problems for which a Pareto optimal solution is not available. Applying the Running Performance metric [33] in Figure 13, it was possible to appreciate how the algorithm improves its performance from the start to the end of 600 generations.

The optimal solution to the optimization problem is given in Table 6. The values of the design variables represent the optimum geometry of the core and low- and high-voltage windings. It can also be seen that the apparent power and reactance, related to the equality constraints, are fulfilled with a minimum difference. Table 6 also gives the 3D FE transformer model simulation results for the optimum solution, where there are close values with respect to the RSM-based optimization with NSGA-III; this demonstrates the validity of the RSM models and the NSGA-III optimization. An accuracy in the apparent power of 0.01% and an accuracy value of 0.01% for the reactance were obtained. There is an overprediction of 7.03 W in the electrical power losses. The copper and core weights are almost identical to those obtained by the RSM-based NSGA-III optimization. The tank has a height (th) of 1243 mm, a width (tw) of 1345 mm, and a depth (td) of 536 mm. The computing time was 2575 s using a desktop computer with a processor i7-7700, 3.6 GHz, and 16 GB of RAM.

Finally, Figure 14 shows the core flux density distribution (Figure 14a) and the core loss distribution (Figure 14b). The 3D FE model solution corresponds to the optimum operating point reported in Table 6. The tangential field is used as a boundary condition to reduce the model by half its mesh size; in all FEA simulations, the FE mesh was constructed using the automatic mesh generator.

6. Conclusions

This paper has presented a multi-objective optimization to find the optimum dimensions of the core and windings of a three-phase power transformer. The cost function has three objectives (electrical power loss, copper, and core weights) and two equality constraints (apparent power and leakage reactance). Firstly, an analytical model was solved using the NSGA-III to find an optimal region; secondly, RSM was applied using 3D FEA to develop polynomial models; and finally, the RSM-based cost function was solved using the NSGA-III. It is concluded that the use of RSM polynomial models allows for faster optimization while maintaining accuracy because FEA was employed to construct the empirical models. The running metric was used to assess the convergence of the optimization process. It was numerically demonstrated that NSGA-III has better convergence performance and can handle more accurately and faster a design problem with three objectives and two equality constraints than the commonly used optimizer, NSGA-II.