Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Distributed Real-time Optimal Power Flow Strategy for DC Microgrid Under Stochastic Communication Networks  PDF

  • Jian Hu
  • Hao Ma
Power Dispatching Control Center, Nanjing Power Supply Company, State Grid Corporation of China, Nanjing, 210019, China; School of Electrical Engineering, Zhejiang University, Hangzhou, China

Updated:2023-09-20

DOI:10.35833/MPCE.2021.000730

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

This paper addresses a distributed real-time optimal power flow (RTOPF) strategy for DC microgrids. In this paper, we focus on the scenarios where local information sharing is conducted in stochastic communication networks subject to random failures. Most existing real-time optimal power flow (OPF) algorithms for the DC microgrid require all controllers to work in concert with a fixed network topology to maintain a zero gap between estimated global constraint violations. Thus, the high reliability of communication is required to ensure their convergence. To address this issue, the proposed RTOPF strategy tolerates stochastic communication failures and can seek the optimum with a constant step size considering the operation limitations of the microgrid. These aspects make the strategy suitable for real-time optimization, particularly when the communication is not reliable. In addition, this strategy does not require information from non-dispatchable devices, thereby reducing the number of sensors and controllers in the system. The convergence of the proposed strategy and the optimal equilibrium points are theoretically proven. Finally, simulations of a 30-bus DC microgrid are performed to validate the effectiveness of the proposed designs.

I. Introduction

MICROGRIDS are generally categorized into AC and DC systems. The DC microgrids are more friendly to some devices such as photovoltaic (PV) panels, batteries, and electric vehicles because they are inherently DC devices. Thus, the DC microgrids can reduce the conversion process, which improves the efficiency and reliability. In addition, the DC microgrids are free of frequency regulation, reactive power control, or AC related power quality problems [

1], which results in less complex control systems. The hierarchical control is frequently utilized in DC microgrids [2], [3]. The most popular solution of the primary control is droop control [4], where load sharing is mainly determined by the droop coefficient. The secondary control is utilized to eliminate the voltage deviation and achieve perfect power sharing through centralized [5] or distributed [6], [7] methods. The tertiary control aims to achieve the optimal operation point by controlling the power flow among distribution systems [8]. Traditionally, the secondary control and tertiary control are implemented in a centralized manner [9]. However, these solutions encounter considerable challenges in large systems, including single-point failures, heavy communication burden, and slow dynamic processes [10]. Therefore, self-organizing distributed approaches, where no control center is required, have recently been used to develop scalable and robust algorithms [11]-[13].

Several algorithms combining real-time coordination and steady-state optimization are developed for both AC systems [

14], [15] and DC systems [16]-[22] to enhance the system performances. These studies are inspiring because they can achieve optimal power flow (OPF) in real-time scenarios. In [16] and [17], the OPF condition is reached by utilizing consensus algorithms to ensure the incremental costs of different converters to be identical. In [18], a sub-gradient method is combined with a consensus approach to enhance the convergence performance of the system. In [19], the OPF problem of a stand-alone DC microgrid is formulated as an exact second-order cone program and solved through the primal-dual decomposition in a distributed manner. In [21], the generation cost minimization and individual bus voltage regulation are obtained in a DC microgrid using a consensus algorithm. In [22] and [23], an event-trigger control algorithm is utilized to reduce the communication burden of the system. However, [23] focuses on the current sharing and voltage regulation, and it cannot achieve the OPF condition.

However, the aforementioned algorithms assume that the communication is ideal, which is not always true. They are vulnerable to some random failures that result in time-varying communication topologies. In addition, the convergence and accuracy of these algorithms would be affected by non-ideal communication, as they require static network topologies. To the best of our knowledge, only a few studies have been devoted to the online optimization of DC microgrids with non-ideal communication. The effects of communication delay are considered in [

24] and [25]. In [26], the economic power dispatch (EPD) problem is discussed through a dynamic communication topology. However, this requires the step size of the system to be diminished, which limits its application in real-time control. Reference [27] solves the EPD problem under uncertainties but its results are limited to static communication topologies, making it inapplicable to time-varying communication topologies. In [28], the package loss of communication is considered and analyzed. In [29], the EPD problem is solved using a gossip algorithm. However, only one communication link can transmit information at each iteration, which limits the convergence speed of the algorithm.

After a careful review of recent literature, we found that the studies that focus on online operations have not considered the communication failures [

16]-[23]. In fact, online algorithms are more sensitive to the communication failures because retransmission is complicated, and this process substantially increases the transmission time. Besides, algorithms that consider the communication failures [26]-[29] have their own limitations. First, these algorithms need to iterate on all variables in the cyber layer until the convergence is reached before the solution is applied to the power system. In general, intermediate results do not satisfy the power flow equations or operation constraints. As a result, these algorithms cannot be directly applied to real-time systems. Second, these algorithms require information from all devices including non-dispatchable devices, e.g., fixed loads, PV panels, and wind turbines, which are not always practical or economical, as the number of non-dispatchable devices is much greater than that of dispatchable ones.

The proposed algorithm is a real-time OPF strategy for DC microgrids, which can work under stochastic communication networks. It is particularly suitable for scenarios where the online optimization is required and the communication topology is time-varying. The proposed algorithm tracks the bus voltage and the current injection of a dispatchable device to calculate their references under a non-ideal communication network. It tends to be robust under uncertainties and disturbances due to the fluctuating loads and volatile renewables. In addition, the proposed algorithm can operate in communication networks with random failures and stochastic topologies, which greatly enhances its reliability in real-time scenarios. Furthermore, the communication with non-dispatchable devices, e.g., fixed loads and renewable energy resources, is not required, as this information can be obtained by local measurement using dispatchable devices. The convergence and optimality of the proposed algorithm are proven in this paper.

The main contributions of this paper are listed as follows.

1) The proposed algorithm can work in communication networks with random failures and stochastic topologies, and its utility can be enhanced when the communication is less reliable. This is a distinctive feature, which is not observed in most existing real-time algorithms.

2) The proposed algorithm can work under dynamic communications, which further enhances its reliability.

3) The proposed algorithm is fully distributed, and only the communications with neighbors are needed. In addition, the information of non-dispatchable devices such as PV or fixed loads is not required. These aspects significantly reduce the number of sensors and controllers.

The remainder of this paper is organized as follows. In Section II, we formulate the model of stochastic communication network and OPF problem. Section III introduces the proposed RTOPF strategy. In Section IV, case studies are presented to verify the accuracy, dynamic performance, robustness to communication failure and delay, and the plug-and-play features of the proposed algorithm. Finally, the conclusion and future work are given in Section V.

II. Model of Stochastic Communication Network and OPF Problem

A. Model of Stochastic Communication Network

The communications between controllers may be unreliable due to errors or failures. In this study, controllers communicate with each other through a stochastic communication network described by graph 𝒢=(𝒱,), where each vertex i𝒱 denotes a controller; each edge eij denotes a communication link with positive weight wij; N is the number of elements in 𝒱; and 𝒩i={j|eij} denotes the neighbors of controller i. In the proposed algorithm, the random communication failures are considered. At each time interval of the algorithm, the communication link eij is active or inactive with probability pij or 1-pij. The weight matrix Wk is used to define the behavior of communication at time interval k, whose element wij(k) is expressed as:

wij(k)=lij(i,j)(k)0(i,j)(k)1-ijwij(k)i=j (1)

where (k) is the set of active communication links at time interval k; and lij is the communication weight between two nodes.

In this paper, the following assumption on Wk is made.

Assumption 1: the weight matrix Wk is drawn independent identically distributed (i.i.d) from probability space such that each Wk fulfills (2) and (3).

1TWk=1TWk1=1 (2)
0<ρ(E(WkTWk-J))<1 (3)

where 1 is the vector containing all ones; and J=11T/N represents the average matrix, which ensures that all elements of a vector are equal to their average value; and ρ() and E() are the spectral radius and expected value, respectively.

Remark 1: this assumption can be easily fulfilled by appropriate tuning of lij. Equation (2) will hold if the two-way communication is utilized, which means lij=lji. Equation (3) will hold if lij=1/(max(Ni)+1), where Ni is the number of elements in 𝒩i.

B. OPF Problem Formulation

Two main objectives need to be addressed with respect to a DC microgrid. One is how to efficiently dispatch the energy from distributed generators (DGs) to meet the demands of all loads economically. The other one is to maintain desirable voltage profiles. Considering a DC microgrid consisting of N DGs, to accommodate the two objectives, a multi-objective optimization goal is set as:

minvi,iiJ=αi=1N(aiPi2+biPi+ci)+βi=1N(vi-vnom)2 (4)
Pi=viii (5)

where ai, bi, and ci are the generation cost coefficients; vi and ii are the bus voltage and current injection, respectively; Pi is the power output of each bus; vnom is the nominal voltage of the system; and the positive parameters α and β are the trade-off weight coefficients balancing the two objectives. In practice, these parameters are selected based on the specific requirements. If the cost is more important, α would be increased. The result is that more emphasis would be placed on the economic aspects of the solution. By contrast, β would be increased if the voltage deviation is more important.

Because of the existence of the bi-linear term (5), the original problem in (4) and (5) is non-convex and it cannot be efficiently solved [

30]. Notice that the voltage of the converter is strictly limited by operation limitations, and the voltage deviation can be ignored when calculating Pi [16], [21], [31]. Combining (4) and (5) by replacing vi with vnom yields a simplified convex problem as:

minvi,iiJ=αi=1N(ai'ii2+bi'ii+ci)+βi=1N(vi-vnom)2 (6)

where ai'=aivnom2 and bi'=bivnom. Based on Ohm’s law, the bus voltage and current injections are related to the conductance matrix:

IIF=GDDGDFGFDGFFVVF (7)

where V, VF and I, IF are the vectors of the bus voltages and current injections connected to dispatchable and non-dispatchable devices, respectively; and GDD, GDF, GFD, and GFF are the conductance matrices that represent the relationships between bus voltages and current outputs. Then, I is solved as:

I=(GDD-GDFGFF-1GFD)V+GDFGFF-1IF (8)

Three constraints in (9)-(11) are considered in the operation. Formulas (9) and (10) are the voltage and current constraints, and (11) represents the physical limitation of the real-time operation.

v̲viv¯        i=1,2,...,N (9)
i̲iiii¯i        i=1,2,...,N (10)
I=GV+F (11)

where v̲, v¯ and i̲i, i¯i are the lower and upper bounds of the bus voltage and current injection, respectively; and G=GDD-GDFGFF-1GFD and F=GDFGFF-1IF are the constant matrices determined by conductance matrices and non-dispatchable devices, respectively.

An assumption is made regarding the OPF problem:

Assumption 2 (Slater’s condition): V and I exist such that constraints (9)-(11) hold [

32]. In other words, the OPF problem is feasible.

III. RTOPF Strategy

The proposed algorithm has a hierarchical structure, as shown in Fig. 1.

Fig. 1  Control diagram of proposed algorithm.

The primary control level uses droop control to maintain the power balance of the system:

vi=vi,ref-mi(ii-ii,ref)    i=1,2,...,N (12)

where vi,ref and ii,ref are the reference bus voltage and current injection, respectively; and mi is the droop coefficient. No communication is required at this level. Thus, the power balance of the system can be obtained in real time. The goal of the secondary control is to solve the constrained optimization problem in a distributed manner to determine the reference values vi,ref and ii,ref. The voltage vector V and current vector I fulfill the operation relationship in (11), because they are real-time physical values measured by the controllers.

Combining (11) and (12) yields:

I=AVref+BIref+C (13)

where Vref and Iref are the vectors of reference voltage and current, respectively; and A=(E+GM)-1G, B=(E+GM)-1GM, and C=(E+GM)-1F are the constant matrices determined by system parameters, M=diag{mi}, and E is the N-order identity matrix.

Note that if ii,ref=ii, we have vi=vi,ref based on (13). Therefore, (11) can be replaced with (14).

I=Iref (14)

The dual method [

33] is applied to develop a distributed iterative approach, in which vi,ref and ii,ref are moved in the direction of minimizing the dual function. The dual function ϕ(λ) is defined as:

ϕ(λ)=infVref𝕍, Iref𝕀L(Vref,Iref,λ)=infVref𝕍,Iref𝕀J(Vref,Iref)+i=1N(λi)TAivi,ref+(Bi-Ei)ii,ref+CN (15)

where 𝕍={vi,ref|v̲vi,refv¯} and 𝕀={ii,ref|i̲iii,refi¯i} are the feasible regions of Vref and Iref to fulfill constraints (9) and (10), respectively; λi is the local estimation of dual vector related to constraint (14), and λ=[λ1,λ2,,λN]T is the matrix consist of estimations in each controller; and Ai, Bi, and Ei are the ith columns of A, B, and E, respectively. As Salter’s condition (Assumption 2) holds, strong duality is achieved [

34]. Therefore, the optimization problem in (6) is equivalent to (16).

maxϕ(λ)s.t.  λi=λj   i,j={1,2,...,N} (16)

In addition, vi,ref and ii,ref can be solved analytically as:

vi,ref=vnom-(λi)TAi2βv̲v¯ii,ref=-2αbi+(λi)T(Bi-Ei)2αaii̲iii¯ (17)

where the projection function [x]ab=max(a,min(x,b)) is defined to limit vi,ref and ii,ref within the feasible region.

Based on Lemma 1 presented in [

34], the gradient of ϕ(λ) is defined as:

ϕ(λ)=AVm+(B-E)Im+1NC1T (18)

where Vm=diag{v1,ref,v2,ref,,vN,ref} and Im=diag{i1,ref,i2,ref, ,iN,ref}. It is impossible to calculate the gradient ϕ(λ) without communication as it requires the information of the entire system. However, the current imbalance g(λ)=I-Iref can be directly measured by each converter, which provides the information gradient ϕ(λ). Combining (12) and (16) yields:

gi(λ)=ii-ii,ref=j=1Nϕ(λ)λij (19)

where gi(λ) and λij are the ith elements of g(λ) and λj, respectively.

The pseudocode of the OPF algorithm is shown in Algorithm 1. In the algorithm, we define yi as the local estimation of g for controller i; y=[y1,y2,,yn]T is the status matrix of the entire system. In addition, yji is the jth element of vector yi; si and qi are the temporary variable vectors exchanged between converters; and γ is the step size of the system.

Algorithm 1  : secondary control level of OPF algorithm

Step 1: initialization. For each converter, set the iteration k=0, λi(0)=0;     initialize vi,ref(0)=vnom and ii,ref(0)=-bi/(2ai); measure the local current   imbalance gi(λ(0))=ii(0)-ii,ref(0); and set yii(0)=Ngi(λ(0)), while other   elements of yi are initialized to be 0.

Step 2: local optimization of λi. For controller i, compute:

si(k)=λi(k)-γyi(k)λi(k+1)=si(k)+jNiwij(k)(sj(k)-si(k))

Step 3: primary problem solution. Calculate primary variables vi,ref(k+1)     and ii,ref(k+1) by:

vi,ref(k+1)=vnom-(λi(k))TAi2βv̲v¯ii,ref(k+1)=-bi+(λi(k))T(Bi-Ei)2αaii̲ii¯i

Step 4: dynamic average consensus of yi. Measure the current imbalance    gi(λ(k+1))=ii(k+1)-ii,ref(k+1) under new references vi,ref(k+1) and     ii,ref(k+1). Then, update yi(k+1) with

qii(k)=yii(k)+Ngi(λ(k+1))-gi(λ(k))qji(k)=yji(k)jiyi(k+1)=qi(k)+jNiwij(k)(qj(k)-qi(k))

Step 5: set k: k+1 and go to Step 2.

In the initial stage, the initial values of the local variables λi, vi,ref, and ii,ref are constants. The current imbalance ii-ii,ref can be measured using a local controller, which is used to set yi. In this paper, the iteration k is only used to distinguish the current and next iterations. In real-time operations, this variable is not recorded. The controller only requires the current status of the local variables, and the historical variables would not be recorded. Any DG connected to the microgrid initializes its local variables based on constants and local measurements. It can then be introduced into the algorithm. This implies that the plug-and-play capacity is not affected by the initialization stage of the proposed algorithm.

Remark 2: in the proposed algorithm, only the generator cost is considered. If other power management algorithms exist such as energy storage management algorithms, the controlled devices are considered non-dispatchable ones. This means that the proposed algorithm can work in parallel with other real-time power management algorithms.

Remark 3: in the proposed algorithm, a quadratic cost function is applied. If another convex function is applied as a cost function, Step 3 would become an optimal problem to solve vi,ref and ii,ref based on local variables. This problem can be solved using a local controller. However, a non-convex cost function cannot be applied because it makes the OPF problem non-convex.

The convergence and optimality of the proposed algorithm are then given in Theorem 1.

Theorem 1: assuming that Assumptions 1 and 2 hold, consider the sequences λ(k) and y(k) generated by the proposed algorithm. Let the vector λ¯(k)=Jλ(k) be an average matrix and λ˜(k)=λ(k)-λ¯(k) is the corresponding disagreement matrix. Then, a positive number γ* exists such that if γγ*, we have limkλ˜(k)=0 and limk ϕ(λ(k))=ϕ*, where ϕ* is the solution to (16).

The proof of Theorem 1 is inspired by [

35]. In [35], an unconstrained distributed optimization method is provided for stochastic communication networks. Thus, the original OPF problem is converted to a dual problem, which aims to solve the dual variable λ. A major obstacle is that ϕ(λ) is coupled between different agents, which differs from the assumption in [35]. To solve this problem, we use g(λ), which can be directly measured by each agent, to obtain the information about ϕ(λ).

The detailed proof is given in Appendix A.

IV. Case Study

A. Test System

To validate the effectiveness of the proposed algorithm, a DC microgrid utilizing the skeleton of the IEEE 30-bus test feeder system is set, as shown in Fig. 2. We choose IEEE 30-bus test system as the skeleton because it is a meshed grid that can verify our algorithm in a microgrid with meshed structure. The nominal voltage of the DC microgrid is set to be 1000 V within the ±5% deviation range. The line resistance values are listed in Table I.

Fig. 2  Test system for proposed algorithm.

Table I  Line Resistance Values
Line No.Resistance (Ω)Line No.Resistance (Ω)Line No.Resistance (Ω)
1-2 0.06 1-3 0.19 2-4 0.17
3-4 0.04 2-5 0.02 2-6 0.18
4-6 0.04 5-7 0.12 6-7 0.08
6-8 0.04 6-9 0.21 6-10 0.56
9-11 0.21 9-10 0.11 4-12 0.26
12-13 0.14 12-14 0.26 12-15 0.13
12-16 0.20 14-15 0.20 16-17 0.19
15-18 0.22 18-19 0.13 19-20 0.07
10-20 0.21 10-17 0.08 10-21 0.07
10-22 0.15 21-22 0.02 15-23 0.20
22-24 0.18 23-24 0.27 24-25 0.33
25-26 0.38 25-27 0.21 28-27 0.40
27-29 0.42 27-30 0.60 29-30 0.45
8-28 0.20 6-28 0.06

The interval time of the secondary control is set to be 0.2 s. In this paper, each communication link is assumed to be subject to random failure following a certain Bernoulli process. In other words, in each iteration, each communication link will be activated with a probability of p or deactivated with a probability of 1-p. Thus, when p=1, the simulation is retrograded to a fixed scenario. In this paper, p is set to be 0.5 in most scenarios except for a special statement.

Nine DGs are modeled as Capstone microturbines, which are connected to DC mircrogird via power electronic converters. The proposed controller is also integrated in converters. The economic parameters of DGs are shown in Table II. Three PV plants are set to generate 100 kW power. We set the droop coefficient m=1 and the step size γ=0.03. The trade-off parameters are set as α=1 and β=0.01.

Table II  Economic Parameters of DGs
ModelDG No.Economic parameter

Imax

(A)

Imin

(A)

ai

(¢/kW2h)

bi

(¢/kWh)

ci

(¢/h)

Capstone 330

(high pressure)

DG1 0.0248 2.366 21.43 28 5

Capstone 330

(liquid fuel)

DG2, DG5 0.0680 1.730 21.46 26 5
Capstone C65 DG3, DG7, DG8 0.0045 3.253 29.51 65 5
Capstone C200 DG4, DG6, DG9 0.0019 2.232 82.33 200 5

B. Accuracy Validation

To validate the accuracy of the proposed algorithm, the CVX tool [

36] in MATLAB is used to solve the original problem given in (5). As shown in Table III, the optimal solutions of problem (5) solved by CVX, the steady-state bus voltages and output currents controlled by proposed algorithm, and the relative errors between the optimal solutions and the steady-state values are listed. It can be derived that relative error is less than 0.6%, which validates the accuracy of the proposed algorithm.

Table III  Accuracy Comparison
DG No.CVXProposed algorithmRelative error (%)
Voltage (V)Current (A)Voltage (V)Current (A)VoltageCurrent
DG1 992.1 9.96 992.1 9.99 0 0.305
DG2 995.9 8.18 995.9 8.19 0 0.065
DG3 995.3 5.00 995.3 5.00 0 0.000
DG4 1001.3 159.54 1001.3 159.52 0 -0.016
DG5 989.4 8.52 989.4 8.51 0 0.077
DG6 1003.7 166.62 1003.7 166.62 0 -0.002
DG7 983.4 5.00 983.4 5.03 0 0.577
DG8 996.6 5.00 996.6 5.00 0 0.090
DG9 1042.1 32.16 1042.1 32.17 0 0.041

C. Dynamic Process

As in a time-varying environment, the available power of PVs fluctuates rapidly over time [

37]-[39]. In this subsection, the impact of the load change and PV fluctuation are analyzed. In this scenario, the power of PV1 drops by 80% of the output power at a speed of 8 kW/s at t=10 s. The power of PV1 recovers at the same speed at t=60 s. These fluctuations resemble those reported in [37]-[39]. In addition, a new load L=70 kW connected to bus 3 begins to drain the power at t=40 s. Dynamic output currents of DG1-DG4 and dynamic bus voltages of DG7 and DG9 are shown in Fig. 3.

Fig. 3  Dynamic output currents of DG1-DG4 and dynamic bus voltages of DG7 and DG9. (a) Dynamic output currents of DG1-DG3. (b) Dynamic output current of DG4. (c) Bus voltages of DG7 and DG9.

DG1-DG4 are four different types of generators for the output current; DG7 and DG9 have the maximum and minimum bus voltages, respectively. After the PV power drops, all DGs begin to increase the output power to maintain the power balance. Simultaneously, the secondary control also adjusts references ii,ref and vi,ref to find the new optimal working point of the system. Finally, both the current injection and bus voltage converge to the optimal solution shortly after the PV power stops decreasing. Similarly, when the PV power recovers, the DGs decrease their power outputs to reach the new optimum. After the loads are connected to the microgrid, the output current ii of each converter immediately increases to satisfy the power balance derived from the primary control. The secondary control then begins to adjust the references ii,ref and vi,ref to find the new optimal working point of the microgrid. Finally, both the current injection and bus voltage converge to an optimal solution within 20 s. It can also be observed that the ranges of the current and voltage are recovered after the transient process. In summary, the proposed algorithm can drive the microgrid to a new optimal state under fluctuating PV and load power.

D. System Performance Under Different Communication Failure Ratios

Next, the dynamic current injection of DG1 using the proposed algorithm under different p is examined in Fig. 4.

Fig. 4  Dynamic current injection of DG1 under different p.

In this scenario, a new load L=70  kW connected to bus 3 starts to drain the power at t=4 s. Simultaneously, the power of PV4 drops to 70 kW. When p is close to 1, e.g., p=0.9, the curve is smooth and resembles the dynamics of the fixed topology case. When p is low, e.g., p=0.1, which means that only 10% of communication is successful, the oscillation occurs and the convergence is slower. However, the algorithm still converges to the optimal solution. The ability of the proposed algorithm to withstand communication failure is thus validated. To investigate the effect of p on the convergence time, a Monte Carlo test is performed. The simulation contains 10 sets with even distributed p from 0.1 to 1. Each set contains 500 cases. A box plot of Monte Carlo test of convergence time is presented in Fig. 5, where the middle-horizon red line represents the median of each set, the blue box represents the middle 50% of all cases, and the vertical blue lines represent the ranges of all cases.

Fig. 5  Box plot of Monte Carlo test of convergence time with different p.

If p=0.1, the convergence is slow. The median of the convergence time is approximately 110 s. In addition, the convergence time is distributed in a long range from 50 s to 200 s, which means that the randomness has a tremendous effect on the convergence. When p increases, the convergence becomes faster and the distribution range becomes narrower. This indicates that with a larger p, the randomness will have a smaller effect on the convergence. When p0.5, the convergence time stops decreasing and nearly all cases converge at the same time. It could be observed that under these cases, the random communication failures have little effect on the convergence and the system behavior is similar to one without random communication failure.

E. System Performance Under Communication Delay

To verify the performance of the proposed algorithm under a communication delay, the dynamic response of the proposed algorithm under random communication is performed. The delay time is set randomly between 0 s and 2 s, and the failure ratio of the communication remains 50%. The dynamic current injection of DG1 under a random communication delay of 0-2 s is illustrated in Fig. 6. It could be concluded that the algorithm converges to an optimal solution under a random delay. However, the oscillation occurs and the algorithm converges more slowly. Again, it could be concluded that the algorithm could operate under a communication delay. However, this delay will degrade the performance of the algorithm.

Fig. 6  Dynamic current injection of DG1 under a random communication delay of 0-2 s.

F. Plug-and-play Capacity

The plug-and-play capacity of the proposed algorithm is illustrated in Fig. 7.

Fig. 7  Plug-and-play capacity of proposed algorithm. (a) Dynamic output currents of DG1-DG3 and DG9. (b) Dynamic output current of DG4.

The dynamic responses of DG1-DG4 and DG9 are provided. When DG9 is disconnected from the microgrid at t=5 s, other DGs increase their outputs to compensate the power imbalance derived from the disconnection. The secondary control then navigates the microgrid to a new optimal point in 10 s. At t=25 s, DG9 is recovered. The outputs of other DGs decrease. The secondary control starts to lower the reference current to recover the original optimal point. Finally, the system returns to its original working point at t=45 s.

G. Comparison of Original and Simplified Models

Finally, to justify the rationality of the simplification of the OPF problem as described in Section II, a Monte Carlo test containing 4000 cases is performed to compare the solutions of the OPF problems (3)-(5). In each case, the parameters of the system are selected with a normal distribution with expectation μ and variance σ. The system setup for Monte Carlo simulation is presented in Table IV.

Table IV  System Setup for Monte Carlo Simulation
ParameterμσParameterμσ
ai 0.03 0.01 Imax (A) 100 33
bi 3.00 1.00 Imin (A) 5 1
ci 50.00 25.00 Load (kW) 100 33

First, all random cases are solved based on the simplified model (5) to obtain a sub-optimal solution. Then, the global optimization method provided by MATLAB, which uses the scatter-search mechanism, is used to find the optimal solution of (5). The relative error Er between the two models is calculated by:

Er=Jsub-JoptJopt (20)

where Jopt and Jsub are the values of the objective function based on the original model (3) and simplified model (5), respectively.

A histogram of relative error Er is shown in Fig. 8. It could be concluded that the relative error is less than 1% for all cases. In addition, in most cases, the relative error is less than 0.3%. Based on the aforementioned discussion, we could conclude that the simplification described in Section II is reasonable.

Fig. 8  Histogram of relative error.

V. Conclusion and Future Work

In this paper, a real-time OPF algorithm for DC microgrids is proposed. We prove that our algorithm converges to the optimal solution, even under a stochastic communication network. This aspect significantly enhances the reliability of the proposed algorithm. Moreover, our algorithm can obtain information from non-dispatchable devices by local measurement of dispatchable devices, which significantly reduces the need for controllers and communication lines. To validate our algorithm, simulations on a IEEE 30-bus DC microrid are adopted including the accuracy, dynamic performance, and plug-and-play capacity.

Our future research will attempt to extend the algorithm to the following scenarios.

1) The proposed algorithm will be extended for OPF control for AC microgrids.

2) Because the proposed algorithm is not suitable for microgrids with varying electrical topologies, developing an algorithm that can solve the OPF problem under these conditions is an interesting topic for future work. The proposed algorithm is validated under random communication delay, but analytical analysis of the performance of the proposed algorithm under communication delay remains an open problem.

Appendix

Appendix A

A. Preliminary

In our proof, we use ||X||F, ||X||1, and ||X||2 for the Frobenius norm, 1-norm, and 2-norm of matrix X, respectively, and A,BF for the Frobenius inner product of two matrices. We define the average vector X¯=JX, where all its elements are equal to the average of the elements of vector X. Several important lemmas are presented before the optimum and convergence of the proposed algorithm are discussed.

Lemma 1: an equation is given as A,B¯F=A¯,B¯F.

Proof 1: based on the definition of the Frobenius inner product, we have:

A,B¯F=tr(ATJB)=tr((JA)T(JB))=A¯,B¯F (A1)

where the third equality is derived from J=J2.

Lemma 2: let λ˜(k)=λ(k)-λ¯(k) and y˜(k)=y(k)-y¯(k) be the corresponding disagreement matrices. Let X(k)=Ei=0k||λ˜(k)||2, Y(k)=Ei=0k||y˜(k)||2, and Z(k)=Ei=0k||y¯(k)||2 be expected energy from 0 to k. Then, we have:

X(k)s1p21-s1s2Z(k)+q1+s1q21-s1s2 (A2)
Y(k)p21-s1s2kZ(k)+q2+s1q11-s1s2 (A3)

where s1=2γη1-η; q1=2E||λ˜(0)||21-η2; s2=2NL(1+η)1-η; p2=2NL(1+η)γ1-η; and q2=2E||y˜(0)||21-η2.

Proof 2: since WkJ=JWk=J, we could obtain:

λ˜(k+1)=Ak(λ˜(k)-γy˜(k)) (A4)

where Ak=Wk-J. Then, based on the inequality of the 2-norm [

40], we have:

E||λ˜(k+1)||2=E||Ak(λ˜(k)-γy˜(k))||2E||Akλ˜(k)||2+E||γAky˜(k)||2ηE||λ˜(k)||2+γηE||y˜(k)||2 (A5)

According to the convergence property of energy of random variables [

35], let v(k)=E||λ˜(k)||2 and w(k)=γηE||y˜(k)||2, and then we have:

X(k)s1Y(k)+q1 (A6)

Similarly, we can yield:

y˜(k+1)=Aky˜(k)+AkΔgm(λ(k))) (A7)

Then, taking the norm expectation of both sides yields:

E||y˜(k+1)||2E||Aky˜(k)||2+E||AkΔgm(λ(k))||2ηE||y˜(k)||2+E|maxg(λ(k+1))-g(λ(k)|ηE||y˜(k)||2+E||ϕ(λ(k+1))-ϕ(λ(k))||1ηE||y˜(k)||2+NE||ϕ(λ(k+1))-ϕ(λ(k))||FηE||y˜(k)||2+NLE||λ˜(k+1)-λ˜(k)-γy¯(k)||FηE||y˜(k)||2+NLE||λ˜(k+1)||2+E||λ˜(k)||2+γE||y¯(k)||2 (A8)

where Δgm(λ(k))=Ndiag{g(λ(k+1))-g(λ(k))}; ① is derived from the definition of Δgm(λ(k)) and the fact that E||Ak||2<1; ② is derived in the next subsection; and ③-⑤ are derived from the inequality between matrix norms [

40]. Combining (A5) and (A8) yields:

E||y˜(k+1)||2ηE||y˜(k)||2+NL(1+η)E||y˜(k)||2+γE||y˜(k)||2 (A9)

Utilizing the convergence of the energy of random variables [

35] and letting v(k)=E||y˜(k)||2 and w(k)=NL(1+η)E||λ˜(k)||2+γ||y˜(k)||2, we obtain:

Y(k)s2X(k)+p2Z(k)+q2 (A10)

Combining (A6) and (A10) completes the proof.

Lemma 3: convergence of random sequence. Let (Ω,) be a probability space and 01...k be a sequence of σ subfields of . Let v(k),a(k), and w(k) be non-negative random variables, the following relationship then holds for k0:

Ev(k+1)(1+a(k))v(k)-u(k)+w(k) (A11)

where k=0a(k) and k=0w(k). Then, v(k) converges to some random variables v, and we further have k=0u(k).

B. Proof of Convergence Theorem

Proof 3: by properly intertwining the optimization and dynamic consensus steps, we can rewrite the proposed algorithm in a compact form as:

λ(k+1)=Wk(λ(k)-γy(k)) (A12)
y(k+1)=Wk(y(k)+Δgm(λ(k))) (A13)

To investigate the convergence properties of the proposed algorithm, we first consider the following auxiliary sequence that runs analogously to (A12) and (A13) in the average space:

λ¯(k+1)=λ¯(k)-γy¯(k) (A14)
y¯(k+1)=y¯(k)+Δg¯m(λ(k)) (A15)

Function ϕ(λ) is a concave second-order Lipschitz continuity with constant L. In other words, for λ1,λ2n×n, we have:

ϕ(λ2)ϕ(λ1)+ϕ(λ1),λ2-λ1F+L2||λ2-λ1||F2 (A16)

In addition, based on (A16) and the conservation property of the averaging matrix [

42], we immediately obtain:

ϕ¯(λ(k))=Jϕ(λ(k))=g¯m(λ(k))=y¯(k) (A17)

Taking conditional expectation on k and plugging λ2=λ¯(k+1) and λ1=λ¯(k) into (A16) yields:

E(ϕ(λ¯(k+1)))ϕ(λ¯(k))-γϕ(λ¯(k)),y¯(k)F+Lγ22||y¯(k)||F2=ϕ(λ¯(k))-γϕ¯(λ(k),y¯(k)F+Lγ22||y¯(k)||F2+γϕ¯(λ(k)-ϕ(λ¯(k))),y¯(k)F=ϕ(λ¯(k))-γ||y¯(k)||F2+Lγ22||y¯(k)||F2+γϕ¯(λ(k)-ϕ(λ¯(k))),y¯(k)F (A18)

where ① and ③ are derived from (A16) and (A17), respectively. In addition, based on the Lipschitz continuity of ϕ(λ), we have:

ϕ¯(λ(k)-ϕ(λ¯(k))),y¯(k)F=ϕ¯(λ(k),y¯(k)F-ϕ(λ¯(k)),y¯(k)F=ϕ(λ(k)),y¯(k)F-ϕ(λ¯(k)),y¯(k)F=ϕ(λ(k))-ϕ(λ¯(k)),y¯(k)F||ϕ(λ(k))-ϕ(λ¯(k))||F||y¯(k)||FL||λ˜(k)||F||y¯(k)||F (A19)

Combining (A18) and (A19) yields:

E(ϕ(λ¯(k+1)))ϕ(λ¯(k))-γ-Lγ22||y¯(k)||F2+γL||λ˜(k)||F||y¯(k)||F (A20)

Let v(k)=ϕ(λ¯(k)) and we have:

E(v(k+1))v(k)-a'||y¯(k)||22+γNL||λ˜(k)||2||y¯(k)||2 (A21)

where a'=Nγ21/γ-L/2.

Taking the total expectation and sum (A21) from 0 to k, we obtain:

E(v(k+1))E(v(0))-a'Ei=0k||y¯(i)||22+γNLEi=0k||λ˜(i)||2||y¯(i)||2E(v(0))-a'Z2(k)+γNLX(k)Z(k) (A22)

Then, the last inequality is derived from the Cauchy-Schwarz inequality.

Invoking Lemma 2, we obtain:

X(k)Z(k)b1Z2(k)+b2Z(k) (A23)

where b1=s1p21-s1s2 and b2=q1+s1q21-s1s2.

Combining (A21) and (A22) yields:

E(v(k+1))E(v(0))-a0Z2(k)+b0Z(k) (A24)

where a0=a'-γNLb1 and b0=γNLb2 are the constants that depend on η and γ. Because E(v(k))0, we have:

-a0Z(k)2+b0Z(k)+E(v(0))0 (A25)

Based on the definition of a0, we can derive that a0>0 when the step size γ is sufficiently small. Because b0>0 and E(v(0))>0, it follows from (A25) that

limkZ(k)< (A26)

Thus, applying Markov’s inequality [

43] to any ϵ>0, we obtain:

i=0P(||y¯(i)||2>ϵ)i=0E||y¯(i)||2ϵ2< (A27)

By the Borel-Cantelli Lemma in [

43] and Proposition 1.2 in [44], we have limk||y¯k||2=0. From Lemma 2, we can derive:

X(k)s1p2+p11-s1s2Z(k)+q1+s1q21-s1s2< (A28)

It shows that limk||λ˜(k)||2=0. Likewise, using Lemma 2, we can deduce that Y(k) is bounded and limk||y˜(k)||2=0. Then, using (a+b)22a2+2b2, we can rewrite (A21) as:

E(v(k+1))E(v(k))-u(k)+w(k) (A29)

where u(k)=a'||y¯(k)||22 and w(k)=2NL(||λ˜(k)||2+||y¯(k)||22). Based on the previous discussion, we have:

k=0w(k)2NLlimkX2(k)+limkZ2(k) (A30)

Then, when Lemma 3 is applied, the sequence v(k) converges to some random variables. Because ϕ(λ) is radically unbounded, λ¯(k) is also most likely bounded.

Based on the inequality between Frobenius norm and 2-norm, we have limk||y¯k||F=0, limk||y˜k||F=0, and limk||λ˜k||F=0.

Because ϕ(λ) is concave and second-order Lipschitz continuity with constant L, we obtain:

ϕ(λ*)-ϕ(λ¯(k))ϕ(λ¯(k)),λ*-λ¯(k)F=ϕ(λ(k)),λ*-λ¯(k)F+ϕ(λ¯(k))-ϕ(λ(k)),λ*-λ¯(k)F=Δg¯m(k),λ*-λ¯(k)F+ϕ(λ¯(k))-ϕ(λ(k)),λ*-λ¯(k)F||y¯(k)||F||λ*-λ¯(k)||F+L||λ˜(k)||F||λ*-λ¯(k)||F (A31)

where λ* is the optimum of problem (16). Because λ¯(k) is bounded, we can derive that ||λ*-λ¯(k)||FM, where M is a constant. Thus, from (A31), and recalling the Lipschitz continuity of ϕ(λ), we can claim that

ϕ(λ*)-limkϕ(λ¯(k))Mlimk(||y¯(k)||F+L||λ˜(k)||F)=0 (A32)

Based on the definition of λ*, we can obtain ϕ(λ*)-ϕ(λ¯(k))0 . Combining it with (A32) yields:

limkϕ(λ¯(k))=ϕ(λ*) (A33)

Moreover, based on Cauchy’s mean value theorem, we have:

ϕ(λ(k))=ϕ(λ¯(k))+ϕ(λ(k)+ξλ˜(k)),λ˜(k)F (A34)

where 0ξ1. Because ||λ(k)+ξλ˜(k)||F||λ(k)||F+ξ||λ˜(k)||F and ||ϕ(λ(k)+ξλ˜(k))||F are both bounded, we have:

limk|ϕ(λ(k))-ϕ(λ¯(k))|limk||ϕ(λ(k)+ξλ˜(k))||F||λ˜(k))||F=0 (A35)

Combining (A33) and (A35) yields limkϕ(λ(k))=ϕ(λ*), which completes the proof.

References

1

T. Dragičević, X. Lu, J. C. Vasquez et al., “DC microgrids–Part I: a review of control strategies and stabilization techniques,” IEEE Transactions on Power Electronics, vol. 31, no. 7, pp. 4876-4891, Jul. 2016. [Baidu Scholar] 

2

F. Gao, R. Kang, J. Cao et al., “Primary and secondary control in DC microgrids: a review,” Journal of Modern Power Systems and Clean Energy, vol. 7, no. 2, pp. 227-242, Mar. 2019. [Baidu Scholar] 

3

Q. Shafiee, T. Dragičević, J. C. Vasquez et al., “Hierarchical control for multiple DC-microgrids clusters,” IEEE Transactions on Energy Conversion, vol. 29, no. 4, pp. 922-933, Dec. 2014. [Baidu Scholar] 

4

A. Maknouninejad, Z. Qu, F. L. Lewis et al., “Optimal, nonlinear, and distributed designs of droop controls for DC microgrids,” IEEE Transactions on Smart Grid, vol. 5, no. 5, pp. 2508-2516, Sept. 2014. [Baidu Scholar] 

5

J. M. Guerrero, J. C. Vasquez, J. Matas et al., “Hierarchical control of droop-controlled AC and DC microgrids—a general approach toward standardization,” IEEE Transactions on Industrial Electronics, vol. 58, no. 1, pp. 158-172, Aug. 2011. [Baidu Scholar] 

6

M. Baharizadeh, M. S. Golsorkhi, M. Shahparasti et al., “A two-layer control scheme based on p-droop characteristic for accurate power sharing and voltage regulation in DC microgrids,” IEEE Transactions on Smart Grid, vol. 12 no. 4, pp. 2776-2787, Jul. 2021. [Baidu Scholar] 

7

Z. Li, Z. Cheng, J. Si et al., “Distributed event-triggered secondary control for average bus voltage regulation and proportional load sharing of DC microgrid,” Journal of Modern Power Systems and Clean Energy, vol. 10 no. 2, pp. 678-688, May 2022. [Baidu Scholar] 

8

P. Wang, W. Wang, N. Meng et al., “Multi-objective energy management system for DC microgrids based on the maximum membership degree principle,” Journal of Modern Power Systems and Clean Energy, vol. 6, no. 4, pp. 668-678, Jul. 2018. [Baidu Scholar] 

9

J. Xiao, P. Wang, and L. Setyawan, “Hierarchical control of hybrid energy storage system in DC microgrids,” IEEE Transactions on Industrial Electronics, vol. 62, no. 8, pp. 4915-4924, Aug. 2015. [Baidu Scholar] 

10

M. Yazdanian and A. Mehrizi-Sani, “Distributed control techniques in microgrids,” IEEE Transactions on Smart Grid, vol. 5, no. 6, pp. 2901-2909, Nov. 2014. [Baidu Scholar] 

11

Z. Zhang and M. -Y. Chow, “Convergence analysis of the incremental cost consensus algorithm under different communication network topologies in a smart grid,” IEEE Transactions on Power Systems, vol. 27, no. 4, pp. 1761-1768, Nov. 2012. [Baidu Scholar] 

12

S. Yang, S. Tan, and J. Xu, “Consensus based approach for economic dispatch problem in a smart grid,” IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 4416-4426, Nov. 2013. [Baidu Scholar] 

13

W. Zhang, Y. Xu, W. Liu et al., “Distributed online optimal energy management for smart grids,” IEEE Transactions on Industrial Informatics, vol. 11, no. 3, pp. 717-727, Jun. 2015. [Baidu Scholar] 

14

F. Dörfler, J. W. Simpson-Porco, and F. Bullo, “Breaking the hierarchy: Distributed control and economic optimality in microgrids,” IEEE Transactions on Control of Network Systems, vol. 3, no. 3, pp. 241-253, Sept. 2016. [Baidu Scholar] 

15

G. Chen, J. Ren, and E. N. Feng, “Distributed finite-time economic dispatch of a network of energy resources,” IEEE Transactions on Smart Grid, vol. 8, no. 2, pp. 822-832, Mar. 2017. [Baidu Scholar] 

16

J. Hu, J. Duan, H. Ma et al., “Distributed adaptive droop control for optimal power dispatch in DC microgrid,” IEEE Transactions on Industrial Electronics, vol. 65, no. 1, pp. 778-789, Jan. 2018. [Baidu Scholar] 

17

Y. Li, Z. Zhang, T. Dragičević et al., “A unified distributed cooperative control of DC microgrids using consensus protocol,” IEEE Transactions on Smart Grid, vol. 12, no. 3, pp. 1880-1892, May 2021. [Baidu Scholar] 

18

Z. Wang, W. Wu, and B. Zhang, “A distributed control method with minimum generation cost for DC microgrids,” IEEE Transactions on Energy Conversion, vol. 31, no. 4, pp. 1462-1470, Dec. 2016. [Baidu Scholar] 

19

Z. Wang, F. Liu, Y. Chen et al., “Unified distributed control of stand-alone DC microgrids,” IEEE Transactions on Smart Grid, vol. 10, no. 1, pp. 1013-1024, Jan. 2019. [Baidu Scholar] 

20

J. Yang, W. Feng, X. Hou et al., “A distributed cooperative control algorithm for optimal power flow and voltage regulation in DC power system,” IEEE Transactions on Power Delivery, vol. 35, no. 2, pp. 892-903, Apr. 2020. [Baidu Scholar] 

21

J. Peng, B. Fan, and W. Liu, “Voltage-based distributed optimal control for generation cost minimization and bounded bus voltage regulation in DC microgrids,” IEEE Transactions on Smart Grid, vol. 12, no. 1, pp. 106-116, Jan. 2020. [Baidu Scholar] 

22

J. Peng, B. Fan, Q. Yang et al., “Distributed event-triggered control of DC microgrids,” IEEE Systems Journal, vol. 15, no. 2, pp. 2504-2514, Jun. 2021. [Baidu Scholar] 

23

L. Xing, Q. Xu, F. Guo et al., “Distributed secondary control for DC microgrid with event-triggered signal transmissions,” IEEE Transactions on Sustainable Energy, vol. 12, no. 3, pp. 1801-1810, Jul. 2021. [Baidu Scholar] 

24

L. Ding, Q. Han, L. Wang et al., “Distributed cooperative optimal control of DC microgrids with communication delays,” IEEE Transactions on Industrial Informatics, vol. 14, no. 9, pp. 3924-3935, Sept. 2018. [Baidu Scholar] 

25

G. Chen and Z. Zhao, “Delay effects on consensus-based distributed economic dispatch algorithm in microgrid,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 602-612, Jan. 2018. [Baidu Scholar] 

26

M. Hamdi, M. Chaoui, L. Idoumghar et al., “Coordinated consensus for smart grid economic environmental power dispatch with dynamic communication network,” IET Generation, Transmission & Distribution, vol. 12, no. 11, pp. 2603-2613, Apr. 2018. [Baidu Scholar] 

27

G. Wen, X. Yu, Z. Liu et al., “Adaptive consensus-based robust strategy for economic dispatch of smart grids subject to communication uncertainties,” IEEE Transactions on Industrial Informatics, vol. 14, no. 6, pp. 2484-2496, Jun. 2018. [Baidu Scholar] 

28

J. Duan and M. Chow, “Robust consensus-based distributed energy man agement for microgrids with packet losses tolerance,” IEEE Transactions on Smart Grid, vol. 11, no. 1, pp. 281-290, Jan. 2020. [Baidu Scholar] 

29

R. Wang, Q. Li, G. Li et al., “A gossip-based distributed algorithm for economic dispatch in smart grids with random communication link failures,” IEEE Transactions on Industrial Electronics, vol. 67, no. 6, pp. 4635-4645, Jun. 2020. [Baidu Scholar] 

30

K. G. Murty and S. N. Kabadi, “Some NP-complete problems in quadratic and nonlinear programming,” Mathematical Programming, vol. 39, no. 2, pp. 117-129, Mar. 1987. [Baidu Scholar] 

31

J. Zhao and F. Dörfler, “Distributed control and optimization in DC microgrids,” Automatica, vol. 61, pp. 18-26, Nov. 2015. [Baidu Scholar] 

32

D. P. Bertsekas. Constrained Optimization and Lagrange Multiplier Methods. Nashua: Athena Scientific, 1996. [Baidu Scholar] 

33

D. P. Palomar and M. Chiang, “A tutorial on decomposition methods for network utility maximization,” IEEE Journal on Selected Areas in Communications, vol. 24, no. 8, pp. 1439-1451, Jul. 2006. [Baidu Scholar] 

34

D. G. Luenberger. Linear and Nonlinear Programming, 4th ed. Berlin: Springer International Publishing, 2016. [Baidu Scholar] 

35

J. Xu, S. Zhu, Y. C. Soh et al., “Augmented distributed gradient methods for multi-agent optimization under uncoordinated constant stepsizes,” in Proceedings of 2015 54th IEEE Conference on Decision and Control (CDC), Osaka, Japan, Dec. 2015, pp. 2055-2060. [Baidu Scholar] 

36

M. Grant and S. Boyd. (Jan. 2014). CVX: MATLAB software for disciplined convex programming, version 2.1. [Online]. Available: http://cvxr.com/cvx [Baidu Scholar] 

37

E. Dallnese and A. Simonetto, “Optimal power flow pursuit,” IEEE Transactions on Smart Grid, vol. 9, no. 2, pp. 942-952, Mar. 2018. [Baidu Scholar] 

38

X. Zhou, E. Dall’Anese, L. Chen et al., “An incentive-based online optimization framework for distribution grids,” IEEE Transactions on Automatic Control, vol. 63, no. 7, pp. 2019-2031, Jul. 2018. [Baidu Scholar] 

39

A. Bernstein and E. Dall’Anese, “Real-time feedback-based optimization of distribution grids: a unified approach,” IEEE Transactions on Control of Network Systems, vol. 6, no. 3, pp. 1197-1209, Sept. 2019. [Baidu Scholar] 

40

G. H. Golub and C. F. Loan, Matrix Computations. Baltimore: JHU Press, 2012. [Baidu Scholar] 

41

H. Robbins and D. Siegmund, “A convergence theorem for non negative almost supermartingales and some applications,” in Herbert Robbins Selected Papers. New York: Springer, 1985. [Baidu Scholar] 

42

M. Zhu and S. Martinez, “Discrete-time dynamic average consensus,” Automatica, vol. 46, no. 2, pp. 322-329, Feb. 2010. [Baidu Scholar] 

43

A. Klenke. Probability Theory: A Comprehensive Course. New York: Springer Science & Business Media, 2007. [Baidu Scholar] 

44

A. Gut. Probability: A Graduate Course. New York: Springer Science & Business Media, 2005. [Baidu Scholar]