Distributed Optimal Economic Dispatch Based on Multi-agent System Framework in Combined Heat and Power Systems

In this paper, a novel distributed method is presented to solve combined heat and power economic dispatch problem, which is formulated as a distributed coupled optimization problem. The optimization goal is achieved by establishing two modified consensus protocols with two corresponding feedback parts while satisfying the electrical and heat supply–demand balance. Moreover, an alternating iterative method is proposed to handle the heat-electrical coupling problem existed in the objective function and the feasible operating regions. In addition, the proposed distributed method is implemented by a multi-agent system framework, which only requires local information exchange among neighboring agents. Simulation results obtained on a 16-bus test system are provided to illustrate the effectiveness of the proposed distributed method.


Introduction
Economic dispatch, as one of fundamental and key problem in power system, is formulated to an optimization problem that aims at minimizing the total operating cost while satisfying the supply-demand balance constraint and some capacity limit. With pressures of energy crisis [1] and environment pollution [2,3], there is an increasing interest in utilizing combined heat and power (CHP) into the power system, to decrease carbon emission, increase energy efficiency, and save operational costs in recent years, which brings many challenges for the economic dispatch problem (EDP).
During these decades, there have been plentiful literature studies concerned with solving the combined heat and power economic dispatch problem (CHPEDP) in a centralized manner, such as the lambda iteration method [4], the gradient method [5], multi-team particle swarm optimization method [6], evolutionary method [7], etc. However, with the rapid development of the CHP system, centralized method has become the bottleneck of solving CHPEDP and cannot satisfy the requirement of practical applications. To be specific, the centralized methods have to rely on a central controller to communication with all the other system components to collect global information, resulting in great implementation cost and high sensitivity to single-point failures [8][9][10]. Moreover, both the physical and communication topology of the future CHP system are prone to have a huge and variable topology with partially-unknown characteristic, which may seriously degrade the effectiveness of centralized methods. In addition, each component in CHP system should possess the plug-and-play feature, which is one of the key characteristics of the future CHP system [11][12][13]. Thus, the centralized method is not very suitable to address the distributed features of the future CHP system.
The distributed methods featuring fast computation, high flexibility and reliability do not relay on the central controller and can provide cost-effective distributed solutions to meet the requirements of It is worth noting that the energy resources discussed above may generate different types of energy (electric and/or heat) power. To facilitate the system modeling, three types of agents are defined in this paper according to different energy power generation types. The three different agents include electricity-only agents (EOAs), co-generation agents (CGAs) and heat-only agents (HOAs). Therein, the EOAs, which can generate electrical power only, consist of various kinds of DGs and EESDs. The CGAs, which can generate not only electrical power but also heat power, consist of various kinds of CHPs. The HOAs, which can generate heat power only, consist of various kinds of DHDs and HESDs. Then, each type of agent will be formulated to build the system model in the next section.

Modeling
The cost function for CGAs is often modeled as the following form [6]: The feasible operating region for electrical and heat power output of a CGA is show in Figure 2. The region of ABCD denotes the bounds of the feasible operating region, which can be represented as a set of linear constraints, i.e., ). Furthermore, the set of linear constraints can be re-written as the following forms [4]:  It is worth noting that the energy resources discussed above may generate different types of energy (electric and/or heat) power. To facilitate the system modeling, three types of agents are defined in this paper according to different energy power generation types. The three different agents include electricity-only agents (EOAs), co-generation agents (CGAs) and heat-only agents (HOAs). Therein, the EOAs, which can generate electrical power only, consist of various kinds of DGs and EESDs. The CGAs, which can generate not only electrical power but also heat power, consist of various kinds of CHPs. The HOAs, which can generate heat power only, consist of various kinds of DHDs and HESDs. Then, each type of agent will be formulated to build the system model in the next section.

Modeling
The cost function for CGAs is often modeled as the following form [6]: The feasible operating region for electrical and heat power output of a CGA is show in Figure 2. The region of ABCD denotes the bounds of the feasible operating region, which can be represented as a set of linear constraints, i.e., d m x pq + e m x qp + f m ≥ 0 (m = 1, 2, 3, 4). Furthermore, the set of linear constraints can be re-written as the following forms [4]: The reasons for obtaining Equations (1b) and (1c) are as follows. For each feasible operation point (x pq , x qp ), we can get x pq ≥ − (e m x qp + f m ) /d m . Note that there always exist one (or several) e m /d m ≥ 0 and one (or several) e m /d m < 0 for different m. Thus, x pq ≥ − (e m x qp + f m ) /d m can be re-written as the form F max (x qp ) ≥ x pq ≥ F min (x qp ). Therein, F min (x qp ) and F max (x qp ), which are determined by x qp , are the upper and lower bound functions of x pq . In this paper, we define x pq min (x qp ) and x pq max (x qp ) as F min (x qp ) and F max (x qp ), respectively. Similarly, we can define x qp min (x pq ) and x qp max (x pq ) as the lower and upper bound functions of x pq , respectively. Therefore, the expression of d m x pq + e m x qp + f m ≥ 0 is equivalent to the expressions of Equations (1b) and (1c).  The cost function of EOAs are modeled as the following form [25]: where i b and i c are both equal to zero for each EESD. Emphatically, of RG i are not constants, i.e., the maximum electrical power output of each wind generator will be changed following the change of wind speed. The cost function of HOAs can also be modeled as the following form: where β i and γ i are both equal to zero for every HESD. Similarly, for RHDs, The CHPEDP is to coordinate all the agents to minimize the system production cost ( ) C x while the electrical and heat power balance constraints, and a set of capacity limits constraints are satisfied. The objective can be mathematically represented as follows: subject to The cost function of EOAs are modeled as the following form [25]: where b i and c i are both equal to zero for each EESD. Emphatically, of RG i are not constants, i.e., the maximum electrical power output of each wind generator will be changed following the change of wind speed.
The cost function of HOAs can also be modeled as the following form: where β i and γ i are both equal to zero for every HESD. Similarly, for RHDs, The CHPEDP is to coordinate all the agents to minimize the system production cost C(x) while the electrical and heat power balance constraints, and a set of capacity limits constraints are satisfied. The objective can be mathematically represented as follows: Appl. Sci. 2016, 6, 308 5 of 20 subject to and Equations (1b), (1c), (2b) and (3b). In this paper, when i = 1, · · · , n 1 , it expresses the EOA. Similarly, when i = n 1 + 1, · · · , n 1 + n 2 and i = n 1 + n 2 + 1, · · · , n 1 + n 2 + n 3 , it expresses the CGA and HOA, respectively.

Centralized Method
The Lagrange multiplier theory can be used to analyze the CHPEDP. Inspired by the Lagrange multiplier theory, a fully distributed method will be introduced in Section 3. The Lagrange function is expressed as following equation: When the inequality constraints are not taken into account, the Lagrangian operator provides the following set of equations: i + β i be the heat incremental cost. According to Equation (7), the necessary conditions of the existing optimal operating point are that the entire electrical incremental cost and the heat incremental cost must be equal to λ p and λ q , respectively. If there exists the optimal operating point, the calculations of the electrical and heat power are given by: Furthermore, when the inequality constraints are taken into account, the necessary conditions to solve the problem can be expanded as: where if = p, λ * ν , x i_min and x i_max denote λ * p , x p i_min and x p i_max , respectively. If = pq, λ * ν , x i_min and i_min and x q i_max , respectively.

Remark 1:
Compared with the distributed methods [8,[17][18][19][20][21][22][23][24], in CHP system, the relationship between electrical and heat power are coupled, which existed in the objective function caused by the nonzero coefficient ξ i of Equation (1a), and in the feasible operating region of CGAs mathematically represented  (9). Thus, the optimal solutions cannot be separately calculated by the power system or heat system.

Basic Knowledge
Therein, the ith agent can receive information from jth agent if the corresponding arc e ij exists. If agent i can get (send) information from (to) agent j, then agent j is called in-neighbors (out-neighbors) of agent i. The set of in-neighbors . It is assumed that every agent has self-loops, which means every agent can send and receive information to itself. An undirected graph is connected if there exists an undirected path between any pairs of distinct nodes. A directed graph is strongly connected if there exists a path between any pair of two nodes with respect to the orientation of arcs. A nonnegative matrix is called row (or column) stochastic if all of its row-sums (or column-sums) are 1. In this paper, graph nodes represent the agents, and the graph arcs represent the communication links among agents.
(2) Basic Definitions: Consider a CHP system G with n 1 EOAs, n 2 CGAs and n 3 HOAs. Let the CHP system include two subgraphs G EC and G HC , where G EC composed of EOAs and CGAs, and G HC composed of HOAs and CGAs, meanwhile the two subgraphs are both strongly connected digraphs. Moreover, define two matrices R, S in subgraph G EC , where R is row stochastic and S is column Similarly, define two matrices R and S in subgraph G HC , where R is row stochastic and S is column (3) Consensus algorithm (or protocol): Based on the definitions discussed above, the first-order consensus protocol is given by: Since R is row stochastic, Equation (10) forms the first-order consensus protocol [26]. Each agent reaches the final common value χ i (∞)= ω T χ (0).

Distributed Method
Inspired by the consensus algorithm, this paper focuses on combing the consensus algorithm and the optimization theory to design a suitable, distributed method to solve the CHPEDP. Therein, there exist strong power-heat-coupling features both in the objective function and constraint limits. To address this issue, from the analysis in Section 2.3, there is not only an optimal electrical incremental cost λ * p but also an optimal heat incremental cost λ * q . More importantly, λ * p and λ * q are global variables. However, there is no centralized controller in this paper. Thus, the basic concept is to establish two different consensus protocols, in which the one is to make the electrical incremental cost run an common value (i.e., the final consensus value) and the other is to make heat incremental cost run an another common value, meanwhile a feedback mechanism is designed and added into the corresponding consensus protocol to obtain the optimal λ * p and λ * q . In addition, to solve the electrical-heating-coupling problem, which has been illustrated in Remark 1, the basic idea is to establish an effective alternating iterative method to solve this problem, where electricity and heat alternately converge to the optimal solution. Then, the updating rules of the proposed distributed method are designed as: (1) In the even-numbered iteration where θ denotes p or pq, ϑ denotes qp or q, and t = 2k or 2k + 1. In addition, represent all the variables relating to heat (electrical) power. Equation (11) shows the updating rules for Θ(2k). In this step, the purpose is to make Θ go to the optimal solution, meanwhile Φ(2k) = Φ(2k − 1) so that they can be seen as unchanged for further updating Θ. The detailed updating processes are as the following statements. Firstly, λ qp i (2k) and λ q i (2k) in Equation (11a) compose to the first consensus protocol, which includes consensus part and feedback part. The consensus part is associated to row stochastic R which drives λ qp i (2k) and λ q i (2k) to a common value. The feedback part includes ηy qp i (2k − 1) and ηy q i (2k − 1) which can let HOAs and CGAs change their incremental cost to satisfy the electrical power supply-demand balance to further obtain the optimal λ * q . That is, each HOA or CGA decrease (or increase) its corresponding heat incremental cost and heat power output when y q i (2k − 1) or y qp i (2k − 1) is negative (or positive) for the common value further adjustment to λ * q . Secondly, according to Equation (8), z qp i (2k) and z q i (2k) in Equation (11b) are used to update the estimated heat power outputs before considering the inequality constraints. Then, according to Equation (9), Equation (11c) defines the projection operators to take account of the inequality constraints, on the basis of the first assumption that x pq i_w_min ≤ x pq i (2k − 1) ≤ x pq i_w_max must be satisfied. Finally, inspired by literature [21], y qp i (2k) and y q i (2k) in Equation (11d) are used to update the estimated the heat power mismatch in a fully distributed fashion, which will be used to adjust the heat power incremental cost in the next cycle.
Similarly, Equation (12) shows the updating rules for Φ(2k + 1) meanwhile Θ(2k + 1) = Θ(2k). λ p i (2k + 1) and λ pq i (2k + 1) in Equation (12a) compose to the second modified consensus protocol which also includes its corresponding consensus part associated to row stochastic R and feedback part (i.e., ηy p i (2k) and ηy pq i (2k)) to obtain the optimal λ * p . z According to the aforementioned statement, we can give the major conclusion as follows. If the studied problem (Equation (4)) is feasible, then the method composed of Equations (11a)−(11c) and (12a)−(12c) is stable. Furthermore, each variable can converge to its optimal value, i.e., λ as t → ∞ (he proof can be seen in Appendix A).
According to Remark 2, for each CGA, the initial values of electrical and heat power outputs, i.e., , are set as the difference between the local electrical (or heat) load demands and the initial electrical (or heat) power outputs. To sum up, each variable is initialized as follows: Remark 2: There are two assumptions which must be satisfied during the updating process. In the even-numbered iteration, we have ) is within the area between L1 and L2 as shown in Figure 2. Due to the calculation of z qp i (2k) without considering the inequality constraints, (z pq i (2k), z qp i (2k)) may be out of the feasible operating region, e.g., the point E shown in Figure 2. However, the upper and lower bounds of z qp i (2k) can be easily obtained, because z pq i (2k) is known. Then, by using the corresponding projection operator (11c), the infeasible point can be mapped into the feasible operating region to further obtain (x pq i (2k), x qp i (2k)), e.g., the point G shown in Figure 2. It is not very difficult to verify that i_w_max . Thus, in the next iteration, the second assumption is satisfied. Then, (z pq i (2k + 1), z qp i (2k + 1)) is within the area between L3 and L4. By using the projection operator 12(c), (x pq i (2k + 1), x qp i (2k + 1)) can be further calculated, which is within the feasible operating region. By such analogy, the two assumptions will be always satisfied if only the initial value of (x  4). Let t = t + 1, turn to step 1.

Remark 3:
The distributed fashion of our method is reflected in that each agent only requires its own information and neighbors' information to calculate the optimal solution. Specifically, each agent only needs to exchange the local information of λ ) with its neighbors. Furthermore, the updating processes for all of variables are implemented through local calculation. On the contrary, centralized methods require a central controller and rely on a two-way communication network structure between the central controller and each agent. The central controller should collect all of the information of each agent to calculate the optimal solution, and then sent the optimal solution to each agent.

Remark 4:
Existing methods to solve the CHPEDP are mainly based on the centralized fashion, such as [4,7], etc. However, they have some obvious demerits as follows. (1) Since the centralized methods rely on a central controller, they are susceptible to single-point failures. To be specific, if the communication line between the central controller and any agent is failure, then the corresponding agent cannot keep working in the optimal operating point. Moreover, if the centralized controller is failure, then the whole system will lose control which may result in serious consequence, e.g., supply-demand unbalance. (2) When a new agent is connected to the system, the communication link between the new agent and the central controller must be established, which cannot be flexible for possible integration and expansions of devices. (3) Centralized methods cannot provide a satisfied plug-and-play function. On the contrary, the proposed distributed method is more reliable, robust, flexible and cost-effective and can address the plug-and-play feature in a better way. To be specific, since each agent can calculate its optimal operating point by using its own and neighbors information, a single (even several) link failures will not affect the optimal performance as long as the communication network maintains strong connection. Furthermore, when a new agent is required to be added into the system, the agent needs only establish the communication links with its neighbor agents, which simplifies the system maintenance and possible expansions. Moreover, the communication structure used in this paper with strong connection is a sparse communication structure, in which each agent only needs its local information and its neighbor's information to calculate the optimal solution. Thus, distributed method used in this paper is less expensive and more cost-effective. In addition, the distributed method provides desired plug-and-play feature for the future CHP system. It is very difficult for centralized method to meet this feature when the system becomes larger and larger.

Remark 5:
According to the definitions in graph theory, it can be seen that strongly connected structure belongs to directed graph while connected structure belongs to undirected graph. This also means, under strongly connected structure, each agent can send information to or receive information from its neighbors through single-way communication network. However, if it is under connected structure in undirected graph, each agent must send information to and receive information from its neighbors in the same time through two-way communication network. In other words, the strongly connected structure requires only a single-way communication network while the connected structure requires a two-way communication network. Thus, the strongly connected topology possesses less conservatism compared with connected graph. In addition, it is well-known that the convergence analysis, in the field of distributed control, is much more challenging under the directed graph setting.

Simulation Results
The first case is used to show the effectiveness of the proposed method via comparing with the well-known Lagrangian relaxation method. In the second case, one-day load demand of Northeastern University is used to test the performance of the proposed method under time-varying demand. The third case shows the requirement of plug and play. The last case considers the change of the maximum peak power tracking point for RGs and RHDs. The physical and communication structure of the test system is shown in Figure 3, in which EOA1, -3 and -5 represent three different FGs, EOA2 represents a EESD, EOA4 and -6 represent two different RGs, CGA1 and -2 represent two different CHPDs, HOA1 and -2 represent two different FHDs, HOA3 represents a RHDs and HOA4 represents a HESD. The parameters of each agent are listed in Table 1. The feasible operating region of CGA1 and CGA2 are shown in Figure 4. In the first three cases, the maximum peak power tracking point for EOA4, EOA6 and HOA3 are 90 (p.u.), 130 (p.u.) and 180 (p.u.), respectively.

Case Study 1: Comparison with Centralized Method
In this case study, the focus is on comparing the optimal values between the proposed method and the Lagrangian relaxation method [4]. The initial values can be obtained by Equation (13). The total electrical load demand is 750 (p.u.) with five electrical loads of 150 (p.u.) and the total heat load demand is 800 (p.u.) with five heat loads of 160 (p.u.). Then, the proposed method is used to calculate the optimal solutions. The estimated electrical and heat power incremental cost, mismatch and  In this case study, the focus is on comparing the optimal values between the proposed method and the Lagrangian relaxation method [4]. The initial values can be obtained by Equation (13). The total electrical load demand is 750 (p.u.) with five electrical loads of 150 (p.u.) and the total heat load demand is 800 (p.u.) with five heat loads of 160 (p.u.). Then, the proposed method is used to calculate the optimal solutions. The estimated electrical and heat power incremental cost, mismatch and outputs are shown in Figure 5. It can be seen that all of the electrical and heat power mismatch converge to zeros, which means the supply-demand balance are satisfied. After about 250 iterations, all of the electrical power incremental cost converge to a common value and all of the heat power incremental cost converge to another common value, meanwhile all agents are running in their feasible operating regions and the electrical and heat power outputs are listed in Table 2. In order to further verify the effectiveness of the proposed method, the Lagrangian relaxation method is used to solve this problem in a centralized fashion. The calculation results are also listed in Table 2. It can be seen that the solutions of the two methods are very close, which means the proposed method can also model the optimal values with distributed fashion. solve this problem in a centralized fashion. The calculation results are also listed in Table 2. It can be seen that the solutions of the two methods are very close, which means the proposed method can also model the optimal values with distributed fashion.

Case Study 2: Time-Varying Demand
In this case study, the focus is on the performance of the proposed method under time-varying loading conditions. The variable load demands of one day, shown in Figure 6a, are taken into consideration. The scheduling cycle is set to one hour. The initial values can be obtained by Equation (13). Then, the proposed method is used to calculate the optimal solution. From the results shown in Figure 6b-d, it can be observed that the system can automatically converge to new solutions at each scheduling cycle. Both the estimated electrical and heat power incremental cost are modified following the new demand, the estimated electrical and heat power mismatches also converge to zeros meeting the supply-demand balance, and all the electrical and heat power outputs are within their corresponding feasible range. Therefore, this case study verifies that the proposed method can automatically adjust to meet the time-varying loading conditions.

Case Study 2: Time-Varying Demand
In this case study, the focus is on the performance of the proposed method under time-varying loading conditions. The variable load demands of one day, shown in Figure 6a, are taken into consideration. The scheduling cycle is set to one hour. The initial values can be obtained by Equation (13). Then, the proposed method is used to calculate the optimal solution. From the results shown in Figure 6b-d, it can be observed that the system can automatically converge to new solutions at each scheduling cycle. Both the estimated electrical and heat power incremental cost are modified following the new demand, the estimated electrical and heat power mismatches also converge to zeros meeting the supply-demand balance, and all the electrical and heat power outputs are within their corresponding feasible range. Therefore, this case study verifies that the proposed method can automatically adjust to meet the time-varying loading conditions.

Case Study 2: Time-Varying Demand
In this case study, the focus is on the performance of the proposed method under time-varying loading conditions. The variable load demands of one day, shown in Figure 6a, are taken into consideration. The scheduling cycle is set to one hour. The initial values can be obtained by Equation (13). Then, the proposed method is used to calculate the optimal solution. From the results shown in Figure 6b-d, it can be observed that the system can automatically converge to new solutions at each scheduling cycle. Both the estimated electrical and heat power incremental cost are modified following the new demand, the estimated electrical and heat power mismatches also converge to zeros meeting the supply-demand balance, and all the electrical and heat power outputs are within their corresponding feasible range. Therefore, this case study verifies that the proposed method can automatically adjust to meet the time-varying loading conditions.

Case Study 3: Plug and Play Test
In this study case, the focus is on the plug and play adaptability, which is one of the most important advantages of the future CHP system. The initial conditions of are the same as in case 1. The electrical and heat load demands maintain unchanged during the simulation. At time step 400 t = , CGA2 is removed from the test system. It is not very difficult to verify the remaining subgraphs EC G and HC G are also strongly connected. The electrical (heat) loads supplied by CG2 are sent fifty-fifty to its out-neighbors in EC G ( HC G ). From the results in Figure 7, both the electrical and heat incremental cost increase and converge to their new optimal values, meanwhile the electrical and heat power mismatch go to zeros. Of course, the remaining agents have to supply more electrical and heat power to compensate for the electrical and heat power previously supplied by the disconnected CGA2. Moreover, at time step 800 t = , the CGA2 is plugged again in the test system, and the output of CGA2 is set to 11 Table 2. It can be seen that they are very close to the solution prior to disconnection. Therefore, this case study can verify the proposed method performs good plug and play capability.  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

Case Study 3: Plug and Play Test
In this study case, the focus is on the plug and play adaptability, which is one of the most important advantages of the future CHP system. The initial conditions of are the same as in case 1. The electrical and heat load demands maintain unchanged during the simulation. At time step t = 400, CGA2 is removed from the test system. It is not very difficult to verify the remaining subgraphs G EC and G HC are also strongly connected. The electrical (heat) loads supplied by CG2 are sent fifty-fifty to its out-neighbors in G EC (G HC ). From the results in Figure 7, both the electrical and heat incremental cost increase and converge to their new optimal values, meanwhile the electrical and heat power mismatch go to zeros. Of course, the remaining agents have to supply more electrical and heat power to compensate for the electrical and heat power previously supplied by the disconnected CGA2. Moreover, at time step t = 800, the CGA2 is plugged again in the test system, and the output of CGA2 is set to x In order to accommodate CGA2, the other agents have to reduce or not change (if the ones have been running on the upper bounders) their electrical and heat power outputs, while CGA2 increases its corresponding outputs. Thus, the electrical and heat incremental cost drop due to the lower average electrical and heat power outputs. In addition, the finally stable values are listed in Table 2. It can be seen that they are very close to the solution prior to disconnection. Therefore, this case study can verify the proposed method performs good plug and play capability.

Case Study 4: Maximum-Power-Output-Varying Condition
In this study case, the focus is on testing the performance of the proposed method when the maximum power output points of RGs and RHDs are changed. In this paper, we let EOA4, EOA6 and HOA3 represent a photovoltaic, a wind generator and a solar water heater, respectively. Their maximum power output profiles are shown in Figure 8a. In this case, the scheduling cycle is set to 20 s, i.e., EOA4, EOA6 and HOA3 measure and reset their corresponding maximum power output point every 20 s. The initial values can be obtained by Equation (13). Then, the proposed method is used to find the optimal operating point during each scheduling cycle. In Figure 8b, it can be observed that, within each scheduling cycle, all of the electrical power incremental costs converge to a common value and all of the heat power incremental costs converge to another common value. That implies the optimization objective is achieved. It is worth noting that the electrical (or heat) incremental cost decreases as the total electrical (or heat) power outputs of EOA4 and EOA6 (or HOA3) increasing. That is because EOA4, EOA6 and HOA3 are renewable energy resources, which have lower cost compared with fuel energy resources. Moreover, Figure 8c shows that both the estimated electrical and heat power mismatches can also converge to zeros within each scheduling cycle, which indicates

Case Study 4: Maximum-Power-Output-Varying Condition
In this study case, the focus is on testing the performance of the proposed method when the maximum power output points of RGs and RHDs are changed. In this paper, we let EOA4, EOA6 and HOA3 represent a photovoltaic, a wind generator and a solar water heater, respectively. Their maximum power output profiles are shown in Figure 8a. In this case, the scheduling cycle is set to 20 s, i.e., EOA4, EOA6 and HOA3 measure and reset their corresponding maximum power output point every 20 s. The initial values can be obtained by Equation (13). Then, the proposed method is used to find the optimal operating point during each scheduling cycle. In Figure 8b, it can be observed that, within each scheduling cycle, all of the electrical power incremental costs converge to a common value and all of the heat power incremental costs converge to another common value. That implies the optimization objective is achieved. It is worth noting that the electrical (or heat) incremental cost decreases as the total electrical (or heat) power outputs of EOA4 and EOA6 (or HOA3) increasing. That is because EOA4, EOA6 and HOA3 are renewable energy resources, which have lower cost compared with fuel energy resources. Moreover, Figure 8c shows that both the estimated electrical and heat power mismatches can also converge to zeros within each scheduling cycle, which indicates that the total electrical and heat supply-demand balance constraints are satisfied. In addition, the electrical and heat power outputs are shown in Figure 8d. It can be seen that, within each scheduling cycle, each agent can find its corresponding operating point, which is limited in its feasible operating region. Based on the discussion above, it can be concluded that, by using the proposed method, the system can respond automatically to converge to new solution as the changes of the maximum power outputs of EOA4, EOA6 and HOA3. Therefore, this case study clearly shows that the proposed method can automatically adjust to meet the maximum-power-output-varying condition.
Appl. Sci. 2016, 6, 308 16 of 23 that the total electrical and heat supply-demand balance constraints are satisfied. In addition, the electrical and heat power outputs are shown in Figure 8d. It can be seen that, within each scheduling cycle, each agent can find its corresponding operating point, which is limited in its feasible operating region. Based on the discussion above, it can be concluded that, by using the proposed method, the system can respond automatically to converge to new solution as the changes of the maximum power outputs of EOA4, EOA6 and HOA3. Therefore, this case study clearly shows that the proposed method can automatically adjust to meet the maximum-power-output-varying condition.

Conclusions
A distributed solution for the CHPEDP, taking into consideration multiple-energy-carriers, has been presented. Various kinds of energy conversion devices have been modeled and studied in the problem formulation. The proposed method reaches consensus on the electrical and heat incremental cost, respectively, using corresponding feedback terms that ensure the two supply-demand balance constraints. Because the relationship between the electrical and heat power are coupled, the concept of alternating iterative method has been presented and used to solve the electrical-heating problem, resulting in better solutions. In addition, several simulations have been presented to demonstrate the effect of the proposed method.