Power Transformer’s Electrostatic Ring Optimization Based on ANSYS Parametric Design Language and Response Surface Methodology

: In this paper, in view of the low efficiency of the traditional finite element method (FEM), which has been widely used in the insulation design of power transformers, the response surface methodology (RSM) is proposed to optimize the insulation structure of a power transformer electrostatic ring. Firstly, the power transformer model was built using the ANSYS parametric design language (APDL) to realize the automatic pre-processing of numerical calculation. Then with the objective of reducing the maximum electric ﬁeld intensity, the Taguchi method was used to select the parameters that have a greater impact on the maximum electric ﬁeld intensity, by which the subsequent optimization process could be e ﬀ ectively simpliﬁed. The test points were constructed by the central composite design (CCD) and a response surface model was established by the mutual calls of MATLAB and ANSYS. Finally, the variance analysis, diagnostic analysis, and signiﬁcance test of regression were carried out to obtain the ﬁnal response surface model. By comparing the result of RSM with that of FEM, we can ﬁnd that the results obtained by the two methods are consistent and the maximum electric ﬁeld strength is obviously reduced. The RSM is more systematic and convincing, which improves the optimization e ﬃ ciency and provides a reliable and fast way for the optimization of power transformers.


Introduction
As the main device of power systems, the power transformer is essential to the safe and reliable operation of the whole electrical network [1].The insulation structure determines the electric field distribution; thus, the optimization of the main insulation structure is of great importance to ensure the safe and stable operation of the power transformer [2,3].With the increase of voltage levels, more and more attention has been paid to the electrical insulation problem of the power transformer. However, the insulation structure of a transformer is particularly complex, which makes the electric field distribution non-uniform. Therefore, to confirm the rationality and reliability of the main insulation structure of the transformer, it is necessary to test, calculate, and analyze the electric field distribution [4].
In engineering applications, electric field analysis methods can be divided into two types, the analytical methods and the numerical methods. A few simple problems may be solved by the analytic method; but in the majority of cases, an analytical method simply does not exist or it is extremely difficult to obtain.As a numerical method, the FEM (finite element method) has been applied

Calculation Model
The calculation model of a 500 kV power transformer is shown in Figure 1 and the values of the initial parameters are listed in Table 1.

Calculation Model
The calculation model of a 500 kV power transformer is shown in Figure 1 and the values of the initial parameters are listed in Table 1. The first boundary condition: 0kV The first boundary condition: 200kV The first boundary condition: 500kV The first boundary condition: 0kV The second homogeneous boundary condition   (1) where ε is the dielectric constant; ϕ is the electric scalar potential of the electrostatic field; Ω is the whole computational domain; Г1 is the high voltage winding boundary; Г2 is the low voltage winding boundary;Г3 is the regulating winding, the upper ,and the right boundary; and Г4,5 are the left and lower boundary, respectively.  The boundary value problem of the electrostatic field can be expressed by Equation (1).
where ε is the dielectric constant; φ is the electric scalar potential of the electrostatic field; Ω is the whole computational domain; Г 1 is the high voltage winding boundary; Г 2 is the low voltage winding boundary; Г 3 is the regulating winding, the upper, and the right boundary; and Г 4,5 are the left and lower boundary, respectively. The APDL is applied to realize automatic pre-processing as follows: (1) The relative dielectric constants of oil and pressboard were 2.2 and 4.5, respectively; (2) the voltage applied on the high voltage winding was 500 kV (Г 1 ), the low voltage winding was 200kV (Г 2 ), and the voltage regulating winding and the upper and the right boundary were 0kV (Г 3 ). These are the first kind of boundary conditions and (3) the left and lower boundary are the second kind of homogeneous boundary conditions (Г 4,5 ). In addition, the boundary conditions of Figure 1b were interpolated based on the results of Figure 1a.
The computer processor we used had an Intel(R) Core(TM) i5-7500 with 8.00 GB of memory. After modeling, meshing, applying boundary conditions, and solving by APDL, the electric field intensity nephogram of the initial structure is shown in Figure 2. By analyzing the initial insulation structure of the power transformer, it can be found that the maximum electric field intensity appeared at the upper right corner of the electrostatic ring. Therefore, we try to reduce the maximum electric field intensity by changing the structure parameters related to the electrostatic ring. The APDL is applied to realize automatic pre-processing as follows: (1) The relative dielectric constants of oil and pressboard were 2.2 and 4.5, respectively;(2) the voltage applied on the high voltage winding was 500 kV (Г1), the low voltage winding was 200kV (Г2),and the voltage regulating winding and the upper and the right boundary were 0kV (Г3).These are the first kind of boundary conditions and (3) the left and lower boundary are the second kind of homogeneous boundary conditions (Г4,5).In addition, the boundary conditions of Figure 1bwere interpolated based on the results of Figure 1a.
The computer processor we used had an Intel(R) Core(TM) i5-7500 with 8.00 GB of memory. After modeling, meshing, applying boundary conditions, and solving by APDL, the electric field intensity nephogram of the initial structure is shown in Figure 2. By analyzing the initial insulation structure of the power transformer, it can be found that the maximum electric field intensity appeared at the upper right corner of the electrostatic ring. Therefore, we try to reduce the maximum electric field intensity by changing the structure parameters related to the electrostatic ring.

Taguchi Method
In the optimization experiments, the selection of the parameter that has no influence on the optimization result as a variable will not only increase the workload, but will also affect the accuracy of the optimization result. Therefore, the Taguchi method was employed to eliminate the irrelevant variables and select the parameters which have great influence on the maximum electric field intensity [40,41]. Then the response surface experiments were carried out with these filtered parameters, by which the optimization efficiency can be significantly improved.
The Taguchi method designs the horizontal combination of each parameter by establishing an orthogonal table, carries on an orthogonal experiment to calculate the proportion of the influence of each parameter on the optimization target, and filters the parameters according to the proportion. Each parameter that needs to be optimized usually takes three values in the range of change, which is called the impact factor.

RSM
The response surface model is also called the response surface function and the regression equation [42,43]. A proper function can make the approximation more precise and the design space more suitable for wider use. The quadratic polynomials can describe the real function and ensure that the model is as simple as possible with fewer undetermined coefficients. Its expression can be written as

Taguchi Method
In the optimization experiments, the selection of the parameter that has no influence on the optimization result as a variable will not only increase the workload, but will also affect the accuracy of the optimization result. Therefore, the Taguchi method was employed to eliminate the irrelevant variables and select the parameters which have great influence on the maximum electric field intensity [40,41]. Then the response surface experiments were carried out with these filtered parameters, by which the optimization efficiency can be significantly improved.
The Taguchi method designs the horizontal combination of each parameter by establishing an orthogonal table, carries on an orthogonal experiment to calculate the proportion of the influence of each parameter on the optimization target, and filters the parameters according to the proportion. Each parameter that needs to be optimized usually takes three values in the range of change, which is called the impact factor.

RSM
The response surface model is also called the response surface function and the regression equation [42,43]. A proper function can make the approximation more precise and the design space more suitable for wider use. The quadratic polynomials can describe the real function and ensure that the model is as simple as possible with fewer undetermined coefficients. Its expression can be written as Appl. Sci. 2019, 9, 4286 5 of 16 where ỹ is the objective function; β 0 , β j , β jj , and β ij are the polynomial coefficients; x j is the jth experimental variable; and n is the number of experimental variables. Assuming that k times tests are needed to obtain the response surface model, Equation (2) can be expressed in the following matrix form: · β 0 β 1 · · · β n β 11 · · · β nn β 12 · · · β nn−1 where µ is the error of y. The polynomial coefficients are obtained according to Equation (5) by using the least squares principle to minimize the sum of the squares of the errors.

Central Composite Design (CCD)
The fitting accuracy of the RSM depends largely on the distribution of sample points in the design space. The CCD method follows specific rules when selecting the sample points [44][45][46]. According to the traditional method, when three variables and five levels are selected for experiments, at least 125(5 3 ) experiments are needed. The CCD can minimize the number of experiments on the premise of good accuracy. Figure 3 shows the distribution of sample points, taking three variables as an example [29]. Where α = 2 n/4 ; (0,0,0) is the center point; (±1,±1,±1) is the cubic point; and (±α,0,0), (0,±α,0), and (0,0,±α) are the axial points. where ỹ is the objective function; β0,βj,βjj,and βij are the polynomial coefficients; xj is the jth experimental variable; and n is the number of experimental variables. Assuming that k times tests are needed to obtain the response surface model, Equation (2) can be expressed in the following matrix form: 2  2  1  11  1  11  1  11 12  1 1 1  2  2  T  2  21  2  21  2  21 22  2 2 1  0  1  11  12  1   2  2  1  1  1 2  1   T  1  2   1  1 1 n n n n n n n n n nn nn where μ is the error of y. The polynomial coefficients are obtained according to Equation (5) by using the least squares principle to minimize the sum of the squares of the errors.

Central Composite Design (CCD)
The fitting accuracy of the RSM depends largely on the distribution of sample points in the design space. The CCD method follows specific rules when selecting the sample points [44][45][46]. According to the traditional method, when three variables and five levels are selected for experiments, at least 125(5 3 ) experiments are needed. The CCD can minimize the number of experiments on the premise of good accuracy. Figure 3 shows the distribution of sample points, taking three variables as an example [29]. Where α=2 n/4 ; (0,0,0) is the center point; (±1,±1,±1) is the cubic point; and (±α,0,0), (0,±α,0), and (0,0,±α) are the axial points.

Test for Significance of Regression
A test which is usually used to verify the overall applicability of the model and the significance of regression is called null hypothesis [21,47], as follows: H0: β1=β2=…=βn=0, H1: βj0 for at least one j. The value H0means that the variables have the same effect on the response value, and the error is only caused by the selection of the sample points. The rejection ofH0 indicates that at least one of the variables has a significant contribution to the model, that is, the model is reasonable. In addition,

Test for Significance of Regression
A test which is usually used to verify the overall applicability of the model and the significance of regression is called null hypothesis [21,47], as follows: H 1 : β j 0 for at least one j.
The value H 0 means that the variables have the same effect on the response value, and the error is only caused by the selection of the sample points. The rejection of H 0 indicates that at least one of the variables has a significant contribution to the model, that is, the model is reasonable. In addition, the p-value approach, which reflects the probability of the H 0 , can be used to test the null hypothesis. If the p-value is less than 0.05, the null hypothesis is rejected, that is, the H 0 is invalid.
The value R 2 is the determination coefficient and R 2 adj is the adjusted R 2 . The closer their values are to 1, the better is the fitting effect of the model.

Experimental Process and Analysis
As can be seen in Section 2.1, the maximum electric field intensity emerged at the upper right corner of the electrostatic ring. Therefore, we can select the parameters related to the electrostatic ring as the variables to optimize. The structure and parameters of the electrostatic ring are shown in Figure 4.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 16 the p-value approach, which reflects the probability of the H0, can be used to test the null hypothesis. If the p-value is less than 0.05, the null hypothesis is rejected, that is, the H0 is invalid. The value R 2 is the determination coefficient and R 2 adj is the adjusted R 2 .The closer their values are to 1, the better is the fitting effect of the model.

Experimental Process and Analysis
As can be seen in subsection2.1, the maximum electric field intensity emerged at the upper right corner of the electrostatic ring. Therefore, we can select the parameters related to the electrostatic ring as the variables to optimize. The structure and parameters of the electrostatic ring are shown in Figure 4. Due to the existence of the pressboard, for easy analysis, H is simplified as h, the distance from the top surface of the electrostatic ring (without insulation layer) to the first pressboard above the electrostatic ring. The relationship between the two arcs of the electrostatic ring can be expressed by Equation (6).
Among them, the coordinates of the center of the two arcs are(x1, y1) and (x2, y2). Values T1 and T2 are the ordinate values of the upper and lower surfaces of the electrostatic ring (excluding the insulating layer). As can be seen from the above equation, R2 and θ2 can be derived from R1 and θ1. Therefore, the four optimization variables R1, θ1, R2, and θ2 can be simplified into two variables, namelyR1 and θ1.
The optimization process is demonstrated in Figure 5. Due to the existence of the pressboard, for easy analysis, H is simplified as h, the distance from the top surface of the electrostatic ring (without insulation layer) to the first pressboard above the electrostatic ring. The relationship between the two arcs of the electrostatic ring can be expressed by Equation (6).
Among them, the coordinates of the center of the two arcs are (x 1 , y 1 ) and (x 2 , y 2 ). Values T 1 and T 2 are the ordinate values of the upper and lower surfaces of the electrostatic ring (excluding the insulating layer). As can be seen from the above equation, R 2 and θ 2 can be derived from R 1 and θ 1 . Therefore, the four optimization variables R 1 , θ 1 , R 2 , and θ 2 can be simplified into two variables, namely R 1 and θ 1 .
The optimization process is demonstrated in Figure 5.

Taguchi Experiment
Through the above analysis, the Taguchi experiment was carried out with R 1 , θ 1 , h, s, and w as variables to be optimized. Taking the median values of the variables on both sides and their adjacent values as the impact factors, the impact factors of the optimized parameters are listed in Table 2.  Figure 5.Flowchart of insulation structure optimization of power transformer based on RSM.

Taguchi Experiment
Through the above analysis, the Taguchi experiment was carried out with R1, θ1, h, s, and w as variables to be optimized. Taking the median values of the variables on both sides and their adjacent values as the impact factors, the impact factors of the optimized parameters are listed in Table 2. The orthogonal tableL9 (3 5 ) is obtained by the Taguchi method, where L represents an orthogonal matrix; 9 stands for the number of trials; 3 is the number of times each variable is taken, and 5 symbolizes the number of optimization parameters. According to the orthogonal table, the maximum electric field intensity can be achieved by the mutual call between MATLAB and ANSYS, as shown in the Table 3. Table 3.Orthogonal matrix and results of L9(3 5 ).  The orthogonal table L 9 (3 5 ) is obtained by the Taguchi method, where L represents an orthogonal matrix; 9 stands for the number of trials; 3 is the number of times each variable is taken, and 5 symbolizes the number of optimization parameters. According to the orthogonal table, the maximum electric field intensity can be achieved by the mutual call between MATLAB and ANSYS, as shown in the Table 3. Table 3. Orthogonal matrix and results of L 9 (3 5 ). From the results obtained by Table 3, the proportion of each parameter is calculated by Equation (7).

Number of Experiments
where x i represents R 1 , θ 1 , h, s, and w; EE is the proportion of each optimization parameter; m xi (E i ) is the average value of E max under the ith impact factor of x; and m(E) is the average value of E max . The proportion of each optimization parameter is listed in Table 4. It can be seen from Table 5 that R 1 , h, s, and w have greater influence on the maximum electric field intensity,while θ 1 can be eliminated from the impact variables. The Taguchi method effectively reduced the workload of the numerical analysis of the subsequent response surface. In the following experiment, for the convenience of analysis, the electrostatic ring was divided into two arcs; both of them were 90 degrees. The values R 1 , h, s, and w were selected as experimental variables. Note: The −α and α levels of R 1 are −1 and 27, respectively, which exceed the range of variation and are appropriately adjusted.

Experimental Data
After determining the experimental variables, the appropriate range of values was selected for the response surface experiment. There are four optimization variables in this experiment, so α=2. The design and combination of the experimental sample points were carried out by the CCD method, as shown in Table 5. Table 6 is the experimental data results of the sample points which were obtained by calling each other between ANSYS and MATLAB.  The last two columns of Table 6 are the number of elements and nodes (degrees of freedom) of the local finite element model in Figure 1b. According to statistics, it took only 61.8 seconds to complete the above 30 CCD experiments. Without the CCD method, 625 RSM experiments would be required. Thus, by using CCD method, the optimization efficiency was greatly improved.
Referring to the above data, the coefficients of the response surface model were obtained by Equation (5). We can get the initial response surface model from Table 7.

Analysis of Variance
The ANOVA table summarizes the results of the variance analysis, as shown in Table 8. As shown in Table 8, the p-value of the regression model was less than 0.0001, which means that the choice of the quadratic polynomial model is reasonable, and the relationship between the maximum electric field intensity value, y, and the regression equation is extremely significant. Thus, the null hypothesis, H 0 , is not valid and can be rejected. The p-value of all linear terms were less than 0.0001, indicating that each variable in R 1 , h, s, and w is an extremely significant factor affecting the optimization objective. Except that the interaction between R 1 and s was significant, the p-value of each of the other two was greater than 0.05, demonstrating the interaction between the two has little effect on the response. The p-values of the quadratic terms R 1 2 and s 2 were less than 0.0001, and the significance was higher than the other two, which is an extremely significant factor. In addition, R 2 adj is 0.9841, which is very close to 1. The model can explain 98.41% of the test changes and better reflect the law of data changes with less experimental error and better fitting. Furthermore, from the confidence interval listed in Table 8, the significant items are also R 1 , h, s, w, R 1 s, R 1 2 , and s 2 , for the upper and lower bounds of their confidence intervals do not contain zero. Note: *** Extremely significant; ** Very significant; * Significant

Diagnostic Analysis
To verify the applicability of the model, the diagnostic analysis was generally made by making the regression analysis and giving some diagnostic plots.
The design of the normal probability plot made the cumulative normal distribution a straight line. In Figure 6a, the points on the normal probability plot were approximately distributed in a straight line, which indicates the accuracy of the data for calculating the electric field intensity.  As can be seen from the residuals versus run plot in Figure 6b, the residuals of runs were distributed within the prescribed appropriate range, which implies that choosing the quadratic polynomial as the response surface function is particularly suitable. In Figure 6c and Table 9, there was a good correlation between the predicted value and the actual value, which indicates that the polynomial model has a good predictive ability. In a word, the response surface model adopted in this study is acceptable.

Run
Actual Values/(V/mm) Predicted Values/(V/mm) As can be seen from the residuals versus run plot in Figure 6b, the residuals of runs were distributed within the prescribed appropriate range, which implies that choosing the quadratic polynomial as the response surface function is particularly suitable. In Figure 6c and Table 9, there was a good correlation between the predicted value and the actual value, which indicates that the polynomial model has a good predictive ability. In a word, the response surface model adopted in this study is acceptable.

3-D Response Surface Analysis
For a better demonstration of the influence of the interaction between two variables, the 3-D plots were adopted. If the interaction has a significant impact on the response surface, the variation of 3-D plot will be accordingly large. Figure 7a shows that when R 1 and s interact, if we fix one of the variables and increase the other variable, the maximum electric field intensity will decrease. As can be seen from the trend of the surface graph, the maximum electric field intensity decreased faster with R 1 than with s, indicating that R 1 has a greater influence on the response value. Comparing Figure 7a with Figure 7b, the response surface of Figure 7b fluctuated slightly, which means that the change of the interaction term hw has no significant effect on the maximum electric field intensity. Similarly, hs and sw are not significant factors either.
As can be seen from Figure 7c,d, compared with s and w, the change of the maximum electric field intensity was larger with the change of R 1 . The maximum electric field intensity increases with the increase of w, which has negative effects on the model. Therefore, in the structural design process, a smaller starting position w should be selected as far as possible within a reasonable range. From the variation range of the two figures, the interactions between R 1 and h, R 1 , and w may be a significant factor affecting the maximum electric field intensity. However, since their p-values in Table 8 are little larger than 0.05, R 1 h and R 1 w can be considered insignificant factors. Figure 7a shows that when R1 and s interact, if we fix one of the variables and increase the other variable, the maximum electric field intensity will decrease. As can be seen from the trend of the surface graph, the maximum electric field intensity decreased faster with R1 than with s, indicating that R1 has a greater influence on the response value. Comparing Figure 7a with Figure  7b, the response surface of Figure 7b fluctuated slightly, which means that the change of the interaction term hw has no significant effect on the maximum electric field intensity. Similarly, hs and sw are not significant factors either. As can be seen from Figure 7cand Figure 7d, compared with s and w, the change of the maximum electric field intensity was larger with the change of R1. The maximum electric field intensity increases with the increase of w, which has negative effects on the model. Therefore, in the structural design process, a smaller starting position w should be selected as far as possible within a reasonable range. From the variation range of the two figures, the interactions between R1 and h, R1, and w may be a significant factor affecting the maximum electric field intensity. However, since their p-values in Table 8 are little larger than 0.05, R1h and R1wcan be considered insignificant factors.
Through the above analysis, the insignificant factors, such as R1h, R1w, hs, hw, sw, h 2 , and w 2 , are removed, and the final response surface equation with a higher prediction ability is shown in Equation (8). Through the above analysis, the insignificant factors, such as R 1 h, R 1 w, hs, hw, sw, h 2 , and w 2 , are removed, and the final response surface equation with a higher prediction ability is shown in Equation (8)

Results and Discussion
It is an important goal of the response surface optimization to get the optimal results with the appropriate values assigned to each variable. Equation (8) is the quadratic response surface model with single objective optimization, which can be solved quickly by quadratic programming. By comparing the maximum electric field intensities on the electrostatic ring before and after optimization in Figure 8, the decrease of the maximum electric field intensity can be observed more intuitively. The result shows that the maximum electric field intensity has been significantly reduced.

Results and Discussion
It is an important goal of the response surface optimization to get the optimal results with the appropriate values assigned to each variable. Equation (8) is the quadratic response surface model with single objective optimization, which can be solved quickly by quadratic programming. By comparing the maximum electric field intensities on the electrostatic ring before and after optimization in Figure 8, the decrease of the maximum electric field intensity can be observed more intuitively. The result shows that the maximum electric field intensity has been significantly reduced. In order to verify the validity of the RSM, we compared the result of RSM with that of traditional FEM. The FEM taking much time and requiring many experiments, to ensure the fairness of the experiment and save the calculation time, thirty experiments were carried out by using the FEM.As can be seen from Table 10, the results obtained by RSM and FEM were consistent, while the time the of RSM experiments was much less, which demonstrates the effectiveness of RSM and also shows the high optimization efficiency and good prediction ability of this method. In order to verify the validity of the RSM, we compared the result of RSM with that of traditional FEM. The FEM taking much time and requiring many experiments, to ensure the fairness of the experiment and save the calculation time, thirty experiments were carried out by using the FEM. As can be seen from Table 10, the results obtained by RSM and FEM were consistent, while the time the of RSM experiments was much less, which demonstrates the effectiveness of RSM and also shows the high optimization efficiency and good prediction ability of this method. In Figure 9, a comparison of the variations of variables predicted by the FEM and those actually obtained by the RSM was made. It can be observed that the variations of variables in this paper are monotonic, so FEM can obtain the optimal results by inference in this special case. However, this method is no longer applicable when the variables are not monotonic. The RSM has a wider range for applications and the result obtained by it is more rigorous and convincing.

Conclusion
In this paper, the response surface methodology, Taguchi method, and ANSYS parametric modeling language were combined to optimize the electrostatic ring of a 500 kV power transformer. The APDL extracted the maximum electric field intensity automatically and the mutual calls between MATLAB and ANSYS solved the problem of repeated manual modeling and meshing in traditional methods. To avoid the increased calculating workload of developing experiments aroused by irrelevant variables, the Taguchi method was utilized to filter the variables before making response surface experiments, by which the optimization process was effectively simplified. The experimental points were constructed by the experimental design method of CCD, and the response surface model was then obtained. In order to ensure the accuracy and predictability of the model, the response surface model was analyzed by ANOVA, diagnostic analysis, and the significance test of regression, and the final response surface equation was determined. Comparing the optimization results with the traditional finite element method, on the premise of ensuring accuracy, the proposed method in this paper has a higher optimization efficiency.
In general, the structure of the electrostatic ring is optimized more systematically and completely through the Taguchi method and response surface experiments. The feasibility and efficiency of this method were verified, which provides a new and more systematic, fast, and efficient way for the optimization of the transformer insulation structure.

Conclusions
In this paper, the response surface methodology, Taguchi method, and ANSYS parametric modeling language were combined to optimize the electrostatic ring of a 500 kV power transformer. The APDL extracted the maximum electric field intensity automatically and the mutual calls between MATLAB and ANSYS solved the problem of repeated manual modeling and meshing in traditional methods. To avoid the increased calculating workload of developing experiments aroused by irrelevant variables, the Taguchi method was utilized to filter the variables before making response surface experiments, by which the optimization process was effectively simplified. The experimental points were constructed by the experimental design method of CCD, and the response surface model was then obtained. In order to ensure the accuracy and predictability of the model, the response surface model was analyzed by ANOVA, diagnostic analysis, and the significance test of regression, and the final response surface equation was determined. Comparing the optimization results with the traditional finite element method, on the premise of ensuring accuracy, the proposed method in this paper has a higher optimization efficiency.
In general, the structure of the electrostatic ring is optimized more systematically and completely through the Taguchi method and response surface experiments. The feasibility and efficiency of this method were verified, which provides a new and more systematic, fast, and efficient way for the optimization of the transformer insulation structure.