Identifying Multiple Local Optima for Optimal Power Flow Based on Nonlinear Dynamics

An optimal power flow (OPF) problem of power systems can have multiple local optimal solutions, which is worthwhile studying both in theory and practice. Based on the existing nonlinear dynamic systems, this paper proposes an efficient deterministic algorithm to solve multiple or all local optimal solutions of OPF, which takes some numerical improving measures to enhance the numerical convergence for integration process of dynamics and adapt to OPF problem. The steps of this algorithm are as follows: 1. The reflected gradient system (RGS) is used to calculate the decomposition points to locate different feasible components. 2. The quotient gradient system (QGS) is used to calculate feasible points in different feasible components, and we numerically integrate projected gradient system (PGS) with these feasible points as initial points forward until the trajectories approach the local optima. 3. Slack variable perturbation method (SVPM) is proposed to help escape from the saddle points to the adjacent local optima when the trajectories fall into saddle points. Compared with the interior point method (IPM) with random initialization, multiple IEEE test cases show that the proposed algorithm can identify much more local optimal solutions, and meanwhile, significantly reduce the calculation time.


I. INTRODUCTION
The optimal power flow (OPF) problem is an important issue in power system operation, planning and control. The purpose of OPF is to optimize the objective function (such as generator operation cost, network loss, voltage offset, etc.) by adjusting the active power output of the generators and other control variables. Meanwhile, the constraints of power balance and safe operation of the power systems should be strictly satisfied.
Because OPF problem usually has many constraints, especially the existence of power flow equality constraints makes the feasible regions non-convex and disconnected, resulting in more than one optimal power flow solution [1]. Therefore, identifying multiple or all local optimal solutions can better find and recognize the global optimal solution, and enhance the confidence of the result as the global optimal solution [2]. The complex behaviors of power systems such as load changes [3], bifurcation [4] and non-convexity [5] The associate editor coordinating the review of this manuscript and approving it for publication was Easter Selvan Suviseshamuthu . of feasible regions, are often accompanied by changing, generating and disappearing of the local optima. Therefore, identifying and analyzing the local optima for OPF is helpful to better understand the operating conditions of power systems, and help us to better understand OPF, which is a complex constrained optimization problem [2]. In addition, many numerical optimization methods have difficulty converging to the local optimal solutions due to the stagnation points such as saddle points and maximum points [6], especially for OPF problem which has many saddles. Hence an effective method is needed to help escape from the saddle points to adjacent local optima when optimization algorithms converge to the saddle points. Therefore, it is of great significance in theory and practice to efficiently identify multiple or all local optimal solutions for OPF problem.
Many optimization algorithms are proposed to solve the OPF problems, including numerical calculation methods such as interior point method [7], sequential quadratic programming method [8], etc. And heuristic algorithms such as particle swarm optimization [9], genetic algorithm [10], etc. These algorithms can guarantee the local optimality of the solutions, but cannot guarantee the global optimality of the solutions, nor can they find all or most of the local optimal solutions deterministically within a valid time. In reference [3], repeated calculations are carried out by many NLP solvers, which is cumbersome. In reference [2], branch tracing method and monotone search strategy are proposed to locate multiple local optimal solutions of the OPF, but the rules are complicated and will result in mass computation, so it is difficult to implement when the scale of nodes of test cases is up to dozens. At present, there are few papers about how to solve the multiple or all solutions for OPF accurately and efficiently.
Nonlinear dynamics is a science for studying the law of changes in variety of systems. Since the 1970s, people began to solve the optimal solutions of nonlinear optimization problems or the root of nonlinear algebraic equations by tracing the solution trajectories of ordinary differential equations [11], which was called the trajectory method. In reference [12]- [15], a series of nonlinear dynamic models have been established to solve nonlinear optimization problems. In reference [12]- [13], a gradient vector system has been established based on KKT condition, and a nonlinear optimization method was initially developed. In reference [14], a quasi-gradient system has been established, and then a method for solving a class of unconstrained optimization problems has been proposed. In reference [15], by alternating between two dynamic models, multiple or all local optimal solutions with constrained optimization can be solved, which was verified in a small numerical example, but large-scale engineering applications still remain to be further investigated. In reference [16]- [18], the nonlinear dynamics methods have been used to solve the optimization problems of power systems. In reference [16], the reactive power optimization problem has been solved, and in reference [17], [18], the active power optimization problem has been solved, but the upper limit of objective function needs to be adjusted many times, resulting in mass computation, and it is also unable to calculate multiple or all local optimal solutions deterministically.
The optimization algorithms based on integration [12]- [15] are different from these based-on iteration represented by the interior point method, but the engineering applications are insufficient. In this paper, this optimization idea and method are applied to identify multiple local optima for OPF, and some measures are adopted to improve the convergence of numerical integration and make it more appropriate for OPF problems, and then, a deterministic algorithm for efficiently calculating multiple or all local optimal solutions for OPF is formed: the reflected gradient system (RGS) is used to calculate decomposition points x d , so as to locate multiple or all feasible components of OPF problem, which is implemented by integrating the quotient gradient system (QGS) with exit points X ex ,so that we can calculate the initial points in different feasible components, with which we could integrate projected gradient system (PGS) [15] numerically forward until the trajectories approach to adjacent local optima, and for more details, see section IV. When the trajectories fall into saddle points, slack variable perturbation method (SVPM) can help escape from saddle points to adjacent local optima. Multiple or all local optima for OPF can be identified by alternating between these procedures. The algorithm process is shown in Figure 1.
The crucial contributions of this paper are presented as follows: 1) Apply the dynamic models proposed in reference [12]- [15] to the OPF problem, the challenging engineering application, and some improvement measures are taken to enhance the numerical convergence property for integration process of dynamics and make it appropriate to the characteristic of OPF problem. 2) Based on the theory of stability region, we combine several important dynamic models organically and propose a deterministic algorithm, which can identify multiple or all solutions for OPF efficiently.
3) The slack variable perturbation method (SVPM) can help optimization algorithms escape from saddle points and reach the local optima again more effective, and it is especially suitable for OPF problems, which have many saddle points.
This paper is organized as follows: In section II and III, the mathematical model of OPF and KKT condition are briefly introduced respectively; in section IV, the crucial dynamic model, projected gradient system is introduced, and several methods to enhance the numerical convergence property of PGS are proposed, and then, the deterministic optimization method is proposed; in section V, the test results of IEEE cases in comparison to interior point method (IPM) are shown and discussed. Finally, the conclusion and research prospects are illustrated in section VI.

II. MATHEMATICAL MODEL OF OPF
In this paper, the objective function is the minimal total generation cost. The OPF problem in polar coordinate can VOLUME 8, 2020 be described as follows: where S nB , S nG , and S nL represent the set of nodes, generation units and transmission lines respectively. P Li and Q Li denote are the active and reactive load demand of the node i. G ij and B ij are the conductance and susceptance of the line between the node i and j. Voltage phase angle of the slack bus is set to 0. θ ij = θ i − θ j , which represents the voltage phase angle difference between node i and j. V i max and V i min are the upper and lower voltage limits of node i. P Gi max , P Gi min , Q Gi max and Q Gi min are the upper and lower limits of the active and reactive output of generator i respectively. S fi and S ti are the apparent power at the head and end of the line i respectively. S fi max and S ti max are the upper limit of the apparent power at the head and end of the line i respectively. By introducing the slack variable vector s for the inequality constraints, (1) is transformed into the following compact form [1]: H : 6nG+4nB+2nL → 4nB+4nG+2nL+1 is the constraint set, where g : 2(nB+nG) → 2nB+1 is the set of equality constraints, h : 2(nB+nG) → 4nG+2nB+2nL is the set of inequality constraints.

III. KKT CONDITIONS FOR DEALING WITH CONSTRAINED OPTIMIZATION PROBLEMS
For the constrained optimization problem (2), lagrangian multiplier vector λ ∈ 4nG+4nB+2nL+1 is added, and then (2) can be transformed into an unconstrained optimization problem as shown in equation (3): where F : 10nG+8nB+4nL+1 → 1 . The gradient of F is: DH(X) represents the Jacobian matrix of H(X). Let (X * , λ * ) be the point satisfying (4), which is the stagnation point (KKT point) of problem (1). KKT points includes the local optima, local maxima and saddle points. It is necessary to further verify the second-order sufficient condition [1] to distinguish the local optima from these stagnation points. Calculate the Hesse matrix of (3): The second-order sufficient condition is simplified by transforming (1) into the compact form (2): for any vector u in the set = u = 0 : DH(X * )u = 0 , which is the basis of the tangent space of the constraint set H(X) at X * , if u T ∇ 2 F(X * )u > 0, ∀u ∈ , then X * is a strict local optimal solution of problem (2).
2) Every local optimal solution of optimization problem (2) corresponds to a stable equilibrium point of system (7).
From these properties, we can search local optima of (2) by integrating system (7) with initial points until the trajectories approach to stable equilibrium point of (7).
According to property 1), the initial value of system (7) can be selected at random, but it is different from the simple mathematical examples mentioned in [12]- [15]. In OPF problems, the dimension of the vector (X, λ) is generally high, and the orders of magnitude of the elements of (X, λ) are significantly different, so the integration trajectories of (7) tend to diverge due to the limitations of numerical integration method without any processing, although it is theoretically globally convergent. In order to promote the numerical convergence, the following methods are adopted:

B. METHODS TO IMPROVE THE NUMERICAL CONVERGENCE OF PGS (7) 1) SELECT THE APPROPRIATE INITIAL VALUES OF INTEGRATION FOR (7)
If the integral initial values of (7) is selected at random, more divergence and lower efficiency will occur in practice.
Selecting the feasible points of (2), which satisfy H(X) = 0, and it can make (7) converge quickly and steadily. The quotient gradient system (QGS) [15] has good robustness and strong search capabilities, which can quickly generate widely distributed feasible points. Therefore, in this paper, the QGS is used to obtain the feasible points as the integral initial points for PGS (7).
For problem (2), the QGS is established with the constraint set: The feasible region (FR) of OPF is often disconnected [1], where FR i is the i th feasible component, and N is the number of all feasible components. Reference [1] proves that every regular stable equilibrium manifold of QGS (8) is corresponding to one feasible component of (2). And (8) has good robustness, which can be integrated forward with any initial points (infeasible) to the stable equilibrium points of (8), that is, the feasible points of (2). With the increase of search times, all the stable equilibrium manifolds of (8) will be gradually found, which is the complete feasible region of (2). Moreover, (8) is globally convergent, and each trajectory will converge to an equilibrium point of (8), which guarantees the efficiency of generating feasible initial points of (7).
In reference [19], the QGS is used to calculate the complete feasible region of OPF. But in this paper, it is not necessary to obtain the complete feasible region, but just need to acquire the integral initial points of (7) lying inside as many disconnected feasible components as possible, which can expand the search space for local optima. To accomplish this, besides relying on the search capability of QGS, it is also necessary to find the decomposition points (DP) of QGS, which connects different feasible components, so as to ensure that feasible points in multiple or all feasible components of (2) can be found.
The type-1 unstable equilibrium points located on the boundaries of the stability domains of the stable equilibrium manifolds of QGS (8) are called the decomposition points (DP), represented by x d in Fig.1. In reference [15], it is proved that the unstable equilibrium manifolds of DP converge to two adjacent stable equilibrium manifolds of QGS (8). From this property, the disconnected feasible components of (2) can be found through the DP of QGS (8). The methods for finding the type-1 unstable equilibrium points include eigenvector tracing method [21], homotopy-based method [22], etc. In this paper, eigenvector tracking method is used to find x d of (8).
The eigenvector tracing method can find the type-1 unstable equilibrium points of QGS (8) by the reflected gradient system (RGS). By establishing RGS (9), the type-1 unstable equilibrium points of (8) are transformed into the stable equilibrium points of RGS (9). For the n-dimensional vector X, let λ 1 (X) ≤ λ 2 (X) . . . ≤ λ n (X) be the eigenvalue of the Jacobian matrix of X at QGS (8), and v i (X) be the eigenvector corre- RGS is defined as the following nonlinear dynamic system: where DQ H (X) represents the Jacobian matrix of Q H (X). The selection method of the integral initial points of (9) is generally empirical [22], which is related to the structure of the dynamic system under study. For simplicity, in this paper, the integral initial values of (9) are chosen randomly. The searching method of type-1 of QGS (8) is as follows: Numerically integrate (9) forward with an random initial point.
1) If the integral trajectory converges to the equilibrium point of QGS (8), then judge whether is the type-1 unstable equilibrium point of the system (8) and filter DP. The judgment method is to calculate the number of positive real part eigenvalues corresponding to eigenvectors of DQ H ( ) on the normal space N x ( ) (orthogonal complement of tangent space T x ( )) [1], and the details are as follows: A, and if J Q has k eigenvalues of positive real parts, then is the typek unstable equilibrium point of (8), when k = 1, it is what we desire.
The method for filtering DP from {x d } [14] is: start from x d and move a small distance ε along the two directions of the unstable eigenvector E u d of DQ H (x d ) to obtain the points x d ± εE u d , and take x d ± εE u d as the initial values to integrate (8),if they converge to different stable equilibrium points, then x d is DP.
When we get the DP of QGS (8), we can search the feasible points in different feasible components through the known feasible points and DP, so that the integral initial values of PGS (7) are widely distributed and the searching space for local optima is vast. The method [24] is as follows: let X ex = X 0 + σ (DP − X 0 ) where X 0 is the point inside the known feasible component, X ex is named exit point [14], and σ is greater than 1 to ensure that X ex falls into the stability domain of the adjacent stable manifold. Integrate (8) forward with X ex until the trajectory approaches the stable equilibrium point X 1 , where X 1 is the new feasible point inside the adjacent feasible component (see Figure 1).

2) IMPROVE THE SINGULARITY OF DH(X)
Slack variable vector s is introduced so as to deal with the inequality constraints, and s is not dependent, which satisfies h(x) + s 2 = 0. If we update s with (7) as if it is independent without any processing method, case study results indicate the poor convergence.
The method is as follows: at each step of integration for (7), if the i th inequality constraint h i < 0, then let VOLUME 8, 2020 −h i ,which can keep the (2nB +1+ i) th element of H(X) satisfying h i + s 2 i = 0, and enhancing the convergence. When h i ≥ 0, let s i = δ, instead of 0, δ is small positive real number, the reasons are as follows: At each integration step of (7), when the row of DH(X) is not full rank, formula (7) will diverge because of the singularity of DH(X)DH(X) T . Besides, the first order derivative of h i +s 2 i with respect to s i is 2s i . When k th inequality constraint violate the limit, Which will not only improve the singularity of DH(X) but also basically satisfy the feasible conditions of (2), so as to enhance the numerical convergence of (7).

3) SLACK VARIABLE PERTURBATION METHOD
Numerically integrate (7) with different feasible points of (2) as the integral values until the trajectories converge to the adjacent local optima or saddle points (SD), and then the integration stops. The saddle points are the unstable equilibrium points of PGS (7), and their stable manifolds form the boundaries of the stability domains of the stable equilibrium points [15]. OPF problem usually has many saddle points [2], which causes much difficulties when Newton method is adopted to solve KKT equations (4), because Newton method is easy to fall into saddle points, especially for largescale cases. According to theory of stability regions [14], the stable manifolds of the saddle points are much smaller in the whole space than that of stable equilibrium points. Therefore, when dynamic method is adopted to solve OPF problem, the trajectories are more likely to converge to the local optima of (2) (stable equilibrium points of (7)) than converge to the saddle points. Even so, some trajectories will still fall into the saddle points, which is meaningless for optimization. Reference [15] proposed the eigenvector perturbation method to escape from saddle points, but it is tedious to calculate multiple unstable eigenvectors of the Jacobian matrix of (7) at SD. In this paper, a simple and effective method to escape from the saddle points is proposed: slack variable perturbation method (SVPM).
At every integration step of (7), s is updated as follows: where ε ∈ (0.99, 1), δ is a small positive number as described in 2). The saddle points (SD) are unstable equilibrium points, of PGS (7), and SVPM is adopted to give the value of (7) at SD a perturbation by changing the value of ε, and then the PGS (7) on the right generates a mismatch, which means max(P(SD)) = 0, and then numerically integrate (7) forward with initial point SD. During the first few integration steps, s is updated according to (11). After a short distance of integration, s is updated according to the method in 2). Due to the global convergence of (7), the mismatch value will eventually be reduced to 0. Meanwhile, the integration trajectory of (7) approaches to one adjacent local optimal solution of OPF problem (1) which is stable equilibrium point of PGS (7). The numerical experiments of many different scale test cases show that trying different ε for 1-2 times can guarantee to escape from saddle points to adjacent local optima. Here, the integration trajectories of (7) falling into saddle points can be replaced by any optimization algorithms' trajectories converging to saddle points which violate the second-order sufficient condition, and then the SVPM can help them escape from the saddle points and reach the local optima again.
Based on the discussion above, X 0 can be selected as the integral initial value of system (7), where: FR is the feasible region of (2). Measures proposed above in this section are taken to enhance the numerical convergence of PGS (7) and help the integral trajectories escape from saddle points. The algorithm flow chart is shown in Figure 2.

V. NUMERICAL EXAMPLES
Multiple IEEE test cases are used to verify the effectiveness of the proposed algorithm, which is compared with the interior point method (IPM). The ode15s integrator in MATAB is adopted for the integration of QGS and RGS, and the hybrid integration method in reference [15] is adopted for the integration of PGS (7), which combines the convergence speed of Newton method and the robustness of the steepest descent method. The numerical measures proposed in B of section IV are taken to enhance convergence and adapt to OPF problem at each integration step of PGS (7). The parameters, wiring diagram and optimal solution information of the cases can be found on the website [23]. For the convenience of description, in the following figures, INI stands for the integral initial points of system (7). X ex stands for the exit points. DP stands for the decomposition points. And LOS stands for the local optimal solutions. SD stands for saddle points. And subscript numbers are used to distinguish each point. FR stands for feasible region.

A. CASE STUDY 1) LMBM3 TEST
The three-bus case contains 5 local optimal solutions and 2 disconnected feasible components. In order to show the relative positions of the local optima and integration trajectories, the complete feasible region (grey area) of this case is characterized in Figure 3 using the method in [19], and the optimization process is also shown in Figure 3. The green curves are the trajectories of searching the decomposition points, and the red curves are the trajectories of searching different feasible components from the exit points, the blue curves are the trajectories of identifying the local optima, and the local optima are marked with red crosses. In the following figures, the axis label P 1 means the active power of the first generator, and so on; and the axis label Q 1 means the reactive power of the first generator, and so on.
The optimization process is as follows in detail: Loop 1: Integrate (8) with random initial point to feasible point INI4, and integrate (7) with INI4 as the initial value to local optimal solution LOS4. Move a distance along the line (dashed line) between INI4 and the decomposition point DP1 (cyan diamond) to the exit point X ex1 . And then integrate (8) with X ex1 as the initial point to the feasible point INI3 (in another feasible component, marked in yellow). Integrate (7) to another local optimal solution LOS3 with INI3 as the initial value. Then extend the line between INI4 and the decomposition point DP2 to the exit point X ex2 and integrate (8) with X ex2 as the initial value to the feasible point INI2, and integrate (7) (7) with INI5 as the initial value to LOS5, and all local optima have been identified.
When the integration trajectories of (7) converge to the saddle points, the slack variable perturbation method can help the trajectories escape from the saddle points and reach the adjacent local optimal solutions. As shown in Figure 4, the red curves are the trajectories of falling into saddle points, and the blue curves are the trajectories of escaping from the saddle points to local optimal solutions. Integrate (7) with INI1 and INI2 and the trajectories fall into saddle points SD1 and SD2 respectively. By SVPM to change ε for perturbation, two trajectories can escape from the saddle points to the adjacent local optimal solution LOS1. The red trajectories in the figure can also be replaced by the trajectories falling into saddle points of other numerical algorithms (such as the interior point method), which can be improved by SVPM.
In order to test the convergence of the algorithm, the method in [19] is adopted to characterize the complete feasible region of the 3-bus case, with a total of 10,000 feasible points as initial value to integrate (7), so as to identify corresponding local optimal solutions. Among them, 112 points FIGURE 4. The trajectories escaping from saddle points to local optima. VOLUME 8, 2020 converge to the saddle points, which all escape to adjacent local optima through SVPM, and the rest all converge to local optima. And then characterize the convergence domain of each local optimal solution with different colors, whose boundary are smooth, as shown in Figure 5. The PGS (7) presents adjacent convergence and global convergence, that is, the feasible points will converge to the adjacent local optimal solutions, and each feasible point will converge to an equilibrium point of system (7), and most points will converge to the stable equilibrium points of (7), which are local optima of (2). This property enables the algorithm to identify multiple or all local optima for OPF efficiently and deterministically.

2) WB5MOD TEST
The 5-bus case has 4 disconnected feasible components and 5 local optimal solutions. The optimization process is shown in Figure 6.
The green curves are the trajectories of identifying the decomposition points, and the red curves are the trajectories of searching different feasible components from the exit points, and the blue curves are the trajectories of identifying local optima. And the gray area is the feasible region.
The optimization process is as follows in detail: Loop 1: Integrate (8) with random initial point to the feasible point INI3, and then integrate (7) with INI3 as the initial value to the local optimal solution LOS3. Integrate (7) with feasible points INI2, INI1 and INI4 (located in different feasible components) as initial points, which are acquired by integrating (8) with corresponding exit points X ex1 , X ex2 ,and X ex3 searched by extending the dashed line between INI3 and DP1, DP2, and DP3, to local optimal solutions LOS2, LOS1 and LOS4 respectively.
Loop 3: Integrate (8) with random initial point to the feasible point IN5, and integrate (7) with INI5 as the initial  value to the local optimal solution LOS5. And all local optima have been identified. Figure 7 shows the process of escaping from the saddle points to adjacent local optima. The complete feasible region is characterized, including 10000 feasible points, and we numerically integrate (7) with them as initial values to identify corresponding local optima. Among them, 108 points converge to saddle points, which all escape to adjacent local optima through SVPM, and the rest converge to local optimal solutions.
The convergence domain of each local optimal solution are characterized with different colors, as shown in Figure 8, which shows the fine global convergence and adjacent convergence of PGS (7).

3) CASE9MOD TEST
The 9-bus case has 3 disconnected feasible components and 4 local optimal solutions. The optimization process is shown in Figure 9.
The green curves are the trajectories of identifying the decomposition points, and the red curves are the trajectories of searching different feasible components from the exit  points, and the blue curves are the trajectories of identifying the local optima. And the gray area is the feasible region. The optimization process is as follows: Loop 1: Integrate (8) with random initial point to the feasible point INI1, and then integrate (7) with INI1 as the initial value to the local optimal solution LOS1. Integrate (7) with feasible points INI2 and INI3 (located in different feasible components) as initial points, which are acquired by integrating (8) with corresponding exit points X ex1 and X ex2 searched by extending the dashed line between INI1and DP1, DP2, to local optimal solutions LOS2, LOS3 respectively. Loop 2: Integrate (8) from random initial point to the feasible point INI5, and then integrate (7) with INI5 as initial point to the local optimal solution LOS2. Integral (7) with INI4 (located in different feasible component), which are acquired by integrating (8) with exit points X ex3 searched by extending the dashed line between INI5 and DP3, to the local optimal solution LOS4. Figure 10 shows the process of esfrom the saddle points to adjacent local optima. The complete feasible region is characterized, including 10000 feasible points, and we numerically integrate (7) with these feasible points as initial points to identify corresponding local optima. Among them, 98 points converge to the saddle points, which all escape to adjacent local optima through SVPM, and the rest all converge to the adjacent local optimal solutions. Figure 11 shows the convergence region of each local optimal solution with different colors.

4) CASE118MOD TEST
It should be noted that in large-scale cases, it is very difficult to calculate the decomposition points [15]. The QGS has good robustness and strong search ability, which can generate widely distributed feasible points, and as the number of searching increases, it can efficiently search for feasible points located in multiple or all feasible components. Therefore, the step of searching the decomposition points can be FIGURE 11. The convergence domain of each local optimal solution. VOLUME 8, 2020 omitted. Choose integral initial points randomly as dispersedly as possible in the whole space, for example, by Latin hypercube sampling [19], with which integrate QGS (8) to the feasible points of (2), and then integrate (7) with these feasible points until the trajectories approach the local optima of (2). It takes a little longer time to implement this method, but the optimization effect can also be guaranteed.
There are two disconnected feasible components and three local optimal solutions in the 118-bus case. The optimization process is shown in Figure 12. After the QGS calculating 39 feasible points of (2) and repeating the process above, all 3 local optimal solutions are identified, and 33 points converge to LOS1, 2 points converge to LOS2, and 3 points converge to LOS3, and 1 point converges to saddle point (cyan). The trajectory falling into the SD is marked in red and the trajectory escaping from the SD to LOS1 by SVPM is marked in blue and amplified in Figure 12.

B. ALGORITHM COMPARISON
This section, we will compare the calculation effect of the proposed algorithm with interior point method (IPM). The IPM solvers adopted the MIPS solver based on the primaldual interior point method and the IPOPT solver based on the prediction-corrected interior point method. The all 18 test cases which are from school of mathematical, university of Edinburgh (see website [23]) are used to verify the effectiveness of the proposed algorithm, which contains many types and is extensive and persuasive. For IPM solvers, the maximum number of experiments for each case is 10000 (including the number of divergence and convergence to saddle points), and if all the local optimal solutions cannot be found within 10000 times or all are found, then experiment will stop. For MIPS and IPOPT, iteration calculation from the random initial point to the stagnation of (4) or to infeasible points(divergence) is considered as an experiment. For the proposed algorithm, integration calculation from the random initial point to one equilibrium of (7) is considered as an experiment, if all local optima are found or the same solution is found for consecutive times, then experiment will stop. The total calculation time includes the time of calculating the decomposition points and integrating QGS and PGS. The results are shown in Table 1. It can be seen from Table 1 that the proposed algorithm can identify all or most of the local optimal solutions of the above cases within the effective time. For some cases, such as WB3, WB5 and case30, the calculation effect of the proposed method is about the same as that of IPM, but for most cases above, especially for the cases with many extreme points and large-scale cases, the proposed method can identify far more local optima (see deep red column in Table 1) within much less time (see pink column) compared with IPM.
The qualitative interpretation is as follows: it can be seen from detailed case analyses in A of section V that the convergence regions of the local optima are continuous, with similar sizes and smooth boundaries, which are different from the fractal structure of Newton method [25], and there is only one local optimal solution in each convergence region [14], which is regular and orderly. The inherent good convergence of PGS (7) and the numerical improvement measures can ensure the trajectories of (7) approach to local optima as long as the initial feasible points fall into the corresponding convergence regions. Besides, the decomposition points can help locate all disconnected feasible components to ensure the initial feasible points are widely distributed so that the search space of local optima can be expanded. Moreover, the good robustness of system (8) can also generate widely distributed feasible points with random initial points and expand the search space. Besides, stability region theory [14]- [15] indicates that the trajectories of system (7) are much more likely to converge to local optima than converge to saddle points and other invalid points, and the SVPM can help escape from saddle points to adjacent local optima, even if the trajectories of (7) converge to saddle points. The good property of related dynamic systems and design of the proposed method are conducive to identifying most or all local optima in a shorter time compared with IPM with random initialization.

VI. CONCLUSION
In this paper, the method based on nonlinear dynamics is applied to identify multiple local optima for OPF problem. Combined with the RGS to locate the disconnected feasible components, the QGS to generating feasible points that are widely distributed, and the PGS (7) which has fine global convergence and adjacent convergence, as well as the numerical improvement measures to enhance the numerical convergence of PGS (7) and adapt to OPF problem, this algorithm can efficiently identify multiple or even all local optima within effective time. The slack variable perturbation method enables the optimization algorithms to escape from the saddle points and reach the local optima again, especially suitable for OPF problem which has many saddle points. Multiple IEEE test cases verify the effectiveness of the proposed method.
Further researches will focus on the following three aspects: Firstly, this method can efficiently handle the optimization problems, which are differentiable. For problems with nondifferentiable objective functions and constraints, the subdifferential or the generalized gradient [26] can be adopted to form the elements of Jacobian matrix DH(X) at the nondifferentiable points, which will enlarge the application range of the proposed algorithm.
Secondly, this paper mainly focused on the design of the optimization, so the OPF model (1) adopts the concise and classical mathematical form, which is relatively preliminary. In future studies, various OPF problems need to be considered. For example, considering reactive power optimization, dynamic economic dispatching, and the optimization of power systems with renewable energies represented by wind power and photovoltaics.
Thirdly, for constrained optimization problems, with the unified form (2), the various swarm algorithms which have good global optimization characteristics can be adopted to search the promising solutions, and on the basis, the proposed method which have good adjacent convergence property can search the higher-quality local optimal solutions through fine tuning. This hybrid design can achieve better optimization effect [24].