Distributed Absorption and Half-Search Approach for Economic Dispatch Problem in Smart Grids

The economic dispatch problem (EDP) is a significant class of optimization issues in the power system, which works on minimizing the total cost when generating a certain amount of power. A novel distributed approach for EDP is proposed in this paper. The presented approach consists of two steps. The first step, named absorption search, is to simplify the network structure through absorption searching. A flooding-based consensus approach is applied in the first step, which can be used to achieve consensus information among nodes. After the first step, only the generation nodes are kept in the network. The data collection can be completed by local computation and communication between neighbors. The first step can be considered as the stage of gathering information. In the second step, a distributed half-search algorithm makes the nodes obtain the final optimal solution in a distributed way. The results on three case studies demonstrate that the proposed approach is highly effective for solving the EDP.


Introduction
For power generation, smart grids exploit intelligent controls and developed technology to control the power generation combination composed of renewable resources [1].Because of advantages such as environmental friendliness and sustainable development, cost savings, electricity market and shared economy, smart grid and its related research topics have attracted much attention of scientific researchers.Ref. [2] surveyed the enabling technologies for Smart Grid.The authors of [3][4][5] analyzed the stability of nonlinear power systems based on backstepping control approaches.The optimal control and management for a large-scale battery energy storage system with wind and photovoltaic power station is introduced in [6][7][8].Ref. [9] discussed optimization of sustainable microgrid considering cost analysis, carbon emission and availability of energy resources.In the wake of development and expansion of the scale of smart grid, energy management systems (EMS) [10,11], especially online EMS [12,13], are becoming important research subjects.As an essential research direction of EMS, the economic dispatching problem (EDP) has been deeply studied due to its significant benefit of economical efficiency for smart grid.The EDP is a resource allocation problem, which minimizes the total generation cost while meeting the load demands.Several classical techniques such as Lambda-iteration method [14], Gradient methods [15], and Newton's method [16] have been developed to solve EDP whose cost function is convex.From the converge procedure point of view, 1. Simplified network structure by absorption search: In the network simplification stage, nodes in the system can share information and communicate with neighbors.Each time information is shared, the generation nodes (nodes in which the generator unit is connected) will update their status information by recording the current information of their neighbors, and the load nodes (nodes in which no generator unit is connected) whose information is recorded by their generation neighbor will disappear from the network at the current stage.2. Flooding-based consensus (FBC) approach for collection information: A flooding-based algorithm is used to gather information, and each node (agent) can obtain consistent information about its neighbors in the power network; this process of collecting data is carried out in a distributed way.3. Distributed half-search algorithm to obtain optimal solution in distributed manner: When only generation nodes are included in the network, each agent in the network can be competent for the result of EDP.A fully decentralized system based on half-search is used to solve EDP.
The remainder of this paper is organized as follows.Section 2 introduces the basic theory used in this article, including graph theory and consensus protocols.In Section 3, the EDP and the optimal solution are discussed.Section 4 presents the key network simplification and distributed half-search algorithm.Section 5 presents the analysis of the numerical simulation of a series of case studies.In Section 6, we conclude the paper.

Preliminary
In this section, basic preliminaries including some concepts in graph theory and consensus algorithm are introduced, which is needed in the following discussion.

Graph Theory
A (directed) graph [40,41] G is composed of vertices and edges, and is denoted by ordered sets G = {V(G), E(G)}, where V(G) = {1, 2, ..., n} and E(G) ⊆ V(G) × V(G) are the finite nonempty sets of the vertices and the edges, respectively.Define a symbolic relationship between ordered vertices in directed graph G, that is, (e) = (µ, ω), which means that (µ, ω) denotes a path from node µ to node ω. µ and ω are the tail and head of e, respectively.An undirected graph G is a special form of the directed graph with a bidirectional path by exchanging the direction of each edge (µ, ω) in G.
Define an adjacency matrix A(G) = (a ij ) n×n that associates with G (directed or undirected), and the index set of generation nodes in the smart grid can be represented by V(G).The directed path (i, j) ∈ E(G) shows that node i can send information to node j, and node j can receive information from node i.In this paper, we assume each node i ∈ V(G) does not belong to the neighbors of itself, that is, there is no self-loop in a directed graph.However, node i is allowed to receive its own information.Let } denote the in-neighbor set and out-neighbor set of ith node, respectively.Physically, it implies that a node i ∈ V(G) can send information to any node j, (j ∈ N − i ) and receive information from any node j, (j represent the in-degree and out-degree of node i, respectively, where | • | denotes the cardinality of a set; it is the number of elements in set "•".If there is a path, (i, j) between any two nodes (i and j) in G, the graph can be treated as strongly connected.The connection between the smart grid nodes constitutes the strong connected graph, distinctly, d ↑ i = 0 and d ↓ i = 0 for each node i ∈ V(G) in a strongly connected graph.Further, considering undirected graphs, it can be concluded that } and the number of neighbors of node i is represented by d i .
Based on the characteristics of power system, the concepts of directed and undirected graphs are used to study power flow and communication between nodes in power networks.Especially in solving consistency based on multi-agent system [42,43], network computing is indispensable.

Consensus Algorithm
For G = (V, E) that is strongly connected, define a non-negative adjacency matrix Q ⊆ R n×n as follows: It is not difficult to prove that [Q] ij > 0 for (j, i) ∈ E, and we can verify that the sum of column entries of Q is 1, i.e., ∑ N i=1 q ij = 1.This shows that the matrix Q is a column stochastic matrix associated with graph G. Similarly, Q T is a row stochastic matrix.According to the properties of the stochastic matrix, we have We consider all nodes in the smart grid work together as a team of N autonomous agents; it is assumed that each agent can define a subjective probability distribution for itself but does not know the subjective probability distribution of other agents.Label them from 1 to n and they can achieve a common goal value with certain constraint conditions.To explore the feasibility of this problem, we can learn from DeGroot's model of consensus [44].In this model, the opinions of each agent are gathered into an opinion pool and each agent adjusts the subjective probability distribution of itself by merging its own subjective probability distribution and learning the subjective probability distribution of other agents.Finally, the subjective probability distribution of all agents can reach a common value, and a certain subjective probability distribution parameter θ is formed on the basis, where θ can be considered as any discretionary variable whose value is a subset of the abstract parameter space Ψ.The subjective probability distribution assigned by each monomer i to parameter θ is represented by Φ i ∈ R, and the subjective probability distributions of all agents is denoted by It is assumed that the subjective probability distributions of all agents in different backgrounds on parameter space θ are denoted by Q i = [q i1 , q i2 , ..., q in ], and ∑ n i=1 q ij = 1, q ij ≥ 0. The stochastic Q associated with G denotes the n × n matrix consisting of the elements q ij , and we assume that it is feasible for each monomer i to update its subjective probability distribution from Φ(t) to Φ(t + 1).The linear iteration process for the update is as follows: The iteration index of the linear iteration is denoted by t = 0, 1, 2, ..., and the initial value is denoted by Φ(0).For the column stochastic matrix Q associated with a undirected graph G, the element Φ(t + 1) shown in Equation (2) can be expressed as: where q ii and q ij are the elements located on the diagonal and the ith row, jth column in matrix Q, respectively.Before studying the asymptotic consensus performance shown in Equation (2), we give two related theorems on the consensus problem.Theorem 1. [45] Given a non-negative matrix Q associated with a strongly connected graph G, the Q is a primitive matrix if G is aperiodic.Theorem 2. [46] If matrix Q ∈ R n×n is a primitive matrix and q ij ≥ 0, which is the element of Q, then, where ρ(Q) is the spectral radius of Q. U and V are the right and left Perron vectors of Q, respectively. where T is the right eigenvector with the character ζ i > 0 corresponding to the matrix Q at the eigenvalue 1.For all nodes i, there exists a properties equation, that is, Considering Theorem 2, let ρ(Q) > 0 denote a specific eigenvalue of Q (spectral radius).Consider the matrix Q, which is nonnegative and stochastic, such that ρ(Q) = 1.U > 0 and V > 0 are the right Perron vector and left Perron vector of Q with the following properties: From Equations ( 2) and ( 4), the consensus algorithm satisfies the following property: For any initial state Φ(0), if there exists Φ * ∈ R, such that lim t→∞ Φ i (t) = Φ * , ∀i = 1, 2, ..., n, the system state (shown in Equation ( 3)) has a stable value after t times iteration.
Based on the above mentioned theorems, we can draw a conclusion: the system (shown as in Equation ( 2)) state of each monomer i can converge to a stable value, and the convergence performance is determined by the initial value Φ(0) and the strongly connected topology Q.

Problem Formulation
In this paper, we only consider that EDP with a quadratic cost function [47].The basic problem of economic operation of power system with quadratic cost function can be considered as a convex optimization problem.It is widely used in the traditional EDP.It can be solved analytically without any approximations.Furthermore, the duality theory can be used to solve this problem.

Economic Dispatch Problem
We first establish a traditional power system model with M buses in which N generation units are included.Considering that the actual power system network may not contain a generation unit on every bus, we can prescribe a limit to M > N. The EDP can be expressed as the following form: where C g (P Gg ) is generation cost function.P Gg and P Dj denote the power generated and load demand by generator g and bus j, (g = 1, 2, ..., N; j = 1, 2, ..., M), respectively.From inequality in Equation ( 7), we can know that the power generated P Gg is restricted in the corresponding maximum P Gg and minimum P Gg bounds.The total load demand on all buses is denoted by P 0 and the generation cost function is denoted by C i (P Gg ), which is approximated as the follow quadratic function: Energies 2019, 12, 1527 6 of 21 where a i > 0, b i ≥ and c i ≥ 0 are the parameters of cost function.To simplify the expression, Equation ( 9) can be rewritten as: where β i > 0, α i ≤ 0. Obviously, the necessary and sufficient condition for the above problems to be feasible is The model of optimization problem shown as above is applied in this paper to solve EDP.

Solution
It can easily be proven that generation cost function C g (P Gg ) is strictly convex, and all the constraints of the EDP are linear, therefore the EDP can be treated as s convex-quadratic optimization problem.In other words, such a strict convex function can guarantee the establishment of the strong duality theorem, and then, we can set up the Lagrange dual problem of the original problem shown as Equations ( 6)-( 8), and the dual optimal solution obtained is also the solution of the original problem.According to the inverse form of the dual problem, the Lagrange dual problem is established as a non-differentiable concave function with the following structure: where λ denote the Lagrangian multiplier and C L g (λ) denotes the the Lagrange dual cost function.The incremental cost rates of the generator unit g is denoted by where λ ∈ R is the Lagrange multiplier and The above-mentioned Lagrangian dual problem can be established depending on whether the following equation or inequality condition can be satisfied, which includes power generated upper and lower bounds constraints and power generated balance constraints.The constraints mentioned above are written as, Energies 2019, 12, 1527 7 of 21 Based on the above analysis, we can define a mapping as follows Based on the above discussion, we can conclude that, under the restriction of power generated and load demand balance, if the primal problem shown as Equations ( 6)-( 8) is solvable, we can easily get λ * as the unique optimal solution, which satisfies According to the property of Lagrange's dual theorem and the strong duality theorem, we can verify that there is no duality gap between the solution of the primal problem and the dual problem, which means that the primal problem in Equations ( 6)-( 8) exists where the unique optimal primal P * Gg and optimal dual −τ g (λ * ) are equivalent in function value.That is, P * Gg = −τ g (λ * ), g = 1, 2, ..., n, i.e., We can discover from the optimal solution of the EDP that the constant coefficient γ i has no effect on the cost increment.From Equation ( 17), we can conclude that, once the value of parameter λ * is determined, the optimal solution of the EDP is determined.

Fully Distributed Solution for EDP
We divide the proposed approach into two stages.The first stage is network simplification, in which the load connected to each bus is incorporated into the nodes with generator units.A distributed approach based on FBC is used as the way of information dissemination.In the second stage, we propose a distributed solution for the optimal solution of EDP.The algorithm is implemented in the platform of Java agent development framework (JADE), which is an effective middleware for the development of multi-agent systems.

Absorption Search
To solve the EDP, we need to collect the load power demand in the first phase.However, it is not easy to calculate the load demand P 0 = ∑ M j=1 P j with a distributed method.We complete data collection with two steps: shrinking the size of network and exchanging information based on FBC.The concept of flooding algorithms [48,49] are applied by FBC, which is used in data networks for broadcasting.In this paper, this communication method is used between adjacent agents.
Each node in the power network may only contain pure generator unit, pure load unit, or both.Assume that an agent is embedded in each node to share information.It means that agents can send (receive) messages to (from) their neighbors in the power system.We create agents in JADE; each agent has a unique identification (ID) that consists of the order number and the node type.There are two kinds of node types in this paper, generating node (the node with generator unit) and load node (pure load connected to the node), and they are represented by (ID) N and (ID) M , respectively, where (ID) is the order number, and subscripts N and M denote generation node and load node, respectively.Assume that each agent knows its neighbor agent's ID and has the sufficient ability of computation.Each node i updates information by collecting and recording the information which includes the power data of the neighbor nodes and the ID of all neighbors of the neighbor.where msg _j is the set of information, ID M is the ID of the current load node, pm id and ID Mneibs are the load demand and the ID set of its neighbors, and [ID * Mneibs ] is alternative neighbor set, whose specific meaning is given below.Load agent will log off after successfully sending his information to the adjacent generation agent.Figure 1 gives the life cycle of load agents.Nneibs ]} where ID N represents the ID of the current generation node, and pm id is a parameter set that includes the load demand P Di , the upper bound (P Gg ) and lower bound (P Gg ) of the generator unit that is connected to the current node.ID M is the ID set of absorbed nodes at current iteration.ID Mneibs is the set of neighbors ID that connect to the absorbed nodes, [ID Nneibs ] is the set of neighbors ID of the generation node, and [ID * Nneibs ] is the set of alternative neighbors, which is designed for the generation agent that contains only one neighbor, and the neighbor corresponds to the type of load.Alternative neighbor sets prevent isolated nodes in the process of network simplification.When a generation agent receives a massage from its load neighbors, it updates its message.Furthermore, define a virtual agent (VA) in JADE, and each generation agent sends a message MSG VA to VA after each iteration step.MSG VA = 1 if the message sender receive new message from his neighbors, otherwise MSG VA = 0.All generation agents can receive feedback information from VA.If the received messages from generation agents are equal to zero in one iteration step, then VA will send message MSG VA = 0 as feedback information to all generation agent.Here, we obtain the final simplified network without load nodes.
The process of simplifying network is shown in Figure 2.
To illustrate the process of simplifying network, we take IEEE 14-bus system as an example, and the multi-agent system based on IEEE 14-bus system is shown in Figure 3.The blue dashed line connecting agents 8N and 4M represents that 4M is an element of 8 * Nneibs .

Distributed Half-Search Algorithm
After the first stage of information collection, the structure of the original communication network G M is reduced to the communication network G N = N , E N ).For every node i ∈ V n , we can get the summation of load demand P i collected by each generation node.Since all load nodes act as the dependency nodes of the bus where the generation node is determined, we now propose a fully distributed algorithm in the communication network G N to solve EDP.For the needs of the algorithm, we define an adjacency matrix A = {a} ij associated with G N .
Now, it is reasonable to assume that the non-negative matrix A can be used as the state matrix of each generation node in the topological G N .For i ∈ V N , a corresponding variable g i (t) is established for each generation node as the carrier variable of the power load carried by the generation node.And the total power collected by each generation node from the neighboring node is taken as the initial value of the variable, that is g i (0) = P i .Then, the following iterative algorithm is applied to reallocate the load power of each generation node.
where g i (t) is the carrier variable of the power load carried by the generation node, and j ∈ N n,i denotes all neighbors of node i in G N .For any node i ∈ V n , the computation in Equation ( 19) applied to each generation node will converge after t-time iterations.That is, g * = lim t→∞ g(t).The following equation is derived from Equation (5). where T is the right eigenvector with the character ζ i > 0 corresponding to the matrix A at the eigenvalue 1, and g * i is the final load information obtained by each generation node in G N .Define the basic variables λ(k), λ ↑ (k), λ ↓ (k) required by the algorithm, where λ(k) represents the Lagrange multiplier for each iteration, and the upper and lower bounds of λ(k) are denoted by λ ↑ (k) and λ ↓ (k), respectively.The step index of the half-search algorithm is denoted by k ≥ 0. The initializations of λ ↑ (k) and λ ↓ (k) may take any values as long as λ ↑ (k) is large enough and λ ↓ (k) is small enough.To make the interval as tight as possible, we initialize with the following equation: This means that the closer are λ ↓ (0) and λ ↑ (0) to the optimal Lagrangian λ * , fewer iterations are needed.There are multiple ways to define and solve for variable λ(k), and for simplification, λ(k) as an approximation of the optimal Lagrangian λ * , that is In the communication network graph G N = (V n , E n ), for all nodes i ∈ V n , the output of each generation node is assigned as follows: Then , for all nodes i ∈ V n , define an auxiliary variable X i (t) as the power generated carried by each generation node and initialized by X i (0) = P Gi (k).Then, the following iterative algorithm is applied to compute X i (t).
where X i (t) is the power generated carried by node.Let us define X * = lim t→∞ X i (t), and then, we get the convergent generators output variable X * corresponding to the current λ(k) according to Equation ( 5): Now, the half-search algorithm is proposed.In the communication topology graph G N , generation node i updates the values of the current upper bound λ ↑ (k + 1) and lower bound λ ↓ (k + 1) of the Lagrange multiplier λ(k) by comparing the magnitude of the local load information g * and the local output X * as follows: From Equations ( 21), ( 25) and ( 26), it is easy to get that λ * = lim k→∞ λ(k), and then, the output of each generation node can get a local optimal solution from Equation (22), that is The distributed half-search algorithm is summarized in Algorithm 1.
Algorithm 1 Distributed half-search algorithm.
Input: P i , i = 1, 2, ..., N:load demand for each generation agent; P min Gi : lower limit for generation unit i; P max Gi : upper limit for generation unit i; a i , b i α i , β i : cost coefficient; k=0.

Algorithm Analysis
We divide the proposed algorithm into two parts for analysis.First, we need to prove that the algorithm is convergent asymptotically; and second, we need to formulate a stop criterion for the algorithm.

Proof of Convergence
Before proving the convergence of the proposed algorithm, the following two remarks are given, which are used in the following proof.

Remark 1.
If there exists an optimal solution to the original EDP (shown as in Equations ( 6)- (8), then there exists a positive k that, when k → ∞, the proposed algorithm can converge to the globally unique optimal solution asymptotically.
Proof.From Equation (13), for any node i ∈ V n , we can get the incremental cost rates of the generation node i u i (P Gi ).λ + = max(u i (P Gi )) and λ − = min(u i (P Gi )) denote the upper and lower bounds of λ, respectively.Since the original EDP is feasible, we can get λ − ≤ λ ≤ λ + .
From Equation ( 16), for any node i ∈ V n , we can easily get that the function According to Remark 2, we can easily get that the distributed algorithm we proposed is convergent.In this paper, the EDP is constructed as astrictly convex optimization problem, therefore the distributed algorithm proposed can converge to a unique global optimal solution.

Stopping Criterion
For the half-search algorithm proposed in this paper, we give the algorithm stop criterion.We know from the discussion of the above algorithm that the lower and upper bounds of λ(k) as the search interval of the algorithm will be reduced by half after each search.Therefore, λ(k) in this algorithm can rapidly converge to the unique optimal solution λ * .At the same time, the optimal solution of the original power system EDP can be directly obtained as P * Gi = −τ i (λ * ) by solving the multiplier optimal solution λ * .Theoretically, as long as the preset step size k value can ensure the convergence of the algorithm, the stopping criterion can be set to stop when the algorithm runs to the preset k value.However, this stopping criterion does not satisfy the flexibility of the algorithm.Now, a stopping criterion with algorithm change is given as follow: where k is the iteration step size, and δ is a preset positive number small enough.By applying this stopping criterion, the algorithm can be stopped after convergence to unique optimal solution.

Simulation Results
Three simulation cases were studied to verify the effectiveness of the proposed algorithm.Firstly, we studied the simulation results of the algorithm based on IEEE14-bus data with generators constraints.Secondly, on the basis of simulation Case 1, we verified the plug-and-play performance of the proposed algorithm by removing load nodes of IEEE14-bus and adding a generation unit.Lastly, a case study based on IEEE57-bus was used to test the performance of the proposed algorithm in a wide range of power system.

Case 1: With Generator Constraints
We first studyied the case with generation constraints based on the IEEE 14-bus system.In the course of the study, all research data were derived from IEEE14-bus data.In our study, buses {1, 2, 3, 6, 8} were chosen as the generation buses and the load buses were {2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14}.Note that bus 7 was selected as neither a generation bus nor a load bus since the power generated and the load demand of bus 7 were equal to zero.The local power load were: P D2 = 21.7 MW, P D3 = 94.2MW, P D4 = 47.8MW, P D5 = 7.6 MW, P D6 = 11.2MW, P D9 = 29.5 MW, P D10 = 9 MW, P D11 = 3.5 MW, P D12 = 6.1 MW, P D13 = 13.5 MW, and P D14 = 14.9 MW.We calculated the summation of the load demand as P 0 = 259 MW, which was unknown to each individual bus node.Five generation nodes were located on buses 1, 2, 3, 6, 8, and the generator parameters of each generation unit are shown in Table 1.We took δ = 0.0001 as the stopping criterion.Each generator unit was initialized with the sum of the load power collected in the first stage: P G1 = 332.4MW, P G2 = 70 MW, P G3 = 100 MW, P G6 = 100 MW, and P G8 = 100 MW.We set the lower and upper bounds of Lagrange multiplier as: λ ↓ (0) =2 $/MW and λ ↑ (0) = 72 $/MW.The numerical simulation results with generators constraints are shown in Figure 5.To clearly analyze the simulation results of this case, we show the Lagrange multiplier, λ(k), in Figure 5 (top), the evolution of ∑ P j (j = 1, 2, 3, 6, 8) in Figure 5 (middle), and the generator output P Gg (g = 1, 2, 3, 6, 8) in Figure 5 (bottom).we subjectively set the iteration step to k = 20, while the stopping condition was satisfied at k = 12.
The Lagrangian multiplier λ(k) gradually approached an optimum stable value λ * = 5.4165 $/MW at iteration step k = 12, which affected the output of each generator unit.As λ(k) was asymptotically stable, the total power generated by all generators units was gradually ).The optimal generators output P * g and the optimal Lagrangian multiplier λ * satisfied the generator constraints.At the same time, the total power generated equaled the total load demand.Since the generator constraints and the power generated cost coefficient of generator units 3, 4, and 5 in the IEEE14-bus data were identical, the output of each generator unit was completely uniform.adapted to the total static load demand.When λ(k) converged asymptotically to the optimal λ * , the combination of generators output of each generator unit at the optimal λ * was optimized, and the output powers of the generator units were: P * G1 = 39.70 MW, P * G2 = 6.84 MW, P * G3 = 70.82MW, P * G4 = 70.82MW, and P * G5 = 70.82MW, and ∑ P j = 259 MW (j = 1, 2, 3, 6, 8).The optimal generators output P * g and the optimal Lagrangian multiplier λ * satisfied the generator constraints.At the same time, the total power generated equaled the total load demand.Since the generator constraints and the power generated cost coefficient of generator units 3, 4, and 5 in the IEEE14-bus data were identical, the output of each generator unit was completely uniform.

Case 2: Plug and Play Capability
The characteristics of plug and play reflect the flexibility of power supply and demand of smart grid.In this case, we slightly changed the IEEE14-bus power system structure in Case 1 to study the plug-and-play performance of the proposed half-search algorithm in power systems.In this case, all research data were the same as in Case 1 until the system was changed.To design a more realistic power supply and demand scenario experiment, we first removed a load power based on Case 1, which caused a sharp reduction in the total load power of the system.The plug-and-play performance of the proposed half-search algorithm was effectively tested when the total load power was sharply reduced.Besides, to verify the plug-and-play performance of the algorithm under different generator conditions, we added a generator unit to change the generator conditions of the system for a period of

Case 2: Plug and Play Capability
The characteristics of plug and play reflect the flexibility of power supply and demand of smart grid.In this case, we slightly changed the IEEE14-bus power system structure in Case 1 to study the plug-and-play performance of the proposed half-search algorithm in power systems.In this case, all research data were the same as in Case 1 until the system was changed.To design a more realistic power supply and demand scenario experiment, we first removed a load power based on Case 1, which caused a sharp reduction in the total load power of the system.The plug-and-play performance of the proposed half-search algorithm was effectively tested when the total load power was sharply reduced.Besides, to verify the plug-and-play performance of the algorithm under different generator conditions, we added a generator unit to change the generator conditions of the system for a period of time after the system converged again.Note that, whether removing a to reduce the total power of the system or adding a generator unit to change the system's generator conditions, the total power generated must equal the total power demand.
At iteration step k = 23, the load demand on bus 4 was removed, and the situation of sharply reduced power demand was simulated.At iteration step k = 43, a generator unit was added between the 13th and 14th buses.The constraints of the generator units were P G6 = 200 MW and P G6 = 50 MW, and the power generated cost coefficient was the same as the power generated coefficient of the second generator unit.That is, the power generated cost coefficient of the added generator unit was: a = 0.25, b = 2, c = 0.The results of the case study are shown in Figure 6.time after the system converged again.Note that, whether removing a load to reduce the total power of the system or adding a generator unit to change the system's generator conditions, the total power generated must equal the total power demand.
At iteration step k = 23, the load demand on bus 4 was removed, and the situation of sharply reduced power demand was simulated.At iteration step k = 43, a generator unit was added between the 13th and 14th buses.The constraints of the generator units were P G6 = 200 MW and P G6 = 50 MW, and the power generated cost coefficient was the same as the power generated coefficient of the second generator unit.That is, the power generated cost coefficient of the added generator unit was: a = 0.25, b = 2, c = 0.The results of the case study are shown in Figure 6.Since this case study was based on Case 1, the simulation data were the same as the simulation data of Case 1 at the beginning.Therefore, before the reduction of the total power demand, the power generated by each generator unit reached a steady state value.At iterative step k = 23, the total load demand decreased by 47.8 MW.After another 13 iterations of the proposed algorithm, the Lagrangian multiplier λ(k) could approach a steady state value λ * = 5.124 $/MW again.At the same time, the total power generated and the total load demand were balanced again, and the power generated by each generator unit was also stabilized at a steady state value  Since this case study was based on Case 1, the simulation data were the same as the simulation data of Case 1 at the beginning.Therefore, before the reduction of the total power demand, the power generated by each generator unit reached a steady state value.At iterative step k = 23, the total load demand decreased by 47.8 MW.After another 13 iterations of the proposed algorithm, the Lagrangian multiplier λ(k) could approach a steady state value λ * = 5.124 $/MW again.At the same time, the total power generated and the total load demand were balanced again, and the power generated by step of absorption search, load agents hold the message shown as follow, msg _j = {ID M , pm id , [ID Mneibs ], [ID * Mneibs ]}

Figure 1 .
Figure 1.Life cycle of load agent.We create generation agent with tuple, msg _i = {[ID N , (pm id ), MSG VA ], [< ID M > (ID Mneibs )], [ID Nneibs ], [ID *Nneibs ]} where ID N represents the ID of the current generation node, and pm id is a parameter set that includes the load demand P Di , the upper bound (P Gg ) and lower bound (P Gg ) of the generator unit that is connected to the current node.ID M is the ID set of absorbed nodes at current iteration.ID Mneibs is the set of neighbors ID that connect to the absorbed nodes, [ID Nneibs ] is the set of neighbors ID of the generation node, and [ID * Nneibs ] is the set of alternative neighbors, which is designed for the generation agent that contains only one neighbor, and the neighbor corresponds to the type of load.Alternative neighbor sets prevent isolated nodes in the process of network simplification.When a generation agent receives a massage from its load neighbors, it updates its message.Furthermore, define a virtual agent (VA) in JADE, and each generation agent sends a message MSG VA to VA after each iteration step.MSG VA = 1 if the message sender receive new message from his neighbors, otherwise MSG VA = 0.All generation agents can receive feedback information from VA.If the received messages from generation agents are equal to zero in one iteration step, then VA will send message MSG VA = 0 as feedback information to all generation agent.Here, we obtain the final simplified network without load nodes.The process of simplifying network is shown in Figure2.To illustrate the process of simplifying network, we take IEEE 14-bus system as an example, and the multi-agent system based on IEEE 14-bus system is shown in Figure3.The blue dashed line connecting agents 8N and 4M represents that 4M is an element of 8 * Nneibs .

Figure 2 .
Figure 2. Flow chart of simplifying network.

Figure 3 .
Figure 3. Multi-agent system for IEEE 14 bus system.
(a) First search step (b) Second search step

Table 1 .
Generator parameters in IEEE 14-bus system.