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

Abstract

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.

1. 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 real-time application in modern power systems [14,15,16]. In fact, various distributed methods proposed only solve the EDP on the power grid side. For example, an incremental cost consensus method was proposed in [8] to solve the EDP in a distributed fashion, where a leader agent is needed to satisfy the electrical power balance constraint. In [17], a fully distributed method was proposed, in which the communication among generators are strongly connected and the estimation of power mismatch does not require a central controller or leader. Furthermore, the transmission losses, demand response and ramping rate limits have been taken into account and discussed in [18,19,20,21,22], respectively. In addition, there are also some researches analyzing the EDP from other perspectives, such as non-convex dispatch model [23], online generation adjustment [24], etc. It is worth noting that the aforementioned distributed methods are mainly focusing on the economic dispatch problem on the power side. However, there exist strong power-heat-coupling features both in the objective function and constraint limits in the CHP system scenario. Hence, the distributed methods above cannot apply to this kind of distributed coupling optimization problem.

Based on above discussions, this paper presents a novel distributed method to solve the CHPEDP. The proposed method is implemented based on two consensus procedure modified their corresponding feedback parts where the one is used to get the optimal electrical incremental cost and the other is used to get the optimal heat incremental cost, resulting in achieving the optimal solutions. Meanwhile, the alternating iterative method is proposed and used to solve the electrical–heat-coupling problem. In addition, the estimation of the electrical and heat power are also through distributed fashion, which are then used as a feedback mechanism to meet the supply–demand balance. The salient features of this paper are listed as follows:

(1) Various kinds of energy conversion devices fed by multiple-energy-carriers are modeled and considered to analyze the electrical and heat power production in CHP system, which can effectively increase the energy efficiency.

(2) The proposed method is fully distributed without a leader or central controller, which means it is more reliable and robust. In addition, the requirement on the communication topology of this paper is a strong connection that possesses less conservatism compared with undirected, connected communication topology.

(3) The CHPEDP has been formulated into a class of coupling optimization problem, which can be effectively solved by the proposed method in a distributed fully fashion. To the best knowledge of the authors, there still exist very few papers concerned about solving the CHPEDP in a distributed manner.

The rest of this paper is organized as follows. In Section 2, typical structure of CHP system, agent definitions and modeling are introduced. Then, the centralized method to solve the CHPEDP is discussed. In Section 3, some basic knowledge is briefly introduced, and the proposed distributed method is presented. In Section 4, several case studies are presented to show the effectiveness of the proposed method. Section 5 concludes the paper.

2. Problem Formulation

2.1. Structure and Agent Definitions

A typical architecture of CHP system is shown in Figure 1. Fed by multiple-energy-carriers, the energy resources of the CHP system consist of available distributed generators (DGs) that contain renewable generators (RGs) and fuel generators (FGs), distributed heating devices (DHDs) that contain renewable heating devices (RHDs) and conventional fuel heating devices (FHDs), distributed energy storage devices (DESDs) that contain electrical energy storage devices (EESDs) and heat energy storage devices (HESDs), CHPs that contain renewable combined heat and power (RCHP) and conventional fuel combined heat and power (FCHP), and loads that contain electrical loads and heat loads.

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.

2.2. Modeling

The cost function for CGAs is often modeled as the following form [6]:

$$C_c^i(x_i^{pq},x_i^{qp})=a_i{(x_i^{pq})}^2+b_i{x}^{pq}+{\alpha }_i{(x_i^{qp})}^2+{\beta }_i{x}_i^{qp}+{\gamma }_i+{\mathsf{\xi }}_i{x}_i^{pq}{x}_i^{qp}$$

(1a)

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\ge 0$ ($m=1,2,3,4$). Furthermore, the set of linear constraints can be re-written as the following forms [4]:

$$x_{i\_\min }^{pq}(x_i^{qp})\le x_i^{pq}\le x_{i\_\max }^{pq}(x_i^{qp})$$

(1b)

$$x_{i\_\min }^{qp}(x_i^{pq})\le x_i^{qp}\le x_{i\_\max }^{qp}(x_i^{pq})$$

(1c)

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}\ge -\left(e_m{x}^{qp}+f_m\right)/d_m$. Note that there always exist one (or several) $e_m/d_m\ge 0$ and one (or several) $e_m/d_m<0$ for different m. Thus, $x^{pq}\ge -\left(e_m{x}^{qp}+f_m\right)/d_m$ can be re-written as the form $F^{\max }(x^{qp})\ge x^{pq}\ge 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_{\min }^{pq}(x^{qp})$ and $x_{\max }^{pq}(x^{qp})$ as $F^{\min }(x^{qp})$ and $F^{\max }(x^{qp})$, respectively. Similarly, we can define $x_{\min }^{qp}(x^{pq})$ and $x_{\max }^{qp}(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\ge 0$ is equivalent to the expressions of Equations (1b) and (1c).

The cost function of EOAs are modeled as the following form [25]:

$$C_e^i(x_i^p)=a_i{(x_i^p)}^2+b_i{x}_i^p+{\mathsf{\gamma }}_i$$

(2a)

$$x_{i\_\min }^p\le x_i^p\le x_{i\_\max }^p$$

(2b)

where $b_i$ and $c_i$ are both equal to zero for each EESD. Emphatically, $x_{i\_\max }^p=\min \{x_{i\_MPPT\_\max }^p,x_{i\_w\_\max }^p\}$ 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:

$$C_h^i(x_i^q)={\mathsf{\alpha }}_i{(x_i^q)}^2+{\mathsf{\beta }}_j{x}_i^q+{\mathsf{\gamma }}_i$$

(3a)

$$x_{i\_\min }^q\le x_i^q\le x_{i\_\max }^q$$

(3b)

where ${\mathsf{\beta }}_i$ and ${\mathsf{\gamma }}_i$ are both equal to zero for every HESD. Similarly, for RHDs, $x_{i\_\max }^q=\min \{x_{i\_MPPT\_\max }^q,x_{i\_w\_\max }^q\}$.

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:

$$C(x)={\displaystyle \sum _{i=1}^{n_1}{C}_e^i}(x_i^p)+{\displaystyle \sum _{i=n_1+1}^{n_1+n_2}{C}_c^i}(x_i^{pq},x_i^{qp})+{\displaystyle \sum _{i=n_1+n_2+1}^{n_1+n_2+n_3}{C}_h^i}(x_i^q)$$

(4)

subject to

$${\displaystyle \sum _{i=1}^{n_1}{x}_i^p}+{\displaystyle \sum _{i=n_1+1}^{n_1+n_2}{x}_i^{pq}}-D_P=0$$

(5a)

$${\displaystyle \sum _{i=n_1+1}^{n_1+n_2}{x}_i^{qp}}+{\displaystyle \sum _{i=n_1+n_2+1}^{n_1+n_2+n_3}{x}_i^q}-D_Q=0$$

(5b)

and Equations (1b), (1c), (2b) and (3b).

In this paper, when $i=1,\cdots ,n_1$, it expresses the EOA. Similarly, when $i=n_1+1,\cdots ,n_1+n_2$ and $i=n_1+n_2+1,\cdots ,n_1+n_2+n_3$, it expresses the CGA and HOA, respectively.

2.3. 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:

$$\begin{array}{c}L=C(x)-{\mathsf{\lambda }}_p({\displaystyle \sum _{i=1}^{n_1}{x}_i^p}+{\displaystyle \sum _{i=n_1+1}^{n_1+n_2}{x}_i^{pq}}-D_P)-{\mathsf{\lambda }}_q({\displaystyle \sum _{i=n_1+1}^{n_1+n_2}{x}_i^{qp}}+{\displaystyle \sum _{i=n_1+n_2+1}^{n_1+n_2+n_3}{x}_i^q}-D_Q)\\ -{\displaystyle \sum {\mathsf{\mu }}_{i\_\min }}(x_i-x_{i\_\min })-{\displaystyle \sum {\mathsf{\mu }}_{i\_\max }}(x_{i\_\max }-x_i)\end{array}$$

(6)

When the inequality constraints are not taken into account, the Lagrangian operator provides the following set of equations:

$$\{\begin{array}{l}\frac{\partial L}{\partial x_i^p}=2a_i{x}_i^p+b_i-{\mathsf{\lambda }}_p=0,\frac{\partial L}{\partial x_i^{pq}}=2a_i{x}_i^{pq}+b_i+{\mathsf{\xi }}_i{x}_i^{qp}-{\mathsf{\lambda }}_p=0\\ \frac{\partial L}{\partial x_i^{qp}}=2{\mathsf{\alpha }}_i{x}_i^{qp}+{\mathsf{\beta }}_i+{\mathsf{\xi }}_i{x}_i^{pq}-{\mathsf{\lambda }}_q=0,\frac{\partial L}{\partial x_i^q}=2{\mathsf{\alpha }}_i{x}_i^q+{\mathsf{\beta }}_i-{\mathsf{\lambda }}_q=0\end{array}$$

(7)

Let $\partial C_e^i(x_i^p)/\partial x_i^p=2a_i{x}_i^p+b_i$ and $\partial C_c^i(x_i^{pq},x_i^{qp})/x_i^{pq}=2a_i{x}_i^{pq}+b_i+{\mathsf{\xi }}_i{x}_i^{qp}$ be the electrical incremental cost, and $\partial C_c^i(x_i^{pq},x_i^{qp})/x_i^{qp}=2{\mathsf{\alpha }}_i{x}_i^{qp}+{\mathsf{\beta }}_i+{\mathsf{\xi }}_i{x}_i^{pq}$ and $\partial C_h^i(x_i^q)/x_i^q=2{\mathsf{\alpha }}_i{x}_i^q+{\mathsf{\beta }}_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 ${\mathsf{\lambda }}_p$ and ${\mathsf{\lambda }}_q$, respectively. If there exists the optimal operating point, the calculations of the electrical and heat power are given by:

$$\{\begin{array}{l}{x}_i^p=({\mathsf{\lambda }}_p^{*}-b_i)/2a_i,x_i^{pq}=({\mathsf{\lambda }}_p^{*}-b_i-{\mathsf{\xi }}_i{x}_i^{qp})/2a_i\\ x_i^{qp}=({\mathsf{\lambda }}_q^{*}-{\mathsf{\beta }}_i-{\mathsf{\xi }}_i{x}_i^{pq})/2{\mathsf{\alpha }}_i,x_i^q=({\mathsf{\lambda }}_q^{*}-{\mathsf{\beta }}_i)/2{\mathsf{\alpha }}_i\end{array}$$

(8)

Furthermore, when the inequality constraints are taken into account, the necessary conditions to solve the problem can be expanded as:

$$\{\begin{array}{ll}\partial C^i/\partial x_i^{\hslash }={\mathsf{\lambda }}_{\nu }^{*},& x_{i\_\min }^{\hslash }<x_i^{\hslash }<x_{i\_\max }^{\hslash }\\ \partial C^i/\partial x_i^{\hslash }\le {\mathsf{\lambda }}_{\nu }^{*},& x_{i\_\max }^{\hslash }=x_i^{\hslash }\\ \partial C^i/\partial x_i^{\hslash }\ge {\mathsf{\lambda }}_{\nu }^{*},& x_{i\_\min }^{\hslash }=x_i^{\hslash }\end{array}$$

(9)

where if $\hslash =p$, ${\mathsf{\lambda }}_{\nu }^{*}$, $x_{i\_\min }^{\hslash }$ and $x_{i\_\max }^{\hslash }$ denote ${\mathsf{\lambda }}_p^{*}$, $x_{i\_\min }^p$ and $x_{i\_\max }^p$, respectively. If $\hslash =pq$, ${\mathsf{\lambda }}_{\nu }^{*}$, $x_{i\_\min }^{\hslash }$ and $x_{i\_\max }^{\hslash }$ denote ${\mathsf{\lambda }}_p^{*}$, $x_{i\_\min }^{pq}(x_i^{qp})$ and $x_{i\_\max }^{pq}(x_i^{qp})$, respectively. If $\hslash =qp$, ${\mathsf{\lambda }}_{\nu }^{*}$, $x_{i\_\min }^{\hslash }$ and $x_{i\_\max }^{\hslash }$ denote ${\mathsf{\lambda }}_q^{*}$, $x_{i\_\min }^{qp}(x_i^{pq})$ and $x_{i\_\max }^{qp}(x_i^{pq})$, respectively. If $\hslash =q$, ${\mathsf{\lambda }}_{\nu }^{*}$, $x_{i\_\min }^{\hslash }$ and $x_{i\_\max }^{\hslash }$ denote ${\mathsf{\lambda }}_q^{*}$, $x_{i\_\min }^q$ and $x_{i\_\max }^q$, 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 ${\mathsf{\xi }}_i$ of Equation (1a), and in the feasible operating region of CGAs mathematically represented as Equations (1b) and (1c). The former results in the calculation of $x_i^{pq}$ ($x_i^{qp}$) depend on not only ${\mathsf{\lambda }}_p^{*}$ (${\mathsf{\lambda }}_q^{*}$) but also $x_i^{qp}$ ($x_i^{pq}$), as shown in Equation (8). The latter expresses that the lower and upper bounds of $x_i^{pq}$ ($x_i^{qp}$) are not two constants but turning to two different functions which collectively depend on $x_i^{pq}$ and $x_i^{qp}$ as shown in Equation (9). Thus, the optimal solutions cannot be separately calculated by the power system or heat system.

3. Distributed Method

3.1. Basic Knowledge

(1) Graph Theory: A weighted directed graph $G=(V,E)$ is composed of a set of nodes $V=(v_1,v_2,\cdots ,v_n)$ and a set of arcs $E=\{e_{ij}=(v_i,v_j)\}\subseteq V\times V$. Therein, the $i$th agent can receive information from $j$th 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 (out-neighbors) is denoted by $N_i^{in}$ ($N_i^{out}$) with cardinality $d_i^{in}$ ($d_i^{out}$). 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 stochastic. Let $r_{i,j}=1/d_i^{in}$ if $j\in N_i^{in}$, $r_{i,j}=0$ if $j\notin N_i^{in}$, and $s_{i,j}=1/d_j^{out}$ if $i\in N_j^{out}$, $s_{i,j}=0$ if $i\notin N_j^{out}$. Similarly, define two matrices $\overline{R}$ and $\overline{S}$ in subgraph $G_{HC}$, where $\overline{R}$ is row stochastic and $\overline{S}$ is column stochastic. Let ${\overline{r}}_{i,j}=1/{\overline{d}}_i^{in}$ if $j\in {\overline{N}}_i^{in}$, ${\overline{r}}_{i,j}=0$ if $j\notin {\overline{N}}_i^{in}$, and ${\overline{s}}_{i,j}=1/{\overline{d}}_j^{out}$ if $i\in {\overline{N}}_j^{out}$, ${\overline{s}}_{i,j}=0$ if $i\notin {\overline{N}}_j^{out}$.

(3) Consensus algorithm (or protocol): Based on the definitions discussed above, the first-order consensus protocol is given by:

$${\chi }_i\left(t+1\right)={\displaystyle \sum _{j\in N^{in}}{r}_{i,j}}{\chi }_j\left(t\right)$$

(10)

Since R is row stochastic, Equation (10) forms the first-order consensus protocol [26]. Each agent reaches the final common value ${\chi }_i\left(\infty \right)$$={\mathit{\omega }}^T\mathit{\chi }\left(0\right)$.

3.2. 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 ${\mathsf{\lambda }}_p^{*}$ but also an optimal heat incremental cost ${\mathsf{\lambda }}_q^{*}$. More importantly, ${\mathsf{\lambda }}_p^{*}$ and ${\mathsf{\lambda }}_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 ${\mathsf{\lambda }}_p^{*}$ and ${\mathsf{\lambda }}_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

$$\{\begin{array}{l}{\mathsf{\lambda }}_i^p(2k)={\mathsf{\lambda }}_i^p(2k-1),{\mathsf{\lambda }}_i^{pq}(2k)={\mathsf{\lambda }}_i^{pq}(2k-1)\\ {\mathsf{\lambda }}_i^{qp}(2k)={\displaystyle \sum _{j\in {\overline{N}}_i^{in}}{\overline{r}}_{i,j}{\mathsf{\lambda }}_j^{\mathsf{\vartheta }}(2k-1)}+\mathsf{\eta }{y}_i^{qp}(2k-1),{\mathsf{\lambda }}_i^q(2k)={\displaystyle \sum _{j\in {\overline{N}}_i^{in}}{\overline{r}}_{i,j}{\mathsf{\lambda }}_j^{\mathsf{\vartheta }}(2k-1)}+\mathsf{\eta }{y}_i^q(2k-1)\end{array}$$

(11a)

$$\{\begin{array}{l}{z}_i^p(2k)=x_i^p(2k-1),z_i^{pq}(2k)=x_i^{pq}(2k-1)\\ z_i^{qp}(2k)=({\mathsf{\lambda }}_i^{qp}(2k)-{\mathsf{\beta }}_i-{\mathsf{\xi }}_i{x}_i^{pq}(2k-1))/2{\mathsf{\alpha }}_i,z_i^q(2k)=({\mathsf{\lambda }}_i^q(2k)-{\mathsf{\beta }}_i)/2{\mathsf{\alpha }}_i\end{array}$$

(11b)

$$\{\begin{array}{l}{x}_i^p(2k)=z_i^p(2k),x_i^{pq}(2k)=z_i^{pq}(2k)\\ x_i^{qp}(2k)=\{\begin{array}{l}{z}_{i\_\min }^{qp}(2k)(z_i^{pq}(2k))\\ if{z}_i^{qp}(2k)<z_{i\_\min }^{qp}(2k)(z_i^{pq}(2k))\\ z_{i\_\max }^{qp}(2k)(z_i^{pq}(2k))\\ if{z}_i^{qp}(2k)>z_{i\_\max }^{qp}(2k)(z_i^{pq}(2k))\\ z_i^{qp}(2k)\mathrm{otherwise}\end{array},x_i^q(2k)=\{\begin{array}{l}{x}_{i\_\min }^{q}if{z}_i^q(2k)<x_{i\_\min }^q\\ x_{i\_\max }^{q}if{z}_i^q(2k)>x_{i\_\max }^q\\ z_i^q(2k)\mathrm{otherwise}\end{array}\end{array}$$

(11c)

$$\{\begin{array}{l}{y}_i^p(2k)=y_i^p(2k-1),y_i^{pq}(2k)=y_i^{pq}(2k-1)\\ y_i^{qp}(2k)={\displaystyle \sum _{j\in {\overline{N}}_i^{out}}{\overline{s}}_{i,j}}{y}_j^{\mathsf{\vartheta }}(2k-1)-x_i^{qp}(2k)+x_i^{qp}(2k-1)\\ y_i^q(2k)={\displaystyle \sum _{j\in {\overline{N}}_i^{out}}{\overline{s}}_{i,j}}{y}_j^{\mathsf{\vartheta }}(2k-1)-x_i^q(2k)+x_i^q(2k-1)\end{array}$$

(11d)

(2) In the odd-numbered iteration

$$\{\begin{array}{l}{\mathsf{\lambda }}_i^p(2k+1)={\displaystyle \sum _{j\in N_i^{in}}{r}_{i,j}{\mathsf{\lambda }}_j^{\mathsf{\theta }}(2k)}+\mathsf{\eta }{y}_i^p(2k),{\mathsf{\lambda }}_i^{pq}(2k+1)={\displaystyle \sum _{j\in N_i^{in}}{r}_{i,j}}{\mathsf{\lambda }}_j^{\mathsf{\theta }}(2k)+\mathsf{\eta }{y}_i^{pq}(2k)\\ {\mathsf{\lambda }}_i^{qp}(2k+1)={\mathsf{\lambda }}_i^{qp}(2k),{\mathsf{\lambda }}_i^q(2k+1)={\mathsf{\lambda }}_i^q(2k)\end{array}$$

(12a)

$$\{\begin{array}{l}{z}_i^p(2k+1)=({\mathsf{\lambda }}_i^p(2k+1)-b_i)/2a_i,z_i^{pq}(2k+1)=({\mathsf{\lambda }}_i^{pq}(2k+1)-b_i-{\mathsf{\xi }}_i{x}_i^{qp}(2k))/2a_i\\ z_i^{qp}(2k+1)=x_i^{qp}(2k),z_i^q(2k+1)=x_i^q(2k)\end{array}$$

(12b)

$$\{\begin{array}{l}{x}_i^p(2k+1)=\{\begin{array}{l}{x}_{i\_\min }^{p}if{z}_i^p(2k+1)<x_{i\_\min }^p\\ x_{i\_\max }^{p}if{z}_i^p(2k+1)>x_{i\_\max }^p\\ z_i^p(2k+1)\mathrm{otherwise}\end{array},x_i^{pq}(2k+1)=\{\begin{array}{l}{z}_{i\_\min }^{pq}(2k+1)(z_i^{qp}(2k+1))\\ if{z}_i^{pq}(2k+1)<\\ z_{i\_\min }^{pq}(2k+1)(z_i^{qp}(2k+1))\\ z_{i\_\max }^{pq}(2k+1)(z_i^{qp}(2k+1))\\ if{z}_i^{pq}(2k+1)>\\ z_{i\_\max }^{pq}(2k+1)(z_i^{qp}(2k+1))\\ z_i^{pq}(2k+1)\mathrm{otherwise}\end{array}\\ x_i^{qp}(2k+1)=z_i^{qp}(2k+1),x_i^q(2k+1)=z_i^q(2k+1)\end{array}$$

(12c)

$$\{\begin{array}{l}{y}_i^p(2k+1)={\displaystyle \sum _{j\in N_i^{out}}{s}_{i,j}}{y}_j^{\theta }(2k)-x_i^p(2k+1)+x_i^p(2k)\\ y_i^{pq}(2k+1)={\displaystyle \sum _{j\in N_i^{out}}{s}_{i,j}}{y}_j^{\theta }(2k)-x_i^{pq}(2k+1)+x_i^{pq}(2k)\\ y_i^{qp}(2k+1)=y_i^{qp}(2k),y_i^q(2k+1)=y_i^q(2k)\end{array}$$

(12d)

where $\mathsf{\theta }$ denotes $p$ or $pq$, $\mathsf{\vartheta }$ denotes $qp$ or $q$, and $t=2k\mathrm{or}2k+1$. In addition, we let $\mathsf{\Phi }=[{\mathsf{\lambda }}_i^{\theta },z_i^{\theta },x_i^{\theta },y_i^{\theta }]$ ($\mathsf{\Theta }=[{\mathsf{\lambda }}_i^{\mathsf{\vartheta }},z_i^{\mathsf{\vartheta }},x_i^{\mathsf{\vartheta }},y_i^{\mathsf{\vartheta }}]$) represent all the variables relating to heat (electrical) power.

Equation (11) shows the updating rules for $\mathsf{\Theta }(2k)$. In this step, the purpose is to make $\mathsf{\Theta }$ go to the optimal solution, meanwhile $\mathsf{\Phi }(2k)=\mathsf{\Phi }(2k-1)$ so that they can be seen as unchanged for further updating $\mathsf{\Theta }$. The detailed updating processes are as the following statements. Firstly, ${\mathsf{\lambda }}_i^{qp}(2k)$ and ${\mathsf{\lambda }}_i^q(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 $\overline{R}$ which drives ${\mathsf{\lambda }}_i^{qp}(2k)$ and ${\mathsf{\lambda }}_i^q(2k)$ to a common value. The feedback part includes $\mathsf{\eta }{y}_i^{qp}(2k-1)$ and $\mathsf{\eta }{y}_i^q(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 ${\mathsf{\lambda }}_q^{*}$. That is, each HOA or CGA decrease (or increase) its corresponding heat incremental cost and heat power output when $y_i^q(2k-1)$ or $y_i^{qp}(2k-1)$ is negative (or positive) for the common value further adjustment to ${\mathsf{\lambda }}_q^{*}$. Secondly, according to Equation (8), $z_i^{qp}(2k)$ and $z_i^q(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_{i\_w\_\min }^{pq}\le x_i^{pq}(2k-1)\le x_{i\_w\_\max }^{pq}$ must be satisfied. Finally, inspired by literature [21], $y_i^{qp}(2k)$ and $y_i^q(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 $\mathsf{\Phi }(2k+1)$ meanwhile $\mathsf{\Theta }(2k+1)=\mathsf{\Theta }(2k)$. ${\mathsf{\lambda }}_i^p(2k+1)$ and ${\mathsf{\lambda }}_i^{pq}(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., $\mathsf{\eta }{y}_i^p(2k)$ and $\mathsf{\eta }{y}_i^{pq}(2k)$) to obtain the optimal ${\mathsf{\lambda }}_p^{*}$. $z_i^p(2k+1)$, $z_i^{pq}(2k+1)$ in Equation (12b) and $x_i^p(2k+1)$, $x_i^{pq}(2k+1)$ in Equation (12c) are used to update the estimated electrical power outputs before and after considering the inequality constraints, respectively, on the basis of the second assumption $x_{i\_w\_\min }^{qp}\le x_i^{qp}(2k)\le x_{i\_w\_\max }^{qp}$. $y_i^p(2k+1)$ and $y_i^{pq}(2k+1)$ in Equation (12d) are used to update the estimated the electrical power mismatch.

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., ${\mathsf{\lambda }}_i^q\to {\mathsf{\lambda }}_q^{*}$, ${\mathsf{\lambda }}_i^{qp}\to {\mathsf{\lambda }}_q^{*}$, ${\mathsf{\lambda }}_i^p\to {\mathsf{\lambda }}_p^{*}$, ${\mathsf{\lambda }}_i^{pq}\to {\mathsf{\lambda }}_p^{*}$, $x_i^q\to x_i^{q*}$, $x_i^{qp}\to x_i^{qp*}$, $x_i^p\to x_i^{p*}$, $x_i^{pq}\to x_i^{pq*}$ as $t\to \infty$ (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., $x_i^{pq}(0)$ and $x_i^{qp}(0)$, can be set to any values in the corresponding feasible operating region. For convenience, we set the value of point B ($x_{i\_B}^{pq}$, $x_{i\_B}^{qp}$) in Figure 2 as the initial value, i.e., $x_i^{pq}(0)=x_{i\_B}^{pq}$ and $x_i^{qp}(0)=x_{i\_B}^{qp}$. In addition, if bus $i$ has electrical or heat load only, then it can be assumed that there is a fictitious EOA (or HOA) with $x_{i\_\min }^p=x_{i\_\max }^p=0$ (or $x_{i\_\min }^q=x_{i\_\max }^q=0$). Then, for each EOA (or HOA), the initial value of electrical (or heat) power output, i.e., $x_i^p(0)$ (or $x_i^q(0)$, can be any value. Furthermore, we set the corresponding lower bound as its initial value. The initial values of electrical and heat incremental costs, i.e., ${\mathsf{\lambda }}_i^p(0)$, ${\mathsf{\lambda }}_i^{pq}(0)$, ${\mathsf{\lambda }}_i^{qp}(0)$ and ${\mathsf{\lambda }}_i^q(0)$, which can also be any values, are set as zeros in this paper. In addition, the initial electrical (or heat) power mismatches, i.e., $y_i^p(0)$ and $y_i^{pq}(0)$ (or $y_i^q(0)$ and $y_i^{qp}(0)$, 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:

$$\{\begin{array}{l}{x}_i^p(0)=x_{i\_\min }^p,x_i^{pq}(0)=x_{i\_B}^{pq},x_i^{qp}(0)=x_{i\_B}^{qp},x_i^q(0)=x_{i\_\min }^q\\ {\mathsf{\lambda }}_i^p(0)={\mathsf{\lambda }}_i^{pq}(0)={\mathsf{\lambda }}_i^{qp}(0)={\mathsf{\lambda }}_i^q(0)=0\\ y_i^p(0)=L_{i\_E}-x_i^p(0),y_i^{pq}(0)=L_{i\_E}-x_i^{pq}(0),y_i^{qp}(0)=L_{i\_H}-x_i^{qp}(0),y_i^q(0)=L_{i\_H}-x_i^q(0)\end{array}$$

(13)

Remark 2:

There are two assumptions which must be satisfied during the updating process. In the even-numbered iteration, we have $x_i^{pq}(2k-1)=z_i^{pq}(2k)=x_i^{pq}(2k)$. Thus, if the first assumption $x_{i\_w\_\min }^{pq}\le x_i^{pq}(2k-1)\le x_{i\_w\_\max }^{pq}$ is satisfied, then $x_{i\_w\_\min }^{pq}\le z_i^{pq}(2k)\le x_{i\_w\_\max }^{pq}$ which indicates $(z_i^{pq}(2k),z_i^{qp}(2k))$ is within the area between L1 and L2 as shown in Figure 2. Due to the calculation of $z_i^{qp}(2k)$ without considering the inequality constraints, $(z_i^{pq}(2k),z_i^{qp}(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_i^{qp}(2k)$ can be easily obtained, because $z_i^{pq}(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_i^{pq}(2k),x_i^{qp}(2k))$, e.g., the point G shown in Figure 2. It is not very difficult to verify that $x_i^{qp}(2k)$ satisfies $x_{i\_w\_\min }^{qp}\le x_i^{qp}(2k)\le x_{i\_w\_\max }^{qp}$. Thus, in the next iteration, the second assumption is satisfied. Then, $(z_i^{pq}(2k+1),z_i^{qp}(2k+1))$ is within the area between L3 and L4. By using the projection operator 12(c), $(x_i^{pq}(2k+1),x_i^{qp}(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_i^{pq}(0)$, $x_i^{qp}(0)$) is in the feasible operating region.

For clarify, the proposed distributed method is summarized as:

Distributed Consensus-Based Method

Initialize: Each agent initializes all of variables based on Equation (13)

Iteration: ($t>0$)

(1).  Each agent exchanges the local information of ${\mathsf{\lambda }}_i^p$, ${\mathsf{\lambda }}_i^{pq}$, ${\mathsf{\lambda }}_i^{qp}$ and ${\mathsf{\lambda }}_i^q$ (and $y_i^p$, $y_i^{pq}$, $y_i^{qp}$ and $y_i^q$) with its neighbors. (2).  If $t=2k$ (i.e., in the even-numbered iteration) (2a).  Update ${\mathsf{\lambda }}_i^p$, ${\mathsf{\lambda }}_i^{pq}$, ${\mathsf{\lambda }}_i^{qp}$ and ${\mathsf{\lambda }}_i^q$ based on Equation (11a) (2b).  Update $z_i^p$, $z_i^{pq}$, $z_i^{qp}$ and $z_i^q$ (i.e., the estimated electrical and heat power outputs before considering inequality constraints) based on Equation (11b) (2c).  Choose projection operations based on Equation (11c) to calculate $x_i^p$, $x_i^{pq}$, $x_i^{qp}$ and $x_i^q$ (i.e., the estimated electrical and heat power outputs after considering inequality constraints) (2d).  Update $y_i^p$, $y_i^{pq}$, $y_i^{qp}$ and $y_i^q$ based on Equation (11d) (3).  Else (i.e., in the odd-numbered iteration) (3a).  Update ${\mathsf{\lambda }}_i^p$, ${\mathsf{\lambda }}_i^{pq}$, ${\mathsf{\lambda }}_i^{qp}$ and ${\mathsf{\lambda }}_i^q$ based on Equation (12a) (3b).  Update $z_i^p$, $z_i^{pq}$, $z_i^{qp}$ and $z_i^q$ based on Equation (12b) (3c).  Choose projection operations based on Equation (12c) to calculate $x_i^p$, $x_i^{pq}$, $x_i^{qp}$ and $x_i^q$ (3d).  Update $y_i^p$, $y_i^{pq}$, $y_i^{qp}$ and $y_i^q$ based on Equation (12d) (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 ${\mathsf{\lambda }}_i^p$, ${\mathsf{\lambda }}_i^{pq}$, ${\mathsf{\lambda }}_i^{qp}$ and ${\mathsf{\lambda }}_i^q$ (and $y_i^p$, $y_i^{pq}$, $y_i^{qp}$ and $y_i^q$) 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.

4. 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. Parameters of the 16-bus system.

Agent	$\mathit{a}$	$\mathit{b}$	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\epsilon }$	Min	Max
EOA1	0.0174	5.5	-	-	-	60	180
EOA2	0.1880	-	-	-	-	−75	75
EOA3	0.0124	6.4	-	-	-	50	150
EOA4	0.0004	0.106	-	-	-	0	140
EOA5	0.0146	5.8	-	-	-	40	120
EOA6	0.0009	0.02	-	-	-	0	180
CGA1	0.0069	2.9	0.0060	0.84	0.0062	-	-
CGA2	0.0087	5.0	0.0054	0.12	0.0022	-	-
HOA1	-	-	0.0102	3.3	-	40	500
HOA2	-	-	0.0146	2.42	-	30	300
HOA3	-	-	0.0001	0.084	-	0	250
HOA4	-	-	0.1660	-	-	−200	200

. 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.

4.1. 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 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. Comparison for the 16-bus system.

Agent	Case Study 1: Proposed Method	Case Study 1: Centralized Method	Case Study 2: Proposed Method
EOA1	64.1355	64.1987	64.1676
EOA2	20.5636	20.5695	20.5666
EOA3	53.7078	53.7950	53.7544
EOA4	90.0000	90.0000	90.0000
EOA5	66.1650	66.2368	66.2062
EOA6	130.0000	130.0000	130.0000
CGA1-E	215.1412	215.0000	215.0399
CGA2-E	110.2749	110.2000	110.2529
CGA1-H	179.2058	180.0000	179.7754
CGA2-H	134.9488	135.6000	135.1400
HOA1	150.9958	150.1772	150.5662
HOA2	135.6291	135.0553	135.32714
HOA3	180.0000	180.0000	180.0000
HOA4	19.2180	19.1675	19.1914
Total cost	5811.3	5806.9	5809.3

. 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.

4.2. 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.

4.3. 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_{11}^{pq}(800)=44$ (p.u.), $x_{11}^{qp}(800)=75$ (p.u.), $y_{11}^{pq}(800)=-44$ (p.u.), $y_{11}^{qp}(800)=-75$ (p.u.), ${\mathsf{\lambda }}_{11}^{pq}(800)={\mathsf{\lambda }}_{11}^{qp}(800)=0$. Then, the system again converges to the new solutions responding to the new topological change. 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.

4.4. 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.

5. 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.