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

Authors: Peijie Li; Ziyi Yang; Shuchen Huang; Jun Zhang
Journal: Complexity
 (2021)

---

## Abstract

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. The Sequential Quadratic Programming with Gradient Sampling (SQP-GS) is implemented to optimize the droop controller parameters for solving the nonsmooth problem. The 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.


---

## Body

## 1. 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. The 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]. The small-signal stability of the microgrid in islanded mode is strongly affected by the droop controller parameters. Therefore, it is essential to optimize the controller parameters to enhance the small-signal stability of an islanded microgrid.Particle swarm optimization (PSO) method is often used to optimize the controller parameters of the microgrid [6–9]. 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. The 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. The 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. The work in [9] proposed an optimization model with a nonlinear time-domain 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 grid-connected 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]. The solutions derived by the heuristic algorithms are not deterministic and reproducible. These 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–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. The 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. Therefore, 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].The 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)
The 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.The remainder of this paper is organized as follows: Section2 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.


## 2. The Small-Signal Stability Model of the Microgrid

The small-signal stability model of the microgrid can be divided into three parts: inverter, network, and load. These small-signal stability models are deduced indq reference frame [10, 15].


### 2.1. The Model of the Inverter

The 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.


#### 2.1.1. The Model of the Power Controller

The real and reactive power can be obtained when the instantaneous active and reactive power, calculated from the measured output voltage and output current indq frame, passes through the low-pass filters. The small-signal stability model can be given by(1)ΔP˙=−ωc△P+ωcIodΔuod+IoqΔuoq+UodΔiod+UoqΔioq,(2)ΔQ˙=−ωcΔQ+ωcIoqΔuod−IodΔuoq−UoqΔiod+UodΔioq,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, uod,uoq,iod,ioq are the output voltage and current in dq frame, respectively, and ωc is the cut-off frequency of the low-pass filters.The frequency and voltage can be set by the droop gain; and the small-signal stability model can be given by(3)Δω=−mpΔP,(4)Δuod∗=−nqΔQ,(5)Δuoq∗=0,where ω is the frequency, mp,nq are the droop gains related to frequency and voltage amplitude, respectively, and uod∗,uoq∗ are the output voltage reference in dq frame, respectively.The firstdq 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. The small-signal stability model is given by(6)Δδ˙=Δω−Δωcom=−mpΔP−Δωcom,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.


#### 2.1.2. The Model of the Voltage Controller

The 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(7)Δϕ˙d=Δuod∗−Δuod,(8)Δϕ˙q=Δuoq∗−Δuoq,(9)Δild∗=KiuΔϕd+KpuΔuod∗−Δuod−ωnCfΔuoq+FΔiod,(10)Δilq∗=KiuΔϕq+KpuΔuoq∗−Δuoq+ωnCfΔuod+FΔioq,where ϕd,ϕq are ∫uod∗−uoddt and ∫uoq∗−uoqdt, respectively, ild∗,ilq∗ are the output current reference in dq frame, respectively, Kiu,Kpu are the integral and proportional gains of the PI regulator in voltage controller, respectively, wn is the nominal frequency, F is the feed-forward gain, and Cf is the per-phase capacitance.


#### 2.1.3. The Model of the Current Controller

The current controller includes the PI regulator; and the small-signal stability model of the current controller can be given by(11)Δγ˙d=Δild∗−Δild,(12)Δγ˙q=Δilq∗−Δilq,(13)Δuid∗=KicΔγd+KpcΔild∗−Δild−ωnLfΔilq,(14)Δuiq∗=KicΔγq+KpcΔilq∗−Δilq+ωnLfΔild,where γd,γq are ∫ild∗−ilddt and ∫ilq∗−ilqdt, respectively, uid∗,uiq∗ are the output inverter voltage reference in dq frame, respectively, ild,ilq are the sampled filter current in dq frame, respectively, Kic,Kpc are the integral and proportional gains in PI regulator of the current controller, respectively, and Lf is the per-phase inductance.


#### 2.1.4. The Model of the LC Filter and Coupling Inductance

The LC filter and coupling inductance can remove the harmonic wave near the switch frequency; and the small-signal stability model of the LC filter and coupling inductance are given by(15)Δi˙ld=−RfLfΔild+ωoΔiiq−1LfΔuod+1LfΔuid+IlqΔω,(16)Δi˙lq=−RfLfΔilq−ωoΔiid−1LfΔuoq+1LfΔuiq−IldΔω,(17)Δu˙od=1CfΔild+ωoΔuoq−1CfΔiod+UoqΔω,(18)Δu˙oq=1CfΔilq−ωoΔuod−1CfΔioq−UodΔω,(19)Δi˙od=1LcΔuod−RcLcΔiod+ωoΔioq−1LcΔubd+IoqΔω,(20)Δi˙oq=1LcΔuoq−RcLcΔioq−ωoΔiod−1LcΔubq−IodΔω,where ubd,ubq are the bus voltage in dq frame, respectively, uid,uiq are the inverter voltage in dq frame, respectively, ω0,Ild,Ilq,Uod,Uoq,Iod,Ioq are the steady state value at the initial operating point, and Rf,Rc,Lf,Lc,Cf are the resistance, inductance, and capacitance in the microgrid, respectively.


#### 2.1.5. The Model of the Individual Inverter

Through the transformation matrix, other inverter variables can transfer to the common frame. The linearization models of the transformation matrix are obtained:(21)ΔioDQi=cosδo−sinδosinδocosδoΔiodqi+−Iodsinδo−IoqcosδoIodcosδo−IoqsinδoΔδi,(22)Δubdqi=cosδosinδo−sinδocosδoΔubDQi+−UbDsinδo+UbQcosδo−UbDcosδo−UbQsinδoΔδi,where ΔioDQi,ΔubDQi 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(23)Δx˙invi=AINViΔxinvi+BINViΔubDQi+BiωcomΔωcom,(24)ΔωiΔioDQi=CINVωiCINVciΔxinvi,where Δxinvi=ΔδiΔPiΔQiΔϕdqiΔγdqiΔildqiΔuodqiΔiodqiT.


### 2.2. The Model of the Network and Load

The linearization models of the network and load are obtained:(25)Δi˙lineDQ=ANETΔilineDQ+B1NETΔubDQ+B2NETΔω,(26)Δi˙loadDQ=ALOADΔiloadDQ+B1LOADΔubDQ+B2LOADΔω,where ΔilineDQ=ΔilineDQ1,ΔilineDQ2,…,ΔilineDQnT, ΔiloadDQ=ΔiloadDQ1,ΔiloadDQ2,…,ΔiloadDQpT, and n,p are the numbers of the lines and loads, respectively.


### 2.3. The Model of the Microgrid

To define the node voltage clearly, the virtual resistorRN is assumed between each node and ground. The small-signal stability model is given by(27)ΔubDQ=RNMINVΔioDQ+MLOADΔiloadDQ+MNETΔilineDQ,where ΔubDQ=ΔubDQ1,ΔubDQ2,…,ΔubDQmT, m is the number of nodes in the microgrid, MINV maps the inverter connection points onto network nodes, MLOAD maps the load connection point onto the network nodes, and MNET maps the connection lines onto nodes.Combined with (23)-(27), the small-signal stability model of the microgrid is given by(28)Δx˙INVΔi˙lineDQΔi˙loadDQ=AmgΔxINVΔilineDQΔiloadDQ,where ΔxINV=xINV1,xINV2,…,xINVsT, s is the number of inverters, and Amg is the state matrix of the microgrid model.


## 2.1. The Model of the Inverter

The 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.


### 2.1.1. The Model of the Power Controller

The real and reactive power can be obtained when the instantaneous active and reactive power, calculated from the measured output voltage and output current indq frame, passes through the low-pass filters. The small-signal stability model can be given by(1)ΔP˙=−ωc△P+ωcIodΔuod+IoqΔuoq+UodΔiod+UoqΔioq,(2)ΔQ˙=−ωcΔQ+ωcIoqΔuod−IodΔuoq−UoqΔiod+UodΔioq,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, uod,uoq,iod,ioq are the output voltage and current in dq frame, respectively, and ωc is the cut-off frequency of the low-pass filters.The frequency and voltage can be set by the droop gain; and the small-signal stability model can be given by(3)Δω=−mpΔP,(4)Δuod∗=−nqΔQ,(5)Δuoq∗=0,where ω is the frequency, mp,nq are the droop gains related to frequency and voltage amplitude, respectively, and uod∗,uoq∗ are the output voltage reference in dq frame, respectively.The firstdq 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. The small-signal stability model is given by(6)Δδ˙=Δω−Δωcom=−mpΔP−Δωcom,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.


### 2.1.2. The Model of the Voltage Controller

The 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(7)Δϕ˙d=Δuod∗−Δuod,(8)Δϕ˙q=Δuoq∗−Δuoq,(9)Δild∗=KiuΔϕd+KpuΔuod∗−Δuod−ωnCfΔuoq+FΔiod,(10)Δilq∗=KiuΔϕq+KpuΔuoq∗−Δuoq+ωnCfΔuod+FΔioq,where ϕd,ϕq are ∫uod∗−uoddt and ∫uoq∗−uoqdt, respectively, ild∗,ilq∗ are the output current reference in dq frame, respectively, Kiu,Kpu are the integral and proportional gains of the PI regulator in voltage controller, respectively, wn is the nominal frequency, F is the feed-forward gain, and Cf is the per-phase capacitance.


### 2.1.3. The Model of the Current Controller

The current controller includes the PI regulator; and the small-signal stability model of the current controller can be given by(11)Δγ˙d=Δild∗−Δild,(12)Δγ˙q=Δilq∗−Δilq,(13)Δuid∗=KicΔγd+KpcΔild∗−Δild−ωnLfΔilq,(14)Δuiq∗=KicΔγq+KpcΔilq∗−Δilq+ωnLfΔild,where γd,γq are ∫ild∗−ilddt and ∫ilq∗−ilqdt, respectively, uid∗,uiq∗ are the output inverter voltage reference in dq frame, respectively, ild,ilq are the sampled filter current in dq frame, respectively, Kic,Kpc are the integral and proportional gains in PI regulator of the current controller, respectively, and Lf is the per-phase inductance.


### 2.1.4. The Model of the LC Filter and Coupling Inductance

The LC filter and coupling inductance can remove the harmonic wave near the switch frequency; and the small-signal stability model of the LC filter and coupling inductance are given by(15)Δi˙ld=−RfLfΔild+ωoΔiiq−1LfΔuod+1LfΔuid+IlqΔω,(16)Δi˙lq=−RfLfΔilq−ωoΔiid−1LfΔuoq+1LfΔuiq−IldΔω,(17)Δu˙od=1CfΔild+ωoΔuoq−1CfΔiod+UoqΔω,(18)Δu˙oq=1CfΔilq−ωoΔuod−1CfΔioq−UodΔω,(19)Δi˙od=1LcΔuod−RcLcΔiod+ωoΔioq−1LcΔubd+IoqΔω,(20)Δi˙oq=1LcΔuoq−RcLcΔioq−ωoΔiod−1LcΔubq−IodΔω,where ubd,ubq are the bus voltage in dq frame, respectively, uid,uiq are the inverter voltage in dq frame, respectively, ω0,Ild,Ilq,Uod,Uoq,Iod,Ioq are the steady state value at the initial operating point, and Rf,Rc,Lf,Lc,Cf are the resistance, inductance, and capacitance in the microgrid, respectively.


### 2.1.5. The Model of the Individual Inverter

Through the transformation matrix, other inverter variables can transfer to the common frame. The linearization models of the transformation matrix are obtained:(21)ΔioDQi=cosδo−sinδosinδocosδoΔiodqi+−Iodsinδo−IoqcosδoIodcosδo−IoqsinδoΔδi,(22)Δubdqi=cosδosinδo−sinδocosδoΔubDQi+−UbDsinδo+UbQcosδo−UbDcosδo−UbQsinδoΔδi,where ΔioDQi,ΔubDQi 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(23)Δx˙invi=AINViΔxinvi+BINViΔubDQi+BiωcomΔωcom,(24)ΔωiΔioDQi=CINVωiCINVciΔxinvi,where Δxinvi=ΔδiΔPiΔQiΔϕdqiΔγdqiΔildqiΔuodqiΔiodqiT.


## 2.1.1. The Model of the Power Controller

The real and reactive power can be obtained when the instantaneous active and reactive power, calculated from the measured output voltage and output current indq frame, passes through the low-pass filters. The small-signal stability model can be given by(1)ΔP˙=−ωc△P+ωcIodΔuod+IoqΔuoq+UodΔiod+UoqΔioq,(2)ΔQ˙=−ωcΔQ+ωcIoqΔuod−IodΔuoq−UoqΔiod+UodΔioq,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, uod,uoq,iod,ioq are the output voltage and current in dq frame, respectively, and ωc is the cut-off frequency of the low-pass filters.The frequency and voltage can be set by the droop gain; and the small-signal stability model can be given by(3)Δω=−mpΔP,(4)Δuod∗=−nqΔQ,(5)Δuoq∗=0,where ω is the frequency, mp,nq are the droop gains related to frequency and voltage amplitude, respectively, and uod∗,uoq∗ are the output voltage reference in dq frame, respectively.The firstdq 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. The small-signal stability model is given by(6)Δδ˙=Δω−Δωcom=−mpΔP−Δωcom,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.


## 2.1.2. The Model of the Voltage Controller

The 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(7)Δϕ˙d=Δuod∗−Δuod,(8)Δϕ˙q=Δuoq∗−Δuoq,(9)Δild∗=KiuΔϕd+KpuΔuod∗−Δuod−ωnCfΔuoq+FΔiod,(10)Δilq∗=KiuΔϕq+KpuΔuoq∗−Δuoq+ωnCfΔuod+FΔioq,where ϕd,ϕq are ∫uod∗−uoddt and ∫uoq∗−uoqdt, respectively, ild∗,ilq∗ are the output current reference in dq frame, respectively, Kiu,Kpu are the integral and proportional gains of the PI regulator in voltage controller, respectively, wn is the nominal frequency, F is the feed-forward gain, and Cf is the per-phase capacitance.


## 2.1.3. The Model of the Current Controller

The current controller includes the PI regulator; and the small-signal stability model of the current controller can be given by(11)Δγ˙d=Δild∗−Δild,(12)Δγ˙q=Δilq∗−Δilq,(13)Δuid∗=KicΔγd+KpcΔild∗−Δild−ωnLfΔilq,(14)Δuiq∗=KicΔγq+KpcΔilq∗−Δilq+ωnLfΔild,where γd,γq are ∫ild∗−ilddt and ∫ilq∗−ilqdt, respectively, uid∗,uiq∗ are the output inverter voltage reference in dq frame, respectively, ild,ilq are the sampled filter current in dq frame, respectively, Kic,Kpc are the integral and proportional gains in PI regulator of the current controller, respectively, and Lf is the per-phase inductance.


## 2.1.4. The Model of the LC Filter and Coupling Inductance

The LC filter and coupling inductance can remove the harmonic wave near the switch frequency; and the small-signal stability model of the LC filter and coupling inductance are given by(15)Δi˙ld=−RfLfΔild+ωoΔiiq−1LfΔuod+1LfΔuid+IlqΔω,(16)Δi˙lq=−RfLfΔilq−ωoΔiid−1LfΔuoq+1LfΔuiq−IldΔω,(17)Δu˙od=1CfΔild+ωoΔuoq−1CfΔiod+UoqΔω,(18)Δu˙oq=1CfΔilq−ωoΔuod−1CfΔioq−UodΔω,(19)Δi˙od=1LcΔuod−RcLcΔiod+ωoΔioq−1LcΔubd+IoqΔω,(20)Δi˙oq=1LcΔuoq−RcLcΔioq−ωoΔiod−1LcΔubq−IodΔω,where ubd,ubq are the bus voltage in dq frame, respectively, uid,uiq are the inverter voltage in dq frame, respectively, ω0,Ild,Ilq,Uod,Uoq,Iod,Ioq are the steady state value at the initial operating point, and Rf,Rc,Lf,Lc,Cf are the resistance, inductance, and capacitance in the microgrid, respectively.


## 2.1.5. The Model of the Individual Inverter

Through the transformation matrix, other inverter variables can transfer to the common frame. The linearization models of the transformation matrix are obtained:(21)ΔioDQi=cosδo−sinδosinδocosδoΔiodqi+−Iodsinδo−IoqcosδoIodcosδo−IoqsinδoΔδi,(22)Δubdqi=cosδosinδo−sinδocosδoΔubDQi+−UbDsinδo+UbQcosδo−UbDcosδo−UbQsinδoΔδi,where ΔioDQi,ΔubDQi 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(23)Δx˙invi=AINViΔxinvi+BINViΔubDQi+BiωcomΔωcom,(24)ΔωiΔioDQi=CINVωiCINVciΔxinvi,where Δxinvi=ΔδiΔPiΔQiΔϕdqiΔγdqiΔildqiΔuodqiΔiodqiT.


## 2.2. The Model of the Network and Load

The linearization models of the network and load are obtained:(25)Δi˙lineDQ=ANETΔilineDQ+B1NETΔubDQ+B2NETΔω,(26)Δi˙loadDQ=ALOADΔiloadDQ+B1LOADΔubDQ+B2LOADΔω,where ΔilineDQ=ΔilineDQ1,ΔilineDQ2,…,ΔilineDQnT, ΔiloadDQ=ΔiloadDQ1,ΔiloadDQ2,…,ΔiloadDQpT, and n,p are the numbers of the lines and loads, respectively.


## 2.3. The Model of the Microgrid

To define the node voltage clearly, the virtual resistorRN is assumed between each node and ground. The small-signal stability model is given by(27)ΔubDQ=RNMINVΔioDQ+MLOADΔiloadDQ+MNETΔilineDQ,where ΔubDQ=ΔubDQ1,ΔubDQ2,…,ΔubDQmT, m is the number of nodes in the microgrid, MINV maps the inverter connection points onto network nodes, MLOAD maps the load connection point onto the network nodes, and MNET maps the connection lines onto nodes.Combined with (23)-(27), the small-signal stability model of the microgrid is given by(28)Δx˙INVΔi˙lineDQΔi˙loadDQ=AmgΔxINVΔilineDQΔiloadDQ,where ΔxINV=xINV1,xINV2,…,xINVsT, s is the number of inverters, and Amg is the state matrix of the microgrid model.


## 3. The Proposed Optimization Model

The dynamic system is stable in small-signal stability sense when all the real parts of the eigenvalues ofAmg are negative. The largest real part of all eigenvalues is called the spectral abscissa η, which stands for the power system’s stability margin. The value of the spectral abscissa is smaller, and the rate of decay is faster when the system suffers a small disturbance. Therefore, the dynamic system is more stable in the small-signal stability sense. Based on the above small-signal stability model of microgrid, the eigenvalues of Amg can be obtained. Therefore, the spectral abscissa η can be given by(29)ηAmg=maxReλ:λ∈λAmg=Reλη,where λAmg represents all eigenvalues of Amg. 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. Therefore, it is crucial to consider multiple operating conditions in the proposed optimization model to achieve better robustness.The 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 mp,nq and the eigenvalues of medium frequency are largely sensitive to Kpu, Kiu, Kpc, and Kic. The value of η, even the system stability, will vary with droop controller parameters.Therefore, the critical droop controller parameters of each inverter are considered as the variablesx to be optimized:(30)x=mp1,nq1,Kpu1,Kiu1,Kpc1,Kic1,…,mp2,nq2,Kpu2,Kiu2,Kpc2,Kic2,…,mpN,nqN,KpuN,KiuN,KpcN,KicN,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:(31)ηi=fmp,nq,Kpu,Kiu,Kpc,Kic,…,mpN,nqN,KpuN,KiuN,KpcN,KicN.The spectral abscissas of the worst operating condition, also called the largest spectral abscissa, can be represented asηmax:(32)ηmax=maxη1,η2,…,ηl.According to (32), the objective function can be given to achieve maximum stability of the microgrid in the worst operating condition. This objective function also ensures the maximum stability of microgrid in other operating conditions:(33)minηmax.The constraints of controller parameters of thekth inverter can be listed as the constraint of optimization model:(34)m¯pk≤mpk≤m¯pk,n¯qk≤nqk≤n¯qk,K¯puk≤Kpuk≤K¯puk,K¯iuk≤Kiuk≤K¯iuk,K¯pck≤Kpck≤K¯pck,K¯ick≤Kick≤K¯ick.However,ηmax is not locally Lipschitz. Fortunately, ηi has been proved to be locally Lipschitz and continuously differentiable on open dense subsets of ℝ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.Through the above discussion, the complete optimization model can be given by(35)minη¯,where η¯ is the upper limit of spectral abscissa; and the constraints of the optimization model can be listed as follows:(36)ηi≤η¯,(37)m¯pk≤mpk≤m¯pk,(38)n¯qk≤nqk≤n¯qk,(39)K¯puk≤Kpuk≤K¯puk,(40)K¯iuk≤Kiuk≤K¯iuk,(41)K¯pck≤Kpck≤K¯pck,(42)K¯ick≤Kick≤K¯ick,where ⋅¯ and ⋅¯ are the upper and lower limits of controller parameters and k=1,2,…,N.The spectral abscissa of multiple operating conditions moves toward the left as far as possible to ensure maximum stability through the optimization model. Thus, 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]. The 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].


## 4. The SQP-GS Method for Solving the Optimization Problem

### 4.1. The SQP-GS Method

According to the Clarke theorem, if any subgradient is zero, which belongs to the subdifferential ofx0, there is a stationary point x0 in this nonsmooth function [18].The 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.The 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 ℝn.The SQP-GS method is developed to solve the optimization model in the following form:(43)minxfxs.t.gx≤0,where the objective function f:ℝn⟶ℝ and the inequality constraint functions g:ℝn⟶ℝm.The SQP-GS method can compute a search directiondk in the kth iteration by solving a QP dual subproblem:(44)mindk,yf,yg−fxk−gxk⊤yfyg+dk⊤,yf⊤,yg⊤Hkdkyfyg,(45)s.t.Hkdk+∇fxfTyf+∇gxgTyg=0,∀xf∈ℬϵ,kf,∀xg∈ℬϵ,kg,(46)eTyf=ρ,(47)eTyg=1,(48)yf≥0,(49)yg≥0,where xk is the vector of variables in the kth iteration, Hk is the approximate Hessian of the Lagrangian optimization model of (43), ρ is a penalty parameter, and yf=yk,0f,yk,1f,…,yk,pf and yg=yk,0gj,yk,1gj,…,yk,pgj are the approximated minimum-norm elements of the subdifferential of fxk and gxk, respectively, xf and xg are the vectors of the sample points of fxk and gxk, respectively, and(50)ℬϵ,kf≔xk,0f,xk,1f,…,xk,pf,wherexk,0f≔xk,(51)ℬϵ,kgj≔xk,0gj,xk,1gj,…,xk,pgj,wherexk,0gj≔xk,are sets of independent and identically distributed random sample points sampled from the closed ball Bϵxk centered at xk, and(52)Bϵxk=x|x−xk≤ε.For reflecting the iteration process, infeasibility value can be defined as(53)σkxk≔maxgxk,0.The following model reductionΔqk is also defined in terms of primal and dual infeasibility:(54)Δqk=ρfxk+σkxk−ρmaxxk∈ℬϵ,kffxk+∇fxk⊤dk−∑j=1mmaxxk∈ℬϵ,kg,kmaxgjxk+∇gjxk⊤dk,0−12dk⊤Hkdk.When the SQP-GS method converges to the optimal solution, these optimality certificates, such asσkxk, Δqk, can be 0. The 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].


### 4.2. Gradient Evaluation of Spectral Abscissa Function

The gradient of the spectral abscissa∇fx 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. Therefore, the gradient ∇fx of spectral abscissa can be represented as spectral abscissa sensitivity. The spectral abscissa sensitivity can be calculated by the numerical differentiation method, which is numerical spectral abscissa sensitivity.The numerical spectral abscissa sensitivity can be obtained by(55)∂η∂x≈ηAε−ηAε,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 ofAε 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(56)∂λi∂x=ψi⊤∂Amg/∂xϕiψi⊤ϕi,where Amg 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 ∂Amg/∂x can be obtained on matrix calculus.The closed-form spectral abscissa sensitivity can be expressed by(57)∂η∂x=Re∂λη∂x.In (56) and (57), ψη⊤, ϕη, and ∂Amg/∂x 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. Therefore, the closed-form spectral abscissa sensitivity is used to compute the gradient of the spectral abscissa.


### 4.3. 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 variablex should be obtained. The constraint functions in (37)-(42) are smooth; and the gradient of these constraints can be obtained easily. Therefore, 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. The state matrix Amg in (28) for each point in ℬϵ,kgj is established. The most critical eigenvalue λη and the corresponding left and right eigenvectors ψi⊤ and ϕi⊤ for each Amg can be obtained in each iteration. According to (56), the most critical eigenvalue sensitivity ∂λη/∂x with respect to the ith variable for each point is calculated. Therefore, the gradient ∂η/∂x with respect to the ith variable in (57) for each point can be obtained. The initial values of droop controller parameters need to be set in the range of controller parameters limit. Otherwise, they may not converge. The Algorithm 1 is given as follows [14].Algorithm 1: SQP-GS method.
(1)
Initialization: setk=1,kmax=200; Set sample size p>0, sample radius ε>0, constraint violation tolerance θ>0, penalty parameter ρ>0, line search constant ϖ∈0,1, backtracking constant γ∈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: whilek<kmax, if max△qk<νs or ε<νs, and σk<νin, output solution and stop.(3)
Gradient sampling: generate sample pointℬε,kf in (50) and obtain gradient of multiple operating condition ∂ηi/∂x from (56) and (57), Generate sample point ℬε,kg from (51), and get the gradient of smooth constraint in (37)-(42).(4)
Search direction calculation: setHk≻0, solve (45)-(49) to get dk,yf,yg.(5)
L-BFGS update: updateHk for next iteration according to Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method [20].(6)
Parameter update: ifΔqk>ν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,γ,γ2,… such that xk+1←xk+βkdk satisfies: ρfxk+1+σxk+1≤ρfxk+σxk−ϖβkΔqk(8)
Iteration increment: setk←k+1 and go to step 2.(9)
end do


## 4.1. The SQP-GS Method

According to the Clarke theorem, if any subgradient is zero, which belongs to the subdifferential ofx0, there is a stationary point x0 in this nonsmooth function [18].The 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.The 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 ℝn.The SQP-GS method is developed to solve the optimization model in the following form:(43)minxfxs.t.gx≤0,where the objective function f:ℝn⟶ℝ and the inequality constraint functions g:ℝn⟶ℝm.The SQP-GS method can compute a search directiondk in the kth iteration by solving a QP dual subproblem:(44)mindk,yf,yg−fxk−gxk⊤yfyg+dk⊤,yf⊤,yg⊤Hkdkyfyg,(45)s.t.Hkdk+∇fxfTyf+∇gxgTyg=0,∀xf∈ℬϵ,kf,∀xg∈ℬϵ,kg,(46)eTyf=ρ,(47)eTyg=1,(48)yf≥0,(49)yg≥0,where xk is the vector of variables in the kth iteration, Hk is the approximate Hessian of the Lagrangian optimization model of (43), ρ is a penalty parameter, and yf=yk,0f,yk,1f,…,yk,pf and yg=yk,0gj,yk,1gj,…,yk,pgj are the approximated minimum-norm elements of the subdifferential of fxk and gxk, respectively, xf and xg are the vectors of the sample points of fxk and gxk, respectively, and(50)ℬϵ,kf≔xk,0f,xk,1f,…,xk,pf,wherexk,0f≔xk,(51)ℬϵ,kgj≔xk,0gj,xk,1gj,…,xk,pgj,wherexk,0gj≔xk,are sets of independent and identically distributed random sample points sampled from the closed ball Bϵxk centered at xk, and(52)Bϵxk=x|x−xk≤ε.For reflecting the iteration process, infeasibility value can be defined as(53)σkxk≔maxgxk,0.The following model reductionΔqk is also defined in terms of primal and dual infeasibility:(54)Δqk=ρfxk+σkxk−ρmaxxk∈ℬϵ,kffxk+∇fxk⊤dk−∑j=1mmaxxk∈ℬϵ,kg,kmaxgjxk+∇gjxk⊤dk,0−12dk⊤Hkdk.When the SQP-GS method converges to the optimal solution, these optimality certificates, such asσkxk, Δqk, can be 0. The 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].


## 4.2. Gradient Evaluation of Spectral Abscissa Function

The gradient of the spectral abscissa∇fx 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. Therefore, the gradient ∇fx of spectral abscissa can be represented as spectral abscissa sensitivity. The spectral abscissa sensitivity can be calculated by the numerical differentiation method, which is numerical spectral abscissa sensitivity.The numerical spectral abscissa sensitivity can be obtained by(55)∂η∂x≈ηAε−ηAε,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 ofAε 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(56)∂λi∂x=ψi⊤∂Amg/∂xϕiψi⊤ϕi,where Amg 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 ∂Amg/∂x can be obtained on matrix calculus.The closed-form spectral abscissa sensitivity can be expressed by(57)∂η∂x=Re∂λη∂x.In (56) and (57), ψη⊤, ϕη, and ∂Amg/∂x 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. Therefore, the closed-form spectral abscissa sensitivity is used to compute the gradient of the spectral abscissa.


## 4.3. 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 variablex should be obtained. The constraint functions in (37)-(42) are smooth; and the gradient of these constraints can be obtained easily. Therefore, 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. The state matrix Amg in (28) for each point in ℬϵ,kgj is established. The most critical eigenvalue λη and the corresponding left and right eigenvectors ψi⊤ and ϕi⊤ for each Amg can be obtained in each iteration. According to (56), the most critical eigenvalue sensitivity ∂λη/∂x with respect to the ith variable for each point is calculated. Therefore, the gradient ∂η/∂x with respect to the ith variable in (57) for each point can be obtained. The initial values of droop controller parameters need to be set in the range of controller parameters limit. Otherwise, they may not converge. The Algorithm 1 is given as follows [14].Algorithm 1: SQP-GS method.
(1)
Initialization: setk=1,kmax=200; Set sample size p>0, sample radius ε>0, constraint violation tolerance θ>0, penalty parameter ρ>0, line search constant ϖ∈0,1, backtracking constant γ∈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: whilek<kmax, if max△qk<νs or ε<νs, and σk<νin, output solution and stop.(3)
Gradient sampling: generate sample pointℬε,kf in (50) and obtain gradient of multiple operating condition ∂ηi/∂x from (56) and (57), Generate sample point ℬε,kg from (51), and get the gradient of smooth constraint in (37)-(42).(4)
Search direction calculation: setHk≻0, solve (45)-(49) to get dk,yf,yg.(5)
L-BFGS update: updateHk for next iteration according to Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method [20].(6)
Parameter update: ifΔqk>ν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,γ,γ2,… such that xk+1←xk+βkdk satisfies: ρfxk+1+σxk+1≤ρfxk+σxk−ϖβkΔqk(8)
Iteration increment: setk←k+1 and go to step 2.(9)
end do


## 5. Test Variation

The following test system in Figure1 is used to demonstrate the effectiveness of the proposed SQP-GS method. The simulations are carried out on a computer with a four-core 3.70 GHz processor and 16 GB RAM. Also, the SQP-GS method is implemented in MATLAB R2012a. As for the SQP-GS method, those parameters are selected as follows: ρ=5e−4, μρ=0.5, ε=0.1, με=0.1, θ=0.1, and μθ=0.8.Figure 1
Microgrid structure.


### 5.1. Microgrid System Parameters and Four Operating Conditions

We design four operating conditions of the microgrid to verify the effectiveness of the SQP-GS method. The test system in Figure1 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.Table 1
The microgrid system parameters.
ParameterValueLf1.35e−3HCf50e−6FRf0.1ΩRline10.23ΩRline20.35ΩLc0.35e−3HRc0.03ΩF0.75Lline13.2e−4HLline21.9e−3HTable 2
The initial condition of the microgrid system.
ParameterValueVod[380.8, 381.8, 380.4]Iod[11.4, 11.4, 11.4]Ild[11.4, 11.4, 11.4]Vbd[379.5, 380.5, 379]w0314Iline1d−3.8Iline2d7.6Voq[0, 0, 0]Ioq[0.4, −1.45, 1.25]Ild[−5.5, −7.3, −4.6]Vbq[−6, −6, −5]δ[0, 1.9e − 3, −0.0113]Iline1q0.4Iline2q−1.3To ensure that the droop controller parameters are robust after optimization, similar to multiple operating conditions in [6], we design different operating conditions. They are shown in Table 3. Rload of 25 Ω and 20 Ω per phase is equal to the resistive load of 5.6kW and 7.3kW [13]. The objection of this paper is to investigate the optimized parameters in droop controllers that can enhance the small-signal stability of all four operating conditions.Table 3
The multiple operating conditions.
CaseRload1ΩLload1HRload3ΩLload3HCase 1252.5e − 3202.2e − 3Case 2252.2e − 3202.5e − 3Case 3202.5e − 3252.2e − 3Case 4202.2e − 3252.5e − 3The limit of controller parameters can be shown in Table4. 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, we can obtain the value domain of the variables from the root locus of the critical eigenvalues with the variable variation.Table 4
The limits of controller parameters.
VariableLower limitUpper limitmp9e − 61e − 3nq1.3e − 59e − 3Kpv0.00030.1Kiv200450Kpc0.0115.5Kic10e316e4


### 5.2. 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 Table5, 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. Therefore, it demonstrated that the SQP-GS method is effective for optimizing the droop controller parameters.Table 5
The spectral abscissa of microgrid.
Spectral abscissaUnoptimizedSQP-GSCase 172.3025−16.1651Case 272.3024−16.1647Case 372.3112−16.1212Case 472.3110−16.1217The original parameters and the optimized parameters of the controller are given in Tables6 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 Kpv, Kiv, Kpc, and Kic are fairly small, these controller parameters stay the same nearly in the iteration process.Table 6
The original parameter of controller.
DGmpnqKpvKivKpcKicDG16.45e − 42.13e − 30.03943656.1913241DG26.45e − 42.13e − 30.03943656.1913241DG36.45e − 42.13e − 30.03943656.1913241Table 7
The optimal parameters of controller.
DGmpnqKpvKivKpcKicDG17.41e − 51.95e − 40.03973656.1913241DG21.45e − 51.41e − 40.03963656.1913241DG31.57e − 47.97e − 40.03983656.1913241From Figure2, 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. The spectral abscissa of Case 1 is −16.1651. It proves that the stability of the system is improved after the optimization.Figure 2
The comparison of eigenvalues of low frequency in Case 1.In Figure3, 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 small-signal stability sense.Figure 3
Spectral abscissa variation in the iteration process.In Figure4, 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. The stability margin of the system is improved.Figure 4
Objective function value variation in the iteration process.Figure5 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.Figure 5
Sampling radius variation in the iteration process.


### 5.3. Comparison with GA Method

The 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. These 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.Table 8
Comparison of spectral abscissa and computing time between GA method and SQP-GS method.
Spectral abscissaComputing time (h)Methodη¯BestMeanWorstMinMeanMaxGA−15.28−16.05−15.33−15.090.1160.3290.709SQP-GS−16.06−16.12−16.12−16.120.0950.0970.102


## 5.1. Microgrid System Parameters and Four Operating Conditions

We design four operating conditions of the microgrid to verify the effectiveness of the SQP-GS method. The test system in Figure1 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.Table 1
The microgrid system parameters.
ParameterValueLf1.35e−3HCf50e−6FRf0.1ΩRline10.23ΩRline20.35ΩLc0.35e−3HRc0.03ΩF0.75Lline13.2e−4HLline21.9e−3HTable 2
The initial condition of the microgrid system.
ParameterValueVod[380.8, 381.8, 380.4]Iod[11.4, 11.4, 11.4]Ild[11.4, 11.4, 11.4]Vbd[379.5, 380.5, 379]w0314Iline1d−3.8Iline2d7.6Voq[0, 0, 0]Ioq[0.4, −1.45, 1.25]Ild[−5.5, −7.3, −4.6]Vbq[−6, −6, −5]δ[0, 1.9e − 3, −0.0113]Iline1q0.4Iline2q−1.3To ensure that the droop controller parameters are robust after optimization, similar to multiple operating conditions in [6], we design different operating conditions. They are shown in Table 3. Rload of 25 Ω and 20 Ω per phase is equal to the resistive load of 5.6kW and 7.3kW [13]. The objection of this paper is to investigate the optimized parameters in droop controllers that can enhance the small-signal stability of all four operating conditions.Table 3
The multiple operating conditions.
CaseRload1ΩLload1HRload3ΩLload3HCase 1252.5e − 3202.2e − 3Case 2252.2e − 3202.5e − 3Case 3202.5e − 3252.2e − 3Case 4202.2e − 3252.5e − 3The limit of controller parameters can be shown in Table4. 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, we can obtain the value domain of the variables from the root locus of the critical eigenvalues with the variable variation.Table 4
The limits of controller parameters.
VariableLower limitUpper limitmp9e − 61e − 3nq1.3e − 59e − 3Kpv0.00030.1Kiv200450Kpc0.0115.5Kic10e316e4


## 5.2. 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 Table5, 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. Therefore, it demonstrated that the SQP-GS method is effective for optimizing the droop controller parameters.Table 5
The spectral abscissa of microgrid.
Spectral abscissaUnoptimizedSQP-GSCase 172.3025−16.1651Case 272.3024−16.1647Case 372.3112−16.1212Case 472.3110−16.1217The original parameters and the optimized parameters of the controller are given in Tables6 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 Kpv, Kiv, Kpc, and Kic are fairly small, these controller parameters stay the same nearly in the iteration process.Table 6
The original parameter of controller.
DGmpnqKpvKivKpcKicDG16.45e − 42.13e − 30.03943656.1913241DG26.45e − 42.13e − 30.03943656.1913241DG36.45e − 42.13e − 30.03943656.1913241Table 7
The optimal parameters of controller.
DGmpnqKpvKivKpcKicDG17.41e − 51.95e − 40.03973656.1913241DG21.45e − 51.41e − 40.03963656.1913241DG31.57e − 47.97e − 40.03983656.1913241From Figure2, 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. The spectral abscissa of Case 1 is −16.1651. It proves that the stability of the system is improved after the optimization.Figure 2
The comparison of eigenvalues of low frequency in Case 1.In Figure3, 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 small-signal stability sense.Figure 3
Spectral abscissa variation in the iteration process.In Figure4, 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. The stability margin of the system is improved.Figure 4
Objective function value variation in the iteration process.Figure5 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.Figure 5
Sampling radius variation in the iteration process.


## 5.3. Comparison with GA Method

The 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. These 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.Table 8
Comparison of spectral abscissa and computing time between GA method and SQP-GS method.
Spectral abscissaComputing time (h)Methodη¯BestMeanWorstMinMeanMaxGA−15.28−16.05−15.33−15.090.1160.3290.709SQP-GS−16.06−16.12−16.12−16.120.0950.0970.102


## 6. Conclusion

This 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.



---

*Source: 5634274-2021-05-25.xml*