Small-Signal Stability Analysis for Multi-Terminal LVDC Distribution Network Based on Distributed Secondary Control Strategy

: This paper presents a detailed and accurate small-signal analysis model for a four-terminal low-voltage direct current (LVDC) distribution network with distributed secondary control strategy. To tackle the contradiction between power sharing accuracy and average voltage recovery ability of voltage source converters (VSCs), a distributed secondary control strategy is adopted, which is based on average consensus algorithm and local voting protocol. Furthermore, to analyze the stability of LVDC distribution network based on the proposed distributed secondary control strategy, a detailed and accurate small-signal analysis model is formulated which is derived from various non-linear state-space sub-models. The time-domain simulation and electromagnetic simulation are carried out to verify the validity and accuracy of the proposed model. Based on this model, the inﬂuence rules of main eigenvalues are summarized with the variation of PI control and DC line parameters. Time-domain simulations conducted in PSCAD are used to validate the operational limits of the secondary controllers and DC line obtained from the small-signal stability analysis.


Introduction
Recently, the DC distribution system has witnessed a vigorous development because the extensive application of power electronics and large-scale access of renewable energy sources [1]. The DC distribution system has larger power supply capacity and lower line losses with simpler control schemes over traditional AC distribution network [1][2][3]. Multi-terminal DC distribution network is mainly composed by DC converter stations, distributed generations (DGs), load and energy storage systems (ESSs). It has become a research hotspot for superiorities of easy DG access and high system efficiency. The coordinated control among DC converter stations is one of the key technologies to ensure stability, reliability and optimization of the system.
The control targets of the DC converter stations usually involve voltage regulation and power optimization. Generally, the DC converter stations adopt the droop method to implement power distribution according to their own capacity without communication, which can improve the operation reliability of the system [4][5][6]. However, even though the droop control can be implemented easily, it has a contradiction between voltage regulation and power distribution due to the existence of line impedance in the DC system [7]. Therefore, a secondary control for the droop-based DC system is necessary. At present, the centralized and decentralized methods are two main ways to implement secondary control. A single central controller is designed to collect information and send reference to each unit in centralized way [8][9][10][11]. However, the reliability of centralized method is poor due to the failure of the center node or link, and the capability in plug and play is weakened by centralized communication. The decentralized method is more practical than the centralized control, because its reliability and flexibility [12]. In [13], a decentralized droop control method based on virtual resistance in DC system is proposed to compensate power distribution deviations caused by line resistances. However, this method cannot eliminate the voltage deviation. To restore the bus voltage, a hierarchical coordination control based on local information is proposed for DC distribution system to realize smooth operation in steady state and secondary voltage recovery in harsh operating condition as well as global power optimization [14]. Although the simplicity and decentralized nature of the droop control are retained by decentralized method, the global power optimization cannot be realized by only the local information.
Compared with the above two methods, the distributed secondary control method that requires only local computation and information exchange among some neighboring units through a sparse communication network has drawn more attentions [15][16][17]. Reference [7] proposes a hierarchical distributed control strategy for the DC distribution network based on the distributed average algorithm. By designing different weights of the voltage control and power control, the effective balance between the two control objectives can be achieved. A distributed control paradigm based on the reverse droop controller is established to achieve voltage regulation and proportionate load current sharing in a DC network [18]. In [19], a distributed voltage control method for distribution network is proposed to compensate the reactive power and provide the stable voltage for distribution network by using the terminal load and its corresponding power factor correction device. Reference [20] designs a two-level distributed control strategy for the wind farms by regarding the active power utilization rate as the consistent variable. To sum up, although the centralized method can achieve more complex control objectives, it has the risk of system failure. The decentralized method has higher reliability and expansibility, the control effect is not satisfactory. The distributed method can be implemented by multiple distributed controllers working in parallel, which increases the reliability, flexibility, and scalability of the control system. Table 1 lists the advantages and disadvantages according to above three methods. Table 1. Summaries on secondary control methods.

Secondary Control Methods Advantages Disadvantages
Centralized control High control accuracy Poor reliability, scalability Decentralized control High reliability, scalability Poor control accuracy Distributed control High reliability, scalability, accuracy Complicated control scheme Therefore, the distributed manners are considered as promising options for the optimal control of the DC distribution network. In this case, a distributed secondary control strategy for enhancing the power distribution accuracy and recovering the average voltage of the DC converter stations is adopted in this paper. However, there are also many problems remained to be solved before the distributed schemes become more practical for DC distribution network. The stability of the system should be considered first, which is the foundation of voltage regulation and power distribution.
The most commonly used way for stability study is the small-signal analysis method [21]. The small-signal stability analysis of conventional power systems has been well established. More and more researches have focused on the stability analysis of low voltage DC microgrids [22][23][24][25][26] and HVDC transmission & distribution networks [21,[27][28][29][30][31][32][33]. The small-signal model of the converter-based microgrid is established to investigate how circuit and control parameters give rise to the oscillatory modes and the damping performance of these modes [22]. Reference [23] establishes a small-signal model of DC microgrid and analyzes the influence laws of generation parameters to ensure the system stability. A reducedorder model of multi-terminal DC microgrid is established in [24], and the small signal stability of the system is closely related to the sag constant of the distributed generation. In [25], a general methodology considering four kinds of terminals for small-signal stability

•
The detailed and accurate small-signal model of the whole multi-terminal LVDC system containing distributed secondary control is built.

•
The stability of the entire LVDC distribution network is analyzed, and the effects of parameters variations on the small-signal stability are investigated by sensitivity analysis.

•
The time-domain simulation by Simulink and the electromagnetic simulation by PSCAD are both conducted for the entire LVDC system to verify the accuracy of the model.

The Primary Control of VSC-Based Converter Station
In a multi-terminal LVDC distribution system, the coordinated control strategy of VSC-based converter station consists of primary control and secondary control. The primary control usually adopts droop method to distribute active power among all converter stations according to their capacity. The structure of a converter station and its primary control is shown in Figure 2. The primary control consists of P-V droop loop and inner loops, and the SPWM signals achieved by these control loops are applied to converter to accomplish the primary control target.
The ì ï -ḯ ï ï í ï ï ï ï î (1) where U dc * is the reference output voltage of converter station, Udc_ref is the rated voltage provided by secondary control level. The measurement unit of U dc * and Udc_ref is V. k is the droop coefficient. Pdc is the actual output power and the unit is kW. i fd * is the d-axis current reference value of inductance Lf, Udc is the actual output DC voltage and the unit is the same as U dc * . kP,U and kI,U are the parameters of the proportional-integral (PI) controller of the inner voltage loop.

The Primary Control of VSC-Based Converter Station
In a multi-terminal LVDC distribution system, the coordinated control strategy of VSCbased converter station consists of primary control and secondary control. The primary control usually adopts droop method to distribute active power among all converter stations according to their capacity. The structure of a converter station and its primary control is shown in Figure 2. The primary control consists of P-V droop loop and inner loops, and the SPWM signals achieved by these control loops are applied to converter to accomplish the primary control target.
where UdcN is the rated voltage of DC network, Udcmin is the minimum voltage, which can be set to 95% of UdcN. The unit of UdcN and Udcmin is V. PNmax is the maximum output power The P-V droop equation and voltage control loop can be expressed where U dc * is the reference output voltage of converter station, U dc_ref is the rated voltage provided by secondary control level. The measurement unit of U dc * and U dc_ref is V. k is the droop coefficient. P dc is the actual output power and the unit is kW. i f d * is the d-axis current reference value of inductance L f , U dc is the actual output DC voltage and the unit is the same as U dc * . k P,U and k I,U are the parameters of the proportional-integral (PI) controller of the inner voltage loop.
The droop coefficient k can be designed as follows where U dcN is the rated voltage of DC network, U dcmin is the minimum voltage, which can be set to 95% of U dcN . The unit of U dcN and U dcmin is V. P Nmax is the maximum output power of converter station and the unit of it is kW. It can be seen from Equations (1) and (2) that the power sharing of different converter stations is in proportion to their maximum output power [34,35]. However, the DC output voltages of different converter stations are different, which causes the output power deviation between each station and decreases the utilization rate of whole stations. On the other hand, it can be concluded that the elimination of output power deviation will enhance the voltage drop level of different stations. This will increase the probability of voltage violation and the power loss in DC distribution network. Therefore, the primary control has a contradiction between voltage regulation and power distribution. In order to improve the power distribution accuracy and recover DC output voltage, a distributed secondary control strategy is implemented based on the primary control.

The Distributed Secondary Control for DC Distribution Network
As discussed above, the distributed control method requires only local computation and information exchange among some neighboring units through a sparse communication network. It can increase the reliability, flexibility, and scalability of the secondary control level. According to the secondary control target, the distributed secondary control structure of converter station is shown in Figure 3. It can be observed that the secondary control includes voltage secondary control and power secondary control. In voltage regulation step, the average DC output voltage of all converter stations is estimated, and the estimated average voltage U ave is compared with rated voltage U dcN to obtain the compensation item ∆V. In the power regulation step, the average per unit output power P aveλ is estimated, and it is compared with the local value P dcλ of each converter station to generate the compensation item ∆U. Finally, ∆V and ∆U are added to generate the rated voltage U dc_ref for the primary level. Then, the primary control in this paper can be rewritten as follows where ∆V i and ∆U i are the voltage compensation items.  As shown in Figure 3, the estimation of average voltage and per unit power are all realized by distributed algorithm. As the communication delay is too small that the sampling control cycle of the secondary level can be designed shorter. Thus, the continuous model is considered in this paper. The dynamic consensus algorithm is often used to estimate the average value of the whole measurements [36]. Let xi represents the measurement of node i, while xavei is the estimation value of global averaging. Then, the dynamic consensus algorithm can be described as follows where aij is the element of adjacency matrix A. A = [aij] is a (0,1)-matrix with zeros on its diagonal, and the off-diagonal entry aij is one if node j∈Ni.
where CE is the control gain, kP,V and kI,V are the PI parameters of voltage secondary control. Uavei is the average voltage estimated by station i. On the other hand, in order to eliminate the deviation of per unit output power, the local voting protocol is used to estimate the average per unit value, thus the power regulation process can be described as follows = , 1,2,..., .
where Pdcλi is the per unit power of station i, and can be calculated by As shown in Figure 3, the estimation of average voltage and per unit power are all realized by distributed algorithm. As the communication delay is too small that the sampling control cycle of the secondary level can be designed shorter. Thus, the continuous model is considered in this paper. The dynamic consensus algorithm is often used to estimate the average value of the whole measurements [36]. Let x i represents the measurement of node i, while x avei is the estimation value of global averaging. Then, the dynamic consensus algorithm can be described as follows where a ij is the element of adjacency matrix A. A = [a ij ] is a (0,1)-matrix with zeros on its diagonal, and the off-diagonal entry a ij is one if node j∈N i . N i represents the set of adjacent nodes of node i, namely, a node j belongs to N i if there exists a communication link between it and node i. It can be said that the nodes have reached a consensus if and only if x avei = x avej , for all i, j. If the communication is bidirectional, then all nodes will converge to the average value of all measurements. For converter station i, the measurements can be the DC output voltage, per unit output power, et al. Based on the dynamic consensus algorithm, the voltage regulation process can be described as follows where C E is the control gain, k P,V and k I,V are the PI parameters of voltage secondary control. U avei is the average voltage estimated by station i. On the other hand, in order to eliminate the deviation of per unit output power, the local voting protocol is used to estimate the average per unit value, thus the power regulation process can be described as follows where P dcλi is the per unit power of station i, and can be calculated by where P Ni is the rated capacity of station i. P aveλi is the average per unit output power estimated by station i. d ij is the voting weight, where d ij > 0 if node j belongs to N i and d ij = 0, otherwise. As the below steady state analysis shows, when the coefficient matrix D = [d ij ] is designed as a double-stochastic matrix, the power deviation between different Electronics 2021, 10, 1575 7 of 30 stations will be eliminated asymptotically. k P,P and k I,P are the PI parameters of power secondary control. The output voltage reference U dc * then can be calculated according to (3), and the current control loop is used to follow the current reference i f d * from the voltage control loop. The bandwidth of this loop is designed to be large enough to maintain the power balance of the converter station while load changes. Within the distributed secondary control, the average DC output voltage can be recovered and the per unit output power deviation of all stations can be eliminated, respectively.

Steady-State Analysis
Laplace Transform Method is often used to verify the steady control effect of the proposed control strategy. The Laplace Transform of Equations (3), (5) and (6) are shown as follows U dc where U dcN (s) = (U dcN /s)1, 1 is a n × 1 vector whose elements are all equal to one, Let L n2 = D − I n , I n is the unit matrix, then L n1 and L n2 are all the so called Laplace matrix. Suppose the system is stable, then U dc * = U dc , and substitute Equation (9) into (8) [ where H(s) = diag{k P,V + k I,V /s}, G(s) = diag{k P,P + k I,P /s}. Multiplying both sides of Equation (10) by s 2 and implementing the Final Value Theorem, (10) can be described as follows where U s dc is the steady DC voltage vector, and P s dcλ is the steady per unit power vector. Since Equation (11) is guaranteed for any k I,V and k I,P , it can be derived that For Equation (12), according to [37], Q = lim s→0 H obs is the averaging matrix, whose elements are all equal to 1/n. Multiplying both sides by 1 T , Equation (12) can be expressed as follows For Equation (13), because L n2 is the Laplace matrix, (13) only has one solution, that is P s dcλ = L * 1, in other words The above analysis shows that the per unit output power of each converter station will converge to L * , and its value is determined by the net load in the DC distribution network. That is, all converter stations will distribute their output power according to their rated capacity accurately. In addition, the average DC output voltage will converge to the rated value U dcN .

State-Space Model of the Individual Converter Station with Primary Control
To obtain the model of the entire DC distribution system, the state-space model of an individual converter station with primary control needed to be first established. In this section, all elements in the converter station system shown in Figure 2 are expressed in terms of mathematical equations as a basis for the small-signal model.

• Phase Locked Loop Equations
In earlier studies, d PLL q PLL coordinate system based on the AC input voltage u o is usually used for the small-signal modelling. However, when the transient process of PLL cannot be neglected, the state space equations should be established (16) where u oq is the input voltage of filter in the q PLL -axis. ω and ω o are the angular frequency of AC grid and rated angular frequency, respectively. The unit of ω o is rad/s. ω PLL is the output angular frequency of PLL. δ PLL is the output angular of PLL, which is the phase deviation between u o and e, as shown in Figure 2. z is an auxiliary state variable representing the integral part of the PI controller in the PLL, and the parameters of PI controller are k P,PLL and k I,PLL .

• Input Filter Equations
A typical LC filter is composed of the filter inductance L f and its equivalent resistance R f , and the filter capacitor, which is given by C f . Based on the d PLL q PLL coordinate system, the state equations of the LC filter can be presented as where i fd and i fq are the filter inductor current in the d PLL q PLL -axis. u ed and u eq are the input voltage of the converter in the d PLL q PLL -axis, respectively. u od is the input voltage of filter in the d PLL -axis. i od and i oq are the input current in the d PLL q PLL -axis. The units of R f , L f and C f are mΩ, mH and µF, respectively.

• AC Connection Line Equations
The converter station is connected to the AC grid through a transmission line, which consists of resistance R c and inductance L c . As the AC grid is considered as an ideal source, the absolute stable d G q G coordinate system based on it is used to establish the state space equations of the transmission line where e d and e q are the AC grid voltage in the d G q G -axis, and they are constant values with e q = 0. The unit of e d and e q is V. i od_G and i oq_G are the input current in the d G q G -axis, respectively. u od_G and u oq_G are the input voltage of LC filter in the d G q G -axis, respectively.

• Converter Equations
With the losses in switching elements (e.g., IGBTs and diodes) neglected, it can be assumed that the input voltage of the converter station will be equal to its reference value, as follows where u ed * and u eq * are the voltage reference of the inverter in the d PLL q PLL -axis.

• DC Connection Line Equations
The converter station is also connected to the DC bus through the transmission line, with the line parameters consisting of resistance R dc and inductance L dc . The state space equations of the DC connection line are given as follows where I dc is the current of transmission line, and U dc is the DC output voltage of capacitor C dc . I g is the output current of converter and U b is the voltage of DC bus. The measurement units of C dc , L dc and R dc are F, mH and mΩ, respectively.

• Power Calculator
The power calculator is applied in this case for calculating the instantaneous active power. With the losses neglected in converter station, the AC active power is equal to the power transmitted to the DC port, which is shown in (21) where ω c is the cut-off frequency of low-pass filter and the unit is rad/s.

• Voltage and Current Loop Control
The voltage and current loop control block diagram are shown in Figure 4 in detail. k P,U and k I,U are the parameters of the PI controller of the outer voltage loop, and k P,C and k I,C are the parameters of the PI controllers of the inner current loop. The DC output voltage reference (U dc * ) of converter station is obtained through droop equation and the secondary control value ∆V and ∆U, as mentioned in Section 2.2.
Electronics 2021, 10, x FOR PEER REVIEW 10 of 32 The voltage and current loop control block diagram are shown in Figure 4 in detail. kP,U and kI,U are the parameters of the PI controller of the outer voltage loop, and kP,C and kI,C are the parameters of the PI controllers of the inner current loop. The DC output voltage reference (U dc * ) of converter station is obtained through droop equation and the secondary control value ΔV and ΔU, as mentioned in Section 2.2.  According to Figure 4, the space state equations of voltage and current loop control can be written where state variables γ and λ are defined to assess the dynamic properties of the PI controllers in the voltage and current loop, and the corresponding state equation in the dPLLqPLL-axis is shown as According to Figure 4, the space state equations of voltage and current loop control can be written where state variables γ and λ are defined to assess the dynamic properties of the PI controllers in the voltage and current loop, and the corresponding state equation in the d PLL q PLL -axis is shown as

State-Space Model of the Distributed Secondary Control
According to the proposed distributed secondary control strategy described by (5) and (6), the space state model is described as follows T is defined to describe the error accumulation between estimated voltage and rated voltage, while Φ p = [Φ p1 , Φ p2 , . . . , Φ pn ] T is defined to describe the error accumulation between estimated per unit power and local per unit power. The corresponding state equation is

State-Space Model of the DC Distribution Line and Load
As different DC network topologies have different state space models, the fourterminal ring topology DC distribution network is taken for the example in this paper to establish the state space model, as shown in Figure 1. On the other hand, the wye-delta transform is adopted to simplify the DC distribution line and load model, thus the equivalent model of the DC distribution network line and inner constant power load (CPL) is shown in Figure 5. where

State-Space Model of the DC Distribution Line and Load
As different DC network topologies have different state space models, the four-terminal ring topology DC distribution network is taken for the example in this paper to establish the state space model, as shown in Figure 1. On the other hand, the wye-delta transform is adopted to simplify the DC distribution line and load model, thus the equivalent model of the DC distribution network line and inner constant power load (CPL) is shown in Figure 5. •

DC Distribution Line Equations
As Figure 5a shows, all the converter stations are connected to each other through DC distribution lines. The state space equations of lines are written as follows

• DC Distribution Line Equations
As Figure 5a shows, all the converter stations are connected to each other through DC distribution lines. The state space equations of lines are written as follows where U n is the node voltage of the equivalent constant power load, I l is the branch current, L l and R l are line inductance and resistance, respectively. C b represents the grounded capacitance. I CPL is the current of CPL. The measurement units of L l , R l and C b are mH, Ω and µF, respectively •

Constant Power Load Equations
As Figure 5b shows, the equivalent constant power load in the DC distribution system is modeled as DC/DC bulk converter. The state space equations are as follows where C n , L D , C D are the parameters of DC/DC bulk converter, respectively. U L and I L are the capacity voltage and inductor current of resistance side. U L_ref is the rated voltage, and R is the equivalent resistance of load. D is the duty ratio of control signal, k P,L and k I,L are the parameters of the PI controller of the converter. The measurement units of L D , C D and C n are mH, mF and mF, respectively.

Small-Signal Model of the DC Distribution Lines and Loads
Based on (26) and (27), the small-signal state space model contains 28 state variables and 4 output variables, which are given by Thus, the small-signal equation and the output equation are as follows where ∆I dc = [∆I dc1 , ∆I dc2 , ∆I dc3 , ∆I dc4 ] T and ∆R = [∆R 1 , ∆R 2 , ∆R 3 , ∆R 4 ] T . Both of them are the input variables. Appendix A shows the coefficient matrix A s0 , B s0 , C s0 , D s0 .

Small-Signal Model of the Multiple Converter Stations
As previously discussed, all the state space equations needed for modeling the converter stations have been introduced. Based on (16)- (23), the state space model of i-th converter station without proposed distributed secondary control contains 14 state variables and 3 output variables (∆I dci , ∆P dci , ∆U dci ) The small-signal equation of i-th converter station then can be described as follows where As there are four converter stations in the DC system, the entire small-signal model is established as follows: The output equations are as follows where
where ∆x s = [∆x s0 T , ∆x s1 T , ∆x s2 T ] T and represents all state variables in the DC distribution network. ∆R represents the small variation of constant power load. ∆ω represents the small frequency variation of AC grid, and ∆E represents the small voltage magnitude variation of AC grid. Appendix C shows the coefficient matrix

Verification of the Established Small-Signal Model
In this section, the results from the time-domain simulation in MATLAB/Simulink and electromagnetic simulation in PSCAD are presented for the verification of the proposed small-signal model. The structure of the LVDC distribution network for case study is shown as Figure 1. There are four converter stations, and the DC sources (PV, WT, DC loads) are integrated as a constant power load for simplicity. Thus, there are four CPLs in the network. The parameters of the DC distribution network and four converter stations are shown in Table 2. The communication topology of the stations is a ring, and its corresponding adjacency matrix A and coefficient matrix D are shown in Figure 6. Table 2. The parameters of DC distribution system.

Parameter
Value Parameter Value

Verification of the Established Small-Signal Model
In this section, the results from the time-domain simulation in MATLAB/Simulink and electromagnetic simulation in PSCAD are presented for the verification of the proposed small-signal model. The structure of the LVDC distribution network for case study is shown as Figure 1. There are four converter stations, and the DC sources (PV, WT, DC loads) are integrated as a constant power load for simplicity. Thus, there are four CPLs in the network. The parameters of the DC distribution network and four converter stations are shown in Table 2. The communication topology of the stations is a ring, and its corresponding adjacency matrix A and coefficient matrix D are shown in Figure 6.    . It is clearly that the average voltage of DC side is 800 V, and the per unit output power of all converters are all 0.75. Thus, the effectiveness of the distributed secondary control is verified. The introduced small disturbance is a 5 kW step change of the CPL1 from 80 kW to 85 kW at 3.5 s. The curves of different state variables (U n1 , P dc1 , U dc1 , I dc1 , u oq1 , u od1 , i fq1 , i fd1 ) obtained using the small-signal model and electromagnetic simulation are compared in Figure 7. It can be seen from Figure 7 that the time-domain responses of different state variables calculated by the established small-signal model using MATLAB coincide with those obtained by PSCAD simulation. This indicates that the accuracy of the established small-signal model is suitable for the stability analysis of the DC distribution network. different state variables (Un1, Pdc1, Udc1, Idc1, uoq1, uod1, ifq1, ifd1) obtained using the small-signal model and electromagnetic simulation are compared in Figure 7. It can be seen from Figure 7 that the time-domain responses of different state variables calculated by the established small-signal model using MATLAB coincide with those obtained by PSCAD simulation. This indicates that the accuracy of the established small-signal model is suitable for the stability analysis of the DC distribution network. 3

Small-Signal Stability Analysis
The system eigenvalues are calculated and listed in Table 3 based on the parameters in Table 2. In the small-signal analysis, it is necessary to quantify the contribution of each state to a specific eigenvalue according to [38]. Thus, the major participants of each eigenvalue are exhibited in Table 3 as well.

Small-Signal Stability Analysis
The system eigenvalues are calculated and listed in Table 3 based on the parameters in Table 2. In the small-signal analysis, it is necessary to quantify the contribution of each state to a specific eigenvalue according to [38]. Thus, the major participants of each eigenvalue are exhibited in Table 3 as well.
As shown in Table 3, all eigenvalues have negative real parts indicating a stable operating condition for the system. It can be also seen that these eigenvalues fall in to three different clusters, namely, the high damped modes, the medium damped modes and low damped modes. The participation factor indicates that the inductive current of LC filter and inner current loop influence greatly the high damped modes λ 1 -λ 16 . These modes are not detrimental to the system stability. The medium damped modes λ 17 -λ 63 are associated with the states of DC network (DC line currents, bus voltages, node voltages) and their effects on the AC side input currents and voltages. Additionally, the low damped modes considered as the dominate modes are associated with the states of secondary control. These dominate modes are greatly influenced by the parameters of secondary control and PLL. Moreover, the DC line parameters can be also selected to optimize the transient behavior of the DC network.
By plotting the eigenvalue loci in the complex plane, the impacts of the main parameters (such as the PI parameters of the secondary control, communication delay, the parameters of PLL, and the impedance parameters of DC distribution) on system dynamics can be observed. The eigenvalue loci in the complex plane with parameters variations are shown in Figures 8, 12, 13 and 15. In Figures 8, 12, 13 and 15, the eigenvalues are marked with stars. The arrows indicate the direction of the parameter variation and black circles indicate the parameters of critical stability. For simple analysis, only the important eigenvalues are plotted. The following subsections emphasize on the influence of these parameters on the system dynamics.

Distributed Secondary Controller Effect
This section presents the impact of the secondary voltage and power controller parameters k P,V , k I,V , k P,P and k I,P on the small-signal stability and the verification of the stable region through PSCAD simulation results. When the individual parameter is considered and changed for plotting eigenvalue loci, the other parameters are all consistent with the value as presented in Table 1. The varying range of individual parameter k P,V , k I,V , k P,P and k I,P is 1 to 50, 1 to 150, 1 to 10,000 and 1 to 300,000, respectively. The eigenvalue loci are shown in Figure 8, and Figure 9 shows the time-domain simulation results through PSCAD.
With a variation of proportional parameter k P,V , which changed from 1 to 50, the parametric eigenvalue analysis is repeated. In the overview of the eigenvalue loci shown in Figure 8a, λ 80-82 move towards the left half of the complex plane along the real axis. The eigenpair λ 64,65 also move towards left, which indicates that the real part of λ 64,65 decrease while the damping ratio increase. Eigenvalue λ 83 moves towards the original point along the real axis, and the eigenpair λ 53,54 move towards the right part of the complex plane with the real part increasing dramatically. Furthermore, when the value of k P,V is larger than 39.8, an obvious vibration appeared, and the λ 53,54 enter into the unstable region of the complex plane that bringing about the instability of the system. Figure 8b exhibits the eigenvalue loci of the integral parameter k I,V changing from 1 to 150, where eigenpair λ 53-54 are neglected due to unobvious changes in their locations and are not discussed in the following. With the increasing of k I,V , it is clear that eigenvalues λ 80-83 all approach the original point along the real axis and the eigenpair λ 64,65 move towards the right half of the complex plane gradually. When the value of k I,V is larger than 118.5, the λ 53,54 cross the imaginary axis, which means that they enter into the unstable region of the complex plane, generating the instability of the system. Figure 8c shows the eigenvalue loci of the parameter k P,P changing from 1 to 10000. For the power secondary control, eigenpairs λ 53,54 , λ 66,67 , λ 68,69 and λ 70,71 are discussed and analyzed here. It can be seen that with the k P,P increasing, λ 70,71 move towards the left direction with the imaginary part becoming zero, which indicates that the oscillation of the output power disappears gradually. Meanwhile, λ 66,67 and λ 68,69 move towards the opposite direction with the imaginary part becoming zero and real part increasing. In addition, λ 53,54 move towards the right part of the complex plane fast and as soon as k P,P is more than 7500, λ 53,54 cross the imaginary axis immediately, resulting in unstable operating state.
influencing the system stability. Figure 8d reveals the eigenvalue loci with kI,P varying from 1 to 300,000. Eigenpair λ53,54 is neglected while λ64,65 is added in the following analysis. In Figure 8d, λ66,67 turn to the left part of complex plane with moving towards the right direction at first. Eigenvalues λ64,65, λ68,69 and λ70,71 have the same changing patterns, that is the imaginary parts move up and down respectively and the absolute value of imaginary parts increase with the parameter adding up, which means the system would experience a slow response. To verify the eigenvalue loci shown in Figure 8, the time-domain responses of the multi-terminal LVDC distribution network simulated in PSCAD are plotted, as presented in Figure 9. Figure 9a shows that when kP,V steps from 2 to 45 at 3.5 s, the state variable Un1 experiences an amplitude oscillation, which corroborates the results obtained in Figure  8a. Similarly, Figure 9b indicates that when kI,V steps from 10 to 2000 at 3.5 s, the DC voltage Udc1 experiences an amplitude oscillation. However, this oscillation happens at 3.7 s and the value of kI,V in this period is far larger than 118.5 shown in Figure 8b. Furthermore, Figure 9c shows the time-domain simulation results when kP,P changes from 200 to 65,000. It can obviously be seen that state variable Un1 undergoes an amplitude oscillation as well. Nevertheless, the value of kP,P is still well over the stable region shown in Figure 8c. It indicates that the accuracy of MATLAB simulation is a little higher than the results simulated in PSCAD. As shown clearly in Figure 9d, when kI,P varys from 2000 to 500,000, the curve of Udc1 stay unchanged, which is in agreement with the results obtained in Figure  8d. According to the above analysis, the design rule of PI parameters of secondary control for the system stability can be concluded as follows: 1. The augment of kP,V and kI,V may decrease the relative damping with low-frequency eigenvalues, which goes against the system stability. Too large a value of kP,V and kI,V can cause low-frequency oscillation in the multi-terminal LVDC distribution network. 2. Within the range of system stability, the gradual increase of kI,P will decrease the damping ratio of the low-frequency eigenvalues. Furthermore, system also has the possibility of losing stability in the condition of large value of kP,P.

Time Delay Effect
This section presents the impact of the time delay on the small-signal stability and the verification by PSCAD simulation. The impact of time delay cannot be neglected when the sampling rate of the secondary control level is very high. In order to study the impact on system stability, the expression of delay e −τs should be simplified firstly. Using the first- The integral parameter of power secondary control k I,P is also the main parameter influencing the system stability. Figure 8d reveals the eigenvalue loci with k I,P varying from 1 to 300,000. Eigenpair λ 53,54 is neglected while λ 64,65 is added in the following analysis. In Figure 8d, λ 66,67 turn to the left part of complex plane with moving towards the right direction at first. Eigenvalues λ 64,65 , λ 68,69 and λ 70,71 have the same changing patterns, that is the imaginary parts move up and down respectively and the absolute value of imaginary parts increase with the parameter adding up, which means the system would experience a slow response.
To verify the eigenvalue loci shown in Figure 8, the time-domain responses of the multi-terminal LVDC distribution network simulated in PSCAD are plotted, as presented in Figure 9. Figure 9a shows that when k P,V steps from 2 to 45 at 3.5 s, the state variable U n1 experiences an amplitude oscillation, which corroborates the results obtained in Figure 8a. Similarly, Figure 9b indicates that when k I,V steps from 10 to 2000 at 3.5 s, the DC voltage U dc1 experiences an amplitude oscillation. However, this oscillation happens at 3.7 s and the value of k I,V in this period is far larger than 118.5 shown in Figure 8b. Furthermore, Figure 9c shows the time-domain simulation results when k P,P changes from 200 to 65,000. It can obviously be seen that state variable U n1 undergoes an amplitude oscillation as well. Nevertheless, the value of k P,P is still well over the stable region shown in Figure 8c. It indicates that the accuracy of MATLAB simulation is a little higher than the results simulated in PSCAD. As shown clearly in Figure 9d, when k I,P varys from 2000 to 500,000, the curve of U dc1 stay unchanged, which is in agreement with the results obtained in Figure 8d.
According to the above analysis, the design rule of PI parameters of secondary control for the system stability can be concluded as follows: 1.
The augment of k P,V and k I,V may decrease the relative damping with low-frequency eigenvalues, which goes against the system stability. Too large a value of k P,V and k I,V can cause low-frequency oscillation in the multi-terminal LVDC distribution network.

2.
Within the range of system stability, the gradual increase of k I,P will decrease the damping ratio of the low-frequency eigenvalues. Furthermore, system also has the possibility of losing stability in the condition of large value of k P,P .

Time Delay Effect
This section presents the impact of the time delay on the small-signal stability and the verification by PSCAD simulation. The impact of time delay cannot be neglected when the sampling rate of the secondary control level is very high. In order to study the impact on system stability, the expression of delay e −τs should be simplified firstly. Using the first-order approximation of the Taylor series expansion, the e −τs can be written as follows When the delay is small, the equation above stands. Obviously, the equation is a first-order lag function, which will bring a state variable. Based on the simplified equation, the distributed secondary control with communication delay τ can be described as follows where, ∆x s is the state variable including the new variable ∆U ave and ∆P dcλ , which are related to time delay. These eigenvalues related to time delay are calculated and listed in Table 4 based on the parameters of DC distribution network and τ = 3 ms.
As shown in Table 4, all new eigenvalues have negative real parts indicating a stable operating condition for the system when τ = 3 ms. That is, the small delay has little impact on the stability of the system. In order to study the impact of large delay, the eigenvalue loci are plotted according to different time delay. The varying range of time delay is 3 ms to 50 ms, and the eigenvalue loci are shown in Figure 10.
With the increasing of τ, it is clear that all eigenvalues approach the original point along the real axis first. Then, all the eigenvalues are equal to −50 at τ = 20 ms. From then on, some eigenvalues' imaginary parts begin to increase, and they move towards the right half of the complex plane gradually. From Figure 10, it can be indicated that the system is still stable during 3 ms to 50 ms, however, the stability of the DC distribution network is weakened. With the increasing of τ, the oscillation will appear on the estimated average voltages, and the damping ratio decrease gradually. Furthermore, the simulation by PSCAD is implemented as shown in Figure 11. From Figure 11, it can be seen that the U ave will appear low-frequency oscillation after a small disturbance when τ = 30 ms, while the responses are smooth when τ = 3 ms. The damped time of the oscillation is about 0.42 s, which is similar with the result in Figure 10. Thus, the time domain simulation verifies the results from eigenvalue analysis in Figure 10. With the increasing of τ, it is clear that all eigenvalues approach the original point along the real axis first. Then, all the eigenvalues are equal to -50 at τ = 20 ms. From then on, some eigenvalues' imaginary parts begin to increase, and they move towards the right half of the complex plane gradually. From Figure 10, it can be indicated that the system is still stable during 3 ms to 50 ms, however, the stability of the DC distribution network is weakened. With the increasing of τ, the oscillation will appear on the estimated average voltages, and the damping ratio decrease gradually. Furthermore, the simulation by PSCAD is implemented as shown in Figure 11. From Figure 11, it can be seen that the U ave will appear low-frequency oscillation after a small disturbance when τ = 30 ms, while the responses are smooth when τ = 3 ms. The damped time of the oscillation is about 0.42 s， which is similar with the result in Figure 10. Thus, the time domain simulation verifies the results from eigenvalue analysis in Figure 10.  Figure 11. Responses of state variable U ave at the condition that τ = 3 ms and 30 ms.

Phase Locked Loop Effect
This part mainly discusses the effects of PI parameters of phase locked loop on system stability. The secondary control parameters are all original parameters presented in   With the increasing of τ, it is clear that all eigenvalues approach the original po along the real axis first. Then, all the eigenvalues are equal to -50 at τ = 20 ms. From th on, some eigenvalues' imaginary parts begin to increase, and they move towards the rig half of the complex plane gradually. From Figure 10, it can be indicated that the system still stable during 3 ms to 50 ms, however, the stability of the DC distribution network weakened. With the increasing of τ, the oscillation will appear on the estimated avera voltages, and the damping ratio decrease gradually. Furthermore, the simulation PSCAD is implemented as shown in Figure 11. From Figure 11, it can be seen that the U will appear low-frequency oscillation after a small disturbance when τ = 30 ms, while responses are smooth when τ = 3 ms. The damped time of the oscillation is about 0.42 s which is similar with the result in Figure 10. Thus, the time domain simulation verifies results from eigenvalue analysis in Figure 10.  Figure 11. Responses of state variable U ave at the condition that τ = 3 ms and 30 ms.

Phase Locked Loop Effect
This part mainly discusses the effects of PI parameters of phase locked loop on s tem stability. The secondary control parameters are all original parameters presented Figure 11. Responses of state variable U ave at the condition that τ = 3 ms and 30 ms.

Phase Locked Loop Effect
This part mainly discusses the effects of PI parameters of phase locked loop on system stability. The secondary control parameters are all original parameters presented in Table 2 when PI parameters of PLL are studied. The changing range of PI parameters are 1 to 520 and 1 to 1000 respectively. The eigenvalues loci are plotted in Figure 12.  For phase locked loop, with k P,PLL and k I,PLL varying, the eigenvalue loci are shown in Figure 12. In Figure 12a, it is clearly found that λ 88,91 move towards the left direction, which means they have little impacts on system stability when k P,PLL changes. Apart from these, eigenvalue λ 81,82,83 as well as eigenpairs λ 25,26 , λ 29,30 , λ 31,32 and λ 37,38 move towards the right part of complex plane with k P,PLL varying gradually. One point to note is that, when the value of k P,PLL is up to 502.5, the stablity of system will face threats due to eigenpairs λ 25,26 , λ 29,30 , λ 31,32 and λ 37,38 cross the imaginary axis. Figure 12b shows the eigenvalue loci with k I,PLL varying from 1 to 1000. Eigenvalues λ 81,82,83 all move towards the left direction along the negative real axis while λ 89,91 move towards the opposite direction. Furthermore, the absolute values of imaginary part of eigenpairs λ 72,73 , λ 74,75 , λ 76,77 and λ 76,77 increase as k I,PLL changed from 1 to 1000, which means that the damping ratio decrease slowly. Although the response speed of the system declines, the stability of it still can be ensured.
According to the above analysis, the design rule of PI parameters of phase locked loop for the system stability can be concluded as follows: 1.
The augment of k P,PLL has an important effect on the medium damped modes. Too large a value of k P,PLL can cause an unstable state of system.

2.
The parameter k I,PLL has something to do with the system oscillation. And the increase of k I,PLL may decrease the damping ratio with high-frequency eigenvalues without affecting the system stability.

DC Distribution Line Effect
According to the analysis above, the DC line parameters also have an impact on the system's stability. Here, the effects of DC inductance L l1 -L l8 and DC resistance R l1 -R l8 are considered and studied shown in Figure 13, where the varying scopes of L l1 -L l8 and R l1 -R l8 are 3.2 × 10 −5 H to 3.2 × 10 −3 H and 0.01 Ω to 1 Ω respectively. Simulations in PSCAD with parameters varying confirm the correctness of stable region as shown in Figure 14.
participants are Un1-Un4. Figure 13b shows the eigenvalue loci in the parametric sweep of Rl1-Rl8 varying in the range of 0.01 Ω to 1 Ω. Eigenvalues λ49-63 all move towards the imaginary axis or along the imaginary axis in a left direction. Especially, when the value of Rl1-Rl8 is less than 0.04 Ω, the eigenpairs λ49,50 exist in the right part of the coordinate system, causing system unstable.  Figure 14a indicates that when Ll1-Ll8 steps from 0.32 mH to 5 mH at 3.5 s, the voltage of CPL1 Un1 experiences an amplitude oscillation. This oscillation is unobvious at the beginning and becomes conspicuous after 3.7 s and it is a divergent oscillation, which demonstrates that there is small-signal instability appearing in the system when inductance parameter is larger than 1.28 mH. Similarly, in Figure 14b, when Rl1-Rl8 step from 0.1 Ω to 0.001 Ω at 3.5 s, the state variable Un1 also experiences an amplitude oscillation apparently. The unstable state after 3.5 s demonstrates that the time-domain simulation results are in agreement with the stable region obtained in Figure 13b and small value of Rl1-Rl8 might cause instability. Moreover, the influence of Rl1-Rl8 is not as obvious as that of Ll1-Ll8.   Figure 14a indicates that when Ll1-Ll8 steps from 0.32 mH to 5 mH at 3.5 s, the voltage of CPL1 Un1 experiences an amplitude oscillation. This oscillation is unobvious at the beginning and becomes conspicuous after 3.7 s and it is a divergent oscillation, which demonstrates that there is small-signal instability appearing in the system when inductance parameter is larger than 1.28 mH. Similarly, in Figure 14b, when Rl1-Rl8 step from 0.1 Ω to 0.001 Ω at 3.5 s, the state variable Un1 also experiences an amplitude oscillation apparently. The unstable state after 3.5 s demonstrates that the time-domain simulation results are in agreement with the stable region obtained in Figure 13b and small value of Rl1-Rl8 might cause instability. Moreover, the influence of Rl1-Rl8 is not as obvious as that of Ll1-Ll8.  In Figure 13a, all eigenvalues including λ 49-63 move towards the right part of complex plane as parameters L l1 -L l8 changed from 3.2 × 10 −5 H to 3.2 × 10 −3 H. When L l1 -L l8 are all set to 1.28 × 10 −3 H, the eigenpairs λ 49,50 , λ 53,54 and λ 55,56 all enter into the right part of complex plane, indicating that system lose stability. It is noted that when the parameters of L l1 -L l8 are 1.28 × 10 −3 H, the first eigenvalues crossing the imaginary axis are λ 49,50 , whose major participants are U n1 -U n4 . Figure 13b shows the eigenvalue loci in the parametric sweep of R l1 -R l8 varying in the range of 0.01 Ω to 1 Ω. Eigenvalues λ 49-63 all move towards the imaginary axis or along the imaginary axis in a left direction. Especially, when the value of R l1 -R l8 is less than 0.04 Ω, the eigenpairs λ 49,50 exist in the right part of the coordinate system, causing system unstable. Figure 14a indicates that when L l1 -L l8 steps from 0.32 mH to 5 mH at 3.5 s, the voltage of CPL1 U n1 experiences an amplitude oscillation. This oscillation is unobvious at the beginning and becomes conspicuous after 3.7 s and it is a divergent oscillation, which demonstrates that there is small-signal instability appearing in the system when inductance parameter is larger than 1.28 mH. Similarly, in Figure 14b, when R l1 -R l8 step from 0.1 Ω to 0.001 Ω at 3.5 s, the state variable U n1 also experiences an amplitude oscillation apparently. The unstable state after 3.5 s demonstrates that the time-domain simulation results are in agreement with the stable region obtained in Figure 13b and small value of R l1 -R l8 might cause instability. Moreover, the influence of R l1 -R l8 is not as obvious as that of L l1 -L l8 .
According to the above analysis, the design rule of parameters L l1 -L l8 and R l1 -R l8 of DC line for the system stability can be concluded as follows: 1.
The increase of L l1 -L l8 has significant influences on the medium damped modes and low damped modes. Too large a value of L l1 -L l8 can cause an unstable state of system. 2.
The augment of R l1 -R l8 may increase the relative damping with low-frequency and medium-frequency eigenvalues. Too small a value of R l1 -R l8 has a tendency of causing system to lose stability.
These rules are significant for parameters design and optimization of the multiterminal LVDC system such as PI parameters of the secondary control, which are the key parameters for VSC control. During the parameters design based on the small-signal model analysis, it must be ensured that all the eigenvalues of the system in every frequency range are provided with enough damping. The parameter design for the multi-terminal LVDC system must combine the stability of the eigenvalues in every frequency range.

Different Load Effect
This part mainly discusses the effects of constant power load (CPL) and resistance load on system stability. The CPL is the load whose electricity power does not change, while the resistance load will change according to the voltage. The eigenvalues of the system are used to explain the influence of different loads on system stability. Further, the simulation in time domain is implemented to illustrate the stability problem. Firstly, the DC distribution network with CPL and the DC distribution network with resistance load are modeled respectively based on the small signal modeling method presented in this paper. Then, the eigenvalues can be calculated based on A s_CPL and A s_RES , where A s_CPL is the state matrix of the CPL system, and A s_RES is the state matrix of the resistance system. The eigenvalues or modes that are related to U dc , I dc , U n and I l are shown in Figure 15. According to the above analysis, the design rule of parameters Ll1-Ll8 and Rl1-Rl8 of DC line for the system stability can be concluded as follows: 1. The increase of Ll1-Ll8 has significant influences on the medium damped modes and low damped modes. Too large a value of Ll1-Ll8 can cause an unstable state of system 2. The augment of Rl1-Rl8 may increase the relative damping with low-frequency and medium-frequency eigenvalues. Too small a value of Rl1-Rl8 has a tendency of causing system to lose stability.
These rules are significant for parameters design and optimization of the multi-terminal LVDC system such as PI parameters of the secondary control, which are the key parameters for VSC control. During the parameters design based on the small-signal model analysis, it must be ensured that all the eigenvalues of the system in every frequency range are provided with enough damping. The parameter design for the multiterminal LVDC system must combine the stability of the eigenvalues in every frequency range.

Different Load Effect
This part mainly discusses the effects of constant power load (CPL) and resistance load on system stability. The CPL is the load whose electricity power does not change while the resistance load will change according to the voltage. The eigenvalues of the system are used to explain the influence of different loads on system stability. Further, the simulation in time domain is implemented to illustrate the stability problem. Firstly, the DC distribution network with CPL and the DC distribution network with resistance load are modeled respectively based on the small signal modeling method presented in this paper. Then, the eigenvalues can be calculated based on As_CPL and As_RES, where As_CPL is the state matrix of the CPL system, and As_RES is the state matrix of the resistance system The eigenvalues or modes that are related to Udc, Idc, Un and Il are shown in Figure 15. From Figure 15, it can be seen that the number of eigenvalues of As_CPL is more than that of As_RES, because the CPL will bring extra modes from DC/DC converter. The both systems are stable, however, the extra modes that the CPL generated have a large imaginary part, and they are closer to real axis. Thus, it can be indicated that the variation of CPL will bring higher oscillation (about 239 Hz), and the damped time will be longer However, the modes of resistance system have a smaller imaginary part, and the oscillation frequency is about 80 Hz. Moreover, the damped speed is faster. From the distribution of eigenvalues of different systems, the resistance load has a less impact on the stability of the system. Additionally, the simulation by PSCAD is implemented as shown in Figure 16. From Figure 16, it can be seen that the Udc, Idc, and Un will appear high-frequency oscillation after a 4 kW CPL is connecting to network. The damped time of the oscillation mode is about 0.046 s，which is similar with the result in Figure 15. In the same condition From Figure 15, it can be seen that the number of eigenvalues of A s_CPL is more than that of A s_RES , because the CPL will bring extra modes from DC/DC converter. The both systems are stable, however, the extra modes that the CPL generated have a large imaginary part, and they are closer to real axis. Thus, it can be indicated that the variation of CPL will bring higher oscillation (about 239 Hz), and the damped time will be longer. However, the modes of resistance system have a smaller imaginary part, and the oscillation frequency is about 80 Hz. Moreover, the damped speed is faster. From the distribution of eigenvalues of different systems, the resistance load has a less impact on the stability of the system. Additionally, the simulation by PSCAD is implemented as shown in Figure 16. From Figure 16, it can be seen that the U dc , I dc , and U n will appear high-frequency oscillation after a 4 kW CPL is connecting to network. The damped time of the oscillation mode is about 0.046 s, which is similar with the result in Figure 15. In the same condition, the voltage and current change gently after the 4 kW resistance load is connecting to network. Thus, the time domain simulation verifies the results from eigenvalue analysis in Figure 15.

Discussion
In this paper, a small-signal modeling procedure of a four-terminal LVDC distribution network with distributed secondary control is presented and analyzed. The main focus is the small-signal analysis where eigenvalue analysis and major participants analysis are adopted. The benefit of modular method is the rapidity and scalability of the smallsignal modeling. In other words, the procedure of modeling is suitable for multi-terminal LVDC system. The presented procedure of modeling considers three main aspects containing individual converter stations, distributed control part and DC distribution line/load.
With the established small-signal model, the rules of system eigenvalues varying caused by the change of some parameters are carried out and analyzed. In Section 5, distributed secondary controller effect, phase locked loop effect and DC distribution line and load effect on system small-signal stability are considered. On the basis of presented results, some important influence laws are displayed in paper which indicate that PI control parameters of distributed secondary control and DC distribution line parameters such as inductance and resistance have a relative significant effect on system voltage stability through eigenvalue loci and simulations on PSCAD. Moreover, the operational limits of studied parameters are given in this study. The time-domain simulation by PSCAD is conducted for the entire LVDC system to verify the operational limits of the model, and the design rules are given for main parameters of the four-terminal LVDC distribution network, which gives a guideline for future practical parameter selection. The presented analysis only considered some important parameters, and the analysis with some other parameters varying can be studied in the future research.

Conclusions
A distributed secondary control method for VSCs and a detailed and accurate smallsignal analysis model for a four-terminal LVDC distribution network are presented in this paper. Compared with traditional control method, the adopted secondary control strategy based on average consensus algorithm and local voting protocol can promote the power distribution accuracy while increasing the average voltage. This paper focuses on the small-signal stability analysis for the multi-terminal LVDC distribution network. Firstly, various non-linear state-space models of individual entities including individual VSC

Discussion
In this paper, a small-signal modeling procedure of a four-terminal LVDC distribution network with distributed secondary control is presented and analyzed. The main focus is the small-signal analysis where eigenvalue analysis and major participants analysis are adopted. The benefit of modular method is the rapidity and scalability of the small-signal modeling. In other words, the procedure of modeling is suitable for multi-terminal LVDC system. The presented procedure of modeling considers three main aspects containing individual converter stations, distributed control part and DC distribution line/load.
With the established small-signal model, the rules of system eigenvalues varying caused by the change of some parameters are carried out and analyzed. In Section 5, distributed secondary controller effect, phase locked loop effect and DC distribution line and load effect on system small-signal stability are considered. On the basis of presented results, some important influence laws are displayed in paper which indicate that PI control parameters of distributed secondary control and DC distribution line parameters such as inductance and resistance have a relative significant effect on system voltage stability through eigenvalue loci and simulations on PSCAD. Moreover, the operational limits of studied parameters are given in this study. The time-domain simulation by PSCAD is conducted for the entire LVDC system to verify the operational limits of the model, and the design rules are given for main parameters of the four-terminal LVDC distribution network, which gives a guideline for future practical parameter selection. The presented analysis only considered some important parameters, and the analysis with some other parameters varying can be studied in the future research.

Conclusions
A distributed secondary control method for VSCs and a detailed and accurate smallsignal analysis model for a four-terminal LVDC distribution network are presented in this paper. Compared with traditional control method, the adopted secondary control strategy based on average consensus algorithm and local voting protocol can promote the power distribution accuracy while increasing the average voltage. This paper focuses on the small-signal stability analysis for the multi-terminal LVDC distribution network. Firstly, various non-linear state-space models of individual entities including individual VSC with primary control, distributed secondary control and DC distribution line/load are built. Secondly, small-signal model of the entire system is established based on individual nonlinear state-space model. The accuracy of the proposed model is further assessed through a comparison between the results from the electromagnetic simulation in PSCAD and smallsignal model in MATLAB. At last, eigenvalue analysis and major participants analysis are adopted to study the effects of main parameters on the stability of the system involving PI control, time delay and DC distribution line or load. However, the presented results are valid for LVDC distribution network only. If the system's voltage is higher than the studied one in this paper, the MMC structure should be adopted for the converter station. Then, the modeling will be more complicated and the small-signal analysis need more consideration. Furthermore, the time-domain simulation results will be supplemented and verified through laboratory simulation for the next research direction. Data Availability Statement: All data, models, and code generated or used during the study appear in the submitted article.

Conflicts of Interest:
The authors declare no conflict of interest.