## 2 Problem Formulation
### 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}+\vec{\xi}_{i}x_{i}^{pq}x_{i}^{qp} \tag{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}\geq 0\) (\(m=1,2,3,4\)). Furthermore, the set of linear constraints can be re-written as the following forms [4]:  
\[x_{i\_{m}\_{m}}^{pq}(x_{i}^{qp}) \leq x_{i}^{pq}\leq x_{i\_{m}\_{m}}^{pq}(x_{i}^{qp}) \tag{1b}\] \[x_{i\_{m}\_{m}}^{qp}(x_{i}^{pq}) \leq x_{i}^{qp}\leq x_{i\_{m}\_{m}}^{qp}(x_{i}^{pq}) \tag{1c}\]  
Figure 1: A structure of combined heat and power (CHP) system.  
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}\geq-\left(e_{m}x^{qp}+f_{m}\right)/d_{m}\). Note that there always exist one (or several) \(e_{m}/d_{m}\geq 0\) and one (or several) \(e_{m}/d_{m}<0\) for different \(m\). Thus, \(x^{pq}\geq-\left(e_{m}x^{qp}+f_{m}\right)/d_{m}\) can be re-written as the form \(F^{\max}(x^{qp})\geq x^{pq}\geq 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}\geq 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}+\gamma_{i} \tag{2a}\] \[x_{i\_\min}^{p}\leq x_{i}^{p}\leq x_{i\_\max}^{p} \tag{2b}\]  
where \(b_{i}\) and \(c_{i}\) are both equal to zero for each EESD. Emphatically, \(x_{i\_\max}^{p}=\min\left\{x_{i\_MPPT\_\max}^{p},x_{i\_\_\max}^{p}\right\}\) 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})=\alpha_{i}{(x_{i}^{q})}^{2}+\beta_{j}x_{i}^{q}+\gamma_{i} \tag{3a}\] \[x_{i\_\_\min}^{q}\leq x_{i}^{q}\leq x_{i\_\_\max}^{q} \tag{3b}\]  
where \(\beta_{i}\) and \(\gamma_{i}\) are both equal to zero for every HESD. Similarly, for RHDs, \(x_{i\_\_\max}^{q}=\min\left\{x_{i\_MPPT\_\max}^{q},x_{i\_\_\max}^{q}\right\}\).  
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)=\sum_{i=1}^{n_{1}}C_{e}^{i}(x_{i}^{p})+\sum_{i=n_{1}+1}^{n_{1}+n_{2}}C_{ e}^{i}(x_{i}^{pq},x_{i}^{qp})+\sum_{i=n_{1}+n_{2}+1}^{n_{1}+n_{2}+n_{3}}C_{h}^{i }(x_{i}^{q}) \tag{4}\]  
Figure 2: Feasible operating region of co-generation agents (CGA).  
subject to  
\[\sum\limits_{i=1}^{n_{1}}x_{i}^{p}+\sum\limits_{i=n_{1}+1}^{n_{1}+n_{2}}x_{i}^{ pq}-D_{P}=0 \tag{5a}\] \[\sum\limits_{i=n_{1}+1}^{n_{1}+n_{2}}x_{i}^{qp}+\sum\limits_{i=n_{1}+n_{2}+1}^{n _{1}+n_{2}+n_{3}}x_{i}^{q}-D_{Q}=0 \tag{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.