Distributed Algorithm for Economic Dispatch Problem with Separable Losses

TL;DR

This paper presents a distributed algorithm for solving the economic dispatch problem in power systems, incorporating power losses, using convex relaxation and dual decomposition, with robustness to operational changes.

Contribution

It introduces a novel distributed algorithm that handles power losses in economic dispatch, without requiring initialization, and analyzes its behavior under infeasible conditions.

Findings

Algorithm effectively solves the non-convex problem via convex relaxation.

Distributed approach is robust to changing operating conditions.

Behavior under infeasibility is thoroughly analyzed.

Abstract

Economic dispatch problem for a networked power system has been considered. The objective is to minimize the total generation cost while meeting the overall supply-demand balance and generation capacity. In particular, a more practical scenario has been studied by considering the power losses. A non-convex optimization problem has been formulated where the non-convexity comes from the nonlinear equality constraint representing the supply-demand balance with the power losses. It is shown that the optimization problem can be solved using convex relaxation and dual decomposition. A simple distributed algorithm is proposed to solve the optimization problem. Specifically, the proposed algorithm does not require any initialization process and hence robust to various changes in operating condition. In addition, the behavior of the proposed algorithm is analyzed when the problem is infeasible.

Figures

Tables

Table 1. TABLE I : Sampling time used for each node.

Bus	Sampling Time (s)	Bus	Sampling Time (s)
Bus $5$	$0.05$	Bus $11$	$0.03$
Bus $16$	$0.07$	Bus $17$	$0.02$
Bus $21$	$0.04$

Equations

$$\underset{x_1,\dots ,x_N}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$x_i\in {\mathcal{X}}_i,\forall i\in \mathcal{N}$$

$$\varphi (x)={\sum }_{i=1}^N{\varphi }_i(x_i),$$

$$\frac{d{\varphi }_i(x_i)}{d{x}_i}<1,\forall i\in \mathcal{N},$$

$${\sum }_{i=1}^N{\underline{x}}_i-{\varphi }_i({\underline{x}}_i)\le {\sum }_{i=1}^N{d}_i\le {\sum }_{i=1}^N{\bar{x}}_i-{\varphi }_i({\bar{x}}_i)$$

$$0\le {\sum }_{i=1}^N{\bar{x}}_i-{\varphi }_i({\bar{x}}_i)-d_i$$

$$\underset{x_1,\dots ,x_N}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$x_i\in {\mathcal{X}}_i,\forall i\in \mathcal{N}$$

$${\mathcal{L}}^r(x,\lambda )=\sum _{i=1}^N{f}_i(x_i)+\lambda \left(d_i-x_i+{\varphi }_i(x_i)\right)=:\sum _{i=1}^N{\mathcal{L}}_i^r(x_i,\lambda )$$

$$g^r(\lambda )=\underset{x\in \mathcal{X}}{\min }{\mathcal{L}}^r(x,\lambda )=\sum _{i=1}^N\underset{x_i\in {\mathcal{X}}_i}{\min }{\mathcal{L}}_i^r(x_i,\lambda )=:\sum _{i=1}^N{g}_i^r(\lambda ).$$

$$v_i(x_i):=\frac{d{f}_i(x_i)}{d{x}_i}\cdot {\left(1-\frac{d{\varphi }_i(x_i)}{d{x}_i}\right)}^{-1}.$$

$${\hat{x}}_i(\lambda ):=\{\begin{array}{ll}{\underline{x}}_i& 0\le \lambda \le v_i({\underline{x}}_i)\\ v_i^{-1}(\lambda )& v_i({\underline{x}}_i)<\lambda <v_i({\bar{x}}_i)\\ {\bar{x}}_i& v_i({\bar{x}}_i)\le \lambda .\end{array}$$

$$v_i(z_1)-v_i(z_2)=\frac{f_i^{\prime }(z_1)\left(1-{\varphi }_i^{\prime }(z_2)\right)-f_i^{\prime }(z_2)\left(1-{\varphi }_i^{\prime }(z_1)\right)}{\left(1-{\varphi }_i^{\prime }(z_1)\right)\cdot \left(1-{\varphi }_i^{\prime }(z_2)\right)}.$$

$$f_i^{\prime }(z_1)\left(1-{\varphi }_i^{\prime }(z_2)\right)-f_i^{\prime }(z_2)\left(1-{\varphi }_i^{\prime }(z_1)\right).$$

$$\frac{f_i^{\prime }(z_1)}{f_i^{\prime }(z_2)}>1\ge \frac{1-{\varphi }_i^{\prime }(z_1)}{1-{\varphi }_i^{\prime }(z_2)}$$

$$\frac{\partial {\mathcal{L}}_i^r(x_i,\lambda )}{\partial x_i}=\frac{d{f}_i(x_i)}{d{x}_i}+\lambda \frac{d{\varphi }_i(x_i)}{d{x}_i}-\lambda =0$$

$$\frac{\partial {\mathcal{L}}_i^r(x_i,\lambda )}{\partial x_i}<f_i^{\prime }({\bar{x}}_i)-\lambda +\lambda {\varphi }_i^{\prime }({\bar{x}}_i)$$

$$\underset{\lambda \in \mathbb{R}}{\max }$$

$$\mathrm{subject}\mathrm{to}$$

$${\hat{x}}_i(\lambda ):=\{\begin{array}{ll}{\underline{x}}_i& \lambda \le v_i({\underline{x}}_i)\\ v_i^{-1}(\lambda )& v_i({\underline{x}}_i)<\lambda <v_i(\bar{x})\\ {\bar{x}}_i& v_i({\bar{x}}_i)\le \lambda \end{array}$$

$$g^m(\lambda )={\sum }_{i=1}^N{\mathcal{L}}_i^r({\hat{x}}_i(\lambda ),\lambda )=:{\sum }_{i=1}^N{g}_i^m(\lambda )$$

$$g^m(\lambda )={\sum }_{i=1}^N{f}_i({\underline{x}}_i)+\lambda (d_i+{\varphi }_i({\underline{x}}_i)-{\underline{x}}_i).$$

$$\underset{\lambda \in \mathbb{R}}{\max }{g}^m(\lambda )={\sum }_{i=1}^N{g}_i^m(\lambda )$$

$$\dot{\lambda }=\frac{d{g}^m(\lambda )}{d\lambda }={\sum }_{i=1}^N\frac{d{g}_i^m(\lambda )}{d\lambda }.$$

$$\frac{d{g}_i^m(\lambda )}{d\lambda }=d_i-{\hat{x}}_i(\lambda )+{\varphi }_i({\hat{x}}_i(\lambda )),\forall i=1,\dots ,N.$$

$$x_i^{\ast }\in {\mathcal{X}}_i,{\lambda }^{\ast }\ge 0,{\sum }_{i=1}^N{d}_i-x_i^{\ast }+{\varphi }_i(x_i^{\ast })\le 0,$$

$${\lambda }^{\ast }\left({\sum }_{i=1}^N{d}_i-x_i^{\ast }+{\varphi }_i(x_i^{\ast })\right)=0,$$

$$x^{\ast }\in \underset{x\in \mathcal{X}}{\mathrm{argmin}}{\mathcal{L}}^r(x,{\lambda }^{\ast }),$$

$${\sum }_{i=1}^N{d}_i-{\hat{x}}_i({\lambda }^{\ast })+{\varphi }_i({\hat{x}}_i({\lambda }^{\ast }))=0.$$

$${\hat{x}}_i({\lambda }^{\ast })={\hat{x}}_i(0)={\mathrm{argmin}}_{x_i\in {\mathcal{X}}_i}{\mathcal{L}}^r(x_i,0).$$

$${\dot{\lambda }}_i(t)$$

$$x_i(t)$$

Full text

I Introduction

One of the fundamental problems which arise in the operation of the power system is to balance the overall energy demand with generation. In particular, finding the optimal generation has been an important problem due to the socioeconomic impacts of the power system in modern society. The problem of finding the optimal generation is termed as the economic dispatch problem (EDP) [1]. The EDP is often formulated as an optimization problem to decide the power generation of each generator subject to various constraints while minimizing the generation cost. Constraints must include the overall supply and demand balance, while additional constraints such as local generation capacity are often imposed.

The EDP has been studied extensively for the past several decades. Early works were focused on developing centralized algorithms for the EDP. For example, numerical methods [2] and Lagrangian relaxation [3] are developed to solve the EDP. However, the power network is growing with the introduction of the smart grid and distributed energy resources. Hence, there have been significant efforts in recent years to develop a distributed algorithm to solve the EDP due to its scalability and resiliency. Various discrete-time algorithms have been proposed in the literature to solve the EDP in a distributed manner [4, 5, 6]. However, most of these works are not suitable for plug-and-play operation due to requiring an initialization process [4, 5] or decaying step sizes [6].

On the other hand, continuous-time algorithms are also studied due to the ease of applying classical results on the stability of nonlinear systems [7, 8, 9, 10]. In [7], authors considered the EDP with power distribution, but it requires an initialization process and had no capacity constraints. Initialization-free algorithms are proposed in [8, 9] which employed dynamic average consensus and in [10] using dual decomposition and strong coupling.

Most of the works for the EDP mentioned so far considered an ideal scenario where there are no losses in the power system. However, there are various forms of losses occurring in the power network which are significant to the operation of the power system. For example, losses of an electrical generator such as copper losses or core losses are up to $10\%$ of generation depending on the operating condition [11]. In addition, transmission and distribution of the power also cause losses, further lowering the overall efficiency. Therefore, it is important to solve the EDP considering the losses of the system.

Solutions to the EDP with power losses mainly have been developed as a centralized algorithm. For instance, numerical methods are proposed in [12, 2]. For the distributed algorithm, there are only a few works which studied the EDP with power losses. Authors of [13] propose a distributed algorithm considering the transmission losses. However, it is assumed that the power losses can be computed at each iteration and the power losses were not dependent on the decision variables. On the other hand, [14] considered the power losses which depend quadratically on the power generation of each generator. However, it requires an initialization process which is not suitable for plug-and-play operation.

In this work, we propose a continuous-time algorithm which solves the EDP with power losses in a distributed manner. The EDP with power losses are formulated as a non-convex optimization problem with the assumption that the power losses are separable. Despite the non-convexity, it is shown that an optimal solution can be recovered using convex relaxation under mild assumptions. Proposed algorithm does not require any initialization process, thus allowing the plug-and-play operation. In particular, the proposed algorithm is robust to changes such as change of demands or network topology.

The trade-off for having a robust algorithm is that the obtained solution is suboptimal. However, it is shown that with sufficiently high gain, an optimal solution is recovered. Finally, behavior of the algorithm is analyzed when the problem is infeasible.

Notation: For vectors $x_{i}\in{\mathbb{R}}^{n}$ with $i=1,\ldots,N$, $[x_{1}^{T},\ldots,x_{N}^{T}]^{T}$ is denoted by $[x_{1};\cdots;x_{N}]\in{\mathbb{R}}^{\bar{n}}$ where $\bar{n}:=Nn$. An undirected graph is defined as ${\mathcal{G}}=({\mathcal{N}},{\mathcal{E}})$ where ${\mathcal{N}}=\{1,\ldots,N\}$ is the node set and ${\mathcal{E}}\subseteq{\mathcal{N}}\times{\mathcal{N}}$ is the edge set. The Laplacian matrix $L=[l_{ij}]\in{\mathbb{R}}^{N\times N}$ is defined as $l_{ij}:=-1$ if $(j,i)\in{\mathcal{E}}$ and $l_{ij}:=0$ otherwise for $i\neq j$, and $l_{ii}:=-\sum_{j\neq i}l_{ij}$. Eigenvalues of $L$ is denoted as $0=\sigma_{1}(L)\leq\ldots\leq\sigma_{N}(L)$. Given a set ${\mathcal{X}}\subset{\mathbb{R}}^{n}$, let $|x|_{\mathcal{X}}:=\inf_{z\in{\mathcal{X}}}|x-z|$. We denote a set of continuously differentiable functions as ${\mathcal{C}}^{1}$. Derivative of a function $f(x)$ is denoted as $df/dx$ or $f^{\prime}$.

II Problem Formulation

Consider the power network with $N$ nodes in the system. Then, the EDP with power losses can be formulated as the following optimization problem:

[TABLE]

where $x_{i}\in{\mathbb{R}}$ is the power generation (before losses) of the node $i$, $d_{i}\in{\mathbb{R}}$ is the power demand, $f_{i}(x_{i}):{\mathbb{R}}\rightarrow{\mathbb{R}}$ is the local cost function, $x:=[x_{1};\ldots;x_{N}]\in{\mathbb{R}}^{N}$, and $\phi(x):{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}$ represents the power losses. The set ${\mathcal{X}}_{i}:=[\underline{x}_{i},\bar{x}_{i}]$ is a nonempty closed interval where $\underline{x}_{i}$ and $\bar{x}_{i}$ are the minimum and maximum generation of the node $i$ respectively. The objective of the problem (1) is to minimize the generation cost (1a) subject to overall supply and demand balance considering the power losses (1b) and generation capacity constraints (1c).

Let ${\mathcal{X}}:={\mathcal{X}}_{1}\times\cdots\times{\mathcal{X}}_{N}$, $\bar{x}:=[\bar{x}_{1};\ldots;\bar{x}_{N}]$, $d:=[d_{1};\ldots;d_{N}]$ and $\underline{x}:=[\underline{x}_{1};\ldots;\underline{x}_{N}]$. It is supposed that the information such as $f_{i},x_{i},d_{i}$, and ${\mathcal{X}}_{i}$ is private to each node and is not shared with its neighbors.

In this work, we suppose that the loss is separable, i.e.,

[TABLE]

where $\phi_{i}(x_{i}):{\mathbb{R}}\rightarrow{\mathbb{R}}$ is a nonlinear function.

The loss given by (2) includes various forms of losses. For example, it models losses of each generator such as copper losses or mechanical losses [15] where the separability assumption is naturally satisfied. In addition, (2) also includes simplified model for the transmission losses. Separable model for the transmission losses are also employed in previous works, e.g., see [12] and [14].

Note that the optimization problem (1) is not a convex optimization problem due to the nonlinear equality constraint (1b).

Assumption 1

The local cost function $f_{i}$ and loss function $\phi_{i}$ are ${\mathcal{C}}^{1}$, $f_{i}$ is strictly convex and $\phi_{i}$ is convex over ${\mathcal{X}}_{i}$ for $i=1,\dots,N$. Moreover, $\phi_{i}$ satisfies

[TABLE]

for all $\underline{x}_{i}\leq x_{i}\leq\bar{x}_{i}$. $\Diamond$

Inequality (3) of Assumption 1 implies that the incremental loss of each node cannot exceed the incremental generation, which is reasonable. In what follows, we give a necessary and sufficient condition for the feasibility of the optimization problem (1).

Lemma 1

Suppose that Assumption 1 holds. Then,

[TABLE]

holds if and only if the optimization problem (1) is feasible. $\Diamond$

Proof:

Let $D(x):=\sum_{i=1}^{N}x_{i}-\phi_{i}(x_{i})-d_{i}$. Suppose that the problem (1) is feasible. Then, there exists a $z^{*}:=[z_{1}^{*};\ldots;z_{N}^{*}]\in{\mathbb{R}}^{N}$ such that $z_{i}^{*}\in{\mathcal{X}}_{i}$ and $D(z^{*})=0$. Moreover, it follows from Assumption 1 that $\partial D(x)/\partial x_{i}=1-d\phi_{i}(x_{i})/dx_{i}>0$ for all $x_{i}\in{\mathcal{X}}_{i}$. Therefore, $D(x)$ is strictly increasing in each argument. First, we will show that

[TABLE]

holds. If $z^{*}=\bar{x}$, then (5) holds with an equality. If $z^{*}\neq\bar{x}$, then there exists an index $i\in{\mathcal{N}}$ such that $\underline{x}_{i}\leq z_{i}^{*}<\bar{x}_{i}$ since $z^{*}$ is a feasible solution. Therefore, $0=D(z^{*})<D(\bar{x})$ holds since $D(\cdot)$ is strictly increasing in each argument, proving (5). In a similar manner, it can be shown that $D(\underline{x})\leq D(z^{*})$ as well.

Conversely, suppose (4) holds. Consider $D(\underline{x}+\alpha\cdot c)$, where $c\in[0,1]$ and $\alpha:=[\bar{x}_{1}-\underline{x}_{1};\ldots;\bar{x}_{N}-\underline{x}_{N}]\in{\mathbb{R}}^{N}$. Then, $D(\underline{x}+\alpha\cdot c)$ is a continuous function. In addition, (4) can be used to obtain $D(\underline{x})\leq 0\leq D(\bar{x}).$ Therefore, it follows from the intermediate value theorem that there exists a $c^{*}\in[0,1]$ such that $D(\underline{x}+\alpha c^{*})=0$. Hence, the problem (1) is feasible with $\underline{x}+\alpha c^{*}$ as a solution. ∎

Remark 1

It follows from the proof of Lemma 1 that (4) is a sufficient condition for feasibility regardless of (3). However, (4) is not a necessary condition if (3) does not hold. For example, suppose $\sum_{i=1}^{N}\underline{x}_{i}-\phi_{i}(\underline{x}_{i})>\sum_{i=1}^{N}d_{i}$ such that (4) does not hold. Nevertheless, the problem (1) may still be feasible if $d\phi_{i}(x_{i})/dx_{i}>1$. In particular, if one loses more power as $x_{i}$ is increased, a feasible solution may exist. However, if such cases are not allowed (e.g., by assuming (3)), then (4) is indeed a necessary and sufficient condition for the feasibility. $\Diamond$

III A Centralized Solution

In order to solve the non-convex optimization problem, the following assumption is made.

Assumption 2

$df_{i}(x_{i})/dx_{i}>0$* for all $x_{i}\in{\mathcal{X}}_{i}$. $\Diamond$*

Assumption 2 is easily satisfied in practical scenarios. For instance, it is common to assume that the cost function is given by a quadratic function $f_{i}(x_{i})=a_{i}+b_{i}x_{i}+c_{i}x_{i}^{2}$ where $b_{i},c_{i}>0$ and $\underline{x}_{i}\geq 0$. In such case, Assumption 2 holds.

The optimization problem (1) will be relaxed into the following convex optimization problem:

[TABLE]

which will be called as the relaxed problem. Note that the relaxed problem (6) is a convex optimization problem since the equality constraint (1b) is relaxed into an inequality constraint as (6b). Nevertheless, it will be shown that an optimal solution of (1) is obtained by solving (6).

For the optimization problem (6), define the Lagrangian function ${\mathcal{L}}^{r}(x,\lambda):{\mathcal{X}}\times{\mathbb{R}}\rightarrow{\mathbb{R}}$ as

[TABLE]

where $\lambda\in{\mathbb{R}}$ is the dual variable. Then, the dual function $g^{r}(\lambda)$ for (6) can be written as

[TABLE]

The following lemma gives the expression for $g_{i}^{r}(\lambda)$.

Lemma 2

Suppose that Assumptions 1 and 2 hold. Let $v_{i}:{\mathcal{X}}_{i}\rightarrow{\mathbb{R}}$ be defined as

[TABLE]

Then, $v_{i}(x_{i})$ is a strictly increasing function for $\underline{x}_{i}\leq x_{i}\leq\bar{x}_{i}$. Moreover, suppose $\lambda\geq 0$ and let $\hat{x}_{i}(\lambda)$ be

[TABLE]

Then, $\hat{x}_{i}(\lambda)$ is the unique minimizer of ${\mathcal{L}}^{r}_{i}(x_{i},\lambda)$, i.e., $\hat{x}_{i}(\lambda)=\operatorname*{argmin}_{x_{i}\in{\mathcal{X}}_{i}}{\mathcal{L}}^{r}_{i}(x_{i},\lambda).$ $\Diamond$

Proof:

In order to show $v_{i}(x_{i})$ is a strictly increasing function for $\underline{x}_{i}\leq x_{i}\leq\bar{x}_{i}$, let $z_{1},z_{2}\in{\mathcal{X}}_{i}$ such that $z_{1}>z_{2}$. Then, it follows from (7) that

[TABLE]

Since $0<1-\phi_{i}^{\prime}(z_{k})$ holds for $k=1,2$ by (3), it is equivalent to show the strict positivity of

[TABLE]

From Assumptions 1 and 2, it holds that

[TABLE]

which proves $v_{i}(x_{i})$ is strictly increasing. Consequently, $v_{i}^{-1}(\lambda)$ is well-defined for $v_{i}(\underline{x}_{i})\leq\lambda\leq v_{i}(\bar{x}_{i})$.

Next, suppose $\lambda\geq 0$. Then, it follows that ${\mathcal{L}}^{r}_{i}(x_{i},\lambda)$ is a strictly convex function in $x_{i}$ for any fixed $\lambda\geq 0$. Hence, it has a unique minimum and the minimum of ${\mathcal{L}}_{i}^{r}(x_{i},\lambda)$ is obtained when

[TABLE]

which is equivalent to $v_{i}(x_{i})=\lambda$. Therefore, $\arg\!\min_{x_{i}\in{\mathcal{X}}_{i}}{\mathcal{L}}^{r}_{i}(x_{i},\lambda)=v_{i}^{-1}(\lambda)$ holds for $\underline{x}_{i}\leq v_{i}^{-1}(\lambda)\leq\bar{x}_{i}$ which becomes $v_{i}(\underline{x}_{i})\leq\lambda\leq v_{i}(\bar{x}_{i})$ since $v_{i}(x_{i})$ is strictly increasing.

Finally, let $\lambda$ be a fixed scalar such that $\lambda>v_{i}(\bar{x}_{i})$ holds. Then for any $\underline{x}_{i}\leq x_{i}<\bar{x}_{i}$, it follows that

[TABLE]

where the first strict inequality follows since $f_{i}$ is strictly convex. Since ${\mathcal{L}}_{i}^{r}(x_{i},\lambda)$ is strictly decreasing, its minimum is obtained at $x_{i}=\bar{x}_{i}$. The case when $0\leq\lambda<v_{i}(\underline{x}_{i})$ can be proven in a similar manner. ∎

Remark 2

If $p_{i}^{*}:=\bar{x}_{i}=\underline{x}_{i}$, (e.g., node $i$ generates a fixed amount of power or no power at all if $p_{i}^{*}=\bar{x}_{i}=\underline{x}_{i}=0$), then $\hat{x}_{i}(\lambda)=p_{i}^{*}$. This is consistent with Lemma 2. In particular, $\hat{x}_{i}(\lambda)=\operatorname*{argmin}_{x_{i}\in{\mathcal{X}}_{i}}{\mathcal{L}}^{r}_{i}(x_{i},\lambda)$ holds. $\Diamond$

The dual problem of the relaxed problem (6) becomes

[TABLE]

where $\lambda$ is constrained to be non-negative since (6b) is an inequality constraint.

In order to solve the constrained optimization problem such as (9) using continuous-time algorithms, discontinuous vector fields are often employed [9, 16] to constrain variables to the feasible set. However, such methods are harder to implement due to the discontinuity. Different from these approaches, we simply extend the domain of $\hat{x}_{i}(\lambda)$ as below:

[TABLE]

where we have defined $\hat{x}_{i}(\lambda)=\underline{x}_{i}$ for all $\lambda<0$. Accordingly, define the modified dual function as $g^{m}(\lambda):=\sum_{i=1}^{N}{\mathcal{L}}_{i}^{r}(\hat{x}_{i}(\lambda),\lambda)$. Then, the following result holds.

Lemma 3

Suppose that Assumptions 1 and 2 hold. Then, the modified dual function given by

[TABLE]

is ${\mathcal{C}}^{1}$ and concave for all $\lambda\in{\mathbb{R}}$. $\Diamond$

Proof:

For $\lambda>0$, it follows from Lemma 2 that $\hat{x}_{i}(\lambda)=\operatorname*{argmin}_{x_{i}\in{\mathcal{X}}_{i}}{\mathcal{L}}^{r}_{i}(x_{i},\lambda).$ Therefore, [17, Prop.7.1.1] states that $g^{r}(\lambda)$ is concave and ${\mathcal{C}}^{1}$. Since $g^{m}(\lambda)=g^{r}(\lambda)$ for $\lambda>0$, $g^{m}(\lambda)$ is also concave and ${\mathcal{C}}^{1}$ for $\lambda>0$. For $\lambda<0$, it follows from the definition of $\hat{x}_{i}(\lambda)$ that

[TABLE]

Thus, it is obvious that $g^{m}(\lambda)\in{\mathcal{C}}^{1}$ for $\lambda<0$. Moreover, it can be also verified that $g^{m}(\lambda)$ is differentiable at $\lambda=0$. It is left to show $g^{m}(\lambda)$ is concave for $\lambda\leq 0$. However, this follows directly since $g^{m}(\lambda)\in{\mathcal{C}}^{1}$ and it is a linear function of $\lambda$ for $\lambda\leq 0$. ∎

Now, instead of (9), consider the unconstrained optimization problem

[TABLE]

with the gradient ascent algorithm given by

[TABLE]

Using the result of [17, Prop.7.1.1], it can be verified that the derivative of $g_{i}^{m}(\lambda)$ becomes

[TABLE]

Note it follows from Assumption 1 and (10) that $dg_{i}^{m}(\lambda)/d\lambda$ (and hence $dg^{m}(\lambda)/d\lambda$) is monotonically decreasing, uniformly bounded and uniformly continuous.

Before presenting the centralized solution, recall the following result customized from [17, Prop. 6.1.5].

Lemma 4

Let $x^{*}:=[x_{1}^{*};\ldots;x_{N}^{*}]\in{\mathbb{R}}^{N}$ and $\lambda^{*}\in{\mathbb{R}}$. Then, the pair $(x^{*},\lambda^{*})$ satisfies

[TABLE]

if and only if $(x^{*},\lambda^{*})$ is an optimal solution-geometric multiplier pair of (6). $\Diamond$

Now, it will be shown that an optimal solution of (1) can be obtained from (12).

Theorem 1

Suppose that Assumptions 1 and 2 hold and that the optimization problem (1) is feasible (i.e., (4) holds). Consider the gradient ascent algorithm given by (12). Then, $\lim_{t\rightarrow\infty}\lambda(t)=\lambda^{*}$ where $\lambda^{*}$ is an optimal solution of (11). Moreover, $\hat{x}_{i}(\lambda^{*})$ is an optimal solution of the optimization problem (1). $\Diamond$

Proof:

From Lemma 3, it follows that (12) is a gradient ascent algorithm for the concave function $g^{m}(\lambda)$. Hence, it can be easily shown that $\lambda(t)$ converges to an optimal solution of (11) using (4) and (10).

It is left to show $\hat{x}(\lambda^{*}):=[\hat{x}_{1}(\lambda^{*});\ldots;\hat{x}_{N}(\lambda^{*})]\in{\mathbb{R}}^{N}$ is an optimal solution to the problem (1). For this, it will be shown that the pair $(\hat{x}(\lambda^{*}),\lambda^{*})$ satisfies the optimality conditions for the problem (6) provided in Lemma 4 while inequality constraint (6b) is satisfied with an equality.

From the first order optimality condition for (11), it follows that

[TABLE]

Moreover, $\hat{x}_{i}(\lambda^{*})\in{\mathcal{X}}_{i}$ by the definition. Therefore, $\hat{x}_{i}(\lambda^{*})$ is a feasible solution to (6). Now, consider the case when $\lambda^{*}\geq 0$. Then, it follows from Lemma 2 that $\hat{x}_{i}(\lambda^{*})=\operatorname*{argmin}_{x_{i}\in{\mathcal{X}}_{i}}{\mathcal{L}}^{r}_{i}(x_{i},\lambda^{*})$. Therefore, $\hat{x}_{i}(\lambda^{*})$ is an optimal solution of (6) due to Lemma 4. If $\lambda^{*}<0$, then it follows from (10) that

[TABLE]

Hence, we can conclude $(\hat{x}(0),0)$ is an optimal solution-geometric multiplier pair of (6) using Lemma 4. Consequently, $\hat{x}_{i}(\lambda^{*})$ is an optimal solution to (6).

Finally, it follows from (14) that $\hat{x}_{i}(\lambda^{*})$ satisfies constraint (6b) with an equality. Therefore, $\hat{x}_{i}(\lambda^{*})$ is an optimal solution of the problem (1). ∎

From Theorem 1, it can be seen that the optimization problem can be solved using (12). In particular, the optimal generation for each node is obtained using (10).

IV A Distributed Solution

In this section, a distributed algorithm for solving the dual problem (11) (and hence the primal problem (1)) is proposed. Suppose that the node $i$ runs

[TABLE]

where $\lambda_{i}\in{\mathbb{R}}$ is the estimate of the dual variable by the node $i$, $x_{i}(t)$ is the power generation of the node $i$ (at time $t$), $k>0$ is the coupling gain and ${\mathcal{N}}_{i}:=\{j\in{\mathcal{N}}\mid(j,i)\in{\mathcal{E}}\}$. The proposed algorithm (15) is an extension of [10] to the EDP with power losses. For the distributed algorithm, we make the following assumption.

Assumption 3

Communication graph is undirected and connected. $\Diamond$

Note (15) is a distributed algorithm as $dg^{m}_{i}/d\lambda$ can be computed by the node $i$ only using the local information. Moreover, only the estimate of the dual variable is communicated between agents and no private information such as $d_{i}$ or $f_{i}(\cdot)$ are exchanged.

Let $\bm{\lambda}:=[\lambda_{1};\ldots;\lambda_{N}]\in{\mathbb{R}}^{N}$ be the stack of $\lambda_{i}$ and $G(\bm{\lambda}):=[dg^{m}_{1}(\lambda_{1})/d\lambda;\ldots;dg^{m}_{N}(\lambda_{N})/d\lambda]$. Then (15a) can be written as

[TABLE]

where $L\in{\mathbb{R}}^{N\times N}$ is a symmetric Laplacian matrix. Denoting $1_{N}:=[1;\ldots;1]\in{\mathbb{R}}^{N}$, it follows that there exists a matrix $W=[(1/N)1_{N}^{T};R^{T}]\in{\mathbb{R}}^{N\times N}$ and $W^{-1}=[1_{N},Q]$ such that $WLW^{-1}=\mathrm{diag}(0,\sigma_{2}(L),\ldots,\sigma_{N}(L))$ where $R\in{\mathbb{R}}^{N\times(N-1)}$ and $Q\in{\mathbb{R}}^{N\times(N-1)}$. Moreover, it can be checked that $|Q|=\sqrt{N}$ and $|R|=1/\sqrt{N}$ [18].

Now apply the following coordinate transformation

[TABLE]

where $\bar{\xi}\in{\mathbb{R}}$ and $\tilde{\xi}\in{\mathbb{R}}^{N-1}$. In addition, it follows that $\bm{\lambda}=W^{-1}\xi$, or $\lambda_{i}=\bar{\xi}+Q_{i}\tilde{\xi}$ where $Q_{i}$ is the $i$-th row of $Q$. Then the system (16) is transformed into

[TABLE]

where $\xi(0)=W\bm{\lambda}(0)$ and $\tilde{g}(\bar{\xi},\tilde{\xi}):=\sum_{i=1}^{N}({dg_{i}^{m}}/{d\lambda})(\bar{\xi}+Q_{i}\tilde{\xi})-({dg_{i}^{m}}/{d\lambda})(\bar{\xi})$. Convergence of the proposed algorithm is stated below.

Theorem 2

Consider the distributed algorithm (15). Suppose that Assumptions 1, 2, and 3 hold. Also assume that the optimization problem (1) is feasible (i.e., (4) holds). Then, for any $k>0$, the solution of (15) converges to a point and satisfies

[TABLE]

for any initial conditions $\lambda_{i}(0)\in{\mathbb{R}}$. $\Diamond$

Proof:

From (18b), it follows that $\tilde{\xi}(t)$ is bounded since $R^{T}LQ$ is positive definite and $G(\cdot)$ is bounded. It can also be verified that $\bar{\xi}(t)$ is bounded using (4) and (10). Since (17) is a linear transformation, it follows that the solution of (15) is bounded.

Now, let $V(\bm{\lambda})=-\sum_{i=1}^{N}g^{m}_{i}(\lambda_{i})+(k/2)\bm{\lambda}^{T}L\bm{\lambda}$ be a candidate function. Then its time derivative becomes

[TABLE]

Thus, LaSalle’s invariance principle can be applied to conclude that $\lambda_{i}(t)$ approaches to the set $E:=\{\bm{\lambda}\mid\dot{V}(\bm{\lambda})=0\}$. However, note that the set $E$ is the set of equilibrium points of (16). Thus, convergence to a point can be obtained by applying Lemma A.3 from [19].

For the feasibility of the converged solution, let $\hat{\lambda}=[\hat{\lambda}_{1};\ldots;\hat{\lambda}_{N}]:=\lim_{t\rightarrow\infty}\bm{\lambda}(t)\in E$. Then, it holds that $G(\hat{\lambda})-kL\hat{\lambda}=0.$ Multiplying $1_{N}^{T}$ from the left, we obtain

[TABLE]

since $1_{N}^{T}L=0$. ∎

Result of Theorem 1 states that for any $k>0$, the algorithm (15) converges to a feasible solution of (1). However, the optimality of the converged solution has not been stated. In what follows, it is shown that the optimality can be recovered using high coupling gain.

Theorem 3

Consider the distributed algorithm (15) and suppose that the assumptions of Theorem 2 hold. Then, for any $\epsilon>0$, there exists $\bar{k}>0$ and a function $T(\bm{\lambda}(0),k)$ such that for all $k>\bar{k}$, it holds that

[TABLE]

where $\hat{x}_{i}(\lambda^{*})$ is an optimal solution of (1). $\Diamond$

Proof:

Let $p_{i}(z):=z-\phi_{i}(z)$ for all $z\in{\mathcal{X}}_{i}$. Then, $p_{i}(z)$ is uniformly continuous and strictly increasing. Therefore, $p_{i}^{-1}(\cdot)$ is also uniformly continuous and strictly increasing. Thus, there exists $\delta_{1}>0$ such that for all $i\in{\mathcal{N}}$,

[TABLE]

and $\delta_{2}>0$ such that

[TABLE]

Also define $M:=\max_{\lambda}|G(\lambda)|$ which exists since $G(\lambda)$ is bounded. Finally, define $\bar{k}:=2M/\sigma_{2}(L)\delta_{2}$.

Let $V(\tilde{\xi})=(1/2)\tilde{\xi}^{T}\tilde{\xi}$ be a candidate Lyapunov function. Then, the time derivative of $V$ along the trajectories of (18b) becomes

[TABLE]

where $\sigma_{2}(L)$ is the second smallest eigenvalue of the $L$. Therefore, for all $k\geq\bar{k}$, it holds that

[TABLE]

for all $t\geq T_{1}(\bm{\lambda}(0),k)$ where

[TABLE]

and $T_{1}(\bm{\lambda}(0),k):=0$ if $|\tilde{\xi}(0)|\leq\delta_{2}/\sqrt{N}$ (or equivalently $|R^{T}\bm{\lambda}(0)|\leq\delta_{2}/\sqrt{N}$). Thus, $|Q_{i}\tilde{\xi}(t)|\leq\sqrt{N}\cdot(\delta_{2}/\sqrt{N})=\delta_{2}.$ Hence, we obtain

[TABLE]

Define $\Lambda_{\delta_{1}}^{*}:=\{w\in{\mathbb{R}}\mid|dg^{m}(w)/d\lambda|\leq 2\delta_{1}/3\}$ and denote it as $\Lambda^{*}_{\delta_{1}}=[\underline{\Lambda},\bar{\Lambda}]$. In fact, since $(dg^{m}/d\lambda)(\lambda)$ is monotonically decreasing, it holds that $(dg^{m}/d\lambda)(\lambda)\geq 2\delta_{1}/3$ for $\lambda\leq\underline{\Lambda}$ and $(dg^{m}/d\lambda)(\lambda)\leq-2\delta_{1}/3$ for $\lambda\geq\bar{\Lambda}$. From here on, suppose that $t\geq T_{1}(\bm{\lambda}(0),k)$. Then, the trajectory of $\bar{\xi}(t)$ approaches to $\Lambda^{*}_{\delta_{1}}$ with the speed of at least $\delta_{1}/(3N)$. For instance, if $\bar{\xi}(t)\geq\bar{\Lambda}$, it follows from (18a) and (19) that

[TABLE]

Hence, $\bar{\xi}(t)$ converges to $\Lambda^{*}_{\delta_{1}}$ in finite time. In order to compute the convergence time, it follows from (4) that

[TABLE]

Therefore, it follows that

[TABLE]

Define $D^{*}$ as D^{*}:=\max\big{(}|\zeta^{*}|_{\Lambda^{*}_{\delta_{1}}},\left|-\zeta^{*}\right|_{\Lambda^{*}_{\delta_{1}}}\big{)}. Then, it holds that $\bar{\xi}(t)\in\Lambda^{*}_{\delta_{1}}$ for all $t\geq T(\bm{\lambda}(0),k)$, where

[TABLE]

Finally, for some fixed $i$ and for all $t\geq T(\bm{\lambda}(0),k)$,

[TABLE]

where the first equality holds since $p_{j}(\hat{x}_{j}(\cdot))$ is an increasing function, and the second equality holds since \sum_{j=1}^{N}p_{j}\big{(}\hat{x}_{j}(\lambda^{*})\big{)}=\sum_{j=1}^{N}d_{j}. Note, we have

[TABLE]

Therefore, (IV) and (IV) implies

[TABLE]

By the definition of $\delta_{1}$, it follows that

[TABLE]

for all $t\geq T(\bm{\lambda}(0),k)$. ∎

An important feature of the proposed algorithm (15) is that it is an initialization-free algorithm and hence allows plug-and-play operation. In particular, it can be seen from Theorem 2 that the proposed algorithm converges to a feasible solution regardless of the initial condition. Therefore, even if some parameters of the optimization problem (1) changes, the solution of (15) converges to a feasible solution of the new problem. Additionally, converged solution is close to an optimal if the coupling gain $k$ is chosen as stated in Theorem 3. To bound the performance uniformly across changes, the coupling gain must be chosen sufficiently large to incorporate all possible cases. For instance, such gain can be found from the worst case scenario by assuming that the network has a known maximum capacity and parameters such as $f_{i}$ and $\phi_{i}$ are from a finite collection. More detailed discussions can be found in [10, Sec. 6.1] or [20].

Remark 3

Using similar arguments as in (IV), it holds that the constraint violation and the objective error $f_{i}(x_{i}(t))-f_{i}(\hat{x}_{i}(\lambda^{*}))$ become small in finite time. In particular, $|\sum_{i=1}^{N}d_{i}-x_{i}(t)+\phi_{i}(x_{i}(t))|\leq\delta_{1}$ and |f_{i}\big{(}x_{i}(\lambda_{i}(t))\big{)}-f_{i}(\hat{x}_{i}(\lambda^{*}))|\leq c_{i} holds for all $t\geq T(\bm{\lambda}(0),k)$ where $c_{i}>0$ is such that $|a-b|\leq\epsilon$ implies $|f_{i}(a)-f_{i}(b)|\leq c_{i}$ for any $a,b\in{\mathcal{X}}_{i}$. Thus, constraint violation and objective error can be made arbitrarily small by reducing $\epsilon$ which leads to higher coupling gain $k$. $\Diamond$

Results of Theorems 2 and Theorem 3 assume that the problem (1) is feasible. The behavior of the proposed algorithm is also analyzed when the problem is infeasible.

Theorem 4

Suppose that Assumptions 1, 2 and 3 hold. Assume that the optimization problem (1) is infeasible. Specifically, suppose that $\textstyle\sum_{i=1}^{N}d_{i}-\bar{x}_{i}+\phi_{i}(\bar{x}_{i})>0$ holds. Then, for all $i\in{\mathcal{N}}$, $\lambda_{i}(t)$ diverges to $+\infty$ and $\lim_{t\rightarrow\infty}\dot{\lambda}_{i}(t)=D_{0},$ where $D_{0}:=(\sum_{i=1}^{N}d_{i}-\bar{x}_{i}+\phi_{i}(\bar{x}_{i}))/N>0$. Similar result also holds in the case of $\textstyle\sum_{i=1}^{N}d_{i}-\underline{x}_{i}+\phi_{i}(\underline{x}_{i})<0.$ $\Diamond$

Proof:

From the proof of Theorem 2, it holds that $\tilde{\xi}(t)$ is bounded. Moreover, we have

[TABLE]

Therefore, it follows that $\lim_{t\rightarrow\infty}\bar{\xi}(t)=+\infty$. Thus,

[TABLE]

and hence $\lambda_{i}(t)$ diverges to $+\infty$.

From (10), it holds that $dg_{i}^{m}(\bar{\xi})/d\lambda$ is a constant for sufficiently large $\bar{\xi}$. Therefore, it can be verified from (18a) and the boundedness of $\tilde{\xi}$ that $\dot{\bar{\xi}}$ satisfies

[TABLE]

In addition, boundedness of $\tilde{\xi}$ and divergence of $\bar{\xi}$ implies

[TABLE]

In particular, $\tilde{\xi}(t)$ converges to a constant value. Therefore,

[TABLE]

which completes the proof. ∎

V Simulation

V-A Robustness to Changes

For the simulation, the continuous-time algorithm (15) is discretized using the forward difference method. In particular, if we denote the sampling period with $T$, then the distributed algorithm (15) can be discretized into

[TABLE]

for $i\in{\mathcal{N}}$, where $q\geq 0$ is integer and $\lambda_{i}^{\rm d}$ is the discretized state. In particular, $\lambda_{i}^{\rm d}(t)$ is a piecewise constant, right continuous signal which is updated every $T$ seconds. Specifically, the value of $\lambda_{i}^{\rm d}(t)$ is held constant until the next update.

Numerical simulation is done with IEEE 30 bus system [21] to verify the proposed algorithm. The local cost function is given by $f_{i}(x_{i})=a_{i}+b_{i}x_{i}+c_{i}x_{i}^{2}$ where $b_{i}$ and $c_{i}$ are strictly positive for $i\in\{1,2,5,8,11,13\}$ which are the buses with a generator. The power demand of each bus satisfies $d_{i}\in[0,94.2]$ and $\sum_{i=1}^{N}d_{i}=283.4$. The loss function is chosen as $\phi_{i}(x_{i})=\alpha_{i}x_{i}^{2}$ where $\alpha_{i}\in[0.0001,0.0007]$ are chosen randomly. It is assumed that each bus is capable of running the proposed algorithm (22) and that two buses connected by a branch can communicate. Coupling gain of $k=40$ is used for the simulation. For the implementation, (22) is used with the sampling time of $0.005$ seconds. We consider the following scenario:

• S1)

Normal operation condition for $0\leq t\leq 10s$.

• S2)

At $t=10s$, demand at bus $5$ is decreased by $20\%$.

• S3)

At $t=20s$, generator at bus $1$ stops generation and leaves the network.

• S4)

At $t=30s$, bus $1$ joins the network again and the maximum generation at bus $8$ increases by $20\%$.

Note in particular that the above scenarios are simulated in one continuous session. Moreover, network topology changes as a node leaves and joins the network during the operation.

Simulation results are shown in Fig. 1. It can be seen that an optimal solution is obtained in a distributed manner despite the changes in operation conditions. At $t=20s$, the problem (1) becomes infeasible due to the lack of generation at bus $1$. Hence, the trajectory of $\lambda_{i}^{\rm d}(t)$ diverges which verifies the result of Theorem 4. As feasibility is recovered at $t=30s$, $\lambda_{i}^{\rm d}(t)$ converges again. Value of the cost function is also shown in Fig. 1(c) and it is seen that the optimal cost is approximately recovered. The trajectory of power mismatch $\sum_{i=1}^{N}d_{i}-x_{i}(t)+\phi_{i}(x_{i}(t))$ is shown in Fig. 1(d). It is observed that the mismatch converges to zero (except for the case when the problem is infeasible) implying supply-demand balance is satisfied.

From the repeated simulations, we also observed that the coupling gain $k$ and the sampling time $T$ have a close relationship for the stability. Specifically, if the gain $k$ is large, the sampling time must be reduced to yield a stable algorithm. For example, it is observed that (22) diverges with $T=0.01$ but converges if $T=0.005$.

V-B Asynchronous Update

The discretized algorithm (22) supposes that all agents use the same sampling time $T$ and that the update is made synchronously across all agents. However, it may be hard to implement the synchronous algorithms in practice due to the distributed nature of the system. Instead, in this section, we consider the case where each agent uses different sampling time $T_{i}$ and that the update is done asynchronously. In particular, let $q\geq 0$ be an integer representing the time index. Then, the distributed algorithm (15) becomes

[TABLE]

where $\lambda_{i}^{\rm d}(t)$ is a piecewise constant, right continuous signal which is updated every $T_{i}$ seconds. Note that the algorithm (V-B) is equivalent to (22) if $T_{i}=T$ for all $i\in{\mathcal{N}}$. However, if sampling time is different between agents, then some agents update more frequently than others.

The same IEEE 30 bus system is used to simulate the algorithm (V-B). For the simulation, gain of $k=20$ is used while sampling time of $0.009s$ is used for all agents except for the ones denoted in Table I. Specifically, we consider the case when some nodes update less frequently. Simulation results are shown in Fig. 2. Trajectories of $\lambda_{i}^{\rm d}(t)$ for selected bus is shown in Fig. 2(a) and Fig. 2(b). It is clearly seen that each variable is updated asynchronously, and that some agents are updated more frequently. Nonetheless, the solution $\lambda_{i}^{\rm d}(t)$ still converges. Additionally, Fig. 2(c) depicts that the converged solution satisfies supply and demand balance.

VI Conclusion and Future Works

The economic dispatch problem with nonlinear, separable power losses has been studied in this paper. Due to the addition of nonlinear loss, the EDP becomes a non-convex optimization problem. However, it has been shown that convex relaxation with dual decomposition can be used to obtain an optimal solution. The distributed algorithm is proposed and it is shown to converge to a feasible solution while an optimal solution is recovered with sufficiently high coupling gain. Specifically, the proposed algorithm does not require any initialization process and converges from any initial condition. Moreover, the behavior of the proposed algorithm is analyzed when the problem is infeasible. Future works include the theoretical analysis of the discretized version of the proposed algorithm.