## 2 Problem Formulation
### 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}{l}L=C(x)-\lambda_{p}(\sum\limits_{i=1}^{n_{1}}x_{i}^{p}+\sum \limits_{i=n_{1}+1}^{n_{1}+n_{2}}x_{i}^{pq}-D_{P})-\lambda_{q}(\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})\\ \qquad\qquad\qquad-\sum\limits_{i=n_{1}\text{\_min}}(x_{i}-x_{i\text{\_min}})- \sum\limits_{i\text{\_max}}(x_{i\text{\_max}}-x_{i})\end{array}  \tag{6}\]  
When the inequality constraints are not taken into account, the Lagrangian operator provides the following set of equations:  
\[\left\{ \begin{array}{l}\frac{\partial L}{\partial x_{i}^{p}}=2a_{i}x_{i}^{p} +b_{i}-\lambda_{p}=0,\qquad\frac{\partial L}{\partial x_{i}^{pq}}=2a_{i}x_{i}^ {pq}+b_{i}+\xi_{i}x_{i}^{qp}-\lambda_{p}=0\\ \frac{\partial L^{i}}{\partial x_{i}^{pq}}=2\alpha_{i}x_{i}^{qp}+\beta_{i}+\xi _{i}x_{i}^{pq}-\lambda_{q}=0,\qquad\frac{\partial L}{\partial x_{i}^{pq}}=2 \alpha_{i}x_{i}^{q}+\beta_{i}-\lambda_{q}=0\end{array} \right. \tag{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}+ \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\alpha_{i}x_{i}^{qp}+ \beta_{i}+\xi_{i}x_{i}^{pq}\) and \(\partial C_{h}^{i}(x_{i}^{q})/x_{i}^{q}=2\alpha_{i}x_{i}^{q}+\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 \(\lambda_{p}\) and \(\lambda_{q}\), respectively. If there exists the optimal operating point, the calculations of the electrical and heat power are given by:  
\[\left\{ \begin{array}{l}x_{i}^{p}=(\lambda_{p}^{*}-b_{i})/2a_{i},x_{i}^{pq}=( \lambda_{p}^{*}-b_{i}-\xi_{i}x_{i}^{qp})/2a_{i}\\ x_{i}^{qp}=(\lambda_{q}^{*}-\beta_{i}-\xi_{i}x_{i}^{pq})/2\alpha_{i},x_{i}^{ q}=(\lambda_{q}^{*}-\beta_{i})/2\alpha_{i}\end{array} \right. \tag{8}\]  
Furthermore, when the inequality constraints are taken into account, the necessary conditions to solve the problem can be expanded as:  
\[\left\{ \begin{array}{l}\partial C^{i}/\partial x_{i}^{h}=\lambda_{v}^{*}, \quad x_{i\text{\_min}}^{h}<x_{i}^{h}<x_{i\text{\_max}}^{h}\\ \partial C^{i}/\partial x_{i}^{h}\leq\lambda_{v}^{*},\qquad\quad x_{i\text{\_ max}}^{h}=x_{i}^{h}\\ \partial C^{i}/\partial x_{i}^{h}\geq\lambda_{v}^{*},\qquad\quad x_{i\text{\_ min}}^{h}=x_{i}^{h}\end{array} \right. \tag{9}\]  
where if \(\hbar=p\), \(\lambda_{v}^{*}\), \(x_{i\text{\_min}}^{h}\) and \(x_{i\text{\_max}}^{h}\) denote \(\lambda_{p}^{*}\), \(x_{i\text{\_min}}^{p}\) and \(x_{i\text{\_max}}^{p}\), respectively. If \(\hbar=pq\), \(\lambda_{v}^{*}\), \(x_{i\text{\_min}}^{h}\) and \(x_{i\text{\_max}}^{h}\) denote \(\lambda_{p}^{*}\), \(x_{i\text{\_min}}^{pq}(x_{i}^{qp})\) and \(x_{i\text{\_max}}^{pq}(x_{i}^{qp})\), respectively. If \(\hbar=qp\), \(\lambda_{v}^{*}\), \(x_{i\text{\_min}}^{h}\) and \(x_{i\text{\_max}}^{h}\) denote \(\lambda_{q}^{*}\), \(x_{i\text{\_min}}^{qp}(x_{i}^{pq})\) and \(x_{i\text{\_max}}^{qp}(x_{i}^{pq})\), respectively. If \(\hbar=q\), \(\lambda_{v}^{*}\), \(x_{i\text{\_min}}^{h}\) and \(x_{i\text{\_max}}^{h}\) denote \(\lambda_{q}^{*}\), \(x_{i\text{\_min}}^{qp}\) and \(x_{i\text{\_max}}^{qp}(x_{i}^{pq})\), respectively. If \(\hbar=q\), \(\lambda_{v}^{*}\), \(x_{i\text{\_min}}^{h}\) and \(x_{i\text{\_max}}^{h}\) denote \(\lambda_{q}^{*}\), \(x_{i\text{\_min}}^{q}\) and \(x_{i\text{\_max}}^{q}\), respectively.