Distributed Economic Dispatch Scheme for Droop-Based Autonomous DC Microgrid

In this paper, a distributed economic dispatch scheme considering power limit is proposed to minimize the total active power generation cost in a droop-based autonomous direct current (DC) microgrid. The economical dispatch of the microgrid is realized through a fully distributed hierarchical control. In the tertiary level, an incremental cost consensus-based algorithm embedded into the economical regulator is utilized to search for the optimal solution. In the secondary level, the voltage regulator estimating the average voltage of the DC microgrid is used to generate the voltage correction item and eliminate the power and voltage oscillation caused by the deviation between different items. Then, the droop controller in the primary level receives the reference values from the upper level to ensure the output power converging to the optimum while recovering the average voltage of the system. Further, the dynamic model is established and the optimal communication network topology minimizing the impact of time delay on the voltage estimation is given in this paper. Finally, a low-voltage DC microgrid simulation platform containing different types of distributed generators is built, and the effectiveness of the proposed scheme and the performance of the optimal topology are verified.


Introduction
Microgrid is a new form of grid consisting of distributed generations (DGs), energy storages, loads and power electronic devices, which is considered as a solution for flexible and reliable integration of renewable energy resources into the grid [1,2]. Compared with alternating current microgrid (AC MG), direct current microgrid (DC MG) has attracted more attentions due to the advantages that it can overcome some disadvantages of AC system, such as transformer inrush current, frequency synchronization, reactive power flow, power quality issues, etc. Therefore, DC microgrid is widely used in data centers and IT facilities, where the load is DC in nature [3,4]. Similar to AC microgrid, in order to improve the reliability of power supply, DC load is usually supplied by multiple distributed generators operated in parallel when in autonomous operation mode. However, unlike the synchronous generators in transmission network, there are different types of DGs in microgrid, such as internal combustion engines, dispatchable clear resources (fuel-cell), and non-dispatchable renewable resources (solar, wind turbines) [5]. Because of different generation costs of these distributed generators, the traditional way of power sharing based on capacity leads to higher operation cost.

1.
Decentralized economic dispatch scheme that requires local information only. This type of optimal control is less expensive because no communication is needed. However, the lack of common information (such as the common bus voltage in DC MG) leads to the challenge in utilizing all available resources in the network. Additionally, the power and voltage oscillation may happen, if every DGs regulate their output power just taking their own generation costs into account.

2.
Distributed economic dispatch scheme that requires only local computation and information exchange among some neighboring units through a sparse communication network. Based on the well-designed distributed schemes, the EDP can be solved by multiple distributed controllers working in parallel, which increases the reliability, flexibility, and scalability of the system. Therefore, the distributed manners are considered as promising options for the optimal control and management of microgrid. However, there are also many problems remained to be solved before the distributed schemes become more practical for microgrid, such as time delay, the capability of plug and play, the choice of communication topology, and so on.
Based on the previous researches, a distributed economic dispatch scheme considering power limit for droop-based autonomous DC microgrid is developed to achieve the minimal operation cost in this paper. The proposed scheme contains three control levels, namely, the droop controller works as the primary control to maintain the power balance with source/load transient and to ensure the voltage stable. The secondary control, embedded with a voltage regulator, uses the discrete consensus algorithm to estimate the average voltage of the system and recover the voltage deviation caused by primary control without oscillation. The tertiary control, embedded with an economical regulator, regulates the output power of different DGs to ensure the optimal operation of microgrid. Based on the consensus algorithm, the scheme works in a fully distributed manner without central controller and the communication only occurs between device and its neighbors. Moreover, the dynamic model of the system is established to guide the design of parameters, and the impact of time delay, which will cause the deviation of voltage estimation, is analyzed. Finally, an optimal communication topology is selected to minimize the impact of delay. Compared with the existing economic dispatch scheme, the key contributions of this paper are: 1.
The proposed distributed scheme can resolve the problems of power sharing accuracy and voltage oscillation caused by the estimation errors, which cannot be effectively solved by decentralized methods.

2.
The economical regulator updates the control signal with the newest information in each iteration step, instead of requiring the achievement of consensus. Therefore, the convergence speed of the proposed scheme is faster than the traditional consensus-based economic dispatch scheme, which is important for large-scale system. 3.
To support the flexible operation of DGs, the proposed scheme also adds the capability of plug and play into cyber layer to promote the reliability and scalability of the distributed scheme. 4.
The impact of time delay on the voltage estimation is analyzed in detail, and the capability of inhibiting oscillation caused by time delay is used to optimize the topology of the sparse communication network.
The rest of paper is organized as follows: Section 2 introduces the preliminaries, including problem formation and discrete consensus algorithm. The proposed distributed economic dispatch scheme is presented in Section 3. Section 4 optimizes the communication topology based on the analysis of time delay. The simulation case study is carried out in Section 5 to validate the performance of the scheme. Finally, Section 6 presents concluding remarks.

Problem Formation of EDP
There are different DGs in microgrid, such as micro-turbines, diesel generators, fuel cells, solar cells, wind turbines, and so on. The DGs are usually classified into three types, namely, combustion engine, schedulable clear resource and non-schedulable renewable resource. The cost functions of these distributed generators are decided by the type of DGs. For example, if the generation is a diesel generator driven by an internal combustion engine or an electric generator driven by microturbine, its generation cost function usually consists of fuel cost, maintenance cost, and emission penalty, which can be described as [5] follows: where K m , K f , K ξ are the coefficients of maintenance cost, fuel cost, and emission penalty. a, b, and c are coefficients of the fuel consumption curve of the DG, which can be achieved by fitting the experiment data. d, e, f, ξ, and ρ are coefficients of generator emission characteristics [30]. However, the emission penalty coefficient K ξ is much smaller than another two, therefore the emission penalty part is neglected in the paper, and the remaining part can be organized as a second-order quadratic function. Additionally, for clean energy source such as fuel cell, its generation cost function only consists of fuel cost and maintenance cost, thus the type of cost function is the same with the above one. Another type of DG is non-schedulable renewable resource, such as photovoltaic generation and wind generation, the cost function of this type usually contains maintenance cost and power loss cost, which can be described as [5] follows: where K m , K l are the coefficients of maintenance cost and loss cost respectively. v, u, and w are related to the characteristic of the converter, which decide the loss of renewable energy. As (2) shows, the cost function of this type can be still written to a quadratic function. Thus, in this paper, the cost functions of different DGs are modelled uniformly as follows: where P i denotes the output of DG i , C i represents the cost function for generator i, and a i , β i , γ i are cost coefficients decided by the type of DG i . Thereafter, the objective of EDP is defined as minimizing the total operation cost, i.e., where P D is the total active load demand in the MG, P loss is the transmission loss. P i min and P i max are the generation limits of DG i , respectively. The EDP formulated in (4) can be solved by the Lagrange multiplier method, and the solution to (4) can be obtained, which is the equal incremental cost principle, and it takes the following forms, Energies 2020, 13, 404 5 of 27 where ∂C i (P i )/∂P i = λ i is the incremental cost of DG i , and λ op is the optimal incremental cost of DG. The principle in (5) denotes that the total active power generation cost of the MG is minimized, if only the equal incremental cost principle is satisfied.

Discrete Consensus Algorithm
Consensus algorithms are usually applied to solve disagreements between agents with different dynamics via peer-to-peer communications. Let x i denote the state variable of node i. Node i only communicates with its neighbors, and x i can represent physical quantities such as voltage, current, power, etc. We can say the nodes have reached a consensus if and only if x i = x j for all i, j. Considering the discrete nature of digital communication, the consensus algorithm can be written in discrete form as [31], where N i represents the set of adjacent nodes, namely, a node j belongs to N i if there exists a communication link between it and node i. x i [k] is the consensus state variable of node i at kth iteration. d ij is the communication weight, where d ij > 0 if node j belongs to N i and d ij = 0, otherwise. The Equation (6) can be rewritten in matrix form as follows: where X k and X k+1 are the vectors of system state variable at the kth and (k + 1)th iterations respectively. D is the communication weight matrix of the system, which consists of d ij . Each node will converge to a common value when D is designed as a stochastic matrix and satisfies some necessary conditions in [32]. The common value will be: where X[0] is the initial state, f is a nonnegative left eigenvector associated with eigenvalue 1 of D that it satisfies the equation: f T 1 = 1, 1 = [1,1, . . . 1] T . Specially, if D is a double-stochastic matrix, then, all nodes will converge to the average value of the system, that is, In actual control systems, the convergence speed of the consensus algorithm is an important index. It is known that the convergence speed is decided by the spectral radius of D defined as follows where esr(D) is the modulus of the second largest eigenvalue of D, which is decided by the communication network topology. In order to achieve suitable convergence speed, the Metropolis method in [16] is often used to design the D, and the entries in D are calculated according to the network topology, where n i and n j are the numbers of nodes in N i and N j respectively, which are decided by the communication network topology as well. Figure 1 illustrates the proposed distributed hierarchical economic dispatch scheme for the converters in an autonomous DC microgrid, which consists of three control levels.

Primary Control
The primary control based on droop scheme is used to keep the DC bus voltage stable and share the output power properly. It consists of inner control loop and droop loop. In autonomous DC microgrid, the droop-based loop is used for coordinating the operation of the converters, which adds the virtual resistance to distribute the power among different converters. The droop scheme can be designed as follows: where, Uref and Uoi * are the system reference voltage and voltage command of converter respectively.
Ioi is the output current of converter. The value of virtual resistance Rdi determines the ratio of power sharing of different converters, which is updated by economical regulator in tertiary control level. Additionally, voltage correction term δVi is used to recover the average voltage of the islanded microgrid, and the droop character is shifted along the voltage axis by addition of δVi from the voltage regulator in secondary control level. As Figure 1 shows, the inner control loop, which consists of two PI controllers, is used to follow the voltage command Uoi * from the droop loop. The bandwidth of this level is designed to be large enough to maintain the power balance of the system while load changes.

Voltage Regulation in Secondary Control
It is known that the droop scheme in primary level will cause voltage drop inevitably. Reference [3] indicates that in order to improve the power distribution accuracy, the output voltage of each converter should be lower. Therefore, a voltage recovery process is necessary.
In this paper, all converters acquire the average voltage Uave of the DC microgrid via estimation by voltage regulators, which is used as the synchronous regulating variable to improve the dynamic response of the regulation process. Later, the voltage regulator will compare this value with the objective average voltage and the compensation value will be added to Uref to recover the average voltage of the DC microgrid, as shown in Figure 2.

Primary Control
The primary control based on droop scheme is used to keep the DC bus voltage stable and share the output power properly. It consists of inner control loop and droop loop. In autonomous DC microgrid, the droop-based loop is used for coordinating the operation of the converters, which adds the virtual resistance to distribute the power among different converters. The droop scheme can be designed as follows: where, U ref and U oi * are the system reference voltage and voltage command of converter respectively.
I oi is the output current of converter. The value of virtual resistance R di determines the ratio of power sharing of different converters, which is updated by economical regulator in tertiary control level. Additionally, voltage correction term δV i is used to recover the average voltage of the islanded microgrid, and the droop character is shifted along the voltage axis by addition of δV i from the voltage regulator in secondary control level. As Figure 1 shows, the inner control loop, which consists of two PI controllers, is used to follow the voltage command U oi * from the droop loop. The bandwidth of this level is designed to be large enough to maintain the power balance of the system while load changes.

Voltage Regulation in Secondary Control
It is known that the droop scheme in primary level will cause voltage drop inevitably. Reference [3] indicates that in order to improve the power distribution accuracy, the output voltage of each converter should be lower. Therefore, a voltage recovery process is necessary.
In this paper, all converters acquire the average voltage U ave of the DC microgrid via estimation by voltage regulators, which is used as the synchronous regulating variable to improve the dynamic response of the regulation process. Later, the voltage regulator will compare this value with the objective average voltage and the compensation value will be added to U ref to recover the average voltage of the DC microgrid, as shown in Figure 2.
The proposed voltage regulator includes a sample and hold (S/H1) module and a PI controller. The sample cycle of S/H1 is T S/H1 , which is decided by a clock installed in the regulator. Additionally, the communication period is T s , and T S/H1 = NT s . It means that S/H1 will not sample the new voltage data until the voltage consensus iteration process finishes. It is noting that the number of iteration N before converging is decided by the topology of the communication network, thus, N can be chosen offline. The voltage regulation procedure then can be described as follows.  The proposed voltage regulator includes a sample and hold (S/H1) module and a PI controller. The sample cycle of S/H1 is TS/H1, which is decided by a clock installed in the regulator. Additionally, the communication period is Ts, and TS/H1 = NTs. It means that S/H1 will not sample the new voltage data until the voltage consensus iteration process finishes. It is noting that the number of iteration N before converging is decided by the topology of the communication network, thus, N can be chosen offline. The voltage regulation procedure then can be described as follows.
When t = mNTs, m = 0, 1, 2, …, the voltage regulator samples local voltage Uoi and the output of S/H1 is Uoi[mNTs] abbreviated as Uoi [m]. Then, Uoi[m] is given to Ũi[0], which is the initial value used to estimate the average voltage of the microgrid. From then on, each converter starts to communicate with its neighbors to estimate the average voltage of the microgrid according to consensus algorithm, and its update cycle is synchronized with communication cycle. The consensus iteration process according to the message received from its neighbors can be described as follows, As aforementioned analysis, the secondary voltage control cycle is NTs. When t = (m + 1)NTs, each voltage regulator starts to sample again and repeats the consensus iteration process (13) and voltage regulation (14). Because of the existence of secondary voltage regulator, the average voltage of the DC microgrid will be recovered to the objective value. However, when some communication events occur in the information interaction period, such as time delay, (13) will not converge to the common value Uave after N iterations. This will lead to the deviation between δVi, i = 1, 2, …, n, which will cause power and voltage oscillation in turn. Thus, the impact of time delay on voltage estimation is analyzed in Section 4. When t = mNT s , m = 0, 1, 2, . . . , the voltage regulator samples local voltage U oi and the output of S/H1 is U oi [mNT s ] abbreviated as U oi [m]. Then, U oi [m] is given to Ũ i [0], which is the initial value used to estimate the average voltage of the microgrid. From then on, each converter starts to communicate with its neighbors to estimate the average voltage of the microgrid according to consensus algorithm, and its update cycle is synchronized with communication cycle. The consensus iteration process according to the message received from its neighbors can be described as follows, where d v_ij is the communication weight used in voltage regulator, Ũ i [k] and Ũ i [k + 1] are the estimation values at the k th and (k + 1)th iterations respectively. After limited iterations N, the process (13) will direct Ũ i to the average voltage of the DC microgrid. That is, each voltage regulator in secondary level will get this value U ave . It is worth noting that the number of iteration N before converging is decided by the topology of the communication network, thus, N can be chosen offline. At last, each voltage regulator compares this average voltage U ave with the objective value U ave_op , and the mismatch is passed through a PI controller to generate the compensation value δV i , as shown below: As aforementioned analysis, the secondary voltage control cycle is NT s . When t = (m + 1)NT s , each voltage regulator starts to sample again and repeats the consensus iteration process (13) and voltage regulation (14). Because of the existence of secondary voltage regulator, the average voltage of the DC microgrid will be recovered to the objective value. However, when some communication events occur in the information interaction period, such as time delay, (13) will not converge to the common value U ave after N iterations. This will lead to the deviation between δV i , i = 1, 2, . . . , n, which will cause power and voltage oscillation in turn. Thus, the impact of time delay on voltage estimation is analyzed in Section 4.

Economical Regulation in Tertiary Control
To overcome the disadvantages of centralized solution indicated by (4), incremental cost consensus algorithm is applied to develop a distributed solution to minimize the total operation cost of DC microgrid. An economical regulator is then used to provide objective virtual resistance R di * for the primary droop control, by solving EDP with the proposed algorithm, as shown in Figure 3. The economical regulator finds the optimal operational point with an iterative process, whose update cycle is controlled by a clock installed in the tertiary level. To overcome the disadvantages of centralized solution indicated by (4), incremental cost consensus algorithm is applied to develop a distributed solution to minimize the total operation cost of DC microgrid. An economical regulator is then used to provide objective virtual resistance Rdi * for the primary droop control, by solving EDP with the proposed algorithm, as shown in Figure 3. The economical regulator finds the optimal operational point with an iterative process, whose update cycle is controlled by a clock installed in the tertiary level. As shown in Figure 3, the economical regulator contains a sample and hold (S/H2) module and some mathematical modules. The sample cycle of S/H2 is TS/H2 = TS. Then, the economic dispatch procedure can be described as follows.
When t = kTS, k = 0, 1, 2, …, the economical regulator samples local voltage Uoi and output current Ioi and calculates the output power Pi [k]. According to equal incremental cost principle, the incremental cost of each DG is chosen as the consensus variable in this proposed distributed iterative approach, which can be calculated by the following equation: where λi[k] is the local incremental cost achieved in kth sampling, I = 1, 2, …, n. αi and βi are the coefficients of cost function respectively. In the next step, economical regulator i sends this message λi to its neighbors. Because of the bi-direction of the communication network, economical regulator i will also receive the incremental cost information from its neighbors, and calculate the objective incremental cost by where dp_ij is the communication weight used in economical regulator, λi * [k] is the objective incremental cost, which is moved in the direction of minimizing the total operation cost of the microgrid within the process (16). Then, the objective dispatch reference Pi * [k] in kth iteration corresponding to λi * [k] is As shown in Figure 3, the economical regulator contains a sample and hold (S/H2) module and some mathematical modules. The sample cycle of S/H2 is T S/H2 = T S . Then, the economic dispatch procedure can be described as follows.
When t = kT S , k = 0, 1, 2, . . . , the economical regulator samples local voltage U oi and output current I oi and calculates the output power P i [k]. According to equal incremental cost principle, the incremental cost of each DG is chosen as the consensus variable in this proposed distributed iterative approach, which can be calculated by the following equation: where λ i [k] is the local incremental cost achieved in kth sampling, I = 1, 2, . . . , n. α i and β i are the coefficients of cost function respectively. In the next step, economical regulator i sends this message λ i to its neighbors. Because of the bi-direction of the communication network, economical regulator i will also receive the incremental cost information from its neighbors, and calculate the objective incremental cost by where d p_ij is the communication weight used in economical regulator, λ i * [k] is the objective incremental cost, which is moved in the direction of minimizing the total operation cost of the microgrid within the process (16). Then, the objective dispatch reference P i Unlike the PQ control in primary level, the output power of droop-based converter cannot follow the power reference directly. Therefore, in order to control the power indirectly, the virtual resistance is considered. According to the droop characteristic, if the objective output power P i then the virtual resistance should be decreased, and vice versa. In this paper, the method adjusting the virtual resistance is designed according to (12) as follows: where R di * [k] is the objective virtual resistance at kth iteration. Note that the method to achieve the objective virtual resistance is not unique. The mismatch between P i [k] and P i * [k] can be also used through a PI controller to generate the objective virtual resistance [19]. Then, the virtual resistance in (12) is updated by R di * [k] through a low pass filter (LPF), as shown in Figure 1. The first-order low pass filter is used to smoothen the transient process of virtual resistance updating, whose time constant T m should be larger than the communication cycle T s . When t = (k + 1)T s , the economical regulator samples local voltage and current again and repeats the process (15)- (18). With limited circulations, the incremental costs in all economical regulators will be the same, which indicates the optimal operation of the DC microgrid. The output power and virtual resistance of each converter will converge to the optimal value respectively. The optimal values satisfy the following equation, The proposed control scheme of economical regulator overcomes the disadvantages of the existed economic droop control, finding the optimal solution of (4) without the information of line impedance and the cost functions of other distributed generators.

The Consideration of Power Limitation
In the practical low-voltage DC microgrid applications, the power supply distance is usually less than several hundreds of meters that the deviation of different bus voltages are small. Therefore, the voltage constraints are not considered in this paper, only the power constraints are taken into account. It is well-known that the output power of each generator should not be more/less than its maximum/lower generation limits, so the generator reaching its limit must keep the output power at this limit and not to participate in the EDP solving. When considering this problem, the paper changes the communication weight to satisfy the equal incremental cost principle described in (5), which can be described as follows: Supposing that converter i reaches its limitation P max , then the communication weights of this economical regulator will be changed to d p_ii = 1, d p_ij = 0, j∈N i . In this way, the incremental cost information received from its neighbors will not influence its objective incremental cost according to (16). Then, the objective dispatch reference P i * [k] keeps to P max , and the objective incremental cost keeps to λ i * = β i + 2α i P max . That is, the converter reaching its limitation will not participate in EDP solving at the next time until the total load decreases. Similarly, in order to make the remaining converters converge to their optimal value, the converter i should send an inform message to its neighbors, and its neighbors will delete converter i from their neighbor sets and change their communication weights again according to (11). Note that the communication weights will not be changed if the converter is not connected with the one reaching its limitation. Thus, the communication weights of converter i and its neighbors g∈N i can be revised as follows: where N gˆi s the new neighbor set of the converter g, and n gˆi s the numbers of converters in N gˆ. Based on the new matrix D pˆ, the remaining converters' incremental costs will converge to the same value, which indicates the optimal operation according to the equal incremental cost principle. The regulation process that the converter reaches its low limit P min is the same with that discussed above. Note that the communication weight d v_ij in voltage regulator is not changed, because the converter reaching limit can still participate in voltage regulation.

The Global Dynamic Model and Parameter Design
Small-signal methods are often used to build the dynamic model of the system and to tune the design parameters. As the voltage regulation cycle is much larger than the economical regulation cycle, the voltage compensation value δV is considered to be a constant value in the dynamic model. Thus, small-signal modelling is considered here, where each variable , where T(·): R n×1 → R n×n is a transformation that maps a vector to a diagonal matrix, namely Then, the small-signal vector of U d [k] can be written as follows: where ∆R d [k] is the small-signal vector of virtual resistance. Thus, the z transform of (12) can be written as follows:  (15) and (16) in the z domain: where ∆λ * and ∆λ are the z transform of ∆λ * [k] and ∆λ[k], respectively. M 1 [2α 1 I o1 q , 2α 2 I o2 q , . . . , 2α n I on q ] T and M 2 [2α 1 U o1 q , 2α 2 U o2 q , . . . , 2α n U on q ] T are constant vector related to the economical parameters of the generator. Additionally, according to (17), the small-signal part of the objective reference current vector ∆I o * [k] can be organized in the z domain as: . . , (λ * − β n )/2α n U on q2 ] T are constant vector related to the economical parameters of the generator. U oi q is the quiescent part of the output voltage in steady-state.
λ * is the optimal incremental cost of the system, which is the same in all economical regulators. On the other hand, the small-signal vector equation of (18) in the z domain can be written as: where ∆R d * is the z transform of the small-signal part of the objective virtual resistance vector. Q 1 [1/I o1 *q , 1/I o2 *q , . . . , 1/I on *q ] T , Q 2 [(U ref1 q + δV 1 − U o1 q )/I o1 *q2 , (U ref2 q + δV 2 − U o2 q )/I o2 *q2 , . . . , (U refn q + δV n − U on q )/I on *q2 ] T . By substituting (24)- (26) in (23) ∆U behavior of any converter with closed-loop voltage control can be expressed according to [19]: where G i c (z) is the closed-loop transfer function of the ith converter in the z domain. Without loss of generality, the electrical network containing s lines and l buses with loads is used to build the global dynamic model. The currents in lines are the state variables of the system, which can be calculated by . where The mapping matrix M bus is of size s × l, which maps the buses onto lines. For example, if kth line is located between bus i and j, the elements M bus (k,i) and M bus (k,j) will be 1 and −1, respectively. Then, the difference equation of (29) is deduced according to [33] using the sampling cycle T s and the small-signal vector equation of the difference equation in the z domain is where ∆I l is the small-signal vector of line current in the z domain, and ∆U b is the small-signal vector of bus voltage in the z domain.
According to the Kirchhoff's Current Law, the supplied current vector can be written as where I s is the vector of supplied current, the mapping matrix M NET is of size l × s, which maps the line currents onto the supplied currents. g load is the load conductance matrix, and g load = diag {[1/R load1 , 1/R load2 , . . . , 1/R loadl ]}. The small-signal portion of the conductance matrix ∆g load models any small-signal changes in the conductance matrix g load caused by load change, thus, the small-signal formulation of (31) in the z domain is where ∆G load is the z transform of ∆g load . On the other hand, the output current and voltage of the converter can be represented by I o = M CON I s , U o = M CON U b , the mapping matrix M CON is of size n × l, which maps the converter connection points onto the network buses. For example, if ith converter is connected at jth bus, the elements M CON (i,j) will be 1 and all the other elements in that row will be 0. Finally, substituting (29)-(32) into (27) provides the global dynamic model of the microgrid  (34) where β [β 1 , β 2 , . . . , β n ] T and α [α 1 , α 2 , . . . , α n ] T . 1 is a n × 1 vector whose elements are all equal to one. Given all other constant values in (33), one can draw all poles of the transfer function extracted from (33). Then, the desired asymptotically stable dynamic response for the microgrid can be designed by optimizing the converter's transfer function matrix G c and communication cycle T s and filter's time constant T m.
Note that the voltage compensation item vector δV q = [δV 1 , δV 2 , . . . , δV n ] is considered as a constant vector in this dynamic model, however, because of the performance of convergence in voltage estimation, δV q is usually not equal to the desired one, that is, δV q δV * = µ1. µ is a positive scalar. This indicates that the voltage regulation is asynchronous, namely, the voltage compensation part of each converter is not the same at every control cycle, which may cause power and voltage oscillation between converters during the regulation period, as shown in Section 5. Therefore, in order to improve the dynamic response, the performance of voltage estimation is also needed to be considered.

The Performance of Voltage Estimation and Topology Optimization
As mentioned above, the performance of voltage estimation is important for voltage regulation. The factors affecting this may include communication delay, loss of data, and the choice of communication weight, etc. In this paper, the iterative model of voltage estimation containing time delay is established, and the impact of time delay on the voltage regulation is analyzed first. Then, an optimal topology for inhibiting oscillation caused by time delay is given.

The Impact of Time Delay on Average Voltage Estimation
Taking the voltage estimation into consideration, (13) with time delay can be rewritten as follows, where τ is time delay, τ is a positive integer. Thus, the actual time delay τ act should satisfy the following equation considering the discrete nature of digital communication (35) can be transformed to matrix form as follows: where Ũ[k + 1] is the vector estimated by voltage regulator in kth iteration, D on = diag(D v ), and D on + D off = D v . To analyze the impact of time delay, a new vector Y is defined in [34], whose dimension is (37) is reformulated as where I is an identity matrix with compatible dimensions. The impact of time delay then can be described by the dynamic response of (38), which is decided by the second largest eigenvalue of Ψ.
The smaller the esr(Ψ) is, the faster the converging speed will be. Take the ring network with 5 nodes for example, the corresponding matrix D v designed by (11) is shown below, According to (39), the second largest eigenvalues of Ψ corresponding to different time delays (τ = 0-100) are drawn in Figure 4.
where Ũ[k + 1] is the vector estimated by voltage regulator in kth iteration, Don = diag(Dv), and Don + Doff = Dv. To analyze the impact of time delay, a new vector Y is defined in [34], whose dimension is where I is an identity matrix with compatible dimensions. The impact of time delay then can be described by the dynamic response of (38), which is decided by the second largest eigenvalue of Ψ.
The smaller the esr(Ψ) is, the faster the converging speed will be. Take the ring network with 5 nodes for example, the corresponding matrix Dv designed by (11) is shown below, According to (39), the second largest eigenvalues of Ψ corresponding to different time delays (τ = 0-100) are drawn in Figure 4. It can be seen from Figure 4 that the eigenvalues will be directed to (1,0) when τ increases gradually. And the esr(Ψ) increases gradually as well, which indicates a longer convergence time. However, the slow convergence speed will cause large deviation between each bus voltage estimation after N iterations, which leads to voltage or power oscillation during voltage recovery process.

Network Topology Optimization
It can be seen from the analysis above that the distribution of eigenvalues is influenced by time delay, which will further impact the voltage regulation. Therefore, an optimal topology helping to optimize the spectral radius and reduce the impact of delay is necessary. The optimal topology should have the characteristics that the corresponding iterative Equation (35) should converge fast It can be seen from Figure 4 that the eigenvalues will be directed to (1,0) when τ increases gradually. And the esr(Ψ) increases gradually as well, which indicates a longer convergence time. However, the slow convergence speed will cause large deviation between each bus voltage estimation after N iterations, which leads to voltage or power oscillation during voltage recovery process.

Network Topology Optimization
It can be seen from the analysis above that the distribution of eigenvalues is influenced by time delay, which will further impact the voltage regulation. Therefore, an optimal topology helping to optimize the spectral radius and reduce the impact of delay is necessary. The optimal topology should have the characteristics that the corresponding iterative Equation (35) should converge fast when there is no delay and the difference between esr(Ψ) and esr(D v ) should be small. Therefore, the following multi-objective optimization model is established: where, F 1 and F 2 represent the spectral radius of matrix D v and Ψ respectively, which are used to optimize the capability of inhibiting oscillation caused by time delay. F 3 represents the number of network links, which is used to optimize the construction cost of network. G f is the heterogeneous graph set with connectivity feature with n nodes. The adjacency matrix A is a n × n matrix corresponding to its topology, where the off-diagonal entry l ij ∈{0,1} while diagonal entries are zeros. The procedure of seeking for optimal network topology is as follows: 1.
The Determination of G f : The total number of network with n nodes is N G = 2 [(n − 1)n]/2 . The Warshall algorithm is used to judge the connectivity of the N G topologies, which is necessary for consensus algorithm. Then, all connective graphics are classified according to the node degree, and the heterogeneous graphics are selected to form the G f .

2.
The Solution of Multi-objective Optimization Model: Searching for the Pareto frontier of multi-objective optimization model (41) subject to the discrete domain G f determined by (1) based on non-dominated sorting method. At last, the optimal network topology will be selected from the frontier by evaluation function.

Case Study Based on the Proposed Scheme
In this section, four cases are studied to analyze the performance of the proposed strategy and the optimal topology.

Simulation Configuration
A low-voltage DC microgrid with a 500 V dc-nominal voltage is shown in Figure 5. The photovoltaic cell, Li-battery, micro-turbine, and fuel cell models are built in Simulink as well as the dynamic model of converters. In the simulation, the cost parameters of each DG are shown in Table 1 [5]. Initially, load 2, 3, 4, 5 are all 5 kW, load 1 is 10 kW. The connection resistances between DC buses are 0.4 Ω. The initial R d is 1.2 Ω, and other parameters of the primary level are shown in Table 2. The parameters of voltage regulator and economical regulator are shown in Table 3.
where, F1 and F2 represent the spectral radius of matrix Dv and Ψ respectively, which are used to optimize the capability of inhibiting oscillation caused by time delay. F3 represents the number of network links, which is used to optimize the construction cost of network. Gf is the heterogeneous graph set with connectivity feature with n nodes. The adjacency matrix A is a n × n matrix corresponding to its topology, where the off-diagonal entry lij∈{0,1} while diagonal entries are zeros. The procedure of seeking for optimal network topology is as follows: 1. The Determination of Gf: The total number of network with n nodes is NG = 2 [(n − 1)n]/2 . The Warshall algorithm is used to judge the connectivity of the NG topologies, which is necessary for consensus algorithm. Then, all connective graphics are classified according to the node degree, and the heterogeneous graphics are selected to form the Gf. 2. The Solution of Multi-objective Optimization Model: Searching for the Pareto frontier of multiobjective optimization model (41) subject to the discrete domain Gf determined by 1) based on non-dominated sorting method. At last, the optimal network topology will be selected from the frontier by evaluation function.

Case Study Based on the Proposed Scheme
In this section, four cases are studied to analyze the performance of the proposed strategy and the optimal topology.

Simulation Configuration
A low-voltage DC microgrid with a 500 V dc-nominal voltage is shown in Figure 5. The photovoltaic cell, Li-battery, micro-turbine, and fuel cell models are built in Simulink as well as the dynamic model of converters. In the simulation, the cost parameters of each DG are shown in Table  1 [5]. Initially, load 2, 3, 4, 5 are all 5 kW, load 1 is 10 kW. The connection resistances between DC buses are 0.4 Ω. The initial Rd is 1.2 Ω, and other parameters of the primary level are shown in Table  2. The parameters of voltage regulator and economical regulator are shown in Table 3. In this paper, the number of converters is 5, thus there are 19 connective graphs in Gf that can be implemented among these converters. In order to meet the awful conditions, the maximal time delay   In this paper, the number of converters is 5, thus there are 19 connective graphs in G f that can be implemented among these converters. In order to meet the awful conditions, the maximal time delay is set to 10 ms, namely, H = 5. Therefore, the Pareto frontier of multi-objective optimization model (41) can be searched by non-dominated sorting method among G f , as shown in Figure 6. is set to 10 ms, namely, H = 5. Therefore, the Pareto frontier of multi-objective optimization model (41) can be searched by non-dominated sorting method among Gf, as shown in Figure 6.    As shown in Figure 6, there are 7 topologies in the Pareto frontier, such as star topology, ring topology, and mesh topology. An evaluation function considering the construction cost and capability of reducing the impact of delay is given to select the optimal one among the 7 topologies. The difference between esr(Ψ) and esr(Dv) is used to depict the impact of delay, which is designed as: As shown in Figure 6, there are 7 topologies in the Pareto frontier, such as star topology, ring topology, and mesh topology. An evaluation function considering the construction cost and capability of reducing the impact of delay is given to select the optimal one among the 7 topologies. The difference between esr(Ψ) and esr(D v ) is used to depict the impact of delay, which is designed as: where F 2i and F 1i are the objective functions of graph i respectively, i = 1, 2, . . . , 7. Rob i is the evaluation parameter, among which the smallest one represents the smallest impact to control scheme. When considering the construction cost, the comprehensive evaluation function is given by: where Rob max and Rob min are the maximum and minimum value of Rob i respectively. F 3max and F 3min are the maximum and minimum values of F 3i . When ξ 1 = ξ 2 = 0.5, the smallest value of E i is 0.2754 and the corresponding optimal topology is shown in Figure 7. Therefore, the ring topology is deployed on the simulation system as shown in Figure 5.
where Robmax and Robmin are the maximum and minimum value of Robi respectively. F3max and F3min are the maximum and minimum values of F3i. When ξ1 = ξ2 = 0.5, the smallest value of Ei is 0.2754 and the corresponding optimal topology is shown in Figure 7. Therefore, the ring topology is deployed on the simulation system as shown in Figure 5.

Case 1: Performance of the Scheme with Changing Load
In this case, an activation process is conducted to validate the proposed scheme, as shown in Figure 8. The total load demand is 30 kW initially and the output of photovoltaic is assumed large enough to support the stable operation of the corresponding converter. Initially, the voltage regulator and economical regulator are not enabled during 0-0.15 s and the virtual resistance of each converter is constant. From Figure 8a, it can be seen that the output power of each converter is far away from its optimum, which can be calculated by Lagrange multiplier method. Additionally, the average bus voltage falls to 489.76 V caused by droop control in the primary level, as shown in Figure 8f.

Case 1: Performance of the Scheme with Changing Load
In this case, an activation process is conducted to validate the proposed scheme, as shown in Figure 8. The total load demand is 30 kW initially and the output of photovoltaic is assumed large enough to support the stable operation of the corresponding converter. Initially, the voltage regulator and economical regulator are not enabled during 0-0.15 s and the virtual resistance of each converter is constant. From Figure 8a, it can be seen that the output power of each converter is far away from its optimum, which can be calculated by Lagrange multiplier method. Additionally, the average bus voltage falls to 489.76 V caused by droop control in the primary level, as shown in Figure 8f.
Then, the both regulators are enabled at 0.15 s, and the interval time between each communication is set to 2 ms. It can be seen from Figure 8a-c that the economical regulators update virtual resistances after each iteration until all incremental costs converge to a common value λ op = 0.305 $/kWh, and the output power of each converter soon converges to its optimum P op = [10.2 kW 3.16 kW 7.0 kW 2.6 kW 8.18 kW]. The equality of incremental costs among converters is assured by the proposed scheme. In other words, the equal incremental cost principle is satisfied, and the minimum active power generation cost of the DC MG is achieved. The estimation process of average voltage of DC bus is shown in Figure 8d. Within about 11 iterations, the average voltage is achieved by each voltage regulator, which is utilized to recover the actual average voltage of the DC microgrid, as shown in Figure 8e. The average bus voltage of the system is shown in Figure 8f. where F2i and F1i are the objective functions of graph i respectively, i = 1, 2, …, 7. Robi is the evaluation parameter, among which the smallest one represents the smallest impact to control scheme. When considering the construction cost, the comprehensive evaluation function is given by: where Robmax and Robmin are the maximum and minimum value of Robi respectively. F3max and F3min are the maximum and minimum values of F3i. When ξ1 = ξ2 = 0.5, the smallest value of Ei is 0.2754 and the corresponding optimal topology is shown in Figure 7. Therefore, the ring topology is deployed on the simulation system as shown in Figure 5.

Case 1: Performance of the Scheme with Changing Load
In this case, an activation process is conducted to validate the proposed scheme, as shown in Figure 8. The total load demand is 30 kW initially and the output of photovoltaic is assumed large enough to support the stable operation of the corresponding converter. Initially, the voltage regulator and economical regulator are not enabled during 0-0.15 s and the virtual resistance of each converter is constant. From Figure 8a, it can be seen that the output power of each converter is far away from its optimum, which can be calculated by Lagrange multiplier method. Additionally, the average bus voltage falls to 489.76 V caused by droop control in the primary level, as shown in Figure 8f.  Similarly, as long as the bus voltages drop, the voltage regulators achieve the average deviation, and compensate it. The regulation process between 0.4 and 0.7 s is shown in Figure 8. At last, load 3 is added 5 kW at 0.7 s, and the total demand becomes 40 kW. The process of regulation is the same with that during 0.4-0.7 s, as shown in Figure 8. Thus, the effect of the proposed scheme is verified. Furthermore, it can be found from the simulation that the proposed distributed scheme does not degrade the stability of the system.

Case 2: The Comparison of Control Scheme in Different Levels
In this case, the primary control proposed in this paper is first compared with the improved droop control in [35] without considering the secondary level and tertiary level. The control scheme in [35] can be described as where Uref and Uoi * are the system reference voltage and voltage command of converter respectively. Uoi and UBi are the virtual control voltage and bus voltage, respectively. Rdi and Ioi are the virtual resistance and output current of converter, respectively. ε is the control parameter, which can determine the stability of the control scheme. The configuration is the same with case 1, the total demand is 30 kW at first, and 5 kW is added to load 1 and load 3 at 0.4 and 0.7 s, respectively. Figure  9a,b shows the results of the proposed droop control, Figure 9c,d shows the results of the control scheme described by (44). Similarly, as long as the bus voltages drop, the voltage regulators achieve the average deviation, and compensate it. The regulation process between 0.4 and 0.7 s is shown in Figure 8. At last, load 3 is added 5 kW at 0.7 s, and the total demand becomes 40 kW. The process of regulation is the same with that during 0.4-0.7 s, as shown in Figure 8. Thus, the effect of the proposed scheme is verified. Furthermore, it can be found from the simulation that the proposed distributed scheme does not degrade the stability of the system.

Case 2: The Comparison of Control Scheme in Different Levels
In this case, the primary control proposed in this paper is first compared with the improved droop control in [35] without considering the secondary level and tertiary level. The control scheme in [35] can be described as where U ref and U oi * are the system reference voltage and voltage command of converter respectively.
U oi and U Bi are the virtual control voltage and bus voltage, respectively. R di and I oi are the virtual resistance and output current of converter, respectively. ε is the control parameter, which can determine the stability of the control scheme. The configuration is the same with case 1, the total demand is 30 kW at first, and 5 kW is added to load 1 and load 3 at 0.4 and 0.7 s, respectively. Figure 9a,b shows the results of the proposed droop control, Figure 9c,d shows the results of the control scheme described by (44). From Figure 9a,c, it can be seen that the difference of dynamic response between the two control schemes is small. Additionally, the steady-state of the output powers are mainly the same. However, the voltage regulation process is not actually the same, as shown in Figure 9b,d. It can be seen that the voltage response of the proposed scheme is much faster than the compared one. The reason for it is that the control scheme described by (44) brings the first derivative element. This virtual control voltage reduces the response speed of the voltage regulation. However, the steady-state voltages are mainly the same. Therefore, both of the schemes can keep the DC bus voltage stable and share the output power properly. The proposed scheme has a better dynamic response while the stability of the compared scheme may be better. From Figure 9a,c, it can be seen that the difference of dynamic response between the two control schemes is small. Additionally, the steady-state of the output powers are mainly the same. However, the voltage regulation process is not actually the same, as shown in Figure 9b,d. It can be seen that the voltage response of the proposed scheme is much faster than the compared one. The reason for it is that the control scheme described by (44) brings the first derivative element. This virtual control voltage reduces the response speed of the voltage regulation. However, the steady-state voltages are mainly the same. Therefore, both of the schemes can keep the DC bus voltage stable and share the output power properly. The proposed scheme has a better dynamic response while the stability of the compared scheme may be better.
Then, the proposed voltage regulation scheme in secondary level is compared against the decentralized control. Figure 10 shows the results of the two control schemes without considering the tertiary level. It can be seen from Figure 10a,b that the bus voltages and output powers soon reach their target values after the voltage regulators are enabled at 0.15 s. Since the regulators in secondary level get the same average voltage of the system through the consensus algorithm, the voltage compensation items are all the same during the regulation period. Therefore, the voltage regulation process is synchronous and smooth. However, the decentralized control often uses the local bus voltage to generate the compensation item to recover the average voltage. The resistances between different buses lead to the deviation of different bus voltages. Accordingly, the decentralized control without communication cannot achieve the same voltage compensation item only by the local information. The deviation of different voltage compensation items causes the power and voltage oscillation in turn, as shown in Figure 10c,d. Then, the proposed voltage regulation scheme in secondary level is compared against the decentralized control. Figure 10 shows the results of the two control schemes without considering the tertiary level. It can be seen from Figure 10a,b that the bus voltages and output powers soon reach their target values after the voltage regulators are enabled at 0.15 s. Since the regulators in secondary level get the same average voltage of the system through the consensus algorithm, the voltage compensation items are all the same during the regulation period. Therefore, the voltage regulation process is synchronous and smooth. However, the decentralized control often uses the local bus voltage to generate the compensation item to recover the average voltage. The resistances between different buses lead to the deviation of different bus voltages. Accordingly, the decentralized control without communication cannot achieve the same voltage compensation item only by the local information. The deviation of different voltage compensation items causes the power and voltage oscillation in turn, as shown in Figure 10c,d.
Finally, the proposed economical regulation in tertiary level is compared against the decentralized non-linear cost-based droop scheme in [12]. The non-linear cost-based droop scheme can be described as follows where U ref and U oi * are the system reference voltage and voltage command of the converter respectively.
C i is the cost of DG i while C NL,i is the cost in no-load situation. χ is a constant decided by the capacity of the converter. From (45), it can be inferred that the costs of DGs will converge to the same value in steady state. It means that the outputs of high-cost generators are less, while that of low-cost generators are more instead of equal sharing. Figure 11 shows the results of the two control schemes. Finally, the proposed economical regulation in tertiary level is compared against the decentralized non-linear cost-based droop scheme in [12]. The non-linear cost-based droop scheme can be described as follows where Uref and Uoi * are the system reference voltage and voltage command of the converter respectively. Ci is the cost of DGi while CNL,i is the cost in no-load situation. χ is a constant decided by the capacity of the converter. From (45), it can be inferred that the costs of DGs will converge to the same value in steady state. It means that the outputs of high-cost generators are less, while that of low-cost generators are more instead of equal sharing. Figure 11 shows the results of the two control schemes.
(a) (b)  Finally, the proposed economical regulation in tertiary level is compared against the decentralized non-linear cost-based droop scheme in [12]. The non-linear cost-based droop scheme can be described as follows where Uref and Uoi * are the system reference voltage and voltage command of the converter respectively. Ci is the cost of DGi while CNL,i is the cost in no-load situation. χ is a constant decided by the capacity of the converter. From (45), it can be inferred that the costs of DGs will converge to the same value in steady state. It means that the outputs of high-cost generators are less, while that of low-cost generators are more instead of equal sharing. Figure 11 shows the results of the two control schemes.  Figure 11c,d that the incremental costs are closer to each other after the non-linear control scheme is enabled. The operation costs are 7.42 $/h, 9.06 $/h, and 10.85 $/h, respectively. It can be concluded that the non-linear cost-based droop control can really reduce the total cost of the system, however, the scheme cannot make the incremental costs consensus inherently. Therefore, the cost-based droop control cannot achieve the global optimum. Moreover, the effect of optimization based on non-linear control scheme deteriorates because of the deviation of different bus voltages. to each other after the non-linear control scheme is enabled. The operation costs are 7.42 $/h, 9.06 $/h, and 10.85 $/h, respectively. It can be concluded that the non-linear cost-based droop control can really reduce the total cost of the system, however, the scheme cannot make the incremental costs consensus inherently. Therefore, the cost-based droop control cannot achieve the global optimum. Moreover, the effect of optimization based on non-linear control scheme deteriorates because of the deviation of different bus voltages.

Case 3: The Comparison between the Proposed Scheme and the Mathematical Programming Method
In order to verify the optimal values achieved by the proposed scheme, the centralized economic dispatch algorithm based on mathematical programming according to [36] is developed, which can be formulated as follows: where P D,i is the load at bus i, and Y is the admittance matrix corresponding to the electrical network. The optimization model (46) can be solved by mathematical programming such as quadratically constraint quadratic programming. In this case, the configuration is the same with case 1, and the model is solved by commercial solver such as Gurobi or programming in the Matlab platform. The results of the centralized optimization based on model (46) and the results achieved by the proposed scheme are listed in Table 4. It can be seen from Table 4 that the results from the mathematical programming and the proposed scheme are roughly the same. The errors between these results may be brought by the loss of converters. Moreover, the incremental costs calculated by the mathematical programming are 0.298 $/kWh, 0.324 $/kWh, and 0.349 $/kWh, respectively, which are the same with those achieved by the distributed control. Therefore, the accuracy of the proposed scheme is verified.

Case 4: Performance of the Scheme Considering Generation Limitation
In this case, the power limitation is considered in the regulation process. The idea of adding generation limitation into consideration is similar to use the Lagrange multiplier method to solve the EDP, as (5) shows. That is, once one of the DG reaches its limit, the maximum output power of that DG is subtracted from the total load demand, and the EDP is solved for the remaining power demand using the rest of the generation units. The regulation process is described in Section 3.4.
The initial configuration of this case is the same as case 1. Initially, the total demand is 30 kW, the power generation limit is set to 15 kW, and regulation process is enabled at 0.15 s. It can be seen from Figure 12a-c that all the generation powers are smaller than P max and soon converge to their optimum, and the incremental costs converge to a common value before 0.4 s. After 0.4 s, a large demand 20 kW is added to the system, as Figure 12a shows. Consequently, the distributed generation PV + BA unit soon reaches its upper limit, whose generation cost is the lowest. Then, according to (17), the power reference in economical regulator of this converter remains P max , and the corresponding optimal incremental cost λ PV + BA = d(C PV + BA (P max ))/dP = 0.4 $/kWh. At the same time, the PV + BA unit sends information to its neighbors and changes its communication weight as well. The MT1 and FC1 unit communicating with PV + BA unit change their communication weights respectively, and the others remain unchanged. According to (20), the transition matrix D associating with a ring topology is changed to Dˆ, and the remaining topology is drawn by a red line, as shown in Figure 13. From Figure 12, the other converters do not reach the power limit, thus the remaining four incremental costs still converge to another common value λ op = 0.407 $/kWh, and the remaining optimal output power is P op = [6.05 kW 11.68 kW 5.2 kW 13.38 kW], as shown in Figure 12b,c respectively. Under this situation, the system achieves the lowest operation cost, additionally, the average voltage of the microgrid is recovered as shown in Figure 12d,e.      In this case, the plug and play capability of scheme is tested. Take the PV + BA unit for example, the output power of photovoltaic is large enough so that the converter can connect to the microgrid when the day is sunny, otherwise, the converter should disconnect from the microgrid when cloudy. Therefore, the capability of plug and play is important for the scheme.
The initial configuration of this case is the same as case 1. Initially, the total demand is 30 kW, and regulation process is enabled at 0.15 s. Five DGs reach the optimal states before the plug out of PV + BA unit at 0.4 s, and the regulation process is the same with case 1. After 0.4 s, the PV + BA unit plugs out of the microgrid, and the corresponding voltage regulator and economical regulator stop Next, a10 kW demand is cut off at 0.8 s, the total load falls to 40 kW, thus all the generation powers are smaller than P max and the transition matrix is changed from Dˆto D again. Through limited iteration, all generation powers soon converge to their optimum P op = [12.7 kW 4.57 kW 9.3 kW 3.86 kW 10.7 kW] with the same incremental cost λ op = 0.356 $/kWh. From Figure 12, it is illustrated that the proposed scheme considering the power limitation can minimize the total operation cost without any extra algorithm added into the regulator.

Case 5: Plug and Play Capability of the Scheme
In this case, the plug and play capability of scheme is tested. Take the PV + BA unit for example, the output power of photovoltaic is large enough so that the converter can connect to the microgrid when the day is sunny, otherwise, the converter should disconnect from the microgrid when cloudy. Therefore, the capability of plug and play is important for the scheme.
The initial configuration of this case is the same as case 1. Initially, the total demand is 30 kW, and regulation process is enabled at 0.15 s. Five DGs reach the optimal states before the plug out of PV + BA unit at 0.4 s, and the regulation process is the same with case 1. After 0.4 s, the PV + BA unit plugs out of the microgrid, and the corresponding voltage regulator and economical regulator stop working. As illustrated in Figure 14a,b, the output of PV + BA unit soon drops to zero, the remaining DGs have to produce more power to compensate for the amount of power previously generated by it. Additionally, the incremental cost of PV + BA unit drops to 0.1 $/kWh, for the incremental cost of PV + BA unit at no-load is 0.1 $/kWh. Fortunately, under the sustaining working of remaining regulators, the other four incremental costs still converge to their common value from 0.305 $/kWh to 0.372 $/kWh, and the output powers also get their optimum value P op = [5.06 kW 10.1 kW 4.3 kW 11.6 kW], as shown in Figure 14a,b. On the other hand, because of the plug out of PV + BA unit, only four units participate to estimate the average voltage of the microgrid, as shown in Figure 14c. However, the average value estimation is independent of the number of participants according to (9), therefore the average voltage is still achieved to recover the voltage of the microgrid as illustrated in Figure 14d,e.
At last, the seamless plug in of PV + BA unit into the DC MG is achieved at 0.8 s, and the other four DGs reduce their output powers to maintain the system power balance. Because of the total demand does not change, thus the optimal incremental cost and output power are the same with that before 0.4 s. Furthermore, the voltage regulator starts to participate the estimation again, and the average voltage is raised. From Figure 14, it can be seen that the plug out and plug in operations of PV + BA unit do not degrade the convergence of incremental costs and voltage estimation. In other words, the capability of the proposed distributed scheme to meet the requirement for plug and play operation is verified. demand does not change, thus the optimal incremental cost and output power are the same with that before 0.4 s. Furthermore, the voltage regulator starts to participate the estimation again, and the average voltage is raised. From Figure 14, it can be seen that the plug out and plug in operations of PV + BA unit do not degrade the convergence of incremental costs and voltage estimation. In other words, the capability of the proposed distributed scheme to meet the requirement for plug and play operation is verified.  In this case, the performance between optimal topology (ring network) and common topology is compared by considering the time delay. The initial configuration of this case is also the same as case 1. The common topology and its matrix D c are shown in Figure 15.  In this case, the performance between optimal topology (ring network) and common topology is compared by considering the time delay. The initial configuration of this case is also the same as case 1. The common topology and its matrix D c are shown in Figure 15.
First, the performance of common topology is tested under different time delay. Initially, the total demand is 30 kW, and regulation process is enabled at 0.2 s. It can be seen from Figure 16 that the incremental costs and average voltage soon converge to the target values without any oscillation when τ act = 0 s. However, when τ act = 4 ms, the equal incremental costs and average voltage achieve their stable state through a large oscillation process. On the other hand, based on (39), when τ act = 4 ms, the spectral radius of the common topology is calculated: esr(Ψ(2)) = 0.8497, the corresponding eigenvalue is −0.362 ± 0.769i. Thus, the reason for the oscillation process is that the large spectral radius leads to the voltage estimation unfinished under 11 iterations, which causes the voltage correction term in each voltage regulator different. Therefore, the voltage deviation among different converters results in power and voltage oscillation in turn, as shown in Figure 16a,b.
(e) In this case, the performance between optimal topology (ring network) and common topology is compared by considering the time delay. The initial configuration of this case is also the same as case 1. The common topology and its matrix D c are shown in Figure 15.  First, the performance of common topology is tested under different time delay. Initially, the total demand is 30 kW, and regulation process is enabled at 0.2 s. It can be seen from Figure 16 that the incremental costs and average voltage soon converge to the target values without any oscillation when τact = 0 s. However, when τact = 4 ms, the equal incremental costs and average voltage achieve their stable state through a large oscillation process. On the other hand, based on (39), when τact = 4 ms, the spectral radius of the common topology is calculated: esr(Ψ(2)) = 0.8497, the corresponding eigenvalue is −0.362 ± 0.769i. Thus, the reason for the oscillation process is that the large spectral radius leads to the voltage estimation unfinished under 11 iterations, which causes the voltage correction term in each voltage regulator different. Therefore, the voltage deviation among different converters results in power and voltage oscillation in turn, as shown in Figure 16a,b.
Similarly, the performance of optimal topology is also tested under different time delay. It can be seen from Figure 17 that the incremental costs and average voltage soon converge to the objective values without any oscillation when τact = 0 s. When τact = 4 ms, the equal incremental costs and average voltage achieve their stable state without much fluctuation. From the calculation of (39) under the optimal topology, the spectral radius of the optimal topology with the same τact is calculated: esr(Ψ(2)) = 0.7975, the corresponding eigenvalue is 0.526 ± 0.6i. The smaller spectral radius causes smaller voltage deviation among different converters. Thus, the voltage regulation and economical regulation time are less than that of the common one, and the maximum overshoot of the oscillation are smaller as well. Thus, it can be concluded that the performance of both topologies will decrease when time delay occurs. However, the spectral radius of the optimal topology is smaller than that of the common topology at the same time delay, therefore the regulation speed under optimal topology is faster. In other words, the optimal topology has a better effect in inhibiting oscillation. It is noting that the optimal ring topology is selected based on the same weight (ξ1 = ξ2 = 0.5). If the construction cost is neglected, then the weight should be: ξ1 = 1, ξ2 = 0, and the corresponding optimal topology will also change.

Conclusions
A distributed economical dispatch scheme containing three control levels is proposed in this paper for the droop-based autonomous DC microgrid. It can be seen from the simulation results that the proposed scheme can achieve the control target effectively without degrading the stability of the Similarly, the performance of optimal topology is also tested under different time delay. It can be seen from Figure 17 that the incremental costs and average voltage soon converge to the objective values without any oscillation when τ act = 0 s. When τ act = 4 ms, the equal incremental costs and average voltage achieve their stable state without much fluctuation. From the calculation of (39) under the optimal topology, the spectral radius of the optimal topology with the same τ act is calculated: esr(Ψ(2)) = 0.7975, the corresponding eigenvalue is 0.526 ± 0.6i. The smaller spectral radius causes smaller voltage deviation among different converters. Thus, the voltage regulation and economical regulation time are less than that of the common one, and the maximum overshoot of the oscillation are smaller as well.
Thus, it can be concluded that the performance of both topologies will decrease when time delay occurs. However, the spectral radius of the optimal topology is smaller than that of the common topology at the same time delay, therefore the regulation speed under optimal topology is faster. In other words, the optimal topology has a better effect in inhibiting oscillation. It is noting that the optimal ring topology is selected based on the same weight (ξ 1 = ξ 2 = 0.5). If the construction cost is neglected, then the weight should be: ξ 1 = 1, ξ 2 = 0, and the corresponding optimal topology will also change. Thus, it can be concluded that the performance of both topologies will decrease when time delay occurs. However, the spectral radius of the optimal topology is smaller than that of the common topology at the same time delay, therefore the regulation speed under optimal topology is faster. In other words, the optimal topology has a better effect in inhibiting oscillation. It is noting that the optimal ring topology is selected based on the same weight (ξ1 = ξ2 = 0.5). If the construction cost is neglected, then the weight should be: ξ1 = 1, ξ2 = 0, and the corresponding optimal topology will also change.

Conclusions
A distributed economical dispatch scheme containing three control levels is proposed in this paper for the droop-based autonomous DC microgrid. It can be seen from the simulation results that the proposed scheme can achieve the control target effectively without degrading the stability of the system. Moreover, the comparison in the primary level and secondary level show that the proposed droop control has a better dynamic response, and the distributed voltage regulation scheme can recover the average voltage of the system smoothly without voltage and power oscillation. The comparison in the tertiary level illustrates that the proposed economical regulation can achieve the global optimum without the centralized controller, and the accuracy of the optimum is verified by the mathematical programming method. Furthermore, the capability of plug and play of the proposed scheme ensures the convergence performance of the system and meets the flexible operation requirement. Finally, the performance of inhibiting oscillation of the optimal topology is verified under different time delays.

Conclusions
A distributed economical dispatch scheme containing three control levels is proposed in this paper for the droop-based autonomous DC microgrid. It can be seen from the simulation results that the proposed scheme can achieve the control target effectively without degrading the stability of the system. Moreover, the comparison in the primary level and secondary level show that the proposed droop control has a better dynamic response, and the distributed voltage regulation scheme can recover the average voltage of the system smoothly without voltage and power oscillation. The comparison in the tertiary level illustrates that the proposed economical regulation can achieve the global optimum without the centralized controller, and the accuracy of the optimum is verified by the mathematical programming method. Furthermore, the capability of plug and play of the proposed scheme ensures the convergence performance of the system and meets the flexible operation requirement. Finally, the performance of inhibiting oscillation of the optimal topology is verified under different time delays.