Parameter Optimization of Droop Controllers for Microgrids in Islanded Mode by the SQP Method with Gradient Sampling

For enhancing the stability of the microgrid operation, this paper proposes an optimization model considering the small-signal stability constraint. Due to the nonsmooth property of the spectral abscissa function, the droop controller parameters’ optimization is a nonsmooth optimization problem. (e Sequential Quadratic Programming with Gradient Sampling (SQP-GS) is implemented to optimize the droop controller parameters for solving the nonsmooth problem. (e SQP-GS method can guarantee the solution of the optimization problem globally and efficiently converges to stationary points with probability of one. In the current iteration, the gradient of the nonsmooth function can be evaluated on a set of randomly generated nearby points by computing closed-form sensitivities. A test on the microgrid system shows that the optimality and the efficiency of the SQP-GS are better than those of the heuristic algorithms.


Introduction
Microgrid is an alternative solution to integrate distributed energy resources (DERs). It concludes a cluster of loads, distributed generation (DG), and energy storage systems, just like a small-scale power supply network [1,2]. Due to the negligible physical inertia of power-electronic converters, the microgrid operation is susceptible to random disturbances, such as variability of DG and random load fluctuations [3], like charging for electric vehicles [4], especially in an islanded mode. e small-signal stability of an islanded microgrid is critical for its reliable operation. In an islanded microgrid, droop control, following a predefined droop coefficient, is widely used to maintain the proportional load sharing between DERs [5]. e small-signal stability of the microgrid in islanded mode is strongly affected by the droop controller parameters. erefore, it is essential to optimize the controller parameters to enhance the small-signal stability of an islanded microgrid.
In [6], the maximum of the real part of eigenvalues, also called the spectral abscissa, was minimized to enhance the stability in grid-connected mode. On the other hand, a nonlinear time-domain simulation-based objective function was employed to minimize the error in measure power in islanded mode. e study in [7] proposed a method for optimizing the controller parameters of dynamic droop control for a microgrid in islanded mode and used a spectral abscissa objective function to consider small-signal stability. e study in [8] tried to optimize proportional-integral (PI) current controller parameters considering three design criteria and multiple operating conditions which include grid-connected mode and different islanded modes by a weighted sum of multiple objective functions. e work in [9] proposed an optimization model with a nonlinear timedomain simulation-based objective function for a microgrid in islanded mode. In addition to PSO methods, the study in [10] used the Genetic Algorithms (GA) method to optimize the controller parameters by minimizing the power difference-integrating during the changing process between gridconnected mode and islanded mode.
All the studies above adopt the heuristic algorithms to optimize the control parameters of the microgrid due to the nonsmooth property of the optimization model. Although those algorithms can obtain fast and feasible solutions, they are hard to deliver an optimal solution to the vast majority of cases [11,12]. e solutions derived by the heuristic algorithms are not deterministic and reproducible. ese disadvantages will make the debugging and correctness checking difficult, impeding the industrial application. Concerning the optimization model, the study in [6] did not consider the small-signal stability in islanded mode. As for considering stability, the small-signal stability analysis was performed in [8][9][10] to obtain the optimum ranges of the controller parameters, while the study in [7] employed spectral abscissa as an objective function to enhance small-signal stability better. e study in [8] considered multiple operating conditions in the optimization model, which may cause the certain index of the design criteria to be poor due to the considering way.
To tackle the weakness mentioned above, this paper proposes an optimization model considering the spectral abscissa in multiple islanded mode operating conditions. Besides, an SQP method combined with GS (SQP-GS) [13] is presented to solve the proposed optimization model. erefore, droop controllers with optimized parameters can enhance the small-signal stability of the microgrid in islanded mode. As for the SQP-GS method, it had been adopted to optimize the parameters of PSSs and PODs [14]. e major contributions can be summarized as follows: (1) Based on SQP-GS, a mathematical programming method is presented to solve the proposed nonsmooth optimization model and obtain an optimal and deterministic solution without any convergence problem. (2) e small-signal stability in different operating conditions is enhanced using a nonsmooth objective function in the optimization model. It can ensure that optimized parameters are robust under multiple operating conditions. e remainder of this paper is organized as follows: Section 2 presents the small-signal stability model of the microgrid. Section 3 proposes the optimization model for optimizing the droop controller parameters of the microgrid. Section 4 introduces the SQP-GS method to solve the optimization model. In Section 5, the test on the microgrid system under multiple operating conditions is presented. Conclusions are drawn in Section 6.

The Small-Signal Stability Model of
the Microgrid e small-signal stability model of the microgrid can be divided into three parts: inverter, network, and load. ese small-signal stability models are deduced in dq reference frame [10,15].

2.1.
e Model of the Inverter. e small-signal stability model of the inverter contains two parts: the droop controller and LC filter and coupling inductor. Moreover, the dynamic model of the droop controller consists of the power controller, voltage controller, and current controller.

e Model of the Power Controller.
e real and reactive power can be obtained when the instantaneous active and reactive power, calculated from the measured output voltage and output current in dq frame, passes through the low-pass filters. e small-signal stability model can be given by where Δ(·) is the state variable with respect to the time, Δ(·) is the state variable or the algebraic variable after linearization, P, Q are the real and reactive power, respectively, u od , u oq , i od , i oq are the output voltage and current in dq frame, respectively, and ω c is the cut-off frequency of the low-pass filters. e frequency and voltage can be set by the droop gain; and the small-signal stability model can be given by where ω is the frequency, m p , n q are the droop gains related to frequency and voltage amplitude, respectively, and u * od , u * oq are the output voltage reference in dq frame, respectively. e first dq frame of the inverter is set as the common frame; and all the variables from other inverter frames can be converted to the common frame. e small-signal stability model is given by where ω com is the frequency of the common frame taken by the first inverter and δ is the angle between an inverter frame and the common frame.

e Model of the Voltage Controller.
e voltage controller contains the PI regulator and the feed-forward terms; and the small-signal stability model of the voltage controller can be given by 2 Complexity where ϕ d , ϕ q are (u * od − u od )dt and (u * oq − u oq )dt, respectively, i * ld , i * lq are the output current reference in dq frame, respectively, K iu , K pu are the integral and proportional gains of the PI regulator in voltage controller, respectively, w n is the nominal frequency, F is the feedforward gain, and C f is the per-phase capacitance.

2.1.3.
e Model of the Current Controller. e current controller includes the PI regulator; and the small-signal stability model of the current controller can be given by where c d , c q are (i * ld − i ld )dt and (i * lq − i lq )dt, respectively, u * id , u * iq are the output inverter voltage reference in dq frame, respectively, i ld , i lq are the sampled filter current in dq frame, respectively, K ic , K pc are the integral and proportional gains in PI regulator of the current controller, respectively, and L f is the per-phase inductance.

e Model of the LC Filter and Coupling Inductance.
e LC filter and coupling inductance can remove the harmonic wave near the switch frequency; and the smallsignal stability model of the LC filter and coupling inductance are given by where u bd , u bq are the bus voltage in dq frame, respectively, u id , u iq are the inverter voltage in dq frame, respectively, ω 0 , I ld , I lq , U od , U oq , I od , I oq are the steady state value at the initial operating point, and R f , R c , L f , L c , C f are the resistance, inductance, and capacitance in the microgrid, respectively.

e Model of the Individual Inverter.
rough the transformation matrix, other inverter variables can transfer to the common frame. e linearization models of the transformation matrix are obtained: where Δi oDQi , Δu bDQi are the output current and bus voltage of the ith inverter converted to the common frame, respectively, and δ 0 is the angle between the inverter frame and the common frame. Combined by the linearization models (1)-(22), the linearization state equation and the linearization output equation of the ith inverter in the common frame are given by

e Model of the Network and Load.
e linearization models of the network and load are obtained: where Δi loadDQp ] T , and n, p are the numbers of the lines and loads, respectively.

e Model of the Microgrid.
To define the node voltage clearly, the virtual resistor R N is assumed between each node and ground. e small-signal stability model is given by where Δu bDQ � [Δu bDQ1 , Δu bDQ2 , . . . , Δu bDQm ] T , m is the number of nodes in the microgrid, M INV maps the inverter connection points onto network nodes, M LOAD maps the load connection point onto the network nodes, and M NET maps the connection lines onto nodes. Combined with (23)-(27), the small-signal stability model of the microgrid is given by where s is the number of inverters, and A mg is the state matrix of the microgrid model.

The Proposed Optimization Model
e dynamic system is stable in small-signal stability sense when all the real parts of the eigenvalues of A mg are negative. e largest real part of all eigenvalues is called the spectral abscissa η, which stands for the power system's stability margin. e value of the spectral abscissa is smaller, and the rate of decay is faster when the system suffers a small disturbance. erefore, the dynamic system is more stable in the small-signal stability sense. Based on the above smallsignal stability model of microgrid, the eigenvalues of A mg can be obtained. erefore, the spectral abscissa η can be given by where λ(A mg ) represents all eigenvalues of A mg . Re(λ) is the real part of all eigenvalues, and λ η is the eigenvalue with largest real part. Since the load model is contained in the small-signal stability model of the microgrid, the variation of the load can affect the spectral abscissa of the microgrid exceedingly. Also, load condition variation, such as random load fluctuation, frequently happens in the microgrid operation. Although the microgrid can work well at one load condition, also called operating conditions, it cannot guarantee the small-signal stability at other operating conditions. erefore, it is crucial to consider multiple operating conditions in the proposed optimization model to achieve better robustness.
e spectral abscissas of different operating condition can be represented as η i , which stands for the spectral abscissa of the ith operating conditions, i � [1, 2, . . . , l], and l is the number of the operating conditions.
In [7,9,15], through the eigenvalue analysis, the dominant eigenvalues near the imaginary axis are largely sensitive to m p , n q and the eigenvalues of medium frequency are largely sensitive to K pu , K iu , K pc , and K ic . e value of η, even the system stability, will vary with droop controller parameters. erefore, the critical droop controller parameters of each inverter are considered as the variables x to be optimized: x � m p1 , n q1 , K pu1 , K iu1 , K pc1 , K ic1 , . . . , m p2 , n q2 , K pu2 , K iu2 , K pc2 , K ic2 , . . . , m pN , n qN , K puN , where N is the total number of inverters. According to (29) and (30), the spectral abscissa is actually an implicit function of the controller parameters: η i � f m p , n q , K pu , K iu , K pc , K ic , . . . , m pN , n qN , K puN , K iuN , K pcN , K icN .
(31) e spectral abscissas of the worst operating condition, also called the largest spectral abscissa, can be represented as η max : According to (32), the objective function can be given to achieve maximum stability of the microgrid in the worst operating condition. is objective function also ensures the maximum stability of microgrid in other operating conditions: e constraints of controller parameters of the kth inverter can be listed as the constraint of optimization model: m pk ≤ m pk ≤ m pk , n qk ≤ n qk ≤ n qk , However, η max is not locally Lipschitz. Fortunately, η i has been proved to be locally Lipschitz and continuously differentiable on open dense subsets of R n [16], which means that it is continually differentiable almost everywhere. Also, its gradient can be obtained by calculating spectral abscissa sensitivities. Besides, the small-signal stability constraint with the spectral abscissa of multiple operating conditions can be added to the optimization model to ensure that the microgrid works well in the multiple operating conditions; and, due to the small-signal stability constraint in the optimization model, minimizing η is equivalent to minimize η max . 4 Complexity rough the above discussion, the complete optimization model can be given by where η is the upper limit of spectral abscissa; and the constraints of the optimization model can be listed as follows: n qk ≤ n qk ≤ n qk , where (·) and (·) are the upper and lower limits of controller parameters and k � [1, 2, . . . , N]. e spectral abscissa of multiple operating conditions moves toward the left as far as possible to ensure maximum stability through the optimization model. us, the microgrid with optimized controller parameters can operate stably in multiple operating conditions. According to (31), the spectral abscissa cannot be deduced as an analytical and nonlinear or linear expression of controller parameters; and the spectral abscissa with respect to the parameters has been proved to be nonsmooth [17]. e optimization model in (35)-(42) is a nonsmooth programming problem, which cannot be tackled by the methods for linear or nonlinear programming. It is hard to solve the nonsmooth programming in mathematical programming until the breakthrough in [13].

e SQP-GS Method.
According to the Clarke theorem, if any subgradient is zero, which belongs to the subdifferential of x 0 , there is a stationary point x 0 in this nonsmooth function [18]. e nonsmooth problem is hard to be solved until a GS method was proposed in 2005 [18], which can approximate the subgradients for minimizing a nonsmooth objective function. In [18], if the nonsmooth function or constraints are locally Lipschitz continuous, a series of the independent and random gradients of sample point near the nonsmooth point x can be evaluated as convex hull to express the subgradients of x. e traditional SQP method can solve a QP subproblem to obtain a descending direction. It only solves the smooth problem but fails for the nonsmooth problem. In 2012, the Sequential Quadratic Programming with Gradient Sampling (SQP-GS) method was proposed in [13]. It can converge globally to stationary points with probability of one when the objective function or constraints are locally Lipschitz and continuously differentiable on open dense subsets of R n . e SQP-GS method is developed to solve the optimization model in the following form: where the objective function f: R n ⟶ R and the inequality constraint functions g: R n ⟶ R m . e SQP-GS method can compute a search direction d k in the kth iteration by solving a QP dual subproblem: where x k is the vector of variables in the kth iteration, H k is the approximate Hessian of the Lagrangian optimization model of (43), ρ is a penalty parameter, and y f � (y f k,0 , y f k,1 , . . . , y f k,p ) and y g � (y g j k,0 , y g j k,1 , . . . , y g j k,p ) are the approximated minimum-norm elements of the subdifferential of f(x k ) and g(x k ), respectively, x f and x g are the vectors of the sample points of f(x k ) and g(x k ), respectively, and are sets of independent and identically distributed random sample points sampled from the closed ball B ϵ (x k ) centered at x k , and For reflecting the iteration process, infeasibility value can be defined as e following model reduction Δq k is also defined in terms of primal and dual infeasibility: (54) When the SQP-GS method converges to the optimal solution, these optimality certificates, such as σ k (x k ), Δq k , can be 0. e sample size p can be set to 0 when the SQP-GS method is solving a smooth function, in which case the SQP-GS method reduces to the general SQP method [14].

Gradient Evaluation of Spectral Abscissa Function.
e gradient of the spectral abscissa ∇f(x) is crucial for implementing the SQP-GS method; and the spectral abscissa with respect to the controller parameters, also called the spectral abscissa sensitivity, reflects the variation of the variable x on the impact of system stability. erefore, the gradient ∇f(x) of spectral abscissa can be represented as spectral abscissa sensitivity. e spectral abscissa sensitivity can be calculated by the numerical differentiation method, which is numerical spectral abscissa sensitivity. e numerical spectral abscissa sensitivity can be obtained by where η(A) is the spectral abscissa of the state matrix at the equilibrium point, ε is the small quantity of the variable variation, x is the variable, and η(A ε ) is the spectral abscissa of perturbed state matrix. As for different variables, the small quantity range is not defined clearly. It causes the accuracy of numerical spectral abscissa sensitivity to be uncertain. When calculating the spectral abscissa sensitivity for each variable, the eigenvalue of A ε needs to be recalculated more than once in each iteration. It decreases the efficiency of the optimization.
Hence, the closed-form eigenvalue sensitivity can be an alternative option [19]. It is a mathematical eigenvalue derived at an equilibrium point; and the eigenvalue sensitivity with respect to variable x can be expressed as where A mg is the state matrix of microgrid at the equilibrium point, ψ ⊤ i and ϕ i are the left and right eigenvectors of λ i at an equilibrium point, respectively, and (zA mg /zx) can be obtained on matrix calculus. e closed-form spectral abscissa sensitivity can be expressed by In (56) and (57), ψ ⊤ η , ϕ η , and (zA mg /zx) are required to calculate only once in each iteration. Besides, the closed-form spectral abscissa sensitivity is more accurate than numerical spectral abscissa sensitivity in each iteration. erefore, the closed-form spectral abscissa sensitivity is used to compute the gradient of the spectral abscissa.

Solving the Droop Controller Parameters Optimization
Problem by SQP-GS. Before employing the SQP-GS method to solve a QP subproblem to get a descending direction, the gradients of objective function and constraints in the optimization model with respect to variable x should be obtained. e constraint functions in (37)-(42) are smooth; and the gradient of these constraints can be obtained easily. erefore, the sample size p reduces to zero when the SQP-GS method solves the smooth function. Objective function (35) and the small-signal stability constraint (36) are nonsmooth. e state matrix A mg in (28) for each point in B g j ϵ,k is established. e most critical eigenvalue λ η and the corresponding left and right eigenvectors ψ ⊤ i and ϕ ⊤ i for each A mg can be obtained in each iteration. According to (56), the most critical eigenvalue sensitivity (zλ η /zx) with respect to the ith variable for each point is calculated. erefore, the gradient (zη/zx) with respect to the ith variable in (57) for each point can be obtained.
e initial values of droop controller parameters need to be set in the range of controller parameters limit. Otherwise, they may not converge. e Algorithm 1 is given as follows [14]. Figure 1 is used to demonstrate the effectiveness of the proposed SQP-GS method.

Microgrid System Parameters and Four Operating
Conditions. We design four operating conditions of the microgrid to verify the effectiveness of the SQP-GS method.
e test system in Figure 1 contains three inverters, two connection lines, and two loads. Load 1 is in bus 1, and the other one is in bus 3; and the microgrid system parameters are given in Table 1. Also, the initial condition of the microgrid system is given in Table 2.
To ensure that the droop controller parameters are robust after optimization, similar to multiple operating conditions in [6], we design different operating conditions. ey are shown in Table 3. R load of 25 Ω and 20 Ω per phase is equal to the resistive load of 5.6 kW and 7.3 kW [13]. e objection of this paper is to investigate the optimized 6 Complexity parameters in droop controllers that can enhance the smallsignal stability of all four operating conditions. e limit of controller parameters can be shown in Table 4. It is selected based on the trace of the eigenvalues varying with the controller parameters variation [9,10]. When the droop controller parameters vary in the limit of the variable, the critical eigenvalue will move within the stable region. Meanwhile, the stability of the microgrid operation can be guaranteed. According to the above rule, (1) Initialization: set k � 1, k max � 200; Set sample size p > 0, sample radius ε > 0, constraint violation tolerance θ > 0, penalty parameter ρ > 0, line search constant ϖ ∈ (0, 1), backtracking constant c ∈ (0, 1), sampling radius reduction factor μ ε ∈ (0, 1), constraint violation tolerance reduction factor μ θ ∈ (0, 1), penalty parameter reduction factor μ ρ ∈ (0, 1), infeasibility tolerance ] in > 0, stationarity tolerance parameter ] s > 0. (2) Termination conditions check: while k < k max , if max(△q k ) < ] s or ε < ] s , and σ k < ] in , output solution and stop. method [20]. (6) Parameter update: if Δq k > ] s ε 2 , go to Step 7. Otherwise, if σ k ≤ θ, set θ←μ θ θ; if σ k > θ, set ρ←μ ρ ρ. And set ε←μ ε ε, β k ←0 and go to step 8. (7) Line search: set β k as the largest value in the sequence 1, c, c 2 , . . . such that x k+1 ←x k + β k d k satisfies:  Parameter Value    we can obtain the value domain of the variables from the root locus of the critical eigenvalues with the variable variation.

Optimization Method Verification.
In the last subsection, the four operating conditions have been designed; and we implement the SQP-GS method in MATLAB to optimize the droop controller parameters in microgrid. From Table 5, we can see the comparison of spectral abscissas between the unoptimized and the optimized values in the four operating conditions. In Case 1 and Case 3, the unoptimized spectral abscissa is 72.3025 and 72.3112, respectively. It stands that the microgrid system is susceptible to encounter oscillation when random load fluctuation occurs. After the optimization, the optimized spectral abscissa is −16.1651 and −16.1212, respectively. It proves that the small-signal stability of the microgrid in four operating conditions is enhanced through SQP-GS method optimization. erefore, it demonstrated that the SQP-GS method is effective for optimizing the droop controller parameters. e original parameters and the optimized parameters of the controller are given in Tables 6 and 7. As for the original parameters, the droop controller parameter is designed by the traditional method [15]. Also, the voltage and current controller parameters are selected based on the small-signal stability analysis [9,10]. For the closed-form spectral abscissa sensitivities of K pv , K iv , K pc , and K ic are fairly small, these controller parameters stay the same nearly in the iteration process.
From Figure 2, there are two eigenvalues with the positive real part in the original system of Case 1. One positive eigenvalue is 72.3025, and the other one is 15.9420. It means that the original system is unstable. After optimizing the controller parameters by the SQP-GS method, all the optimized eigenvalues lie on the left plane. e spectral abscissa of Case 1 is −16.1651. It proves that the stability of the system is improved after the optimization.
In Figure 3, the spectral abscissa of Case 1 decreases in the iteration process with the SQP-GS method. After 83 iterations, the spectral abscissa of Case 1 varies from 72.3025 to −16.1651. It means that the microgrid system with optimized droop controller parameters is stable in the smallsignal stability sense.
In Figure 4, the upper limit of the spectral abscissa is varied from −15 to −16.061 in the iteration process. It means that the spectral abscissas of all the cases are less than −16.061. e stability margin of the system is improved. Figure 5 depicts the sampling radius variation, which is one of the convergence indicators in SQP-GS method. After 83 iterations, the sampling radius decreases from 10e − 1 to 10e − 7; and it proves that the SQP-GS method arrives at convergence.

Comparison with GA Method.
e SQP-GS and GA methods are called nine times, respectively, for optimizing the droop controller parameters; and the mean value of η, spectral abscissa in Case 4, and the computing time are listed in Table 8. As the table shows, the mean value of η optimized by the GA method is −15.31. Besides, η optimized by the SQP-GS method is −16.06. In the term of the spectral abscissa, the best value optimized by the GA method is −16.05, while the value optimized by the SQP-GS method is −16.12. As for the computing time, the longest time for the   GA method to converge is 0.709 hours, while the shortest time is 0.116 hours. However, the max computing time of the SQP-GS method is 0.102 hours. ese results demonstrate that the optimality and efficiency of the SQP-GS method are better than those of the GA method. Also, Table 8 shows that the SQP-GS can be globally convergent to stationary points with probability of one.

Conclusion
is paper proposes the SQP-GS method, which is a mathematical programming method to optimize the spectral abscissa of the microgrid. Combined with the proposed optimization model and SQP-GS method, the microgrid system is better robust in different operating conditions. Moreover, the SQP-GS method's effectiveness and efficiency are better than those of the GA method.
Data Availability e microgrid system parameters and the initial condition of microgrid system which support the findings of this study are available in [15].

Conflicts of Interest
e authors declare that they have no conflicts of interest.