Fully Distributed Optimal Economic Dispatch for Microgrids under Directed Communication Networks Considering Time Delays

Abstract

Distributed generation and demand-side management are expected to play a more prominent role in future power systems. However, the increased number of generations and load demands pose new challenges to optimal energy management in a microgrid. In this paper, an economic dispatch model for microgrids considering Traditional Generators (TGs), energy storage units, wind turbines (WTs), and flexible loads is established. To tackle the Economic Dispatch Problem (EDP) over directed communication networks, a fully distributed algorithm developed by leveraging a two-step state information exchange mechanism is proposed. In addition, by employing a fixed stepsize, the proposed algorithm demonstrates rapid convergence. Furthermore, our algorithm is well-suited for nonquadratic convex cost functions. Subsequently, we extend our algorithm to address imperfect communication scenarios. Even in the presence of arbitrarily large yet bounded time delays, our algorithm exhibits robustness. Finally, several numerical examples are given to verify the correctness and effectiveness of the developed results.

1. Introduction

Microgrids are important parts supporting the operation of modern power systems [1,2,3,4,5]. The increasing distributed generations and fluctuating load demands have made economic dispatch for microgrids a great challenge. In order to address the characteristics of microgrids, which are flexibility, scalability, and resilience, many researchers have focused on a multi-agent-based distributed framework [6,7,8,9,10]. In recent years, consensus-based algorithms have been widely used to solve the Economic Dispatch Problem (EDP) for microgrids. In [11], a consensus + innovation algorithm with model prediction is designed to realize the dynamic economic dispatch. In [12], based on peer-to-peer communication, a consensus algorithm is used to solve the EDP involving flexible loads. A distributed bisection method is presented in [13] to solve the EDP with general convex cost functions. In [14], a distributed algorithm based on leader–follower consensus is designed to solve the EDP without the need for initialization. In [15], a distributed algorithm is proposed to obtain optimal dispatch results within finite-step iterations.

However, the mentioned methods are only applicable under the undirected graph, which requires that information be exchanged bidirectionally among agents. In an actual communication network, one link may be damaged or data packet loss may occur [16,17]; the condition of an undirected graph is hard to satisfy. In addition, directed networks offer lower cost, better privacy preservation, and improved robustness [18,19,20,21,22]. For the EDP under a directed graph (digraph), a push–pull algorithm is proposed in [23]. A distributed algorithm using the weight matrix of row stochastically is presented in [24]. Based on the push-sum algorithm, a modified algorithm in [25] is proposed to solve the EDP over time-varying digraphs. A robust algorithm considering the packet drops is provided in [26]. Most existing methods use decaying stepsize, resulting in slow convergence. A fixed stepsize can accelerate convergence, but it leads the results to oscillate in an interval. In the power system application, achieving exact optimal results is necessary. In [27], a novel consensus-based algorithm with fixed stepsize is proposed to solve the social welfare maximization problem. However, it is only applicable to the undirected graph. Inspired by the method proposed in [28], we design a distributed consensus-based algorithm using the fixed stepsize to solve the EDP under the unbalanced digraphs.

In modern power systems, the imperative integration of renewable energy generators has emerged as a strategic response to pressing challenges in the realms of energy security and environmental sustainability. Among various kinds of renewable energy generators, wind turbines (WTs) have garnered particular attention, owing to the advantages such as the free availability and environmental friendliness of wind energy, as well as the maturity of turbine techniques [29,30,31,32]. With the increasing penetration of stochastic wind power into microgrids, it is essential to reformulate the Economic Dispatch (ED) problem, incorporating not only Traditional Generators (TGs) but also accounting for the presence of renewable energy generators like WTs. However, the intermittent and stochastic nature of wind speed leads to uncertainty in wind power output, which poses a considerable challenge in designing effective distributed economic dispatch strategies. Due to the integration of energy storage units into wind farms, the gap between the scheduled wind power and the available wind power can be eliminated, and then random wind power can be scheduled like TGs. In [33], a deterministic model was adopted to describe the cost function of wind power, where the underestimation and overestimation costs of available wind power are included. Based on the dispatch model, a distributed strategy using the projected gradient algorithm is proposed in [34]. In [35], a prime–dual consensus algorithm is proposed to solve the EDP considering random wind power. However, the cost functions of wind power are not quadratic anymore, as described in [36]. Hence, some of the distributed methods mentioned above may fail to work. It is necessary to devise a distributed method for solving the EDP with more generalized cost functions.

All the aforementioned distributed algorithms are implemented in a perfect communication network. However, in an actual communication network, the information received from neighbors is always subject to time delays induced by the physical distance between them. Therefore, it is desirable to analyze the effects of nonideal communication conditions on the existing algorithms. In [37], the simulation results show that delays may cause instability when using the consensus algorithm to solve the EDP. In [38], a theoretical upper bound of time delays is obtained, ensuring that distributed EDP can be solved. In [39], both communication delays and measurement delays are considered in designing a distributed algorithm for EDP. Nonetheless, it should be noted that these findings are based on the utilization of quadratic cost functions for modeling the EDP. The effectiveness of the proposed algorithm extends to scenarios where the cost function is nonquadratic. Furthermore, an extension of our algorithm is introduced to enhance its robustness in the presence of communication delays. The main contributions of this paper are as follows:

(1)

An economic dispatch model for microgrids considering random wind power and flexible loads is established. A fully distributed consensus-based algorithm using fixed stepsize is proposed, which is suitable for unbalanced directed communication networks. Different from the traditional distributed methods, where the update of the state information only depends on the information of the previous instant, the latest two-step information is used to iterate, ensuring convergence.

(2)

Our algorithm is initialization-free so that it can converge to the optimal results with arbitrary initial power values. In addition, for nonquadratic convex cost functions, such as the nonlinear function including integral terms that are used to describe the expected cost of wind power, our algorithm is also effective.

(3)

The proposed algorithm can solve the EDP with arbitrarily large but bounded time delays if the upper bound of time delays is preknown. Due to its fast convergence even under delays, our algorithm can be applied to the practical operations of microgrids.

The rest of this paper is organized as follows: Section 2 describes the EDP model for microgrids. A fully distributed algorithm for EDP under the directed networks with time delays is proposed in Section 3. Case studies are carried out in Section 4 to demonstrate the effectiveness of the proposed method. Finally, Section 5 concludes this paper.

2. Problem Formulation

2.1. Cyber–Physical System Architecture for Microgrids

A microgrid can be effectively represented as a cyber–physical system (CPS). Leveraging advanced communication and computing equipment, the implementation of distributed economic dispatch within a microgrid becomes more feasible. The CPS architecture, as shown in Figure 1, comprises both a physical layer and a cyber layer. Each generator or load node corresponds to an agent equipped with communication and computing capabilities. The communication links among these agents exhibit flexibility and may not necessarily mirror the physical connections.

In microgrids, the security and stability of system operation rely on the coordination and cooperation of each CPS unit. The structure of a CPS unit is illustrated in Figure 2. The peer-to-peer communication mode is utilized to establish the connection between different CPS units. Within this framework, communication units are responsible for transmitting and receiving information from their neighboring units. The pivotal role in achieving global control is played by the information processing unit. These units, on the one hand, continuously monitor and sample the operational states of physical devices. On the other hand, they analyze information from neighboring units, compute outgoing information, and generate control instructions. In the following sections, we use this framework as the foundation for our discussion on distributed economic dispatch.

2.2. Directed Communication Networks

Usually, the undirected communication topology is considered for the sake of algorithm design. In our work, the unbalanced directed topology is considered. A typical example is shown in Figure 3. A directed graph is denoted by $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with a set of nodes $\mathcal{V}=\{1,\dots ,N\}$ and edges $\mathcal{E}$. $(i,j)\in \mathcal{E}$ indicates there exists a directed path from the node i to j. The set of nodes that can receive information from node i are called out-neighbors of node i, denoted by ${\mathcal{N}}_i^{\mathrm{out}}\left(t\right)=\{j\mid (i,j)\in \mathcal{E}\}$. ${\mathcal{N}}_i^{\mathrm{in}}\left(t\right)=\{j\mid (j,i)\in \mathcal{E}\}$ is the set of in-neighbors of node i.

2.3. Economic Dispatch Model in Microgrids

Generally, microgrids consist of generators, storage devices, and load units. In this paper, we consider an ED model which involves random wind power. The EDP for microgrids is minimizing the operation cost of generation sides and maximizing the benefit of demand sides. In some issues, the carbon emission cost is also considered in the EDP [40,41]. All these cost functions or benefit functions can be unified as a general function $f_i(·)$ for each agent i. Here, we assume that $f_i(·)$ is strictly convex and continuously differentiable.

For conventional generators, the cost function is usually modeled as the following quadratic function:

$$f_i\left(P_{G,i}\right)=\frac{1}{2}{a}_i{P}_{G,i}^2+b_i{P}_{G,i}+c_i(P_{G,i}^{min}\le P_{G,i}\le P_{G,i}^{max})$$

(1)

where $P_{G,i}$ denotes the output power of the ith generator. $P_{G,i}^{min}$ and $P_{G,i}^{max}$ are the lower and upper limits of output power, respectively. $a_i$, $b_i$, and $c_i$ are the corresponding cost coefficients.

In this paper, we define the utility benefit function for users based on the aggregation of diverse tasks, rather than focusing on individual appliance power consumption. It is important to note that, from a user perspective, the ability to accomplish more tasks often necessitates a higher power consumption. Therefore, it is reasonable to postulate that the utility benefit function exhibits a nondecreasing trend and eventually reaches a state of saturation, a concept well supported by prior studies [42,43]. User demand requests typically comprise two components, schedulable and nonschedulable loads, which can be effectively represented by appropriately defining the demand lower bound.

$$f_j\left(P_{L,j}\right)=\frac{1}{2}{a}_j{P}_{L,j}^2+b_j{P}_{L,j}(P_{L,j}^{min}\le P_{L,j}\le P_{L,j}^{max})$$

(2)

where $P_{L,j}$ denotes the load demand of the jth user. $P_{L,j}^{min}$ and $P_{L,j}^{max}$ are the lower and upper bounds, respectively. $a_j$, $b_j$ are the corresponding benefit coefficients.

For the energy storage units, the cost function is modeled in the following form [11]:

$$f_k\left(P_{B,k}\right)=\frac{1}{2}{a}_k{P}_{B,k}^2+b_k{P}_{B,k}(P_{B,k}^{min}\le P_{B,k}\le P_{B,k}^{max})$$

(3)

where $P_{B,k}$ denotes the power value of the kth energy storage unit. $P_{B,k}^{min}$ and $P_{B,k}^{max}$ are the lower and upper bounds, respectively. $a_k$ and $b_k$ are the corresponding cost coefficients.

For TGs, the scheduled and actual power outputs are always the same. However, owing to the stochastic nature of wind speed, the generated power $P_{W,l,av}$ at the lth WT is a random variable, deviating potentially from the scheduled power $P_{W,l}$. The integration of Energy Storage Systems (ESSs) into the WT units offers a solution to ensure that the total WT output aligns with the scheduled value. For instance, if the scheduled power output $P_{W,l}$ exceeds $P_{W,l,av}$, the ESSs can compensate for the disparity. Conversely, if the scheduled power output $P_{W,l}$ is less than $P_{W,l,av}$, the WT should limit its output to $P_{W,l}$, and the surplus wind power can be utilized to charge the ESSs. To quantify the cost associated with the WT, a methodology involving overestimation and underestimation costs has been proposed in [36]. Therefore, the overall cost for the lth WT can be expressed as

$$f_l\left(P_{W,l}\right)=d_l{P}_{W,l}+{\eta }_l^{ue}\mathbf{E}\left(W_l^{ue}\right)+{\eta }_l^{oe}\mathbf{E}\left(W_l^{oe}\right)(0\le P_{W,l}\le P_{W,l}^r)$$

(4)

where the first term is the direct cost of wind power with the cost coefficient $d_l$. The second and third terms are the underestimation and overestimation costs. ${\eta }_l^{ue}$ and ${\eta }_l^{oe}$ are the corresponding coefficients. $P_{W,l}^r$ is the rated wind power of the lth WT. The definitions of $\mathbf{E}\left(W_l^{ue}\right)$ and $\mathbf{E}\left(W_l^{oe}\right)$ are described as follows:

$$\mathbf{E}\left(W_l^{ue}\right)={\int }_{P_{W,l}}^{P_{W,l}^r}\left(P_{W,l,av}-P_{W,l}\right)f_W\left(P_{W,l,av}\right)d{P}_{W,l,av}$$

(5)

$$\mathbf{E}\left(W_l^{oe}\right)={\int }_0^{P_{W,l}}\left(P_{W,l}-P_{W,l,av}\right)f_W\left(P_{W,l,av}\right)d{P}_{W,l,av}$$

(6)

where $f_W\left(P_{W,l,av}\right)$ denotes the probability distribution function of the available wind power $P_{W,l,av}$, which is described by Weibull distribution in this paper. The closed forms of $\mathbf{E}\left(W_l^{ue}\right)$ and $\mathbf{E}\left(W_l^{oe}\right)$ are given as follows:

$$\begin{array}{cc}\hfill \mathbf{E}\left(W_l^{ue}\right)=& \left(P_{W,l}^r-P_{W,l}\right)\left(exp\left(-\frac{v_r^{\kappa }}{c^{\kappa }}\right)-exp\left(-\frac{v_{out}^{\kappa }}{c^{\kappa }}\right)\right)\hfill \\ \hfill & +\left(\frac{P_{W,l}^r{v}_{in}}{v_r-v_{in}}+P_{W,l}\right)\left(exp\left(-\frac{v_r^{\kappa }}{c^{\kappa }}\right)-exp\left(-\frac{v_l^{\kappa }}{c^{\kappa }}\right)\right)\hfill \\ \hfill & +\frac{P_{W,l}^{r}c}{v_r-v_{in}}\left(\Gamma \left(1+\frac{1}{\kappa },{\left(\frac{v_l}{c}\right)}^{\kappa }\right)-\Gamma \left(1+\frac{1}{\kappa },{\left(\frac{v_r}{c}\right)}^{\kappa }\right)\right)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathbf{E}\left(W_l^{oe}\right)& =P_{W,l}\left(1-exp\left(-\frac{v_{in}^{\kappa }}{c^{\kappa }}\right)+exp\left(-\frac{v_{out}^{\kappa }}{c^{\kappa }}\right)\right)\hfill \\ \hfill & +\left(\frac{P_{W,l}^r{v}_{in}}{v_r-v_{in}}+P_{W,l}\right)\left(exp\left(-\frac{v_{in}^{\kappa }}{c^{\kappa }}\right)-exp\left(-\frac{v_l^{\kappa }}{c^{\kappa }}\right)\right)\hfill \\ \hfill & +\frac{P_{W,l}^{r}c}{v_r-v_{in}}\left(\Gamma \left(1+\frac{1}{\kappa },{\left(\frac{v_l}{c}\right)}^{\kappa }\right)-\Gamma \left(1+\frac{1}{\kappa },{\left(\frac{v_{in}}{c}\right)}^{\kappa }\right)\right)\hfill \end{array}$$

(8)

where $v_r$, $v_{in}$, and $v_{out}$ are the rated, cut-in, and cut-out wind speeds. $\kappa$ and c denote the scale factor and shape factor of the Weibull distribution, respectively. $\Gamma (•)$ is a standard incomplete gamma function. $v_l$ is an intermediary variable, which is defined as

$$v_j=v_{in}+\frac{\left(v_r-v_{in}\right)P_{W,l}}{P_{W,l}^r}$$

(9)

The aim of the optimal dispatch for microgrids is to coordinate all the participants to maintain the active power balance and maximize social welfare. Therefore, the EDP is modeled as

$$\begin{array}{cc}\mathrm{Max}& -\sum _{i\in {\Omega }_G}{f}_i\left(P_{G,i}\right)-\sum _{j\in {\Omega }_L}{f}_j\left(P_{L,j}\right)-\sum _{k\in {\Omega }_B}{f}_k\left(P_{B,k}\right)-\sum _{l\in {\Omega }_W}{f}_l\left(P_{W,l}\right)\hfill \\ \mathrm{s}.\mathrm{t}.& \sum _{i\in {\Omega }_G}{P}_{G,i}+\sum _{j\in {\Omega }_L}{P}_{L,j}+\sum _{k\in {\Omega }_B}{P}_{B,k}+\sum _{k\in {\Omega }_W}{P}_{W,l}=0,\hfill \\ & P_{G,i}^{min}\le P_{G,i}\le P_{G,i}^{max},\hfill \\ & P_{L,j}^{min}\le P_{L,j}\le P_{L,j}^{max},\hfill \\ & P_{B,k}^{min}\le P_{B,k}\le P_{B,k}^{max},\hfill \\ & 0\le P_{W,l}\le P_{W,l}^r.\hfill \end{array}$$

(10)

where ${\Omega }_G$, ${\Omega }_L$, ${\Omega }_B$, ${\Omega }_W$ denote the sets of the TG units, load units, the energy storage units, and the WT units, respectively.

In order to describe the optimization model in a unified form, let $P_i$ denote the power of agent i, which denotes the TG unit the load unit, the storage unit, or the WT unit, i.e., $i\in {\Omega }_G\bigcup {\Omega }_L\bigcup {\Omega }_B\bigcup {\Omega }_W$. Assume that N denotes the number of all participants in a microgrid. For TG agents, $P_i>0$. For load agents, $P_i<0$. For energy storage agents, $P_i>0$ when discharging and $P_i<0$ when charging. For WT agents, $P_i>0$. Mathematically, the EDP for microgrids is modeled as a unified form,

$$min\sum _{i=1}^N{f}_i\left(P_i\right)=\sum _{i\in {\Omega }_G}{f}_i\left(P_{G,i}\right)+\sum _{j\in {\Omega }_L}{f}_j\left(P_{L,j}\right)+\sum _{k\in {\Omega }_B}{f}_k\left(P_{B,k}\right)+\sum _{l\in {\Omega }_W}{f}_l\left(P_{W,l}\right)$$

(11)

subject to two constraints:

(1)

The capacity constraint

$$P_i^{min}\le P_i\le P_i^{max},i=1,\dots ,N$$

(12)

where $P_i^{min}$ and $P_i^{max}$ represent the lower and upper limits, which are presented in Equations (1)–(4). The capacity constraint of agent i can also be written as a box constraint $P_i\in {\Omega }_i=\left\{P_i\mid P_i^{min}\le P_i\le P_i^{max}\right\}$.

(2)

The power balance constraint

$$\sum _{i\in {\Omega }_G}{P}_{G,i}+\sum _{j\in {\Omega }_L}{P}_{L,j}+\sum _{k\in {\Omega }_B}{P}_{B,k}+\sum _{l\in {\Omega }_W}{P}_{W,l}=\sum _{i=1}^N{P}_i=0$$

(13)

Then, the Lagrangian function is defined as follows:

$$L(P,\lambda )=\sum _{i=1}^N{f}_i\left(P_i\right)-\lambda \sum _{i=1}^N{P}_i$$

(14)

where $\lambda$ is the Lagrange multiplier, that is, the incremental cost. The dual problem of the EDP is

$$\underset{\lambda \in {\mathbb{R}}^{+}}{max}minL(P,\lambda )=\underset{\lambda \in {\mathbb{R}}^{+}}{max}\sum _i^N{\Theta }_i\left(\lambda \right)$$

(15)

where

$${\Theta }_i\left(\lambda \right)=\underset{P_i\in {\Omega }_i}{min}{f}_i\left(P_i\right)-\lambda P_i$$

(16)

The subproblem (16) is minimized when satisfying $P_i\left(\lambda \right)={\phi }_i\left(\nabla f_i^{-1}\left(\lambda \right)\right)$, where

$${\phi }_i\left(\nabla f_i^{-1}\left(\lambda \right)\right)=\left\{\begin{array}{cc}{P}_i^{min},\hfill & \nabla f_i^{-1}\left(\lambda \right)<P_i^{min}\hfill \\ \nabla f_i^{-1}\left(\lambda \right),\hfill & P_i^{min}\le \nabla f_i^{-1}\left(\lambda \right)\le P_i^{max}\hfill \\ P_i^{max},\hfill & \nabla f_i^{-1}\left(\lambda \right)>P_i^{max}\hfill \end{array}\right.$$

(17)

where $\nabla f_i^{-1}\left(\lambda \right)$ is the inverse function of $\nabla f_i$, noting that $\nabla f_i^{-1}$ is continuous and strictly monotonic. The gradient of ${\Theta }_i\left(\lambda \right)$ is $\nabla {\Theta }_i\left(\lambda \right)=-P_i\left(\lambda \right)$. Obviously, the gradient $\nabla {\Theta }_i$ is bounded and Lipschitz-continuous. The gradient ascent method (18) is applied to solve the dual problem.

$$\lambda (k+1)=\lambda \left(k\right)+\rho \left(k\right)\sum _{i=1}^N-P_i\left(\lambda \left(k\right)\right)$$

(18)

where $\rho \left(k\right)$ is the stepsize at iteration k. Assume that ${\lambda }^{*}$ is the optimal solution to the dual problem. Accordingly, the optimal solution of the primal EDP is $P_i^{*}={\phi }_i\left(\nabla f_i^{-1}\left({\lambda }^{*}\right)\right)$.

3. Distributed Algorithm for EDP under the Unbalanced Directed Graph with Time Delays

In this section, we propose an agent-based distributed algorithm to solve the EDP under the unbalanced digraph with random but bounded time delays.

3.1. Distributed Algorithm for EDP Using Two-Step Information

Traditional distributed algorithms only use one-step information for iterations. This means the update of the state information only depends on the information at the previous instant. A fully distributed algorithm using the latest two-step information with a fixed stepsize is proposed. The update protocol of agent i is presented as follows:

$$\begin{array}{}& \hfill r_i(k+1)& =r_i\left(k\right)+{\displaystyle \sum _{j=1}^N}{w}_{ij}{r}_j\left(k\right)-{\displaystyle \sum _{j=1}^N}{\tilde{w}}_{ij}{r}_j(k-1)\hfill \mathrm{(19a)}& & -\rho \left(P_i\left(k\right)-P_i(k-1)\right)\hfill \mathrm{(19b)}& \hfill y_i(k+1)& ={\displaystyle \sum _{j=1}^N}{w}_{ij}{y}_j\left(k\right)\hfill \mathrm{(19c)}& \hfill {\lambda }_i(k+1)& =\frac{r_i(k+1)}{y_i(k+1)}\hfill \mathrm{(19d)}& \hfill P_i(k+1)& ={\phi }_i\left(\nabla f_i^{-1}\left[{\lambda }_i(k+1)\right]\right)\hfill \end{array}$$

Then, we can rewrite Equations (19a) and (19b) into the following compact form:

$$\begin{array}{}\mathrm{(20a)}& \mathbf{r}(k+1)=\mathbf{r}\left(k\right)+W\mathbf{r}\left(k\right)-\tilde{W}\mathbf{r}(k-1)-\rho \left(\mathbf{P}\left(k\right)-\mathbf{P}(k-1)\right)\mathrm{(20b)}& \mathbf{y}(k+1)=W\mathbf{y}\left(k\right)\end{array}$$

where $\mathbf{r}={[r_1,\dots ,r_N]}^T$ are the auxiliary variables to obtain the optimal solutions. $\mathbf{y}={[y_1,\dots ,y_N]}^T$ is utilized to address the unbalancedness of the digraphs. $\rho$ is the fixed stepsize within a certain interval. $W=\left\{w_{ij}\right\}$ is the weight matrix, which is designed as follows:

$$w_{ij}=\left\{\begin{array}{cc}1-\epsilon \left|{\mathcal{N}}_i^{out}\right|,\hfill & i=j\hfill \\ \epsilon ,\hfill & i\ne j,i\in {\mathcal{N}}_j^{out}\hfill \end{array}\right.$$

(21)

where $\epsilon$ is a positive constant close to 0. Only the out-degree of each agent i is needed. If agent i is aware of the out-degrees of its in-neighbors, the weight $w_{ij}$ can also be designed as

$$w_{ij}=\left\{\begin{array}{cc}1/\left({\mathcal{N}}_j^{\mathrm{out}}\mid +1\right),\hfill & \mathrm{if}i\in {\mathcal{N}}_j^{\mathrm{out}}\cup \left\{i\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(22)

The $\tilde{W}=\left\{{\tilde{w}}_{ij}\right\}$ is also a weight matrix related to the penultimate state, which is designed by

$${\tilde{w}}_{ij}=\left\{\begin{array}{cc}\mu +(1-\mu )w_{ij},\hfill & i=j\hfill \\ (1-\mu )w_{ij},\hfill & i\ne j\hfill \end{array}\right.$$

(23)

where $\mu \in (0,\frac{1}{2}]$, and the two weight matrices satisfy that $\tilde{W}=\mu I_N+(1-\mu )W$.

In the initialization of our proposed algorithm, $r_i\left(0\right)$ and $P_i\left(0\right)$ can be arbitrary, and $y_i\left(0\right)=1$ for $i=1,\dots ,N$. At the first iteration, let $r_i(-1)=0$ and $P_i(-1)=0$. The iteration is given by

$$\begin{array}{}\mathrm{(24a)}& \hfill {\lambda }_i\left(0\right)& =\frac{r_i\left(0\right)}{y_i\left(0\right)}\hfill \mathrm{(24b)}& \hfill r_i\left(1\right)& =\sum _{j=1}^N{w}_{ij}{r}_j\left(0\right)-\rho \nabla f_i^{-1}\left[{\lambda }_i\left(0\right)\right]\hfill \mathrm{(24c)}& \hfill y_i\left(1\right)& =\sum _{j=1}^N{w}_{ij}{y}_j\left(0\right)\hfill \end{array}$$

In our algorithm, there is no need for strict initialization. Due to the fixed stepsize, the proposed algorithm has an excellent convergence rate. According to definitions, the weight matrices W and $\tilde{W}$ are row-stochastic, satisfying ${\mathbf{1}}^{T}W={\mathbf{0}}^T$ and ${\mathbf{1}}^T\tilde{W}={\mathbf{0}}^T$. In Equation (20), let $g_i\left({\lambda }_i\left(k\right)\right)=P_i\left(k\right)$. Summing the iterations in Equation (20) from 1 to $k+1$, we can obtain $\mathbf{r}(k+1)=W\mathbf{r}\left(k\right)-\rho \mathbf{P}\left(k\right)+{\sum }_{t=0}^{k-1}(W-\tilde{W})\mathbf{r}\left(t\right)$. When $k\to \infty$, ${\lim }_{k\to \infty }\mathbf{r}\left(k\right)={\mathbf{r}}^{*}=W{\mathbf{r}}^{*}$, and ${\lambda }_i$ reaches optimal consensus ${\lambda }^{*}$ due to the role of $\mathbf{y}$. Then, we can obtain $\rho \mathbf{P}(\infty )={\sum }_{t=0}^{\infty }(W-\tilde{W})\mathbf{r}\left(t\right)$. Noting that ${\mathbf{1}}^T(W-\tilde{W})={\mathbf{0}}^T$, we have ${\mathbf{1}}^T\mathbf{P}(\infty )=0$. Thus, the power balance constraint is satisfied. The iteration process can be terminated as the error value satisfies ${\parallel \mathbf{P}(k+1)-\mathbf{P}\left(k\right)\parallel }_2<ϵ$, where $ϵ>0$ is an error threshold.

Remark 1.

The effectiveness of our proposed algorithm is ensured under two conditions: the objective function is strongly convex, and the directed communication network is connected. Firstly, the objective function being strongly convex implies that it has a unique global minimum. This ensures that there is only one optimal solution for the problem at hand. Consequently, any iterative algorithm designed to minimize this objective function will converge to this unique optimal solution. Secondly, the directed communication network being connected means that there is a path of communication between any two nodes in the network. This is crucial for the algorithm to function effectively, as information needs to be exchanged between different nodes to update their values iteratively. Without a connected network, certain nodes may be isolated and not receive or transmit important information, leading to suboptimal or incorrect results. By satisfying these conditions, our proposed algorithm guarantees convergence to the optimal results, regardless of the initial values chosen. This means that the algorithm is robust and reliable in solving the given EDP for microgrids.

3.2. Robustness to Communication Time Delays

In the digraph, ${\tau }_{ji}$ denotes the time delays during the information transmission from the agent j to the agent i. We assume time delays are bounded by

$$0\le {\tau }_{ji}\le \overline{\tau },j\in {\mathcal{N}}_i^{\mathrm{i}n}$$

(25)

Specifically, self-delay ${\tau }_{ii}$ is assumed to be 0, which means that the state update process depending on local information has no delays. Therefore, our algorithm under time delays is executed by

$$\begin{array}{}& \hfill r_i(k+1)=r_i\left(k\right)& +\sum _{j=1}^N{w}_{ij}{r}_j(k-{\tau }_{ji})-\sum _{j=1}^N{\tilde{w}}_{ij}{r}_j(k-1-{\tau }_{ji})\hfill \mathrm{(26a)}& & -\rho \left(P_i\left(k\right)-P_i(k-1)\right)\hfill \mathrm{(26b)}& \hfill y_i(k+1)& =\sum _{j=1}^N{w}_{ij}{y}_j(k-{\tau }_{ji})\hfill \mathrm{(26c)}& \hfill {\lambda }_i(k+1)& =\frac{r_i(k+1)}{y_i(k+1)}\hfill \mathrm{(26d)}& \hfill P_i(k+1)& ={\phi }_i\left(\nabla f_i^{-1}\left[{\lambda }_i(k+1)\right]\right)\hfill \end{array}$$

An augmented digraph method is utilized to model the communication time delays. For each node, there are $\overline{\tau }$ augmented virtual nodes to store the delayed information. As shown in Figure 4, the gray nodes are the virtual nodes. Information sent from the neighboring nodes first reaches the virtual nodes. The number of virtual nodes connected to the real nodes depends on the delay bound. Thus, for each link $(j,i)\in \mathcal{E}$, there are $|{\tau }_{ji}|$ augmented links. The information exchange under the original network with time delays is equivalent to the information exchange without delays under the augmented digraph. The modified update protocol in compact form as

$$\begin{array}{}\mathrm{(27a)}& \hfill \mathbf{r}(k+1)=\mathbf{r}\left(k\right)& +W^{\tau }\mathbf{r}\left(k\right)-{\tilde{W}}^{\tau }\mathbf{r}(k-1)-\rho \left(\mathbf{P}\left(k\right)-\mathbf{P}(k-1)\right)\hfill \mathrm{(27b)}& \hfill \mathbf{y}(k+1)& =W^{\tau }\mathbf{y}\left(k\right)\hfill \mathrm{(27c)}& \hfill {\lambda }_i(k+1)& =\frac{r_i(k+1)}{y_i(k+1)}\hfill \mathrm{(27d)}& \hfill P_i(k+1)& ={\phi }_i\left(\nabla f_i^{-1}\left[{\lambda }_i(k+1)\right]\right)\hfill \end{array}$$

where

$$W^{\tau }=\left[\begin{array}{ccccc}{W}^{\left(0\right)}& I_N& \mathbf{0}& \cdots & \mathbf{0}\\ W^{\left(1\right)}& \mathbf{0}& I_N& \dots & \mathbf{0}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ W^{(\overline{\tau }-1)}& \mathbf{0}& \mathbf{0}& \dots & I_N\\ W^{\left(\overline{\tau }\right)}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{0}\end{array}\right]$$

(28)

where $W^{\left(0\right)}$, $W^{\left(1\right)}$, ⋯, $W^{\left(\overline{\tau }\right)}$ are the matrices relying on the values of time delays, which are defined as

$$W_{ij}^{\left(m\right)}=\left\{\begin{array}{cc}{W}_{ij},\hfill & \mathrm{if}{\tau }_{ji}=m,(j,i)\in \mathcal{E}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(29)

Accordingly, ${\tilde{W}}^{\tau }$ is modified by $W^{\tau }$ using (23). Note that ${\tilde{W}}^{\tau }$ and $W^{\tau }$ are still row-stochastic. Under the algorithm (17), the EDP can be solved over the augmented digraph. Therefore, our proposed algorithm is also efficient under time delays.

4. Simulation Results

In this section, several case studies are presented in order to illustrate and validate the proposed algorithm. The case studies are simulated in the MATLAB R2021a environment on a laptop with Intel Core i5-7300U CPU @2.60 GHz and 8 GB RAM. The test systems are based on the IEEE 14-bus system and the IEEE 118-bus system, considering both without and with time-delay scenarios.

4.1. Tests on the IEEE 14-Bus System

Herein, the simulation cases are based on the IEEE 14-bus system, including three TG units, one energy storage unit, one WT unit, and nine load units. In the test power network, the corresponding generator parameters, energy storage parameters, and load parameters are provided in

Table 1. Parameters of the units in IEEE 14-bus system.

Unit	${\mathit{a}}_{\mathit{i}}$ ($/kW2h)	${\mathit{b}}_{\mathit{i}}$ ($/kWh)	${\mathit{P}}_{\mathit{i}}^{\mathit{min}}$ (kW)	${\mathit{P}}_{\mathit{i}}^{\mathit{max}}$ (kW)	Unit	${\mathit{a}}_{\mathit{i}}$ ($/kW2h)	${\mathit{b}}_{\mathit{i}}$ ($/kWh)	${\mathit{P}}_{\mathit{i}}^{\mathit{min}}$ (kW)	${\mathit{P}}_{\mathit{i}}^{\mathit{max}}$ (kW)
G1	0.080	2.25	30	60	B8	0.035	0.01	−25	25
G2	0.062	4.20	20	50	L9	0.060	8.05	−15	−35
G3	0.075	3.25	30	55	L10	0.078	8.45	−15	−35
L4	0.072	8.25	−10	−30	L11	0.080	8.75	−15	−40
L5	0.066	7.20	−5	−20	L12	0.085	9.00	−15	−40
W6	-	-	-	-	L13	0.069	7.05	0	−15
L7	0.075	7.80	−10	−30	L14	0.077	8.15	−10	−35

[44]. The parameters of wind and WT are presented in

Table 2. Parameters of Wind and WT.

Wind	$v_{in}$ 5	$v_{out}$ 45	$v_r$ 15	$(c,\kappa )$ (8, 2)
Wind Turbine	$d_l$ 5	${\eta }_l^{ue}$ 3.1	${\eta }_l^{oe}$ 3.1	$P_{W,l}^r$ 50

[34].

The communication network among agents is modeled as a directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, with the edge set $\mathcal{E}=\left\{\right(i,i+1)$, $(i,i+2)\mid 1\le i\le 12\}$∪$\left\{\right(13,14)$, $(13,1)$, $(14,1)$, $(3,2)$, $(3,9)$, $(4,10)$, $(5,11)$, $(6,12)\}$. The TG agents are $\{1,2,3\}$. The energy storage agent is $\left\{8\right\}$. The WT agent is $\left\{6\right\}$, and the load agents are $\{4,5,7,9,10,11,12,13,14\}$. The corresponding directed communication network based on the IEEE 14-bus system is illustrated in Figure 5. The stepsize $\rho$ is set as 0.018, and $\mu$ is set as 0.2.

4.1.1. Without Time Delays

In this case study, time delays are not considered. The proposed distributed algorithm is implemented in an ideal communication environment. Figure 6, Figure 7 and Figure 8 show the simulation results of the optimal incremental cost, optimal power, and net power. After about 90 iterations, the optimal solutions are obtained. The incremental costs converge to the optimal value: ${\lambda }^{*}=6.5912$. The optimal power results are listed in

Table 3. The optimal power values.

Unit	Optimal Power (kW)	Unit	Optimal Power (MW)
G1	54.2653	B8	18.8035
G2	38.5681	L9	−24.3129
G3	44.5496	L10	−23.8304
L4	−23.0385	L11	−26.9847
L5	−9.2239	L12	−28.3385
W6	22.5521	L13	−6.6489
L7	−16.1170	L14	−20.2438

. In addition, the convergence performance is compared with the algorithms in [24,25], which are shown in

Table 4. Algorithm comparison results.

	λ* ($/kWh)	Net Power (kW)	Number of Iterations
Algorithm in [24]	6.5866	−0.01	8951
Algorithm in [25]	6.5916	0.0978	2265
Our algorithm	6.5912	8.3393 × 10−4	88

. The algorithms proposed in [24,25] adopt the decaying stepsize to make sure the results can converge to the optimum accurately. Our algorithm adopts a fixed stepsize. However, although a fixed stepsize can accelerate the convergence rate, it only ensures that the results converge within the neighborhood of the optimal value. Therefore, we propose a two-step state information exchange mechanism to address this problem. Due to the fixed stepsize rather than the decaying stepsize, our algorithm has a greatly faster convergence rate than the existing algorithms for EDP under digraph.

4.1.2. Considering Nonquadratic Functions

In this case study, which is different from the traditional quadratic function forms, the cost functions of TG agent 1 and TG agent 3 are replaced by the general convex functions as follows:

$$\begin{array}{cc}\hfill f_1\left(P_1\right)& =\frac{{\left(P_1+25\right)}^2}{25}+50exp\left(\frac{P_1+40}{100}\right)\hfill \\ \hfill f_3\left(P_3\right)& =\frac{{\left(P_3+57.14\right)}^2}{28.58}+7\times {10}^{-6}{P}_3^4\hfill \end{array}$$

As shown in Figure 9, Figure 10 and Figure 11, after about 100 iterations, the optimal results can be obtained. The forms of the cost functions do not influence the convergence rate of our proposed algorithm. In a real economic dispatch model, the cost functions of battery storage and the benefit functions of demand sides are usually modeled as complex functions. If the functions are convex, or if they can be converted into convex functions, the proposed algorithm is still efficient. Therefore, our algorithm is more suitable for practical applications.

4.1.3. With Time Delays

Herein, the bounded communication delays are considered for the EDP. Assume that for each link in the digraph, random time delays bounded by $\overline{\tau }=3$ exist. As shown in Figure 12, Figure 13 and Figure 14, even under delays, our algorithm can converge to the same optimal results after about 2000 iterations. Due to the delay effects of the communication network, the iteration steps increase compared with those in a perfect communication environment. Then, a test with larger time delays is investigated. The delay bound is set as $\overline{\tau }=7$. As shown in Figure 15, Figure 16 and Figure 17, after about 4000 iterations, the optimal results can be obtained. It can be seen that the larger the delay bound, the more iterations the proposed algorithm requires. Fortunately, if the time delays are preknown in a bounded range, our algorithm can be designed for accurate convergence. Therefore, our algorithm is more suitable for practical applications, even when the communication environment is imperfect.

4.2. Tests on the IEEE 118-Bus Test System

In this section, we further apply our proposed algorithm to a larger test system, the IEEE 118-bus test system, which is composed of 50 conventional generators, 4 wind turbines, and 41 load units. Without loss of generality, the parameters of the wind turbines are chosen from Table 2. The cost coefficients of all conventional generators are chosen from the following ranges:

$$\left\{\begin{array}{cc}{a}_i\in [0.07,0.10],\hfill & i=1,\dots ,50\hfill \\ b_i\in [2.5,5.0],\hfill & i=1,\dots ,50\hfill \end{array}\right.$$

The benefit coefficients of all load units are chosen by the following range:

$$\left\{\begin{array}{cc}{a}_i\in [0.06,0.09],\hfill & i=51,\dots ,95\hfill \\ b_i\in [7.0,9.0],\hfill & i=51,\dots ,95\hfill \end{array}\right.$$

4.2.1. Without Time Delays

In this case study, time delays are not considered. The simulation results are shown in Figure 18, Figure 19 and Figure 20. After about 5000 iterations, the optimal results can be obtained. The incremental costs converge to 5.6968. This verifies that our proposed algorithm can be applied in the actual large power system.

4.2.2. With Time Delays

In this case study, time delays are considered. The delay bound is set as $\overline{\tau }=3$. The simulation results of incremental costs, optimal power, and net power are shown in Figure 21, Figure 22 and Figure 23. As the system scales up, and considering the impact of communication delay, the number of iterations required for algorithm convergence also increases accordingly. After about 70,000 iterations, the optimal results can be obtained. The incremental costs converge to 5.6848. This verifies that our algorithm has good robustness to time delays, and it can adapt well to a nonideal communication environment in a large system.

5. Conclusions

In this paper, an economic dispatch model considering random wind power and flexible loads is established, which is described as a general convex cost function. A fully distributed consensus-based algorithm with a fixed stepsize using a two-step information exchange mechanism is proposed to deal with the EDP for microgrids over directed communication networks. Our algorithm is initialization-free and can converge the optimal results at a faster speed than the existing distributed approaches that adopt the decaying stepsize. Communication time delays existing during the data transmission process cannot affect our algorithm’s effectiveness if the upper bound of the time delays is preknown.

Furthermore, we will focus on the resilient consensus-based algorithm for the EDP considering packet loss and channel noises. What is more, due to the need for information transmission among different agents, the privacy of distributed algorithms is compromised. Additionally, distributed algorithms are vulnerable to cyber-attacks, which pose certain risks. Therefore, our next focus of research will be on resilient distributed algorithms that can withstand cyber-attacks. In addition, the economic dispatch model considering the uncertainties of distributed generations and loads in microgrids is another focus in our future work.