FEM-based transformer design optimization technique with evolutionary algorithms and geometric programming

Abstract

Cost optimization of power transformers is a complex engineering task. The design engineers have to solve this multidisciplinary, non-linear optimization problem in a very short-time. Due to the importance of the topic, numerous optimization methods have been developed in the industry to solve this problem. These methods generally use different, simplified two-winding transformer models, where the windings are generally modeled by their copper filling factors. These models cannot consider the eddy current losses and the temperature gradients of the windings properly, because the calculation of these quantities needs knowledge of the conductor dimensions. This paper proposes a novel method, which uses a general geometric programming search based sub-problem to determine the optimal conductor dimensions for the optimal winding shape. The proposed method considers the temperature gradients of the windings and uses FEM to determine the eddy losses of the windings. The whole optimization process is made by an evolutionary algorithm (NSGA-II)-based search. The paper presents this novel transformer optimization methodology and then illustrates it with a practical example.

Keywords

Optimization evolutionary computation FEA power transformers

1. Introduction

Power transformer design and optimization is still a challenging task, where a lot of engineering and economic aspects have to be considered simultaneously. The importance of the topic is supported by the fact that in the past 35 years more than 500 studies were published on this topic [1,2]. This paper deals with calculation of the cost-optimal key-design parameters in the first preliminary or tendering design stage (Fig. 1a). The cost optimization of power transformers is object of research since about 1900 and belongs to the most general branches of the non-linear mixed integer mathematical optimization problems [1]. Therefore, various reduced model-based transformer optimization methods have been proposed in the industry, to replace the judgement and experience of the engineers and increase the competitiveness of the factory [1]. There are generally four different quantities extremized during the optimization: the mass of the transformer, the efficiency, the manufacturing cost of the transformer or the total cost of the ownership. This paper deals with the minimization of the total cost of the ownership (TC) [3].

Fig. 1.

(a) Schematic view of transformer design process. (b) Key parameters of 3-phase core-type transformer. (c) Radial and (d) axial components of eddy loss density in transformer windings.

The majority of the published cost optimization methods consists in modelling the transformer only by the key-design parameters of its core and windings (Fig. 1b) [1]. These methods use a copper filling factor-based winding models to estimate the load loss of the transformer during the optimum search. This analytical formulation approximates the total amount of the load losses well, however, an accurate calculation of the eddy losses and temperature gradients is impossible from knowledge of the copper filling factor and the outline of the winding. Moreover, modern transformer windings are made of more than one winding material. It can be seen in Fig. 1c and Fig. 1d that the radial and axial distributions of the eddy current losses are non-homogeneous in the windings. The most critical part is the top of the HV winding, because it is heated by the radial component of the magnetic flux density and circulating oil is hot near to the top of the tank. To enhance the thermal and mechanical durability of the windings, these winding ends are usually made of different winding materials, and contain more parallel branches to minimize the perpendicular size of the conductors to the radial flux density [4].

This paper proposes a novel method, which uses a combination of two optimization methods: evolution algorithm (NSGA II) and generalized geometric programming (GP). The GP-based part of the optimization algorithm can determine the optimal conductor sizes automatically, from the shape, copper filling factor and other key-design parameters of the winding. Knowledge of the optimal conductor size allows applying a FEM calculation to calculate the eddy current losses, temperature gradients and short-circuit impedance more precisely during the optimization process. Moreover, the model can handle axially segmented windings.

2. Transformer optimization method

The transformer optimization methodology has been realized in the Ārtap framework [5], which is an open-source tool for robust design optimization. In the realized methodology, every possible transformer design is represented as an individual. These individuals are defined by their seven independent parameters: current density of the windings (j_s, j_p, j_r), magnetic flux density (B), core diameter (D_c), height of the low voltage winding (h_s), and main gap distance (g) (Fig. 1b). Every other parameters of the applied transformer model can be determined by knowledge of the independent values and specification. The independent parameters are generated and optimized by the applied evolutionary algorithm (NSGA-II) as shown in Fig. 1a. The calculation of the dependent parameters is made in every iteration step, by the redefined evaluating function of the optimization framework. This calculation consists of three main calculation steps (Fig. 2a): the solution of the analytical model, FEM calculation and embedded geometric programming-based model. The analytical calculations consist of the determination of the geometrical parameters and calculation of the core losses [6]. The role of the FEM solution and geometric programming-based winding model is to determine the load losses, short-circuit impedance and temperature gradient for the new winding model. The paper extended with the temperature gradient calculation of the previously published FEM-based model [7]. The following section shows the objective function of the optimization and details of the new winding model.

2.1. Objective function – Total cost of ownership

The objective function is the total cost of ownership. This function contains the manufacturing cost of the active part and the cost of the calculated losses [6]: $\begin{eqnarray}\displaystyle TC=K_{1}\cdot P_{\text{nll}}+K_{2}\cdot P_{\text{ll}}+C_{0}\cdot M_{\text{C}}+\mathop{\sum }_{k=1}^{n}C_{k}\cdot M_{k}, & & \displaystyle\end{eqnarray}$ (1) where TC is the total cost of ownership of the active part in € and also the objective function of this optimization method. K₁ is the capitalized cost of the no-load loss and K₂ is the load loss capitalization cost in €/kW. P_nll and P_ll are the no-load and load losses of the transformer in kW. M_k is the mass of the k-th part of the model (k represents the core, low voltage (LV), high voltage (HV) and regulating (Reg) windings) in kg and C_k represents the specific cost of the transformer part in €/kg.

Fig. 2.

(a) Schematic view of the realized transformer design optimization algorithm. (b) Temperature notations around the winding-oil contact at approximately the half of winding height [8].

3. The extended winding model

Transformer models in the pre-conceptional design phase model the geometry of the windings by their outer boundary and their copper filling factors [1]. So, the solution of these simplified pre-conceptional transformer models provides the optimal current density, load loss and optimal shape of the winding (see Fig. 1). This data is enough to calculate the magnetic field in the working window by FEM. The magnetic field in the working window of a transformer that is generally nonlinear may be described by the phasor A of magnetic vector potential $\vec{A}$ and the corresponding equation reads $\begin{eqnarray}\displaystyle \text{curl}(\text{curl}\text{}\underline{\boldsymbol{A}})+\text{j}\cdot {\mu}{\omega}{\gamma}\text{}\underline{\boldsymbol{A}}={\mu}\text{}\underline{\boldsymbol{J}}_{\text{ext}}, & & \displaystyle\end{eqnarray}$ (2) where μ denotes the magnetic permeability and 𝛾 is the electric conductivity. Symbol $\vec{J}_{\text{ext}}$ denotes the density of field currents in the windings and ω = 2πf denotes the angular frequency. The boundary condition along a sufficiently distant boundary is of the Dirichlet type. By the solution of this problem, the value of the short-circuit reactance can be calculated from the total magnetic energy, evaluated at the peak current [4,6]: $\begin{eqnarray}\displaystyle x_{\text{pu}}={\displaystyle \frac{{\pi}^{2}f}{VA}}\int _{V}JA\text{d}V & & \displaystyle\end{eqnarray}$ (3) where VA = ∑ VA₊ represents the sum of the amperturns in the defined winding segments. The ratings are the same for both windings and the plus sign indicates the positive sign current.

However, to calculate the eddy losses and the temperature gradients in the winding segments accurately needs to know the exact conductor dimensions and winding layout. If the optimal width (d^∗) and height (h^∗) of the conductor are known, the axial and radial parts of eddy current losses can be expressed by the following formulas [4]: $\begin{eqnarray}\displaystyle P_{\text{ax}}={\displaystyle \frac{1}{24{\rho}}}({\omega}d^{\ast }B_{z})^{2},\quad P_{\text{rad}}={\displaystyle \frac{1}{24{\rho}}}({\omega}h^{\ast }B_{r})^{2}, & & \displaystyle\end{eqnarray}$ (4) where B_r and B_z represent the radial and the axial components of the magnetic flux density. This calculation of eddy current losses in the winding segments assumes that the eddy currents do not modify the magnetic field around the winding segments. From knowledge of the load losses of different winding segments and dimensions of the optimal conductors, the optimal values of temperature gradients can be calculated.

To calculate an optimal shape of the conductors, a geometric programming (GP) sub-problem is defined inside the evolutionary optimization sub-problem. GP is a type of the non-linear mathematical optimization problem characterized by the objective and constraint functions given in the following special form [9]: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \qquad \min f_{0}(x)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \text{s.t. }∼∼f_{i}(x)\leq 1,\quad i=1,\ldots ,m\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \qquad g_{j}(x)=1,\quad j=1,\ldots ,n\end{eqnarray}$ (5) where x = (x₁, x₂, x₃) is a vector containing the variables to be optimized, f₀, …, f_m are the posynomial functions, and g₀, …, g_n are the monomial functions. All elements of x must be positive. The monomial function g (x) is a power product that can be expressed in the following form: $\begin{eqnarray}\displaystyle g_{j}(x)=c_{g}\cdot x_{1}^{{\alpha}_{1}}\cdot x_{1}^{{\alpha}_{2}}\cdot \ldots \cdot x_{n}^{{\alpha}_{n}}, & & \displaystyle\end{eqnarray}$ (6) where c_g is the coefficient of the monomial and $c_{\text{g}}\in \mathbb{R}^{+}$ . Symbol 𝛼_i is the exponent of the variable where ${\alpha}_{i}\in \mathbb{R}$ . As an example, $g(x)=3\cdot x_{1}^{2}\cdot x_{2}^{0.24}\cdot x_{3}^{-1.12}$ is a monomial function of the variables x = (x₁, x₂, x₃). The posynomial function is the sum of monomials: $\begin{eqnarray}\displaystyle f_{i}(x)=\displaystyle \mathop{\sum }_{k=1}^{l}g_{k}(x)=\displaystyle \mathop{\sum }_{k=1}^{K}c_{k}\cdot x_{1}^{{\alpha}_{1k}}\cdot x_{1}^{{\alpha}_{2k}}\cdot \ldots \cdot x_{b}^{{\alpha}_{nk}}, & & \displaystyle\end{eqnarray}$ (7) where c_k >0. The application of this formalism guarantees that the GP solver finds the global optimum of the problem. Solution of this type of problem is really fast, using a cutting-edge solver (CVXOPT).

The objective function of this embedded geometric program to minimize the loss of the winding: $\begin{eqnarray}\displaystyle P_{\text{loss}}=P_{\text{dc}}+P_{\text{ax}}+P_{\text{rad}},\quad P_{\text{dc}}={\rho}\cdot {\displaystyle \frac{2\cdot {\pi}\cdot r_{\text{m}}}{A_{\text{Cu}}}}\cdot I^{2}, & & \displaystyle\end{eqnarray}$ (8) where the load loss of the winding consists of the dc loss (p_dc) and the axial and radial components of the eddy losses (p_ax, p_rad), 𝜌 represents the specific conductivity of the conductor, r_m represents the mean radius, and I is the phase current of the winding. This method uses the B_rad and B_ax values of the FEM calculation and the posynomial form of P_ax + P_rad + P_dc is used to calculate the loss in each segment. The result of this optimization are the optimal conductor sizes (h^∗, w^∗). The following constraints describe the conductor geometry in the segment: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle n=n_{\text{ax}}\cdot n_{\text{rad}},\quad n_{\text{ax}}\cdot h+n_{\text{ax}}\cdot t_{\text{ax}}\leq h_{\text{w}},\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle t_{\text{hor}}\cdot n_{\text{c}}\cdot n_{\text{rad}}+w\cdot n_{\text{c}}\cdot n_{\text{rad}}\leq t_{w},\quad {\lambda}_{\text{ff}}\leq {\displaystyle \frac{n\cdot A_{\text{Cu}}}{h_{\text{w}}\cdot t_{\text{w}}}},\end{eqnarray}$ (10) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle A_{\text{Cu}}=n_{\text{rad}}\cdot n_{\text{c}}\cdot w^{\ast }\cdot h^{\ast },\quad V_{\text{Cu}}=2\cdot {\pi}\cdot r_{m}\cdot n\cdot A_{\text{Cu}},\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle w^{\ast }\leq w_{\max },\quad h^{\ast }\leq h_{\max },\quad n_{c}\leq 1,n_{\text{rad}}\leq 1,\end{eqnarray}$ (12) where n_ax and n_rad represents the radial and the axial number of the conductors, n_c is the total number of the conductors in a winding, t_hor and t_ver are the horizontal and the vertical length of the insulation. A_Cu and V_Cu represents the copper are and the copper volume in a winding. The height of the winding is denoted by h_w and t_w, while w^∗ and h^∗ represent the conductor dimensions in the winding. Finally, w_max and h_max are the upper bounds for the conductor height and the conductor width.

The applied temperature gradient calculation uses an electrical analogy of heat exchange [8] to solve the 2D-Poisson equation of the steady-state heat flow, where the heat flow possible only in the axial or the radial direction (Fig. 2b) [8].

The heat-flow density value can be defined using the following formula, in case of disc-type windings, for the temperature rise calculation [8]: $\begin{eqnarray}\displaystyle P_{o}=\frac{P_{\text{loss}}}{n_{\text{ax}}\cdot 2\cdot {\pi}\cdot r_{m}}, & & \displaystyle\end{eqnarray}$ (13) where P_o is the heat flow density, P_loss is the total loss of the segment and r_m is the mean radius of the segment, $\begin{eqnarray}\displaystyle p_{\text{hor}}=0.5\cdot {\displaystyle \frac{p_{0}\cdot {\beta}_{\text{hor}}}{{\alpha}\cdot w_{\text{t}}}},\quad p_{\text{ver}}=0.5\cdot {\displaystyle \frac{p_{0}\cdot {\beta}_{\text{ver}}}{{\alpha}\cdot h}} & & \displaystyle\end{eqnarray}$ (14) where p_hor and p_ver are the horizontal heat flows. The temperature gradient are calculated from these values and the larger one is selected as the winding gradient $\begin{eqnarray}\displaystyle {\theta}_{p,\text{hor}}=\left({\displaystyle \frac{t_{\text{c}}}{2\cdot k_{c}}}\right)\cdot p_{\text{hor}}\quad {\theta}_{p,\text{ver}}=\bigg(\frac{{\alpha}\cdot n_{\text{rad}}\cdot m_{C}\cdot t_{c}}{2\cdot k_{\text{c}}}\bigg)\cdot {\displaystyle \frac{p_{\text{ver}}}{3}}, & & \displaystyle\end{eqnarray}$ (15) where k_c represents the thermal heat conductivity of the insulation material and t_c is the equivalent length of the insulation.

Temperature gradient through the oil boundary layer, in the case of ON or OF cooling: $\begin{eqnarray}\displaystyle {\theta}_{s}=C\cdot a\cdot \sqrt{p_{\text{hor}}},\quad {\theta}_{s}=C\cdot \sqrt{p_{\text{ver}}} & & \displaystyle\end{eqnarray}$ (16) where the auxiliary parameter a = 1 in the case of zig-zag cooling. Regardless of the used formula, the incorporated constant C depends on the established cooling mechanism. Potential interpretations are shown in [8].

The gradient is defined by the following generalized posynomial: $\begin{eqnarray}\displaystyle \min ({\theta}_{p,\text{hor}}+{\theta}_{s,\text{hor}},{\theta}_{p,\text{ver}}+{\theta}_{s,\text{ver}})\leq {\theta}. & & \displaystyle\end{eqnarray}$ (17) The min function is not a valid generalized posynomial inequality, but it can be handled by a replacement [9].

4. Validation of the transformer model

The physical correctness of the transformer model has been validated by an existing 10 MVA, 33/6.9 kV, star/star connected transformer ([4], page 173–75). The detailed list of the transformer parameters presented in Table 1. The main dimensions of the working window (FEM model) and the resulting flux distribution are depicted in Fig. 3. The resulting value of the short-circuit impedance (SCI) from the finite element analysis is 7.52%, which is very close to the measured value (7.34%). The LV winding is made from 12 parallel conductors, whose thickness is 2.3 mm and the insulating paper thickness is 0.5 mm. The eddy current loss of this winding is calculated by the Rabins and Cogbill method and its value is 489.55 W [4]. The resulting winding dimensions after the GGP optimization differ from the manufactured transformer dimensions, the optimized conductor width is 1.7 mm. The dc loss of the offered winding system is 21.6 kW instead of the original 21.8 kW, and the eddy current loss after the calculation is 583 W. The GGP algorithm calculated 14 K for the temperature gradient.

Fig. 3.

Table 1

Parameters of the optimized transformer

Parameter		Dimension	Value
Nominal power		MVA	10
Frequency		Hz	50
Connection group			YNyn1
Number of phases		#	3
Short circuit impedance		%	7.5
Main gap		mm	50
Sum of the end insulation		mm	80
Phase distance		mm	50
Core-Inner winding distance		mm	20
	Number of legs	#	3
	Flux density limit in columns	T	1.7
Core	Filling factor	%	88
	Material type		M1H
	Material price	€/kg	3.5
	Line voltage	kV	6.9
	Phase voltage	kV	4.0
Low voltage	BIL	kV	125
Winding	AC	kV	50
	Copper filling factor	%	60
	Material and manufacturing price	€/kg	10
	Line voltage	kV	33
	Phase voltage	kV	19.05
High voltage	BIL	kV	125
Winding	AC	kV	50
	Copper filling factor	%	60
	Material and manufacturing price	€/kg	8.5

Table 2

Resulted losses and key-design parameters

Quantity	Dimension	Manufactured	Optimization
Load loss	kW	49.6	51.8
Flux density	T	1.7	1.5
Core mass	t	7.3	7.3
Core loss	kW	7.7	5.4
Turn voltage	V	46.8	39.1
Current density
Inner winding	A/mm²	2.6	2.3
Outer winding	A/mm²	2.5	2.17
Objective function	€	404518	341869

4.1. Illustrative example for optimization

The transformer optimization was made for a solar power plant transformer design to illustrate the optimization method. The economic scenario and the capitalization factors selected from [10], where the capitalized cost of the load losses is k₁ =1000 €/kW and the capitalized loss of the no-load losses is k₂ = 7000 €/kW. The same technical data sheet, which was presented in Table 1, was used for this optimization. The NSGA-II algorithm generated 100 individuals for every 100 generations. The resultant geometry is depicted in Fig. 2 and the resulting main electrical parameters are listed in Table 2. The optimal transformer’s turn voltage is significantly smaller, than in the manufactured case. The smaller turn voltage means that the optimal design contains more steel and smaller amount of copper. This corresponds with the practice, that in the case of high ratio of k₁∕k₂ the transformer design with high load losses and small core losses leads to the economic optimum.

5. Conclusions

This paper proposed a novel transformer design optimization algorithm, which can determine the optimal conductor shapes for the transformer windings in the first design process. By the knowledge of the optimal conductor sizes, the thermal properties of a transformer can be considered at the beginning of the design for different winding layouts. The results of the transformer model were validated on a data of a manufactured 10 MVA transformer. Further research in the domain can use this method to examine the impact of the new biodegradable transformer oils on the key design parameters of large power transformers.

References

Orosz

, Evolution and modern approaches of the power transformer cost optimization methods, Periodica Polytechnica Electrical Engineering and Computer Science (2019).

Mendes

H.F.

Advanced modeling and experimental validation of an optimized power transformer tank, in: ASME 2018, American Society of Mechanical Engineers, 2018, pp. V013T05A043–V013T05A043.

Commission

I.E.

, IEC Temperature rise, in: IEC Std 60076-2, 2006.

Kulkarni

S.V.

and Khaparde

, Transformer Engineering: Design and Practice, vol. 25, CRC Press, 2004.

Pánek

Artap: Robust design optimization framework for engineering applications, in: Proceedings of the ICDS2019, 2019, p. 4, to appear.

Georgilakis

P.S.

, Spotlight on Modern Transformer Design, Springer Science & Business Media, 2009.

Orosz

Pánek

and Karban

, Fem based preliminary design optimization in case of large power transformers, 2019.

Orosz

Nagy

and Tamus

Z.Á.

, A generalized geometric programming sub-problem of transformer design optimization, in: Technological Innovation for Smart Systems, Springer, 2017, pp. 373–381.

Boyd

and Vandenberghe

, Convex Optimization, Cambridge University Press, 2004.

10.

Orosz

, Analysis of the green power transition on optimal power transformer designs, Periodica Polytechnica Electrical Engineering and Computer Science59(3) (2015), 125–131.