Multi-Objective Optimal Control Approach for Static Voltage Stability of Power System Considering Interval Uncertainty of the Wind Farm Output

The static voltage stability margin (SVSM) of a power system considering the uncertain fluctuation range of wind farm (WF) output can be described as an interval value called the SVSM interval. A multi-objective optimal control model for SVS of a power system considering the interval uncertainty of WF output is proposed. The objective functions of the model are to increase the central value and reduce the fluctuation range of the SVSM interval, and the decision variables are the active power output and terminal voltage of generators and the switching capacity of shunt capacitors. Thus, it is a multi-layer optimization model. A parametric approximation (PA) method is used to obtain the approximate functional relationship between the optimal objective function values and the decision variables of the inner-layer and mid-layer optimization models and convert the optimization model into a single-layer bi-objective optimization model. A method for obtaining the continuous Pareto frontier of the bi-objective optimization model is proposed based on the normalized normal constraint and PA methods, and the compromise optimal solution calculated from the continuous Pareto frontier is used as the optimal control scheme. Two methods for improving the calculation efficiency of the PA method are also proposed. Finally, results from experimentation on the IEEE 39-bus system and an actual provincial power grid demonstrate the effectiveness of the proposed method.


P Gi0
active power output of generator i V Gi ref reference value of terminal voltage of generator i Q ci reactive power output of the shunt capacitors i P Li0, Q Li0 initial active and reactive load power of bus i G ij, B ij real and imaginary parts of the i-th line and the j-th column element of the node admittance matrix N g

Number of generators
The associate editor coordinating the review of this manuscript and approving it for publication was Diego Oliva .

I. INTRODUCTION
Affected by the uncertain fluctuation of wind speed, the active-power output of wind farms (WFs) has large uncertainty and volatility, which brings great challenges to the secure operation of power systems [1], [2]. The static voltage stability margin (SVSM) is a common voltage stability evaluation index in power system operation, representing the maximum load increase that the system can bear in the current operating state [3]. Uncertain fluctuation of the WF output will cause uncertain fluctuation of the SVSM of the power system. When an interval is used to describe the uncertain fluctuation of WF output, the SVSM is also correspondingly within an interval instead of being a fixed value [4]. When the fluctuation range of WF output is large, the fluctuation range of SVSM may be resulted large. If the values of the SVSM range are small, the lower bound of the SVSM range is very small, and the system is prone to voltage collapse with the increasing of load. Therefore, when considering the optimal SVS control of a power system with uncertain fluctuation of WF output, it is necessary to increase the central value and reduce the radius of the SVSM interval at the same time, which is a bi-objective optimization problem. Since the calculation of the upper and lower bounds of the SVSM interval is a bi-layer optimization model [5], the optimal SVS control model of a power system considering the uncertain fluctuation of WF output is a bi-objective multi-layer optimization model. How to effectively solve this model is a challenging problem. Many studies have been published on the analysis and control of the SVS of power systems with WFs. The authors of [2] proposed an SVS analysis method suitable for large-scale power systems with high wind-power penetration and verified that the use of doubly-fed induction generator wind turbines can improve the SVS. To evaluate the operating status of large power systems with WFs, the authors of [6] proposed a method based on the static voltage security region, which helps to evaluate the SVSM and security status online. In [7] the authors investigated the feasibility of utilizing the reactive power of grid-connected variable-speed wind generators to enhance the steady voltage stability margin of the system. In [8] data-driven methods were used to evaluate several approaches to fitting the PV curve of WFs and determine the best fitting method based on a variety of indicators, which provided useful information for SVS evaluation. However, the influence of the uncertain fluctuation of WF output on the analysis and control of SVS was not considered in any of these studies [2], [6]- [8]. Methods for considering the influence of the uncertain WF output fluctuation on the analysis and control of SVS include the Monte Carlo method [9], [10], probability method [11]- [13], and interval method [4], [5], [14]. In [9] the authors combined the Monte Carlo method and a neural network-based algorithm to propose a new method to evaluate SVS. Authors of [10] included unstable states caused by insolvability and voltage controllability loss in the SVS probabilistic assessment based on the Monte Carlo method. However, since the Monte Carlo method requires numerous repeated sampling calculations, the calculation time is often long. In [11] authors analyzed the influence of uncertain injected power of WFs on the SVSM and used the stochastic response surface method to evaluate the probability density of SVSM. In [12] a two-stage method to evaluate the impact of uncertain power injections on the SVSM was proposed. Based on the stochastic programming method, a method for evaluating and enhancing the SVS of power systems with uncertain wind-power output was proposed in [13]. However, the probability method usually requires statistical analysis of numerous historical data to establish an accurate probability model of WF output, which is difficult to obtain in practical applications. The authors of [14] present a novel method based on affine arithmetic (AA) for SVS assessment of power systems considering uncertainties of WF output and load power. Authors of [4] and [5] used continuous power flow and optimal power flow methods, respectively, to calculate the SVSM interval of a power system considering the uncertain fluctuation intervals of WF output. The interval method only requires the fluctuation range of WF output, which is easy to obtain from historical data. Thus, the interval method has a great application value. However, there is still little research on the optimal SVSM control method considering the uncertain fluctuation intervals of WF output.
The parametric method refers to a method of obtaining the relationship between certain dependent variables and independent variables in a power system. The parametric method is mainly divided into a sampling method and an approximation method. The sampling method refers to substituting different values of parameters into the model for solving and then using these discrete data points for analysis. A typical application of the sampling method is the Monte Carlo method in uncertainty analysis. The sampling method requires a large number of sampling calculations, so the computational load is large, and the sampling method cannot directly obtain the analytical expression, which greatly limits the flexibility of the analysis [15], [16]. The approximation method is based on the theory of function approximation and uses explicit analytical expressions for the parameters to approximate the relationships between variables. To facilitate calculations, functions with simpler forms are usually used as basis functions. In the general parametric problems of power systems, algebraic polynomials are often selected as basis functions because of their advantages of retaining nonlinear information and facilitating calculations [17]. At present, the parametric approximation (PA) method has been widely used in engineering, computer science, and other fields [18], [19], and it has also been applied in power systems. In terms of SVS, the authors of [20] took the saddle-node bifurcation point condition as the judgement equation of the SVS region of power systems and used the Galerkin method in the PA method to obtain the polynomial approximate expression of the SVS region boundary, which improved the accuracy of the original method. Based on [20], the bound of the reactive power output of generators was further considered in [21], and a more accurate SVS region boundary was obtained. The PA method has also been applied in other power system analysis fields [22]- [24]. Because the relationships between the lower and upper bounds of the SVSM interval and the control variables are relatively complicated, obtaining the simplified approximate expression through the PA method can provide an effective way to solve the multilayer optimal SVSM control problem. However, since largescale power systems have many parameters, the traditional PA method involves a large number of calculations and great computational complexity. Therefore, how to reduce the computational load of the PA method is crucial in applying the method to effectively solve the optimal SVSM control problem of actual large-scale power systems considering the uncertain fluctuation intervals of WF output.
The main contributions of this study are: (1) A bi-objective multi-layer optimization model for the optimal SVS control of a power system considering the uncertain fluctuation intervals of WF output is established. The PA method is used to obtain the approximate functional relationship between the optimal objective function values and the decision variables of the inner-layer and mid-layer optimization models and convert the optimization model into a single-layer bi-objective optimization model, which reduces the computational complexity of solving the problem. (2) Combining the PA method and the normalized normal constraint (NNC) method, a method for obtaining the continuous Pareto frontier (PF) of the bi-objective optimization model is proposed, and it can obtain a better compromise optimal solution from the continuous PF than from the discrete PF. (3) Two methods to reduce the computational load of the PA method using high-order mixed partial derivatives to remove partial cross terms and orthogonal basis decoupling are proposed, and they can effectively improve the computational efficiency.
The rest of this study is organized as follows: Section II proposes a bi-objective multi-layer optimization model for SVS control of a power system considering the uncertain fluctuation intervals of WF output; Section III converts the multi-layer optimization model into a single-layer optimization model by the PA method and proposes a method to solve the continuous PF of the bi-objective optimization problem; Section IV introduces two methods for reducing the computational load of the PA method; Section V discusses case studies carried out on the IEEE 39-bus system and an actual 964-bus provincial power grid; and Section VI presents the conclusions.

II. PROBLEM FORMULATION
To ensure that the whole SVSM interval of a power system considering the uncertain fluctuation interval of WF output can meet the requirements of secure operation, the following bi-objective multi-layer optimal control model for the SVSM of a power system is established.

A. OBJECTIVE FUNCTIONS
The objective functions of the optimal SVS control model are as follows: 1) Increase the central value of the SVSM interval. The interval central value is (λ 1 + λ 2 )/2. Since the coefficient 1/2 does not affect the optimal control result, the objective function is written as follows: 2) Reduce the radius of the SVSM interval. The interval radius is (λ 2 − λ 1 )/2. Since the coefficient 1/2 does not affect the optimal control result, the objective function is written as follows: (2) VOLUME 8, 2020 B. CONSTRAINS 1) Power flow equation at current operating point: 2) The upper and lower limits of bus voltage: 3) Generator output limit: 4) Power output limit of shunt RPC devices:

5) Lower and upper bounds of the SVSM interval
Based on the interval optimization method, the lower and upper bounds of the SVSM interval considering the interval uncertainty of WF output can be calculated by the following two bi-level optimization models [5], (10) and (11), as shown at the bottom of the page, where b p Li = P Li0 and b Q Li = Q Li0 ; and b P Gi is calculated as follows: The proposed optimal SVSM control model (1)-(11) considering the uncertain fluctuation interval of WF output is a bi-objective multi-layer optimization model. The objective functions of the outer-layer optimization problem are the sum or difference of the optimal solutions of two bi-layer optimization models, which is difficult to solve directly using traditional optimization methods. Therefore, this study uses the PA method to simplify the optimization model.

III. MODEL SOLUTION BASED ON PA METHOD
The PA method refers to the quantitative analysis of the functional relationships between multiple variables. In this problem, it refers to the functional relationships between the upper and lower bounds of the SVSM interval and the control variables. If these functional relationships can be obtained, the entire multi-layer optimization model can be transformed into a single-layer optimization model and solved directly. The PA method can obtain an explicit approximate analytical expression between multiple variables by the theory of function approximation, such as the Galerkin method [23].

A. PRINCIPLE OF THE GALERKIN METHOD
The Galerkin method essentially projects the system equations on the normed space and minimizes the residual error to obtain the optimal approximation of the target expression. The principle is described as follows: Assume that the function to be approximated is x t = x t (q). This function is determined by the implicit function F(q, x) = 0, where q is the vector of parameter variables, x is the vector of state variables, and x t is the t-th component of x. Assume that the basis functions used for approximation are ϕ(q), and the approximation of the target expression is as follows: The inner product of two functions h 1 (·) and h 2 (·) on the function space is defined as follows: Substitute the target expression (13) into the implicit function F(q, x) = 0 and use the inner product defined by (14) to project it onto each selected basis function. The projection of implicit function F(q, x) = 0 on each basis can be obtained as follows: By letting this projection equal zero, a series of projection equations can be obtained. It can be seen that the total number of equations determined by (15) is M * N , and the total number of variables c ts to be sought is also M * N . Therefore, the approximate expression coefficient c ts of the target expression (13) to be sought can be obtained by solving these equations.
In this paper, the polynomial function is taken as the basis function. In the Galerkin method, all of the control variables are continuous variables. Thus, for the discrete control variables in the optimal SVSM control model, they can be first considered as continuous variables to find an approximate functional relationship between them and the state variables. Then, when solving the optimal control model, the values of these control variables can be constrained to the corresponding discrete gears.

B. APPLYING GALERKIN METHOD TO SOLVE THE MODEL
For simplicity, in this section, when discussing the proposed multi-layer optimization model in Section II, only the objective function is preserved, and the constraints are omitted. At this time, the proposed optimization model in Section II can be written as follows:

1) PA OF THE INNER-LAYER OPTIMIZATION PROBLEM
The inner-layer optimization problem refers to the problem max x λ in (16), which is the SVSM calculation problem in the current operating state of a power system. For each value of (P G0 , U Gref , Q c , P w ), the solution λ max = max x λ of the inner-layer optimization problem is a certain value. To find the approximate parametric expression λ max (P G0 , U Gref , Q c , P w ) of the max x λ problem, the equality and inequality constraints of the inner-layer optimization problem are written as E(q, x,λ) and H(q, x,λ), and the inner-layer optimization problem is as follows: By the KKT condition, the optimal solution to the innerlayer optimization problem (17) can be transformed into solving the following nonlinear equations: where L = diag(l 1 ,. . . , l r ), U = diag(u 1 ,. . . ,u r ), Z = diag(z 1 ,. . . ,z r ), W = diag(w 1 ,. . . ,w r ). In (18), there are two inequalities that cannot be solved directly as equations. Therefore, the following equations can be added to replace the original inequalities: At this time, the equations in (18) become: As long as a group of basis functions ϕ s (q) are selected, the equations in (20) can be projected onto each basis function by using (15) to obtain the projection equations, and the VOLUME 8, 2020 coefficients of the approximate target expression λ max (P G0 , U Gref , Q c , P w ) can be solved and obtained.

2) PA OF THE MID-LEVEL OPTIMIZATION PROBLEMS
After PA of the inner-layer optimization problem, λ max (P G0 , U Gref , Q c , P w ) can be obtained, and then the entire optimization problem (16) can be transformed into two bi-layer optimization problems as in (21), shown at the bottom of the page.
Two mid-level optimization problems are as follows: Similarly, KKT conditions can be used to transform the mid-layer optimization problems (22) and (23) into nonlinear equations, and then the nonlinear equations can be projected onto each basis function. By solving the projection equations, the approximate target expressions λ 1 (P G0 , U Gref , Q c ) and λ 2 (P G0 , U Gref , Q c ) can be obtained.

3) PA OF THE OUTER-LAYER OPTIMIZATION PROBLEM AND OBTAINING THE CONTINUOUS PARETO FRONTIER
After PA of the middle-layer optimization problems, λ 1 (P G0 , U Gref , Q c ) and λ 2 (P G0 , U Gref , Q c ) are obtained. The entire optimization problem (16) can be transformed into a bi-objective programming problem as follows: Since the NNC method can obtain the uniformly distributed points on the Pareto front and supply relatively complete information for multi-objective optimization decision, the NNC method is selected in this paper. However, methods, such as NNC, can only obtain discrete PF, so some information for further decision may be lost. Therefore, based on NNC and PA methods, this paper proposes a method to obtain the approximate continuous PF of the bi-objective optimization problem.
As shown in Figure. 1, according to the NNC method, a discrete point A on the PF of the bi-objective optimization problem (24) can be obtained by solving the single-objective optimization problem shown in (25), as shown at the bottom of the page, [25], [26].
Each w corresponds to a point on the utopian line J p (w, 1-w), and each point J p corresponds to a point (J 1 , J 2 ) on the PF. By evenly taking the value of w in the interval (0, 1) to obtain a series of evenly distributed points on the utopian line, the evenly distributed discrete PF can be obtained by solving the optimization problem in (25). Therefore, as long as w is selected as a parameter, KKT conditions can be used to convert the optimization problem (25) into nonlinear equations, and these equations can be projected onto the polynomial basis functions. Then the parametric equations (J 1 ,J 2 ) on the PF and variable w can be obtained. If the direct functional relationship between J 1 and J 2 is needed, can also be obtained. Furthermore, the relationship betweenJ 1 andJ 2 can be obtained: Equation (26) is the approximate analytical expression of the continuous PF of the bi-objective optimization problem in (24).

4) DETERMINATION OF COMPROMISE OPTIMAL SOLUTION
After obtaining the continuous PF, a COS is determined from the continuous PF as the final optimal control scheme. Combined with (26), the approximate analytical expression of the continuous PF, the point that is closest to the origin is selected as the COS, as shown in Figure. 2. The COS can be determined by solving the following optimization problem: The continuous PF is obviously more complete than the discrete PF, and it provides more complete information about the coordinated optimization of the two objectives. Therefore, compared with the traditional method for determining the COS based on the discrete PF, the proposed method for determining the COS from the continuous PF can obtain a better COS.

IV. WAYS TO REDUCE COMPUTING SCALE
If the Galerkin method mentioned in Section III is applied to an actual large-scale power system, the following two problems will be encountered, resulting in complex calculations that are very difficult to solve: 1) Too many variables are required to solve. If polynomials are used as the basis function, with the increase of system parameters, the number of polynomial cross-terms will increase significantly, resulting in a large increase of the variables and equations to be solved.
2) The variables to be solved are tightly coupled. The projection equation to be solved is determined by (15). The expansion coefficients of the approximate polynomial of all system variables are tightly coupled together, so all the equations need to be solved simultaneously, which requires much RAM occupation and calculation time.
Therefore, in view of the above two problems, this paper proposes two methods (discussed in Sections IV A and B) to reduce the computing scale of the Galerkin method so that the method can be applied to practical large-scale power systems.

A. REMOVAL OF POLYNOMIAL CROSS-TERMS BY HIGHER ORDER MIXED PARTIAL DERIVATIVES
In the approximate polynomial, the non-cross terms represent the effect of an independent variable on a dependent variable, and the cross-terms represent the joint effect of multiple independent variables on a dependent variable. In the optimal SVS control problem, the joint effect of the output of two generators that are far away from each other on the SVSM may be very small, and it can be ignored. Hence, the corresponding cross-term can be deleted from the approximate target expression. Therefore, if all of the variables in a certain cross-term have a very small joint effect on the SVSM, this cross-term can be removed.
The effect of the independent variable y on the dependent variable λ o can be described by the first-order partial derivative ∂λ/∂y, and the joint effect of the independent variables y 1 and y 2 on λ o can be described by the second-order mixed partial derivative ∂ 2 λ o /∂y 1 ∂y 2 . If the second-order mixed partial derivative is zero, it can be concluded that the independent variables y 1 and y 2 have independent effects on the dependent variable λ o , and the coefficient of the cross-term y 1 y 2 in the approximate polynomial is zero. The same approach can be applied to higher order situations. Assume that the dependent variables before and after removing the cross-term are λ o and λ o , and the set of removed cross-terms is S. VOLUME 8, 2020 Then there is: From (28), as long as the minimum acceptable error d = λ o − λ o is selected, the values ∂ 2 λ o ∂y i ∂y j y i max y j max corresponding to each cross-term are ranked from small to large and summed in turn until the sum is equal to d. At that time, the remaining items are preserved. The total contribution of the cross-terms being summed is less than the minimum error d, so it can be ignored.

B. ORTHOGONAL BASIS DECOUPLING 1) PRINCIPLE OF ORTHOGONAL BASIS DECOUPLING
It can be seen from the projection equations (15) that each equation is equivalent to the integral of the product of the implicit function equations (20) and each basis function. When the rectangular form of the power flow equation is selected, the system of equations (20) represents the quadratic equations of system variables. That is, in the equations F (q, x * (q)) = 0, there are only first-order and second-order terms about x. Suppose there is: The projection equations of the first-order and quadraticorder terms on a basis are expressed as shown in (30) and (31), respectively. Whether it is the first-order term or the second-order term, its projection on a basis is related to the coefficients of all bases in its approximate expression. The coefficients of all bases are tightly coupled, resulting in the large scale of the equations to be solved. The orthogonal basis can be used to decouple the coefficients of all bases: For the first-order term, if the selected bases satisfy the orthogonal relation, then At this time, the projection of the first-order term on the orthogonal basis is as shown in (33). It is only related to the coefficient of this basis in the approximate expression, and the coefficients of other bases all equal zero due to the orthogonal relation: For the quadratic-order term, the selected bases are required to satisfy the following conditions: At this time, the projection of the quadratic-order term on these bases is as shown in (35). It is only related to the coefficient of this basis in the approximate expression: In summary, as long as the selected orthogonal bases satisfy (32) and (34), the projection of the equation system F (q, x * (q)) = 0 on the i-th basis is only related to the coefficients of the i-th basis in the approximate expressions of state variables. Thus, the (M * N )-dimension nonlinear equation system in which the coefficients of all the bases are completely coupled are decoupled into NM -dimension small-scale equation systems, that is, each basis corresponds to a small-scale equation system. Each small-scale equation system can be solved independently, and all of the coefficients can be obtained. The size of each small-scale equation system is much smaller than the original equation system. Hence, the decoupling process can greatly reduce the computational complexity and reduce the solution time. Since each small-scale equation system can be solved independently, the computational efficiency can be further improved by parallel computation.

2) SOLVING THE ORTHOGONAL BASIS
If the linear combination of the original polynomial basis is selected as the orthogonal basis, its degree of freedom cannot satisfy the requirements of (32) and (34). For equation (32), the number of constraints to be satisfied is C 2 n = n(n − 1) 2, and for equation (34), the number of constraints to be satisfied is 2 C 2 n +(n−2)/3 = n 2 +(n−2)/3. The degree of freedom the original polynomial basis can provide is only n 2 . Therefore, linear combinations of the original bases and their squares are used as the orthogonal bases. First, this can increase the degree of freedom enough to make the selected bases satisfy all requirements; second, it adds square terms, which can also improve the accuracy of PA.
Let the orthogonal basis to be solved be as follows: From (32) and (34), the equations that need to be satisfied are as follows: Since (39) corresponds to a large number of equations, the following method is used to satisfy this constraint. If the sum of some bases is equal to a constant (not zero), that is, where S c is a set of certain subscripts, then The number of all s j k dQ(s = j = k) terms is C 3 n , and each equation represented by (38) corresponds to C 2 n equations expressed as shown in (39). Assuming that the number of equations (38) to be supplemented is k, then k needs to satisfy the following conditions: At this time, the coefficients of all the selected orthogonal bases can be solved and obtained.

V. CASE STUDIES
The proposed method is verified by testing on the IEEE-39 bus system and an actual provincial power grid. The computer used for the calculation is a DELL OptiPlex 3010 workstation with a 3.60 GHz Intel Xeon processor E3-1270 and 32 GB RAM.

A. IEEE-39 BUS SYSTEM
The IEEE-39 bus system with 10 generators is shown in Figure. 3. It is assumed that a WF is connected to Bus-14, and its power output fluctuation interval is [300, 500] MW. Five shunt capacitors are installed at Buses 21-25, and their control parameters are shown in Table 1. Control parameters of the active power output and terminal voltage of all generators, except the swing bus generator, are shown in Table 2.

1) PERFORMANCE OF THE REDUCING COMPUTING SCALE METHOD
First, according to the method in Section IV A, the high-order derivatives of λ to control variables are calculated at the current operating point to remove some of the cross-terms. According to (28), the given error d is set to be 10 −4 , and the ∂ 2 λ o ∂y i ∂y j y i max y j max items are calculated. These items are ranked from small to large and then summed one by one starting from the first term until the sum is greater than the given error. It can be seen from the calculation results that when the selected polynomial order is three, the system has 253 quadratic cross-terms and 1771 tertiary cross-terms, and the total number of cross-terms is 2024. The sum of the first 1825 cross-items with smaller values is less than 10 −4 , and these cross-items can be ignored in the PA. Thus, the proposed downscaling method can effectively reduce the number of bases by 90.16% and greatly reduce the number of subsequent calculations.
Then, the remaining cross-terms and all the first-order and constant terms are used as the bases ϕ in the PA. With these bases and the orthogonal basis decoupling method proposed in Section IV B, the orthogonal bases that satisfy (37)-(39) are determined. There are a total of 199 remaining bases. According to Equation (42), the number of equations to be supplemented is calculated to be 67. Table 3 demonstrates the verification of the effect of the orthogonal basis. It can be seen that in the calculated orthogonal bases, the terms need to equal zero are close to zero, and the terms need to be reserved are obviously not zero. Therefore, it is shown that the method for solving the orthogonal basis proposed in Section IV B is effective, and the obtained orthogonal basis can satisfy the requirements of subsequent calculations.

2) PA RESULTS OF INNER-LAYER OPTIMIZATION PROBLEM
In the PA of the inner -layer optimization problem, the third-order polynomial is selected as the basis function, and its approximate expression is as follows: where s (P G0 , U Gref , Q c , P w ) is the orthogonal basis obtained above.
One thousand points of (P G0 , U Gref , Q c ,P w ) are randomly selected to verify the accuracy of the PA equation (43). For each point, the SVSM is calculated by the deterministic inner-layer optimization model. The errors between the SVSM calculated by the PA equation and the SVSM calculated by the inner-layer optimization model are calculated. The PA results of the inner-layer optimization problem when selecting the orthogonal and non-orthogonal bases are shown in Table 4. It can be seen that the PA results using orthogonal basis or non-orthogonal basis have a relatively high accuracy when the polynomial basis order is two or three. As the polynomial basis order increases, the calculation accuracy becomes higher. Although the errors of using orthogonal bases are a little higher than those of using non-orthogonal basis, the CPU times of using orthogonal bases are much shorter than those of using non-orthogonal bases. Hence, the orthogonal decoupling method can effectively improve the calculation efficiency.
The relationship between SVSM and WF output with a certain value of (P G0 , Q c , U Gref ) is shown in Figure. 4, where the value of the original model is used to calculate the SVSM corresponding to the given WF output value by the deterministic inner-layer optimization model. It can be seen that the result obtained using third-order polynomial basis functions for the PA calculation is very close to the result obtained by the deterministic optimization model, and it has high accuracy.

3) PA RESULTS OF MID-LEVEL OPTIMIZATION PROBLEMS
In the PA of the two mid -layer optimization problems, third-order polynomials are selected as the basis functions, and their approximate expressions are as follows: One thousand points of (P G0 , U Gref , Q c ) are randomly selected to verify the accuracy of the PA equations (44). The PA results of the lower bound λ 1 and upper bound λ 2 of two mid-layer optimization models are respectively shown in Tables 5 and 6. It can be seen that, whether it is the PA of λ 1 or λ 2 , after using orthogonal basis decoupling, the CPU time is significantly reduced compared with using the non-orthogonal basis. The approximate accuracy of each basis is relatively high when the polynomial basis order is two or three, indicating the feasibility of the method. The relationship between λ 1 /λ 2 and the generator terminal voltage of Bus-32 with a certain value of (P G0 , Q c ) is shown in Figure. 5. It can also be seen that the results obtained by using  third-order polynomial basis functions for the PA calculation of two mid-layer optimization problems have high calculation accuracy.

4) PA RESULTS OF OUTER-LAYER OPTIMIZATION PROBLEM
After PA of the inner-layer and mid-layer optimization problems, the entire problem is transformed into a single-layer bi-objective optimization problem, which can be solved by conventional optimization methods. Verification of the accuracy of using the PA models to obtain the optimal control scheme is executed. By solving two single-objective optimization models with objective functions (1) and (2) after PA, the optimal control schemes (P G0 , U Gref , Q c ) and the lower and upper bounds of the SVSM interval corresponding to the optimal solutions are obtained. Then the two optimal control schemes (P G0 , Q c , U Gref ) are substituted into the original models (10) and (11) for calculating the lower and upper bounds of the SVSM interval. A comparison of the calculation results is shown in Table 7. It can be seen that for two single-objective optimization problems after PA, the obtained SVSM interval value corresponding to the optimal solution is close to the SVSM interval value obtained after substituting the optimal control scheme into the original model. The errors of the lower and upper bounds of the SVSM interval are all in the level of 10 −4 . Therefore, using the PA models to replace the original model to obtain the optimal control scheme has high accuracy.
To solve the PF, the utopian line is equally divided into 20 segments (i.e., w = {0, 0.05, 0.1. . . 1}), and the discrete PF is obtained by repeatedly solving the formula of the NNC method (25) with the w value of each division point. Meanwhile, the proposed PA method is used to calculate the continuous PF, and the comparative results are shown in Figure. 6. It can be seen that the continuous PF obtained by the proposed PA method is close to the discrete PF obtained by the traditional NNC method, which verifies the correctness of the proposed PA method. Furthermore, the continuous PF obtained by the proposed PA method can provide more complete information for determining the COS than the discrete PF, which is beneficial to obtain a better COS. The method mentioned in Section III is used to determine a COS from the continuous PF as the final optimal control scheme. A comparison of the COS and the two single-objective optimal solutions and the results before optimal control is shown in Table 8. It can be seen that the distance between the COS obtained from the continuous PF and the origin is 0.1385, which is smaller than the distance 0.1438 between the COS obtained from the discrete PF and the origin. The COS obtained from the continuous PF is closer to the origin. This shows that the proposed method can obtain a better COS for bi-objective optimization problems from the continuous PF, and it has good value for practical application.
It can also be observed that the distance between the COS and the origin is much smaller than that between the optimal solutions of the two single-objective optimization problems and the origin. For the single-objective optimal solution of maximizing λ 1 +λ 2 , the value of λ 1 +λ 2 increases from 1.6956 before optimal control to 1.7838, but the value of λ 2 − λ 1 is 0.0892, which is relatively large. At this time, although the overall SVSM level of the system is relatively high, the fluctuation range of SVSM is larger as the WF output fluctuates, which brings many difficulties to the operational dispatch and security monitoring of the system. Similarly, for the single-objective optimal solution of minimizing λ 2 − λ 1 , the value of λ 2 − λ 1 reduces from 0.1312 before optimal control to 0.0425, but the value of λ 1 + λ 2 is 1.7521, which is relatively small. At this time, although the fluctuation range of SVSM is small as the WF output fluctuates, the overall SVSM level of the system is relatively low, which makes it easy for voltage collapse to occur. Therefore, the optimal solutions obtained by the two single-objective optimization problems have not achieved comprehensive optimization for the optimal SVSM interval control. On the contrary, the COS obtained by the proposed bi-objective optimization problem can optimize both objectives simultaneously: the value of λ 1 + λ 2 increases from 1.6956 to 1.7804, and the value of λ 2 − λ 1 reduces from 0.1312 to 0.0466. At this time, the overall SVSM level of the system is high, and the fluctuation range of SVSM is small as the WF output fluctuates, which is desirable since these are more suitable conditions for the actual operating point of the power system. Therefore, the COS obtained by solving the proposed multi-objective optimal control model for SVS considering the interval uncertainty of WF output can effectively increase the center value and reduce the radius of the SVSM interval. This approach can achieve a better optimal control scheme than single-objective optimization, offering good value for practical application.

5) ANALYSIS OF THE OPTIMAL CONTROL SCHEME
The optimal control scheme corresponding to the COS is shown in Tables 9 and 10. Compared with the operation state before control, the reactive power output of shunt capacitors and the terminal voltage of most generators in the optimal control scheme have increased. These control regulations can effectively supply more reactive power injection into the system and increase the overall value of the SVSM interval. In addition, the active power outputs of the generators at  Buses 32 and 34, which are close to the WF, are significantly reduced in the optimal control scheme, whereas the active power outputs of generators far away from the WF are increased. The active power outputs of generators close to the WF are reduced to supply more spinning reserves to balance the uncertain fluctuations of the WF output, which effectively reduces the fluctuation range of the SVSM interval.
From the above analysis, we observe that the proposed PA method can effectively solve the proposed bi-objective multi-layer optimization problem. Regarding the calculation accuracy, the results of the optimization problems of each layer after PA are close to the results of the original optimization problems. The error can be basically controlled in the level of 10 −4 , which shows that the PA method has high calculation accuracy. Moreover, as the polynomial basis order increases, the accuracy of the PA also increases. Hence, the error of the PA can be controlled by selecting a suitable order for the polynomial basis. Regarding the calculation speed, after removing polynomial cross-terms by high-order mixed partial derivatives and orthogonal basis decoupling, the overall calculation time for obtaining an optimal control scheme of COS is 4090 s, and it can be applied to actual power systems. For some large-scale problems, a parallel computing method can be used to solve the projection functions after orthogonal basis decoupling independently to further improve the calculation speed. Therefore, the PA method can effectively solve the proposed bi-objective multi-layer optimal SVS control problem, and it has good value for practical application.

B. ACTUAL 964-BUS PROVINCIAL POWER GRID
The actual provincial power grid contains 964 buses, 1026 branches, and 139 generators. There are three WFs with access to the MMZ21 bus, ZJC21 bus, and YJZ21 bus;  Since the approximate accuracy of the first-order polynomial is poor, only the quadratic-and third-order polynomial approximations are used for the PA, and only the orthogonal basis decoupling with high calculation efficiency is used. For the PA of the inner-layer and mid-layer optimization problems, 5000 points are randomly selected to verify the accuracy. Table 11 shows the calculation time and calculation errors between the results of the PA models and the results of the original models. In Table 11, the parallel computation was conducted using a blade cluster composed of eight DELL PowerEdge M620 computing blades, where each computing blade was composed of two 2.60 GHz Intel Xeon processors E5-2650 v2 (8 cores) and 8 GB of RAM. It can be seen that in actual large-scale power systems, the calculation accuracy of PA models is similar to that of the IEEE-39 Bus system, and the maximum errors and the average errors are also in the level of 10 −4 . This shows that the calculation accuracy of the proposed PA method has little relation with the scale of the power system, and it is mainly affected by the order of the selected polynomial basis. For the calculation time, when the selected polynomial basis order is three, after the introduction of the parameterization method and the use of two downscaling methods, the overall calculation time to obtain the COS is about 13.6 h. Since the projection functions of the inner-layer, upper bound of mid-layer, and upper bound of mid-layer optimization problems after orthogonal basis decoupling are respectively 2641 5959dimension, 2493 2268-dimension and 2377 2268-dimension small-scale equations, after using multiple threads parallel computing, 241, 250, and 238 threads are respectively used to simultaneously solve the corresponding functions, and the overall calculation time is reduced to 343.45 s. The relatively long calculation time of serial computing introduces great difficulty into practical applications. After using parallel computing, the overall computing time decreases significantly, making this method feasible to be practically applied in largescale power grids.
The PA models of the inner-layer and mid-layer optimization problems are used to solve two single-objective optimization models with objective functions (1) and (2), respectively, and the optimal control schemes (P G0 , U Gref , Q c ) and the lower and upper bounds of the SVSM interval corresponding to the optimal solutions are obtained. Two obtained optimal control schemes (P G0 , U Gref , Q c ) are respectively substituted into the original models (10) and (11) to calculate the lower and upper bounds of the SVSM interval. A comparison of the calculation results is shown in Table 12. It can be seen that for the actual large-scale power grid, the SVSM intervals of the optimal solutions obtained by the PA models are very close to the SVSM intervals obtained by substituting the optimal control scheme into the original model. The errors are in the level of 10 −4 . Therefore, the proposed PA method is still feasible for application in actual large-scale power grids, and the PA model can replace the original model for optimal control calculation of SVSM interval. Figure 8 shows the comparison of continuous PF and discrete PF. Table 13 shows the comparison of SVSM intervals at some decision points. The following two conclusions can be drawn from the calculation results. First, the proposed method can effectively obtain an optimal control scheme in the actual large-scale power grid. In the optimal control scheme of COS, the value of λ 1 +λ 2 increases from 1.0443 to 1.0700, and the value of λ 2 − λ 1 decreases from 0.0963 to 0.0795. This shows that after optimization, the overall SVSM  level of the large-scale power grid rises, and the fluctuation range of SVSM decreases as the WF output fluctuates. Hence, the system operates in a more secure state. Second, the COS obtained by the continuous PF is closer to the origin than the COS obtained by the discrete PF. Hence, the proposed method can obtain a better optimal control scheme in a large-scale power grid, and it has good value in practical application.

VI. CONCLUSION
In this study, a bi-objective multi-layer optimization model for SVS control of a power system considering interval uncertainty of WF output is proposed. A PA method is used to first obtain the approximate functional relationship between the optimal objective function values and decision variables of the inner-layer and mid-layer optimization models, and then convert the original optimization model into a single-layer bi-objective optimization model that can be directly solved. A method to obtain the continuous PF of the bi-objective optimization model is proposed, combining the traditional NNC method and the PA method, and the obtained continuous PF can supply more complete information for the optimal control decision and obtain a better COS to be used as the optimal control scheme. To reduce the computational load and improve the calculation efficiency, two methods for reducing the scale of the PA method are proposed, which can effectively reduce the calculation time and make the proposed method more suitable for practical application.
The calculation results of case studies show that the obtained optimal control scheme can effectively increase the overall SVSM level and reduce the fluctuation range of the SVSM of the power system as the WF output fluctuates. Hence, the power system can operate in a more secure state. The obtained approximate expression of the proposed PA method has good accuracy. Furthermore, a better COS can be obtained from the continuous PF, and the two methods for reducing the calculation scale of the PA method can effectively improve the calculation efficiency of the proposed method. MINGBO LIU (Member, IEEE) received the B.S. degree from the Huazhong University of Science and Technology, in 1985, the M.S. degree from the Harbin Institute of Technology, in 1988, and the Ph.D. degree from Tsinghua University, in 1992. He is currently a Professor with the South China University of Technology. He has authored and coauthored four monographs, two standards, and more than 280 articles. His research interests include energy management and operation control of power systems. VOLUME 8, 2020