Double-consensus based distributed optimal energy management for multiple energy hubs

: This paper presents a novel distributed double-consensus algorithm to solve the optimal energy management problem for multiple energy hubs interconnected with each other. The objective is achieved by establishing two interactive and paralleled consensus procedures modiﬁed by their corresponding feedback terms. Meanwhile, a novel projection operation method is proposed to map the infeasible values into the feasible operating region. The proposed algorithm can effectively handle the coupled variables problem existing in the objective function and constraint limits. Moreover, the optimality and convergence analysis are performed strictly under strong connectivity conditions only. Simulations performed on standard test cases are provided to illustrate the effectiveness of the proposed distributed algorithm.


Introduction
The energy crisis and environmental pollution are regarded as a major challenge all over the world.To address these issues, multiple energy systems, aiming at integrating conventional independent energy service networks such as power and natural gas systems, could provide a promising solution [1,2].As one type of multiple energy system, the concept of the energy hub has gained significant attention recently [3,4].An energy hub focuses on feeding different energy loads through multiple energy inputs and outputs, acting as an interface between various energy producers and consumers.It could provide the desired effects of reconstructing traditional decoupled energy supplies, adding operational flexibility of energy services, improving overall efficiency and reliability, reducing pollutant emissions, etc.
The EMP is a fundamental and key topic in different energy and building systems.L. Tronchin et al. [5] made an outstanding contribution in reviewing and analyzing the works related to energy management, from models to technologies, such as the effects of different components [6], cost and performance analysis [7,8], technology evolution [9], etc.In this paper, we focus on investigating EMP within the context of energy hub-based multiple energy systems, which is typically formulated as an optimization problem.Several research works have been reported on this topic, which can be classified into two categories.The first mainly aims at achieving the internal resource management of individual energy hubs, with the objective of maximizing the energy allocation profits or minimizing the energy cost [10][11][12].To further enhance the system robustness and energy efficiency, the other tries to investigate the cooperative energy management problem of multiple energy hubs, where the energy hubs interact with each other to determine the global optimal operations.For instance, a multiagent genetic algorithm was proposed in [13] to solve the economic dispatch problem for multi-hubs with the consideration of uncertain renewable energy resources.In [14], a probabilistic optimization approach was introduced for renewable-based residential energy hubs, in which the interactions between the electricity and natural gas were taken into consideration.In [15], a self-adaptive learning with time varying acceleration coefficient-gravitational search algorithm was proposed to address multi-objective economic dispatch problem for energy hubs.Moreover, the underlying energy flow optimization [16,17], and the reliability of energy demands [18] and criteria [19] have been further taken into consideration and studied in the energy management problem.It is worth noting that the approaches mentioned above for solving the energy management of energy hubs are implemented centrally.However, as the scale of the system expands, centralized approaches may face several challenges, which can be summarized as follows: (1) The centralized approaches rely on a powerful central controller to process a huge mass of data and require two-way and high-bandwidth communication to work on the information gathered for the whole system [20,21].Thus, solutions based on centralized approaches are costly to implement and are prone to single-point failures and modeling errors.(2) Despite the plug-and-play and peer-to-peer nature of distributed energy resources, both the physical and associated communication topologies tend to be subject to topology variabilities, which may undermine the efficacy of the centralized approaches [22,23].
As an alternative, the distributed approaches, which have some advantages compared to centralized approaches, e.g., enhanced reliability and robustness [24], fast computation and reduction in communication [25], are more suitable for accommodating a mass of units.Not surprisingly, only a small number of results have been documented.In [26], the authors proposed the concept of smart energy hubs and formulated the interaction among smart energy hubs as a noncooperative game.Therein, a distributed algorithm was introduced to find the Nash equilibrium.Furthermore, the interaction among smart energy hubs was further modeled as an ordinal potential game in [27], where a distributed algorithm was proposed to solve this problem.Considering the autonomy and self-interest among energy hubs, their cooperative interaction was formulated as a bargaining cooperative game in [28], where an alternating direction method of multipliers was applied to find the optimal solutions.It is worth noting that the existing game-theoretic-based distributed approaches can achieve a partial shift from a centralized manner to a distributed manner.However, they also require a cloud computing center, a central price coordinator, or a virtual coordinator to obtain the optimal operations, which may not be regarded as fully distributed methods.Since there are some global supply-demand balance constraints and strong coupling relationship among variables, it is also a challenge to design an algorithm to solve the EMP of multiple energy hubs in a fully distributed manner.
To address these issues, a novel distributed algorithm, referred as the DDC algorithm, is presented.By introducing the consensus concept and designing a new projection operation method, our proposed approach can make the global computation process divided into each energy hub.As a result, each energy hub can locally calculate its optimal operation.Compared to the existing literature, the major contributions of the presented distributed algorithm are summarized as follows: (1) By making use of some change of variables, the strongly coupled form between variables in global equality constraints are transformed to the local function and inequality constraints.With this effort, the EMP is formulated as a kind of distributed coupled optimization problem.(2) The proposed DDC algorithm, which does not require a price coordinator or a leader to collect global parameters, can solve the EMP in a fully distributed fashion.In contrast to centralized algorithms, the proposed algorithm is more flexible, reliable and robust, etc.To the best knowledge of the authors, there still exist very few papers concerned about solving the EMP for multiple energy hubs in a fully distributed fashion.(3) The strong coupling problem, existing in the objective function and constraint limits, can be effectively solved by implementing the proposed DDC algorithm.Moreover, the rigorous optimality and convergence analysis are presented only under strong connectivity conditions, which is more reasonable and general, and possesses less conservatism.
The rest of this paper is organized as follows.In Section 2, the system model and optimal conditions are briefly introduced.In Section 3, some basic knowledge is firstly introduced.Then, the proposed DDC algorithm and the projection operation method are presented.In Section 4, several case studies are presented to show the effectiveness of the proposed DDC algorithm.Section 5 concludes the paper.

Energy Hub Model
The energy hub is a concept coupling various energy carriers, such as electricity, gas and heat, etc., which can be converted to various forms of energy to provide diversified loads.
Here, for each hub, we consider electricity and natural gas as energy inputs, while electricity and heat are considered as energy outputs.The converter devices consist of a transformer, boiler and CHP unit.The conversion relationship between inputs and outputs is defined as the coupling matrix form given in Appendix A-(a).

EMP of Multiple Energy Hubs
In this paper, we consider the EMP of interconnected energy hubs to achieve the global optimal energy configuration.Here, the anticipated structure of the multiple energy hubs system is shown in Figure 1.Compared to existing multiple energy hubs systems, it has one important difference, which is the paradigm shift from a traditional centralized manner to a distributed manner.Under this distributed structure, the work of each agent that does not rely on a centralized controller is based only on local communication and calculation to find the optimal operating point and dispatch factor.The considered EMP can be formulated as a distributed optimization problem with a global objective function and a set of constraints.
(1) Objective: The total objective function is to cooperatively minimize the aggregated energy costs (i.e., the electricity and natural gas consumption costs associated with the individual energy hub) and the emissions penalty from electricity and heat power generation.The mathematical expression of this objective function can be seen in Appendix A-(b).(2) Constraints: Energy hubs collectively make decisions on their own subject to global and local constraints.Accordingly, the set of local constraints for each energy hub is listed in Appendix A-(c).Moreover, the global constraints imply the supply-demand balance constraints.Specifically, the total electricity and heat power generation, i.e., ∑ n i = 1 E e i,out and ∑ n i = 1 E h i,out , should be equal to the total respective electricity and heat load demands, i.e., ∑ n i = 1 l e i and ∑ n i = 1 l h i , respectively.In terms of the defined coupling matrix, we can further establish the relationship between the input variables and load demands as follows: It is worth noting that Equations (1a) and (1b) are non-affine constraints caused by the coupled product of variables ρ i and E g i,in , which is not beneficial to the design of the distributed optimization algorithm.To address this issue, we employ the following change of variables.
where E g,ρ i and E g,1−ρ i are the natural gas consumed by CHP unit and boiler, respectively.In this novel coordinate system, the above optimization can be re-written in the following form by some operations.
subject to where Appl.Sci.2018, 8, x FOR PEER REVIEW 4 of 18 ( ) e  For convenience, Equations (5b) and (5c) can equally be expressed in the following form ( ) Next, the Lagrange multiplier theory can be used to analyze the problem in Equations ( 3)- (6).Inspired by the Lagrange multiplier theory, a fully distributed algorithm will be introduced in Section 3. The Lagrange function and the corresponding KKT conditions are shown in Appendix A-(d).For convenience, Equations (5b) and (5c) can equally be expressed in the following form for the local objective function.As a result, the non-affine constraints (1a and 1b) are converged to affine ones, i.e., Equations (4a) and (4b).
Next, the Lagrange multiplier theory can be used to analyze the problem in Equations ( 3)- (6).Inspired by the Lagrange multiplier theory, a fully distributed algorithm will be introduced in Section 3. The Lagrange function and the corresponding KKT conditions are shown in Appendix A-(d).

Basic Knowledge
(1) Graph theory: We consider a directed graph G = (V, E, A) with a set of vertices The in-neighbors and out-neighbors of ith node are denoted by . Each node can receive information from its in-neighbors and send information to its out-neighbors; meanwhile, it can be also assumed to communicate with itself.
(2) Basic definitions: With regard to problem Equations ( 3)-( 6), there are three variables that need to be calculated.To this aim, we define two strongly regular graphs G 1 and G 2 .The first consists of n nodes, where each node represents the statement related to E g,1−ρ i . The second is built by dividing each node of graph G 1 into two nodes which represent the statements related to E e i and E g,ρ i , respectively.Furthermore, we define two matrices R = r i,j and S = s i,j associated with G 1 ; meanwhile, let r i,j = 1/d It is not difficult to verify that R is row stochastic and S is column stochastic.Similarly, we define two matrices R = r i,j and S = s i,j for G 2 ; meanwhile, let (3) Consensus algorithm: Based on the definitions discussed above, we consider two different discrete-time systems-shown in Appendix B-(a)-which will be employed in the design of our proposed algorithm.

DDC Algorithm without Inequality Constraints
In this section, none of the inequality constraints are taken into consideration; the consideration of inequality constraints will be discussed in the next section.Based on the KKT conditions, we can obtain the relationship between the optimal solutions and the Lagrange multipliers (see Appendix B-(b) for details).This indicates that the optimal E e i , E g,ρ i and E g,1−ρ i can be obtained if the optimal λ e and λ g can be calculated.Please note that λ e and λ g are both global decision variables.In contrast to centralized approaches, there is no central controller to correct the global information to calculate the optimal λ e and λ g and send the solution to each component.On the contrary, the challenge of this paper is to obtain the optimal λ e and λ g by using only local information and messages exchanged with neighbors.Meanwhile, λ e and λ g are interactive with each other.Inspired by the feature of (P1), the basic concept is to establish two different consensus protocols, in which the one is to make all λ i,e,p (k) and λ i,e,c (k) run a common value and the other is to make all λ i,g,c (k) run another common value.Meanwhile a feedback term is added into the corresponding consensus protocol to obtain the optimal λ e and λ g .Therein, the feedback term uses local estimation of the electricity and heat power mismatches to ensure power and heat supply-demand balances.Then, the designed DDC algorithm, mainly consisting of three parts, is given by: (1) The united updating rules for energy hubs' coordination to estimate the electricity and heat multipliers are designed for where T and y = [y T e,p , y T e,c , y T g,c ] T .
Appl.Sci.2018, 8, 1412 6 of 17 (2) Based on the relationship between the optimal solutions and Lagrange multipliers, the updating rules for energy hubs to estimate the optimal E e i , E g,ρ i and E g,1−ρ i are designed for where (3) Please note that the electricity and heat power supply-demand balance constraints, i.e., Equations (4a) and (4b), are global constraints.Similar to the solving method for λ e and λ g , the concept of local estimation of the mismatch between demand and total generations, obtained by using only neighbors' information, is used to solve the two global constraints in a distributed manner.Inspired by the features of (P2), the updating rules for the energy hubs' coordination to estimate the electricity and heat power mismatches are designed as follows: where S = S 0 0 S , S = S pp S pc S cp S cc .
With regard to the determination of initializations, they are shown in detail in Appendix B-(c).
Remark 2. The coordinative updating rules of Equations ( 7) and (10) show that each component only requires the information received from its in-neighbors to update its corresponding λ and y.In addition, the updating of E is implemented by each energy hub's own information.Therefore, the DDC algorithm is fully distributed, only requiring local communication and calculation among neighbors, without a central controller.
Lemma 1. (Perron-Frobenius, [29]) Let P be a primitive nonnegative matrix with left and right eigenvectors w and v, respectively, satisfying Pv = v, w T P = w T , and v T w = 1.Then lim k→∞ P k = vw T .
Theorem 1.If the problem of Equations ( 3) and ( 4) is feasible, there exists a positive value ς such that for all values of 0 < η < ς, the fixed points of the DDC algorithm determined by Equations ( 7)-( 10) satisfy the KKT conditions of optimality.
The proof of Theorem 1 is provided in Appendix C.

DDC Algorithm with Inequality Constraints
For variable E e i , the corresponding inequality constraints in Equation (5a) can be taken into consideration by introducing the following well-known projection operations.
Different from E e i , E g,ρ i and E g,1−ρ i are coupled by a set of linear inequality constraints, which greatly increases the solving difficulty.To address this issue, the basic concept is to establish a projection operation method that can map the solution without considering inequality constraints onto one point of the feasible operating region determined by Equation ( 6); meanwhile, the values of this point must satisfy the corresponding KKT conditions.Please note that the new solution may be inside the feasible operating region corresponding to no active constraints, on boundaries corresponding to one active constraint, or at the corner point of boundaries corresponding to two active constraints.Then, according to optimality conditions, we establish the identification conditions for the above three cases, which can be seen in Appendix B-(d).
Furthermore, the projection operations for E g,ρ i and E g,1−ρ i , to take Equation ( 6) into account, are given by , if (26)satisfied (13) Theorem 1 can be expanded to solve the problem Equations ( 3)-( 6) by only replacing Equation ( 9), in Theorem 1, with Equations ( 11)- (13).For clarity, the flowchart of the proposed DDC algorithm is summarized in Figure 2. this point must satisfy the corresponding KKT conditions.Please note that the new solution may be inside the feasible operating region corresponding to no active constraints, on boundaries corresponding to one active constraint, or at the corner point of boundaries corresponding to two active constraints.Then, according to optimality conditions, we establish the identification conditions for the above three cases, which can be seen in Appendix B-(d).Furthermore, the projection operations for , to take Equation 6 into account, are given by ( ) Theorem 1 can be expanded to solve the problem Equations ( 3)-( 6) by only replacing Equation ( 9), in Theorem 1, with Equations ( 11)- (13).For clarity, the flowchart of the proposed DDC algorithm is summarized in Figure 2.  Remark 3. The coupled variables problem can be effectively solved by implementing the proposed DDC algorithm.Firstly, we express E g,ρ i and E g,1−ρ i without inequality constraints as the linear combination of λ e and λ g .Then, the direct calculation of E g,ρ i and E g,1−ρ i can be transformed into the indirect calculation of λ p and λ h , which can effectively handle the coupled variables existing in the objective function.Furthermore, along with the KKT conditions and the solution without inequality constraints, the proposed projection operation method maps infeasible values into the feasible operation region, which can further effectively handle the coupled variables existing within the inequality constraint limits.

Application Examples
In this section, several simulations have been studied to verify the effectiveness of the proposed algorithm.The first case is used to show the effectiveness of the proposed DDC algorithm via comparison to a centralized algorithm.The second case shows the performance of the proposed algorithm under time-varying demand.The third case demonstrates the plug and play properties of the proposed algorithm.The configuration and communication structures of the test system with five energy hubs are shown in Figure 3.The energy transformer efficiencies η ee , η e,chp , η h,chp and η boil are selected as 0.98, 0.35, 0.4 and 0.9, respectively.The cost parameters and constraints of each hub are listed in Table 1.Furthermore, along with the KKT conditions and the solution without inequality constraints, the proposed projection operation method maps infeasible values into the feasible operation region, which can further effectively handle the coupled variables existing within the inequality constraint limits.

Application Examples
In this section, several simulations have been studied to verify the effectiveness of the proposed algorithm.The first case is used to show the effectiveness of the proposed DDC algorithm via comparison to a centralized algorithm.The second case shows the performance of the proposed algorithm under time-varying demand.The third case demonstrates the plug and play properties of the proposed algorithm.The configuration and communication structures of the test system with five energy hubs are shown in Figure 3.The energy transformer efficiencies η ee , η e,chp , η h,chp and η boil are selected as 0.98, 0.35, 0.4 and 0.9, respectively.The cost parameters and constraints of each hub are listed in Table 1.

Comparison with Centralized Algorithm
This case study focuses on verifying the effectiveness of the proposed DDC algorithm by comparing it to the centralized algorithm [30].The power and heat load demands for each power and heat load are initialized at 150 kW and 140 kW, respectively.By running the DDC algorithm, the estimated electricity and heat Lagrangian multipliers, mismatches and energy inputs are shown in Figure 4.After the algorithm converges, it can be observed that the electricity multipliers converge to a common value, λ p = 27.0354cents/kWh, and the heat multipliers converge to another value, λ h = 19.3652cents/kWh; meanwhile, the electricity and heat mismatches converge to zero, which means that the power and heat supply-demand equality constraints are satisfied.In addition, the final electricity and gas inputs of each hub are within their corresponding inequality constraints, as shown in Table 2. Therefore, the optimization goal has been fulfilled.Furthermore, the performance of the proposed DDC algorithm is compared with a centralized algorithm.The simulation results provided by the two approaches are shown in Table 2.It can be seen that all of the solutions are very similar, which also verifies the effectiveness of the proposed algorithm.
This case study focuses on verifying the effectiveness of the proposed DDC algorithm by comparing it to the centralized algorithm [30].The power and heat load demands for each power and heat load are at 150 kW and 140 kW, respectively.By running the DDC algorithm, the estimated electricity and heat Lagrangian multipliers, mismatches and energy inputs are shown in Figure 4.After the algorithm converges, it can be observed that the electricity multipliers converge to a common value, 27.0354 λ = p cents/kWh, and the heat multipliers converge to another value, 19.3652 λ = h cents/kWh; meanwhile, the electricity and heat mismatches converge to zero, which means that the power and heat supply-demand equality constraints are satisfied.In addition, the final electricity and gas inputs of each hub are within their corresponding inequality constraints, as shown in Table 2. Therefore, the optimization goal has been fulfilled.Furthermore, the performance of the proposed DDC algorithm is compared with a centralized algorithm.The simulation results provided by the two approaches are shown in Table 2.It can be seen that all of the solutions are very similar, which also verifies the effectiveness of the proposed algorithm.

Time-Varying Demand
This case study focuses on testing the capability of the proposed DDC algorithm to handle changes in load demand.The initial conditions are the same as those in case study 1, and the load demand changes 2 times during the simulation.At t 1 (k = 1000), the total power and heat load demands decrease by 20%, while at t 2 (k = 2000), they are increased by 20%.The simulation results are shown in Figure 5.It can be observed from Figure 5a,b that neither the electricity nor heat mismatch can be maintained at a zero state after t 1 , because of the reduced power and heat load demands.By continually running the proposed DDC algorithm, each energy hub gradually reduces its electricity and gas inputs to meet the new supply-demand balance, as seen in Figure 5e,f; meanwhile, the electricity and heat multipliers decrease in response to the change seen in Figure 5c,d.After about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches converge to zero again, and not only the electricity Lagrangian multipliers, but also the heat Lagrangian multipliers, converge to new common values.This implies that the new supply-demand balance and the optimization goal have been fulfilled.Furthermore, the balance is broken at t 2 due to the increased load demands.From Figure 5a-f, it can be seen that after about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches go to zero, both the electricity and heat multipliers increase and converge to new optimal values, and each hub increases its energy input to compensate for part of the increased load, finally converging to the new solutions.Based on the aforementioned discussion, it can be concluded that the proposed DDC algorithm can respond automatically to converge to new optimal solutions after each load change.This also means that the proposed DDC algorithm works properly under time-varying demand.

Time-Varying Demand
This case study focuses on testing the capability of the proposed DDC algorithm to handle changes in load demand.The initial conditions are the same as those in case study 1, and the load demand changes 2 times during the simulation.At t1 (k = 1000), the total power and heat load demands decrease by 20%, while at t2 (k = 2000), they are increased by 20%.The simulation results are shown in Figure 5.It can be observed from Figure 5a,b that neither the electricity nor heat mismatch can be maintained at a zero state after t1, because of the reduced power and heat load demands.By continually running the proposed DDC algorithm, each energy hub gradually reduces its electricity and gas inputs to meet the new supply-demand balance, as seen in Figure 5e,f; meanwhile, the electricity and heat multipliers decrease in response to the change seen in Figure 5c,d.After about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches converge to zero again, and not only the electricity Lagrangian multipliers, but also the heat Lagrangian multipliers, converge to new common values.This implies that the new supply-demand balance and the optimization goal have been fulfilled.Furthermore, the balance is broken at t2 due to the increased load demands.From Figure 5a-f, it can be seen that after about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches go to zero, both the electricity and heat multipliers increase and converge to new optimal values, and each hub increases its energy input to compensate for part of the increased load, finally converging to the new solutions.Based on the aforementioned discussion, it can be concluded that the proposed DDC algorithm can respond automatically to converge to new optimal solutions after each load change.This also means that the proposed DDC algorithm works properly under time-varying demand.

Plug and Play Capability
This case study focuses on testing the plug and play performance of the proposed algorithm.The initial conditions are the same as those in case study 1.At t 1 (k = 1000), hub3 is removed from the test system and the variables related to them are set to zero.As a result of this change, the original balance is broken, and the electricity and heat mismatches go to a non-zero state, as seen in Figure 6a,b.By continually running the proposed DDC algorithm, it can be observed from Figure 6c-f that the remaining energy hubs have to increase their electricity and gas inputs to generate more electricity and heat to further compensate for the amount of electricity and heat previously generated by hub3; the corresponding electricity and heat Lagrangian multipliers increase, as well, to meet this change.After about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches converge to zero again; meanwhile, the electricity and heat Lagrangian multipliers, as well as the electricity and heat inputs, converge to their corresponding optimal values.As a result, the new supply-demand balance can be satisfied after removing hub3.Then, at t 2 (k = 2000), hub3 is plugged back into the system again.Of course, hub3 gradually increases its electricity and heat generation, while the others decrease their generations to accommodate for hub3.After about 300 iterations, i.e., at k = 2300, it can be seen from Figure 5a-f that the algorithm converges to the new optimal solution in response to the new topological change.Additionally, the final convergence values for all variables are the same as those prior to disconnection.Therefore, the results clearly show that our proposed algorithm can provide good plug and play capability.

Plug and Play Capability
This case study focuses on testing the plug and play performance of the proposed algorithm.The initial conditions are the same as those in case study 1.At t1 (k = 1000), hub3 is removed from the test system and the variables related to them are set to zero.As a result of this change, the original balance is broken, and the electricity and heat mismatches go to a non-zero state, as seen in Figure 6a,b.By continually running the proposed DDC algorithm, it can be observed from Figure 6c-f that the remaining energy hubs have to increase their electricity and gas inputs to generate more electricity and heat to further compensate for the amount of electricity and heat previously generated by hub3; the corresponding electricity and heat Lagrangian multipliers increase, as well, to meet this change.After about 300 iterations, i.e., at k = 1300, the electricity and heat mismatches converge to zero again; meanwhile, the electricity and heat Lagrangian multipliers, as well as the electricity and heat inputs, converge to their corresponding optimal values.As a result, the new supply-demand balance can be satisfied after removing hub3.Then, at t2 (k = 2000), hub3 is plugged back into the system again.Of course, hub3 gradually increases its electricity and heat generation, while the others decrease their generations to accommodate for hub3.After about 300 iterations, i.e., at k = 2300, it can be seen from Figure 5a-f that the algorithm converges to the new optimal solution in response to the new topological change.Additionally, the final convergence values for all variables are the same as those prior to disconnection.Therefore, the results clearly show that our proposed algorithm can provide good plug and play capability.

Conclusions
In this paper, we have presented a novel DDC algorithm to solve the EMP for multiple energy hubs, which is formulated as a class of distributed coupled optimization problem by properly converting some system coordinates.The consensus concept in multiagent systems has been explored and embedded in the design of our distributed algorithm.With this effort, the total computation process can be completely assigned to each individual energy hub without the requirements of a

Conclusions
In this paper, we have presented a novel DDC algorithm to solve the EMP for multiple energy hubs, which is formulated as a class of distributed coupled optimization problem by properly converting some system coordinates.The consensus concept in multiagent systems has been explored and embedded in the design of our distributed algorithm.With this effort, the total computation process can be completely assigned to each individual energy hub without the requirements of a cloud computing center, a central price coordinator or a virtual coordinator.Thus, our proposed algorithm can improve the existing game-theoretic-based distributed approaches, which further solves the EMP for multiple energy hubs in a fully distributed fashion.In addition, a novel projection operation method has also been proposed to address the issue of coupled variables existing in the objective function and constraint limits.Furthermore, we have proved that the proposed algorithm can converge to the optimal solution under strong connectivity conditions only.Column stack vectors form of λ i,e,p , λ i,e,c , λ i,g,c , y i,e,p , y i,e,c and y i,g,c x e , x g,ρ , x g,1−ρ , E e , E g,ρ , E g, where the dispatch factor ρ ∈ [0, 1] is introduced to specify the share of natural gas by the CHP unit and boiler.(b) The total objective function of the EMP for multiple energy hubs is modeled as the following form, The Lagrange function for problem (3)-( 6) is expressed as the following equation: The optimal operating point is determined by the following KKT conditions: and Equation (4a) along with Equation (4b). ( Since R is row stochastic, P1 forms the first-order consensus protocol.Each node reaches the final common value χ i (∞) = w T χ(0).Since S is column stochastic, the summation of all state variables is a constant, i.e., ∑ χ i (k) = ∑ χ i (0).(b) We let u e,min i = u e,max i = u g i,m = 0.Then, according to the first three equations of Equation ( 12), the optimal solution of E e i , E g,ρ i and E g,1−ρ i of each energy hub can be determined by the unique Lagrange multipliers λ e and λ e , given by where (0) can be designed as any admissible values, and the other variables are initialized as To map the infeasible values into the feasible operating region, the cases for no active constraints, one active constraint and two active constraints are separately discussed to further give the corresponding identification conditions as follows: (1) None of the constraints are active constraints, i.e., the optimal solutions before and after considering the inequality constraints have the same values.The identification condition of this case is if then there are no active constraints; otherwise, the optimal solution belongs to another case.(2) Only one inequality constraint is an active constraint.Let represent any one inequality equation.
The identification condition of this case is that if there exist solutions x g,ρ i, , , x g,1−ρ i, , and then is the only active constraint; otherwise, the optimal solution belongs to another case.
(3) Two inequality constraints are active constraints.In this case, let and represent two adjacent constraints (e.g., 2 and 3) and the corner point determined by them be defined as (x g,ρ i, , , x g,1−ρ i, , ).The identification condition of this case is that if there exist solutions u g i, ≥ 0 and u g i, ≥ 0 satisfying The eigenvalue perturbation approach is employed to analyze the convergence properties of the presented algorithm.Along with Equations ( 7)-( 10), we get where Θ = [λ T e,p , λ T e,c , λ T g,c , y T e,p , y T e,c , y T g,c ] T .
Due to the definition of R, R, S and S, it is not difficult to verify that R, R, S and S are all primitive with the maximum eigenvalue 1.Based on Lemma 1, there exist four vectors w T , v, w T and v, satisfying w T R = w T , Sv = v, w T R = w T and Sv = v meanwhile w T 1 = 1, 1 T v = 1, w T 1 = 1 and 1 T v = 1.Therein, 1 denotes a vector with all its elements being 1.To facilitate the analysis, w T and v are expressed in the corresponding block vectors form, i.e., w T = w T p , w T c and v = v T p , v T c T . Furthermore, M is a lower block triangular matrix whose eigenvalues are the union of R, R, S and S, so M has four maximum eigenvalues ψ 1 = ψ 2 = ψ 3 = ψ 4 = 1 while the other eigenvalues lie in the open unit disk on the complex plane.Then, the eigenvalue perturbation approach is used to analyze the eigenvalue changing of M after perturbing by η∆.Let κ pp = ∑ a p + ∑ a p,c , κ pc = ∑ o p,c , κ cp = ∑ o h,c v and κ hh = ∑ a h,c .On this basis, V and W T are defined as the right and left eigenvectors of M, respectively.In this way, W T V = I.Then, along with ∆, V and W T , the matrix W T ∆V can be mathematically represented as follows: Since the object function is a convex function, it is not difficult to verify κ pc κ cp < κ hh κ pp .According to Q, we have the four eigenvalues of W T ∆V satisfying dψ 1 /dη = dψ 2 /dη = 0, dψ 3 /dη < 0 and dψ 4 /dη < 0, which means ψ 1 and ψ 2 do not change against η, and ψ 3 and ψ 4 become smaller when η > 0. Thus, there exists an upper boundary ς such that 0 < η < ς.Except for ψ 1 = ψ 2 = 1, the remaining eigenvalues lie in the open unit disk.Moreover, it can be verified that 1 T , 1 T , 0 T , 0 T , 0 T , 0 T , 0 T T and 0 T , 0 T , 1 T , 1 T , 0 T , 0 T T are the two independent eigenvectors of M corresponding to ψ 1 and ψ 2 , respectively.That means that, as k → ∞ , all of y i,e,p (k), y i,e,c (k) and y i,g,c (k) converge to 0, i.e., the electricity and heat power supply-demand balance constraints are satisfied.At the same time, λ i,e,p (k), λ i,e,c (k) converge to a common value, while λ i,g,c (k) converges to another common value.Therefore, the optimality conditions (KKT conditions of optimality) are satisfied.

Figure 1 .
Figure 1.Structure of a multiple energy hubs system.

Figure 1 .
Figure 1.Structure of a multiple energy hubs system.

Figure 2 .
Figure 2. Flowchart of the proposed DDC algorithm.

Figure 2 .
Figure 2. Flowchart of the proposed DDC algorithm.

Figure 3 .
Figure 3. Structure of multiple energy hub system.

Figure 3 .
Figure 3. Structure of multiple energy hub system.
are the two active constraints; otherwise, the optimal solution belongs to another case.p(I − R pp ) −A p R pc 0 S pp S pc 0 −A pc R cp A pc (I − R cc ) −O pc (I − R) S cp S cc 0 O hc R cp −O hc (I − R cc ) A hc (I − R)

Agent 1 Agent 1 Agent 2 Agent 3 Distributed Fashion Agent n Centralized Fashion Centralized controller Agent 1 Agent 2 Agent 4 Agent 3 Agent n Paradigm Shift Plug-and-play Enhanced robustness Improved flexibility Cost-effectiveness Main supply
) ) can be transformed into the indirect calculation of p λ and h λ , which can effectively handle the coupled variables existing in the objective function. −

Table 1 .
Parameters of the test system.

Table 1 .
Parameters of the test system.

Table 2 .
Comparison with centralized algorithm.

Table 2 .
Comparison with centralized algorithm.
In the simulation part, several examples (NSFC) under Grant No. 61433004, the NSFC under Grant No. 61621004, the NSFC under Grant No. 61603085, and the National Natural Science Foundation of Liaoning Province under Grant No. 201601020.The authors declare no conflict of interest., λ e,c , λ g,c , y e,p , y e,c , y g,c c) The local constraints arise from the limitations of the energy hubs' capability and the dispatch factors, which