Distributed Algorithm for Economic Dispatch Problem with Separable Losses
DDistributed Algorithm for Economic Dispatch Problemwith Separable Losses
Seungjoon Lee and Hyungbo Shim
Abstract — Economic dispatch problem for a net-worked power system has been considered. The ob-jective is to minimize the total generation cost whilemeeting the overall supply-demand balance and gen-eration capacity. In particular, a more practical sce-nario has been studied by considering the powerlosses. A non-convex optimization problem has beenformulated where the non-convexity comes from thenonlinear equality constraint representing the supply-demand balance with the power losses. It is shownthat the optimization problem can be solved usingconvex relaxation and dual decomposition. A sim-ple distributed algorithm is proposed to solve theoptimization problem. Specifically, the proposed al-gorithm does not require any initialization processand hence robust to various changes in operatingcondition. In addition, the behavior of the proposedalgorithm is analyzed when the problem is infeasible.
I. Introduction
One of the fundamental problems which arise in theoperation of the power system is to balance the overallenergy demand with generation. In particular, findingthe optimal generation has been an important problemdue to the socioeconomic impacts of the power systemin modern society. The problem of finding the optimalgeneration is termed as the economic dispatch problem(EDP) [1]. The EDP is often formulated as an opti-mization problem to decide the power generation of eachgenerator subject to various constraints while minimizingthe generation cost. Constraints must include the overallsupply and demand balance, while additional constraintssuch as local generation capacity are often imposed.The EDP has been studied extensively for the pastseveral decades. Early works were focused on develop-ing centralized algorithms for the EDP. For example,numerical methods [2] and Lagrangian relaxation [3] aredeveloped to solve the EDP. However, the power networkis growing with the introduction of the smart grid anddistributed energy resources. Hence, there have beensignificant efforts in recent years to develop a distributedalgorithm to solve the EDP due to its scalability andresiliency. Various discrete-time algorithms have beenproposed in the literature to solve the EDP in a dis-tributed manner [4]–[6]. However, most of these works are
This work was supported by the National Research Foundationof Korea (NRF) grant funded by the Korea government (Ministryof Science and ICT) (No. NRF-2017R1E1A1A03070342).S. Lee and H. Shim are with ASRI, Department of Electrical andComputer Engineering, Seoul National University, Seoul, Korea. [email protected], [email protected] not suitable for plug-and-play operation due to requiringan initialization process [4], [5] or decaying step sizes [6].On the other hand, continuous-time algorithms arealso studied due to the ease of applying classical resultson the stability of nonlinear systems [7]–[10]. In [7],authors considered the EDP with power distribution, butit requires an initialization process and had no capacityconstraints. Initialization-free algorithms are proposed in[8], [9] which employed dynamic average consensus andin [10] using dual decomposition and strong coupling.Most of the works for the EDP mentioned so far con-sidered an ideal scenario where there are no losses in thepower system. However, there are various forms of lossesoccurring in the power network which are significant tothe operation of the power system. For example, losses ofan electrical generator such as copper losses or core lossesare up to 10% of generation depending on the operatingcondition [11]. In addition, transmission and distributionof the power also cause losses, further lowering the overallefficiency. Therefore, it is important to solve the EDPconsidering the losses of the system.Solutions to the EDP with power losses mainly havebeen developed as a centralized algorithm. For instance,numerical methods are proposed in [2], [12]. For thedistributed algorithm, there are only a few works whichstudied the EDP with power losses. Authors of [13] pro-pose a distributed algorithm considering the transmissionlosses. However, it is assumed that the power losses canbe computed at each iteration and the power losses werenot dependent on the decision variables. On the otherhand, [14] considered the power losses which dependquadratically on the power generation of each generator.However, it requires an initialization process which is notsuitable for plug-and-play operation.In this work, we propose a continuous-time algorithmwhich solves the EDP with power losses in a distributedmanner. The EDP with power losses are formulated as anon-convex optimization problem with the assumptionthat the power losses are separable. Despite the non-convexity, it is shown that an optimal solution can be re-covered using convex relaxation under mild assumptions.Proposed algorithm does not require any initializationprocess, thus allowing the plug-and-play operation. Inparticular, the proposed algorithm is robust to changessuch as change of demands or network topology. Thetrade-off for having a robust algorithm is that theobtained solution is suboptimal. However, it is shownthat with sufficiently high gain, an optimal solution is a r X i v : . [ c s . S Y ] A p r ecovered. Finally, behavior of the algorithm is analyzedwhen the problem is infeasible. Notation:
For vectors x i ∈ R n with i = 1 , . . . , N ,[ x T , . . . , x TN ] T is denoted by [ x ; · · · ; x N ] ∈ R ¯ n where¯ n := N n . An undirected graph is defined as G = ( N , E )where N = { , . . . , N } is the node set and E ⊆ N × N isthe edge set. The Laplacian matrix L = [ l ij ] ∈ R N × N isdefined as l ij := − j, i ) ∈ E and l ij := 0 otherwise for i = j , and l ii := − P j = i l ij . Eigenvalues of L is denotedas 0 = σ ( L ) ≤ . . . ≤ σ N ( L ). Given a set X ⊂ R n , let | x | X := inf z ∈X | x − z | . We denote a set of continuouslydifferentiable functions as C . Derivative of a function f ( x ) is denoted as df /dx or f . II. Problem Formulation
Consider the power network with N nodes in the sys-tem. Then, the EDP with power losses can be formulatedas the following optimization problem:min x ,...,x N P Ni =1 f i ( x i ) (1a)subject to P Ni =1 d i = P Ni =1 x i − φ ( x ) (1b) x i ∈ X i , ∀ i ∈ N (1c)where x i ∈ R is the power generation (before losses) ofthe node i , d i ∈ R is the power demand, f i ( x i ) : R → R is the local cost function, x := [ x ; . . . ; x N ] ∈ R N ,and φ ( x ) : R N → R represents the power losses. Theset X i := [ x i , ¯ x i ] is a nonempty closed interval where x i and ¯ x i are the minimum and maximum generation of thenode i respectively. The objective of the problem (1) isto minimize the generation cost (1a) subject to overallsupply and demand balance considering the power losses(1b) and generation capacity constraints (1c).Let X := X × · · · × X N , ¯ x := [¯ x ; . . . ; ¯ x N ], d :=[ d ; . . . ; d N ] and x := [ x ; . . . ; x N ]. It is supposed thatthe information such as f i , x i , d i , and X i is private toeach node and is not shared with its neighbors.In this work, we suppose that the loss is separable , i.e., φ ( x ) = P Ni =1 φ i ( x i ) , (2)where φ i ( x i ) : R → R is a nonlinear function. The lossgiven by (2) includes various forms of losses. For example,it models losses of each generator such as copper lossesor mechanical losses [15] where the separability assump-tion is naturally satisfied. In addition, (2) also includessimplified model for the transmission losses. Separablemodel for the transmission losses are also employed inprevious works, e.g., see [12] and [14].Note that the optimization problem (1) is not a con-vex optimization problem due to the nonlinear equalityconstraint (1b). Assumption 1:
The local cost function f i and lossfunction φ i are C , f i is strictly convex and φ i is convexover X i for i = 1 , . . . , N . Moreover, φ i satisfies dφ i ( x i ) dx i < , ∀ i ∈ N , (3) for all x i ≤ x i ≤ ¯ x i . ♦ Inequality (3) of Assumption 1 implies that the incre-mental loss of each node cannot exceed the incrementalgeneration, which is reasonable. In what follows, we givea necessary and sufficient condition for the feasibility ofthe optimization problem (1).
Lemma 1:
Suppose that Assumption 1 holds. Then, P Ni =1 x i − φ i ( x i ) ≤ P Ni =1 d i ≤ P Ni =1 ¯ x i − φ i (¯ x i ) (4)holds if and only if the optimization problem (1) isfeasible. ♦ Proof:
Let D ( x ) := P Ni =1 x i − φ i ( x i ) − d i . Supposethat the problem (1) is feasible. Then, there exists a z ∗ := [ z ∗ ; . . . ; z ∗ N ] ∈ R N such that z ∗ i ∈ X i and D ( z ∗ ) = 0. Moreover, it follows from Assumption 1 that ∂D ( x ) /∂x i = 1 − dφ i ( x i ) /dx i > x i ∈ X i .Therefore, D ( x ) is strictly increasing in each argument.First, we will show that0 ≤ P Ni =1 ¯ x i − φ i (¯ x i ) − d i (5)holds. If z ∗ = ¯ x , then (5) holds with an equality. If z ∗ = ¯ x , then there exists an index i ∈ N such that x i ≤ z ∗ i < ¯ x i since z ∗ is a feasible solution. Therefore,0 = D ( z ∗ ) < D (¯ x ) holds since D ( · ) is strictly increasingin each argument, proving (5). In a similar manner, itcan be shown that D ( x ) ≤ D ( z ∗ ) as well.Conversely, suppose (4) holds. Consider D ( x + α · c ),where c ∈ [0 ,
1] and α := [¯ x − x ; . . . ; ¯ x N − x N ] ∈ R N .Then, D ( x + α · c ) is a continuous function. In addition,(4) can be used to obtain D ( x ) ≤ ≤ D (¯ x ) . Therefore,it follows from the intermediate value theorem that thereexists a c ∗ ∈ [0 ,
1] such that D ( x + αc ∗ ) = 0. Hence, theproblem (1) is feasible with x + αc ∗ as a solution. Remark 1:
It follows from the proof of Lemma 1that (4) is a sufficient condition for feasibility regardlessof (3). However, (4) is not a necessary condition if (3)does not hold. For example, suppose P Ni =1 x i − φ i ( x i ) > P Ni =1 d i such that (4) does not hold. Nevertheless, theproblem (1) may still be feasible if dφ i ( x i ) /dx i >
1. Inparticular, if one loses more power as x i is increased, afeasible solution may exist. However, if such cases arenot allowed (e.g., by assuming (3)), then (4) is indeed anecessary and sufficient condition for the feasibility. ♦ 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 > x i ∈ X i . ♦ Assumption 2 is easily satisfied in practical scenarios.For instance, it is common to assume that the costfunction is given by a quadratic function f i ( x i ) = a i + b i x i + c i x i where b i , c i > x i ≥
0. In such case,Assumption 2 holds.he optimization problem (1) will be relaxed into thefollowing convex optimization problem:min x ,...,x N P Ni =1 f i ( x i ) (6a)subject to P Ni =1 d i − x i + φ i ( x i ) ≤ x i ∈ X i , ∀ i ∈ N (6c)which will be called as the relaxed problem . Note that therelaxed problem (6) is a convex optimization problemsince the equality constraint (1b) is relaxed into aninequality constraint as (6b). Nevertheless, it will beshown that an optimal solution of (1) is obtained bysolving (6).For the optimization problem (6), define the La-grangian function L r ( x, λ ) : X × R → R as L r ( x, λ )= N X i =1 f i ( x i ) + λ ( d i − x i + φ i ( x i )) =: N X i =1 L ri ( x i , λ )where λ ∈ R is the dual variable. Then, the dual function g r ( λ ) for (6) can be written as g r ( λ ) = min x ∈X L r ( x, λ ) = N X i =1 min x i ∈X i L ri ( x i , λ ) =: N X i =1 g ri ( λ ) . The following lemma gives the expression for g ri ( λ ). Lemma 2:
Suppose that Assumptions 1 and 2 hold.Let v i : X i → R be defined as v i ( x i ) := df i ( x i ) dx i · (cid:18) − dφ i ( x i ) dx i (cid:19) − . (7)Then, v i ( x i ) is a strictly increasing function for x i ≤ x i ≤ ¯ x i . Moreover, suppose λ ≥ x i ( λ ) beˆ x i ( λ ) := x i ≤ λ ≤ v i ( x i ) v − i ( λ ) v i ( x i ) < λ < v i (¯ x i )¯ x i v i (¯ x i ) ≤ λ. Then, ˆ x i ( λ ) is the unique minimizer of L ri ( x i , λ ), i.e.,ˆ x i ( λ ) = argmin x i ∈X i L ri ( x i , λ ) . ♦ Proof:
In order to show v i ( x i ) is a strictly increasingfunction for x i ≤ x i ≤ ¯ x i , let z , z ∈ X i such that z >z . Then, it follows from (7) that v i ( z ) − v i ( z ) = f i ( z ) (1 − φ i ( z )) − f i ( z ) (1 − φ i ( z ))(1 − φ i ( z )) · (1 − φ i ( z )) . Since 0 < − φ i ( z k ) holds for k = 1 , f i ( z ) (1 − φ i ( z )) − f i ( z ) (1 − φ i ( z )) . (8)From Assumptions 1 and 2, it holds that f i ( z ) f i ( z ) > ≥ − φ i ( z )1 − φ i ( z )which proves v i ( x i ) is strictly increasing. Consequently, v − i ( λ ) is well-defined for v i ( x i ) ≤ λ ≤ v i (¯ x i ).Next, suppose λ ≥
0. Then, it follows that L ri ( x i , λ )is a strictly convex function in x i for any fixed λ ≥
0. Hence, it has a unique minimum and the minimum of L ri ( x i , λ ) is obtained when ∂ L ri ( x i , λ ) ∂x i = df i ( x i ) dx i + λ dφ i ( x i ) dx i − λ = 0which is equivalent to v i ( x i ) = λ . Therefore,argmin x i ∈X i L ri ( x i , λ ) = v − i ( λ ) holds for x i ≤ v − i ( λ ) ≤ ¯ x i which becomes v i ( x i ) ≤ λ ≤ v i (¯ x i ) since v i ( x i ) isstrictly increasing.Finally, let λ be a fixed scalar such that λ > v i (¯ x i )holds. Then for any x i ≤ x i < ¯ x i , it follows that ∂ L ri ( x i , λ ) ∂x i < f i (¯ x i ) − λ + λφ i (¯ x i ) < f i is strictlyconvex. Since L ri ( x i , λ ) is strictly decreasing, its mini-mum is obtained at x i = ¯ x i . The case when 0 ≤ λ 0. Ac-cordingly, define the modified dual function as g m ( λ ) := P Ni =1 L ri (ˆ x i ( λ ) , λ ). Then, the following result holds. Lemma 3: Suppose that Assumptions 1 and 2 hold.Then, the modified dual function given by g m ( λ ) = P Ni =1 L ri (ˆ x i ( λ ) , λ ) =: P Ni =1 g mi ( λ )is C and concave for all λ ∈ R . ♦ Proof: For λ > 0, it follows from Lemma 2that ˆ x i ( λ ) = argmin x i ∈X i L ri ( x i , λ ) . Therefore, [17,Prop.7.1.1] states that g r ( λ ) is concave and C . Since g m ( λ ) = g r ( λ ) for λ > g m ( λ ) is also concave and C for λ > 0. For λ < 0, it follows from the definition ofˆ x i ( λ ) that g m ( λ ) = P Ni =1 f i ( x i ) + λ ( d i + φ i ( x i ) − x i ) . hus, it is obvious that g m ( λ ) ∈ C for λ < 0. Moreover,it can be also verified that g m ( λ ) is differentiable at λ =0. It is left to show g m ( λ ) is concave for λ ≤ 0. However,this follows directly since g m ( λ ) ∈ C and it is a linearfunction of λ for λ ≤ λ ∈ R g m ( λ ) = P Ni =1 g mi ( λ ) (11)with the gradient ascent algorithm given by˙ λ = dg m ( λ ) dλ = X Ni =1 dg mi ( λ ) dλ . (12)Using the result of [17, Prop.7.1.1], it can be verified thatthe derivative of g mi ( λ ) becomes dg mi ( λ ) dλ = d i − ˆ x i ( λ ) + φ i (ˆ x i ( λ )) , ∀ i = 1 , . . . , N. Note it follows from Assumption 1 and (10) that dg mi ( λ ) /dλ (and hence dg m ( λ ) /dλ ) is monotonically de-creasing, uniformly bounded and uniformly continuous.Before presenting the centralized solution, recall thefollowing result customized from [17, Prop. 6.1.5]. Lemma 4: Let x ∗ := [ x ∗ ; . . . ; x ∗ N ] ∈ R N and λ ∗ ∈ R .Then, the pair ( x ∗ , λ ∗ ) satisfies x ∗ i ∈ X i , λ ∗ ≥ , P Ni =1 d i − x ∗ i + φ i ( x ∗ i ) ≤ ,λ ∗ (cid:16)P Ni =1 d i − x ∗ i + φ i ( x ∗ i ) (cid:17) = 0 ,x ∗ ∈ argmin x ∈X L r ( x, λ ∗ ) , if and only if ( x ∗ , λ ∗ ) is an optimal solution-geometricmultiplier pair of (6). ♦ Now, it will be shown that an optimal solution of (1)can be obtained from (12). Theorem 1: Suppose that Assumptions 1 and 2 holdand that the optimization problem (1) is feasible (i.e.,(4) holds). Consider the gradient ascent algorithm givenby (12). Then, lim t →∞ λ ( t ) = λ ∗ where λ ∗ is an optimalsolution of (11). Moreover, ˆ x i ( λ ∗ ) is an optimal solutionof the optimization problem (1). ♦ Proof: From Lemma 3, it follows that (12) is agradient ascent algorithm for the concave function g m ( λ ).Hence, it can be easily shown that λ ( t ) converges to anoptimal solution of (11) using (4) and (10).It is left to show ˆ x ( λ ∗ ) := [ˆ x ( λ ∗ ); . . . ; ˆ x N ( λ ∗ )] ∈ R N is an optimal solution to the problem (1). For this,it will be shown that the pair (ˆ x ( λ ∗ ) , λ ∗ ) satisfies theoptimality conditions for the problem (6) provided inLemma 4 while inequality constraint (6b) is satisfied withan equality.From the first order optimality condition for (11), itfollows that P Ni =1 d i − ˆ x i ( λ ∗ ) + φ i (ˆ x i ( λ ∗ )) = 0 . (14)Moreover, ˆ x i ( λ ∗ ) ∈ X i by the definition. Therefore,ˆ x i ( λ ∗ ) is a feasible solution to (6). Now, consider the case when λ ∗ ≥ 0. Then, it follows from Lemma 2 thatˆ x i ( λ ∗ ) = argmin x i ∈X i L ri ( x i , λ ∗ ). Therefore, ˆ x i ( λ ∗ ) is anoptimal solution of (6) due to Lemma 4. If λ ∗ < 0, thenit follows from (10) thatˆ x i ( λ ∗ ) = ˆ x i (0) = argmin x i ∈X i L r ( x i , . Hence, we can conclude (ˆ x (0) , 0) is an optimal solution-geometric multiplier pair of (6) using Lemma 4. Conse-quently, ˆ x i ( λ ∗ ) is an optimal solution to (6).Finally, it follows from (14) that ˆ x i ( λ ∗ ) satisfies con-straint (6b) with an equality. Therefore, ˆ x i ( λ ∗ ) is anoptimal solution of the problem (1).From Theorem 1, it can be seen that the optimizationproblem can be solved using (12). In particular, theoptimal generation for each node is obtained using (10). IV. A Distributed Solution In this section, a distributed algorithm for solving thedual problem (11) (and hence the primal problem (1)) isproposed. Suppose that the node i runs˙ λ i ( t ) = dg mi dλ ( λ i ( t )) + k X j ∈N i ( λ j ( t ) − λ i ( t )) (15a) x i ( t ) = ˆ x i ( λ i ( t )) (15b)where λ i ∈ R is the estimate of the dual variable by thenode i , x i ( t ) is the power generation of the node i (attime t ), k > N i := { j ∈ N | ( j, i ) ∈ E} . The proposed algorithm (15) is an extensionof [10] to the EDP with power losses. For the distributedalgorithm, we make the following assumption. Assumption 3: Communication graph is undirectedand connected. ♦ Note (15) is a distributed algorithm as dg mi /dλ can becomputed by the node i only using the local information.Moreover, only the estimate of the dual variable is com-municated between agents and no private informationsuch as d i or f i ( · ) are exchanged.Let λ := [ λ ; . . . ; λ N ] ∈ R N be the stack of λ i and G ( λ ) := [ dg m ( λ ) /dλ ; . . . ; dg mN ( λ N ) /dλ ]. Then (15a) canbe written as ˙ λ = G ( λ ) − kL λ (16)where L ∈ R N × N is a symmetric Laplacian matrix.Denoting 1 N := [1; . . . ; 1] ∈ R N , it follows that thereexists a matrix W = [(1 /N )1 TN ; R T ] ∈ R N × N and W − =[1 N , Q ] such that W LW − = diag(0 , σ ( L ) , . . . , σ N ( L ))where R ∈ R N × ( N − and Q ∈ R N × ( N − . Moreover, itcan be checked that | Q | = √ N and | R | = 1 / √ N [18].Now apply the following coordinate transformation ξ := (cid:20) ¯ ξ ˜ ξ (cid:21) = W λ = (cid:20) N TN R T (cid:21) λ (17)where ¯ ξ ∈ R and ˜ ξ ∈ R N − . In addition, it follows that λ = W − ξ , or λ i = ¯ ξ + Q i ˜ ξ where Q i is the i -th row of . Then the system (16) is transformed into˙¯ ξ = 1 N N X i =1 dg mi dλ ( ¯ ξ ) + 1 N ˜ g (cid:0) ¯ ξ, ˜ ξ (cid:1) , (18a)˙˜ ξ = − kR T LQ ˜ ξ + R T G (cid:0) N ¯ ξ + Q ˜ ξ (cid:1) , (18b)where ξ (0) = W λ (0) and ˜ g ( ¯ ξ, ˜ ξ ) := P Ni =1 ( dg mi /dλ )( ¯ ξ + Q i ˜ ξ ) − ( dg mi /dλ )( ¯ ξ ). Convergence of the proposed algo-rithm is stated below. Theorem 2: Consider the distributed algorithm (15).Suppose that Assumptions 1, 2, and 3 hold. Also assumethat the optimization problem (1) is feasible (i.e., (4)holds). Then, for any k > 0, the solution of (15)converges to a point and satisfieslim t →∞ P Ni =1 x i ( λ i ( t )) − φ i (cid:0) x i ( λ i ( t )) (cid:1) = P Ni =1 d i for any initial conditions λ i (0) ∈ R . ♦ Proof: From (18b), it follows that ˜ ξ ( t ) is boundedsince R T LQ is positive definite and G ( · ) is bounded. Itcan also be verified that ¯ ξ ( t ) is bounded using (4) and(10). Since (17) is a linear transformation, it follows thatthe solution of (15) is bounded.Now, let V ( λ ) = − P Ni =1 g mi ( λ i ) + ( k/ λ T L λ be acandidate function. Then its time derivative becomes˙ V = ( − G ( λ ) + kL λ ) T ˙ λ = − | G ( λ ) − kL λ | ≤ . Thus, LaSalle’s invariance principle can be applied toconclude that λ i ( t ) approaches to the set E := { λ | ˙ V ( λ ) = 0 } . However, note that the set E is the set ofequilibrium points of (16). Thus, convergence to a pointcan be obtained by applying Lemma A.3 from [19].For the feasibility of the converged solution, let ˆ λ =[ˆ λ ; . . . ; ˆ λ N ] := lim t →∞ λ ( t ) ∈ E . Then, it holds that G (ˆ λ ) − kL ˆ λ = 0 . Multiplying 1 TN from the left, we obtain1 TN G (ˆ λ ) = P Ni =1 d i − ˆ x i (ˆ λ i ) + φ i (cid:0) ˆ x i (ˆ λ i ) (cid:1) = 0since 1 TN L = 0.Result of Theorem 1 states that for any k > 0, thealgorithm (15) converges to a feasible solution of (1).However, the optimality of the converged solution hasnot been stated. In what follows, it is shown that theoptimality 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 (cid:15) > 0, there exists ¯ k > T ( λ (0) , k ) such that for all k > ¯ k , it holds that | ˆ x i ( λ i ( t )) − ˆ x i ( λ ∗ ) | ≤ (cid:15), ∀ i ∈ N , ∀ t ≥ T ( λ (0) , k ) , where ˆ x i ( λ ∗ ) is an optimal solution of (1). ♦ Proof: Let p i ( z ) := z − φ i ( z ) for all z ∈ X i .Then, p i ( z ) is uniformly continuous and strictly increas-ing. Therefore, p − i ( · ) is also uniformly continuous andstrictly increasing. Thus, there exists δ > i ∈ N , | p i ( a ) − p i ( b ) | ≤ δ = ⇒ | a − b | ≤ (cid:15) and δ > | a − b | ≤ δ = ⇒ (cid:12)(cid:12)(cid:12)(cid:12) dg mi ( a ) dλ − dg mi ( b ) dλ (cid:12)(cid:12)(cid:12)(cid:12) ≤ δ N . Also define M := max λ | G ( λ ) | which exists since G ( λ ) isbounded. Finally, define ¯ k := 2 M/σ ( L ) δ .Let V ( ˜ ξ ) = (1 / 2) ˜ ξ T ˜ ξ be a candidate Lyapunov func-tion. Then, the time derivative of V along the trajectoriesof (18b) becomes˙ V ≤ − kσ ( L ) | ˜ ξ | + | R T || G (1 N ¯ ξ + Q ˜ ξ ) || ˜ ξ |≤ − kσ ( L )2 | ˜ ξ | , ∀| ˜ ξ | ≥ Mkσ ( L ) √ N where σ ( L ) is the second smallest eigenvalue of the L .Therefore, for all k ≥ ¯ k , it holds that | ˜ ξ ( t ) | ≤ Mkσ ( L ) √ N ≤ δ √ N for all t ≥ T ( λ (0) , k ) where T ( λ (0) , k ) = ln √ N | R T λ (0) | δ ! · kσ ( L ) , and T ( λ (0) , k ) := 0 if | ˜ ξ (0) | ≤ δ / √ N (or equivalently | R T λ (0) | ≤ δ / √ N ). Thus, | Q i ˜ ξ ( t ) | ≤ √ N · ( δ / √ N ) = δ . Hence, we obtain (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N X i =1 dg mi ( ¯ ξ + Q i ˜ ξ ) dλ − dg mi ( ¯ ξ ) dλ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ N · δ N = δ . (19)Define Λ ∗ δ := { w ∈ R | | dg m ( w ) /dλ | ≤ δ / } anddenote it as Λ ∗ δ = [Λ , ¯Λ]. In fact, since ( dg m /dλ )( λ ) ismonotonically decreasing, it holds that ( dg m /dλ )( λ ) ≥ δ / λ ≤ Λ and ( dg m /dλ )( λ ) ≤ − δ / λ ≥ ¯Λ.From here on, suppose that t ≥ T ( λ (0) , k ). Then, thetrajectory of ¯ ξ ( t ) approaches to Λ ∗ δ with the speed of atleast δ / (3 N ). For instance, if ¯ ξ ( t ) ≥ ¯Λ, it follows from(18a) and (19) that˙¯ ξ ≤ − δ N + δ N = − δ N . Hence, ¯ ξ ( t ) converges to Λ ∗ δ in finite time. In order tocompute the convergence time, it follows from (4) that | ˙¯ ξ ( t ) | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N N X i =1 d i − x i ( λ i ( t )) + φ i (cid:0) x i ( λ i ( t )) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ N N X i =1 ¯ x i − φ i (¯ x i ) − x i + φ i ( x i ) =: ∆ N . Therefore, it follows that | ¯ ξ (cid:0) T ( λ (0) , k ) (cid:1) | ≤ | ¯ ξ (0) | + ∆ N T ( λ (0) , k ) =: ζ ∗ . Define D ∗ as D ∗ := max (cid:0) | ζ ∗ | Λ ∗ δ , |− ζ ∗ | Λ ∗ δ (cid:1) . Then, itholds that ¯ ξ ( t ) ∈ Λ ∗ δ for all t ≥ T ( λ (0) , k ), where T ( λ (0) , k ) := 3 N D ∗ δ + T ( λ (0) , k ) . inally, for some fixed i and for all t ≥ T ( λ (0) , k ), (cid:12)(cid:12) p i (cid:0) ˆ x i ( λ i ( t )) (cid:1) − p i (cid:0) ˆ x i ( λ ∗ ) (cid:1)(cid:12)(cid:12) ≤ N X j =1 (cid:12)(cid:12) p j (cid:0) ˆ x j ( λ i ( t ) (cid:1) − p j (cid:0) ˆ x j ( λ ∗ ) (cid:1)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N X j =1 p j (cid:0) ˆ x j ( λ i ( t )) (cid:1) − p j (cid:0) ˆ x j ( λ ∗ ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) − dg m dλ ( λ i ) (cid:12)(cid:12)(cid:12)(cid:12) (20)where the first equality holds since p j (ˆ x j ( · )) is an in-creasing function, and the second equality holds since P Nj =1 p j (cid:0) ˆ x j ( λ ∗ ) (cid:1) = P Nj =1 d j . Note, we have (cid:12)(cid:12)(cid:12)(cid:12) dg m dλ ( λ i ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12) dg m dλ ( ¯ ξ ) (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) dg m dλ ( ¯ ξ + Q i ˜ ξ ) − dg m dλ ( ¯ ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ δ δ δ . (21)Therefore, (20) and (21) implies (cid:12)(cid:12) p i (cid:0) ˆ x i ( λ i ( t )) (cid:1) − p i (ˆ x i ( λ ∗ )) (cid:12)(cid:12) ≤ δ . By the definition of δ , it follows that | ˆ x i ( λ i ( t )) − ˆ x i ( λ ∗ ) | ≤ (cid:15) for all t ≥ T ( λ (0) , k ).An important feature of the proposed algorithm (15) isthat it is an initialization-free algorithm and hence allowsplug-and-play operation. In particular, it can be seenfrom Theorem 2 that the proposed algorithm convergesto a feasible solution regardless of the initial condition.Therefore, even if some parameters of the optimizationproblem (1) changes, the solution of (15) converges toa feasible solution of the new problem. Additionally,converged solution is close to an optimal if the couplinggain k is chosen as stated in Theorem 3. To boundthe performance uniformly across changes, the couplinggain must be chosen sufficiently large to incorporate allpossible cases. For instance, such gain can be found fromthe worst case scenario by assuming that the networkhas a known maximum capacity and parameters suchas f i and φ i are from a finite collection. More detaileddiscussions can be found in [10, Sec. 6.1] or [20]. Remark 3: Using similar arguments as in (21), itholds that the constraint violation and the objectiveerror f i ( x i ( t )) − f i (ˆ x i ( λ ∗ )) become small in finite time.In particular, | P Ni =1 d i − x i ( t ) + φ i ( x i ( t )) | ≤ δ and | f i (cid:0) x i ( λ i ( t )) (cid:1) − f i (ˆ x i ( λ ∗ )) | ≤ c i holds for all t ≥ T ( λ (0) , k ) where c i > | a − b | ≤ (cid:15) implies | f i ( a ) − f i ( b ) | ≤ c i for any a, b ∈ X i . Thus, constraintviolation and objective error can be made arbitrarilysmall by reducing (cid:15) which leads to higher coupling gain k . ♦ Results of Theorems 2 and Theorem 3 assume thatthe problem (1) is feasible. The behavior of the proposedalgorithm is also analyzed when the problem is infeasible. Theorem 4: Suppose that Assumptions 1, 2 and 3hold. Assume that the optimization problem (1) is infea-sible. Specifically, suppose that P Ni =1 d i − ¯ x i + φ i (¯ x i ) > i ∈ N , λ i ( t ) diverges to + ∞ andlim t →∞ ˙ λ i ( t ) = D , where D := ( P Ni =1 d i − ¯ x i + φ i (¯ x i )) /N > 0. Similar result also holds in the case of P Ni =1 d i − x i + φ i ( x i ) < . ♦ Proof: From the proof of Theorem 2, it holds that˜ ξ ( t ) is bounded. Moreover, we have˙¯ ξ = 1 N N X i =1 d i − ˆ x i ( ¯ ξ + Q i ˜ ξ ) + φ i (cid:0) ˆ x i ( ¯ ξ + Q i ˜ ξ ) (cid:1) ≥ N N X i =1 d i − ¯ x i + φ i (¯ x i ) = D > . Therefore, it follows that lim t →∞ ¯ ξ ( t ) = + ∞ . Thus,lim t →∞ λ i ( t ) = lim t →∞ ¯ ξ ( t ) + Q i ˜ ξ ( t ) = + ∞ and hence λ i ( t ) diverges to + ∞ .From (10), it holds that dg mi ( ¯ ξ ) /dλ is a constant forsufficiently large ¯ ξ . Therefore, it can be verified from(18a) and the boundedness of ˜ ξ that ˙¯ ξ satisfieslim t →∞ ˙¯ ξ ( t ) = D . In addition, boundedness of ˜ ξ and divergence of ¯ ξ implieslim t →∞ R T G (1 N ¯ ξ + Q ˜ ξ ) = R T ( d − ¯ x + φ (¯ x )) . In particular, ˜ ξ ( t ) converges to a constant value. There-fore, lim t →∞ ˙ λ i ( t ) = ˙¯ ξ ( t ) + Q ˙˜ ξ = D which completes the proof. V. Simulation A. Robustness to Changes For the simulation, the continuous-time algorithm (15)is discretized using the forward difference method. Inparticular, if we denote the sampling period with T , thenthe distributed algorithm (15) can be discretized into λ d i (( q + 1) T ) = λ d i ( qT ) + T dg mi dλ ( λ d i ( qT ))+ T k X j ∈N i ( λ d j ( qT ) − λ d i ( qT )) , (22)for i ∈ N , where q ≥ λ d i is the discretizedstate. In particular, λ d i ( t ) is a piecewise constant, rightcontinuous signal which is updated every T seconds.Specifically, the value of λ d i ( t ) is held constant until thenext update.Numerical simulation is done with IEEE 30 bus system[21] to verify the proposed algorithm. The local costfunction is given by f i ( x i ) = a i + b i x i + c i x i where b i and c i are strictly positive for i ∈ { , , , , , } whichare the buses with a generator. The power demand ofeach bus satisfies d i ∈ [0 , . 2] and P Ni =1 d i = 283 . φ i ( x i ) = α i x i where α i ∈ [0 . , . k = 40 is used for thesimulation. For the implementation, (22) is used with thesampling time of 0 . 005 seconds. We consider the followingscenario:S1) Normal operation condition for 0 ≤ t ≤ s .S2) At t = 10 s , demand at bus 5 is decreased by 20%.S3) At t = 20 s , generator at bus 1 stops generation andleaves the network.S4) At t = 30 s , bus 1 joins the network again and themaximum generation at bus 8 increases by 20%.Note in particular that the above scenarios are simulatedin one continuous session. Moreover, network topologychanges as a node leaves and joins the network duringthe operation.Simulation results are shown in Fig. 1. It can be seenthat an optimal solution is obtained in a distributedmanner despite the changes in operation conditions. At t = 20 s , the problem (1) becomes infeasible due to thelack of generation at bus 1. Hence, the trajectory of λ d i ( t ) diverges which verifies the result of Theorem 4.As feasibility is recovered at t = 30 s , λ d i ( t ) convergesagain. Value of the cost function is also shown in Fig.1(c) and it is seen that the optimal cost is approximatelyrecovered. The trajectory of power mismatch P Ni =1 d i − x i ( t ) + φ i ( x i ( t )) is shown in Fig. 1(d). It is observedthat the mismatch converges to zero (except for the casewhen the problem is infeasible) implying supply-demandbalance is satisfied.From the repeated simulations, we also observed thatthe coupling gain k and the sampling time T have a closerelationship for the stability. Specifically, if the gain k islarge, the sampling time must be reduced to yield a stablealgorithm. For example, it is observed that (22) divergeswith T = 0 . 01 but converges if T = 0 . B. Asynchronous Update The discretized algorithm (22) supposes that all agentsuse the same sampling time T and that the update ismade synchronously across all agents. However, it maybe hard to implement the synchronous algorithms inpractice due to the distributed nature of the system.Instead, in this section, we consider the case whereeach agent uses different sampling time T i and that theupdate is done asynchronously. In particular, let q ≥ λ d i (( q + 1) T i ) = λ d i ( qT i ) + T i dg mi dλ ( λ d i ( qT i ))+ T i k X j ∈N i (cid:0) λ d j ( qT i ) − λ d i ( qT i ) (cid:1) (23)where λ d i ( t ) is a piecewise constant, right continuoussignal which is updated every T i seconds. Note that thealgorithm (23) is equivalent to (22) if T i = T for all i ∈ N . However, if sampling time is different between Time[s] λ d i ( t ) (a) Trajectories of λ d i ( t ) Time[s] C o s t Optimal CostProposed (b) Trajectory of P Ni =1 f i (ˆ x i ( t )) and optimal cost Time[s] P o w e r [ M W ] Bus 1 Bus 5 Bus 8 (c) Trajectories of ˆ x i ( t ) for selected buses Time[s] M i s m a t c h [ M W ] (d) Trajectory of power mismatch P Ni =1 d i − x i ( t ) + φ i ( x i ( t )) Fig. 1: Simulation results Bus Sampling Time (s) Bus Sampling Time (s)Bus 5 0 . 05 Bus 11 0 . . 07 Bus 17 0 . . TABLE I: Sampling time used for each node.agents, then some agents update more frequently thanothers.The same IEEE 30 bus system is used to simulatethe algorithm (23). For the simulation, gain of k = 20is used while sampling time of 0 . s is used for allagents except for the ones denoted in Table I. Specif-ically, we consider the case when some nodes updateless frequently. Simulation results are shown in Fig. 2.Trajectories of λ d i ( t ) for selected bus is shown in Fig.2(a) and Fig. 2(b). It is clearly seen that each variableis updated asynchronously, and that some agents areupdated more frequently. Nonetheless, the solution λ d i ( t )still converges. Additionally, Fig. 2(c) depicts that theconverged solution satisfies supply and demand balance. Time[s] − − λ d i ( t ) Bus 1Bus 5 Bus 11Bus 16 Bus 17Bus 21 (a) Trajectories of λ d i ( t ) Time[s] − − λ d i ( t ) Bus 1Bus 5 Bus 11Bus 16 Bus 17Bus 21 (b) Trajectories of λ d i ( t ) for 0 ≤ t ≤ . Time[s] M i s m a t c h [ M W ] (c) Trajectories of power mismatch P Ni =1 d i − x i ( t )+ φ i ( x i ( t )) Fig. 2: Simulation results for asynchronous update. VI. Conclusion and Future Works The economic dispatch problem with nonlinear, sepa-rable power losses has been studied in this paper. Dueto the addition of nonlinear loss, the EDP becomesa non-convex optimization problem. However, it hasbeen shown that convex relaxation with dual decompo-sition can be used to obtain an optimal solution. Thedistributed algorithm is proposed and it is shown toconverge to a feasible solution while an optimal solu-tion is recovered with sufficiently high coupling gain.Specifically, the proposed algorithm does not require anyinitialization process and converges from any initial con-dition. Moreover, the behavior of the proposed algorithmis analyzed when the problem is infeasible. Future worksinclude the theoretical analysis of the discretized version of the proposed algorithm. References [1] A. J. Wood and B. F. Wollenberg, Power generation, opera-tion, and control . John Wiley & Sons, 2012.[2] Z. L. Gaing, “Particle swarm optimization to solving the eco-nomic dispatch considering the generator constraints,” IEEETrans. Power Syst. , vol. 18, no. 3, pp. 1187–1195, 2003.[3] T. Guo, M. I. Henwood, and M. van Ooijen, “An algorithm forcombined heat and power economic dispatch,” IEEE Trans.Power Syst. , vol. 11, no. 4, pp. 1778–1784, 1996.[4] W. T. Elsayed and E. F. El-Saadany, “A fully decentralizedapproach for solving the economic dispatch problem,” IEEETrans. Power Syst. , vol. 30, no. 4, pp. 2179–2189, 2015.[5] S. Yang, S. Tan, and J. X. Xu, “Consensus based approachfor economic dispatch problem in a smart grid,” IEEE Trans.Power Syst. , vol. 28, no. 4, pp. 4416–4426, 2013.[6] S. Kar, G. Hug, J. Mohammadi, and J. M. Moura, “Dis-tributed state estimation and energy management in smartgrids: A consensus+innovations approach,” IEEE J. Sel. Top-ics Signal Process. , vol. 8, no. 6, pp. 1022–1038, 2014.[7] H. S. Ahn, B. Y. Kim, Y. H. Lim, B. H. Lee, and K. K. Oh,“Distributed coordination for optimal energy generation anddistribution in cyber-physical energy networks,” IEEE Trans.Cybern. , vol. 48, no. 3, pp. 941–954, 2018.[8] A. Cherukuri and J. Cort´es, “Initialization-free distributedcoordination for economic dispatch under varying loads andgenerator commitment,” Automatica , vol. 74, pp. 183–193,2016.[9] P. Yi, Y. Hong, and F. Liu, “Initialization-free distributedalgorithms for optimal resource allocation with feasibilityconstraints and application to economic dispatch of powersystems,” Automatica , vol. 74, pp. 259–269, 2016.[10] H. Yun, H. Shim, and H. S. Ahn, “Initialization-free privacy-guaranteed distributed algorithm for economic dispatch prob-lem,” Automatica , vol. 102, pp. 86–93, 2019.[11] R. Dutta, L. Chong, and F. Rahman, “Analysis and experi-mental verification of losses in a concentrated wound interiorpermanent magnet machine,” Progress In ElectromagneticsResearch , vol. 48, pp. 221–248, 2013.[12] T. Yalcinoz and M. Short, “Neural networks approach forsolving economic dispatch problem with transmission capacityconstraints,” IEEE Trans. Power Syst. , vol. 13, no. 2, pp. 307–313, 1998.[13] G. Binetti, A. Davoudi, F. L. Lewis, D. Naso, and B. Turchi-ano, “Distributed consensus-based economic dispatch withtransmission losses,” IEEE Trans. Power Syst. , vol. 29, no. 4,pp. 1711–1720, 2014.[14] C. Zhao, J. He, P. Cheng, and J. Chen, “Consensus-basedenergy management in smart grid with transmission losses anddirected communication,” IEEE Trans. Smart Grid , vol. 8,no. 5, pp. 2049–2061, 2017.[15] A. Grauers, Design of direct-driven permanent-magnet gener-ators for wind turbines . PhD thesis, Chalmers University ofTechnology, 1996.[16] A. Cherukuri, E. Mallada, and J. Cort´es, “Asymptotic con-vergence of constrained primal–dual dynamics,” Systems &Control Letters , vol. 87, pp. 10–15, 2016.[17] D. Bertsekas, Nonlinear Programming . Athena Scientific,2016.[18] J. Kim, J. Yang, H. Shim, J. S. Kim, and J. H. Seo, “Ro-bustness of synchronization of heterogeneous agents by strongcoupling and a large number of agents,” IEEE Trans. Autom.Control , vol. 61, no. 10, pp. 3096–3102, 2016.[19] A. Cherukuri, B. Gharesifard, and J. Cort´es, “Saddle-point dy-namics: Conditions for asymptotic stability of saddle points,” SIAM Journal on Control and Optimization , vol. 55, no. 1,pp. 486–511, 2017.[20] J. Kim, H. Shim, and J. Wu, “On distributed optimal kalman-bucy filtering by averaging dynamics of heterogeneous agents,”in Proc. 55th IEEE Conf. Decision and Control , pp. 6309–6314, Dec. 2016.[21] O. Alsac and B. Stott, “Optimal load flow with steady-statesecurity,”