## Nomenclature  
\begin{tabular}{l l} CHP & Combined heat and power \\ EDP & Economic dispatch problem \\ DGs & Distributed generators \\ RGs & Renewable generators \\ FGs & Conventional fuel generators \\ DHDs & Distributed heating devices \\ RHDs & Renewable heating devices \\ FHDs & Conventional fuel heating devices \\ DESDs & Distributed energy storage devices \\ EESDs & Electrical energy storage devices \\ HESDs & Heat energy storage devices \\ RCHPs & Renewable combined heat and power \\ FCHPs & Fuel combined heat and power \\ EOAs & Electricity-only agents \\ CGAs & Co-generation agents \\ HOAs & Heat-only agents \\ \(i,j\) & Agent indices \\ \(n_{1}\), \(n_{2}\), \(n_{3}\) & The number of EOAs, CGAs and HOAs, respectively \\ \(t,\eta\) & Iteration-step and step-size \\ \(a_{i}\), \(b_{i}\), \(\alpha_{i}\), \(\beta_{i}\), \(\gamma_{i}\) and \(\xi_{i}\) & Cost coefficients \\ \(x_{i}^{pq}\), \(x_{i}^{qp}\) & Electrical power and heat power outputs of CGA \(i\) \\ \(x_{i}^{p}\) & Electrical power output of EOA \(i\) \\ \(x_{i}^{q}\) & Heat power output of HOA \(i\) \\ \(x_{\min}^{pq}(x^{qp})\), \(x_{\max}^{pq}(x^{qp})\) & Upper and lower bound functions of \(x^{pq}\) \\ \(x_{i\min}^{p}(x^{pq})\), \(x_{\max}^{qp}(x^{pq})\) & Upper and lower bound functions of \(x^{pq}\) \\ \(x_{i\min}^{p}\), \(x_{i\max}^{p}\) & Lower and upper bounds of \(x_{i}^{p}\) \\ \(d_{m}\), \(\varepsilon_{m}\), \(f_{m}\) & Coefficients of the set of linear constraints of CGA \\ \(x_{i\min}^{q}\), \(x_{i\max}^{q}\) & Lower and upper bounds of \(x_{i}^{q}\) \\ \(x_{w\min}^{pq}\), \(x_{w\max}^{pq}\) & Whole minimum and maximum values for all of \(x^{pq}\) \\ \end{tabular}  
where \(\rho=\sum\limits_{i=n_{1}+1}^{n_{1}+n_{2}+n_{3}}1/2\alpha_{i}\). In this way, we can get \(WV=I\).  
Furthermore, the change of \(\hbar\) and \(\lambda\)can be quantified by the eigenvalues of \(M\) as follows  
\[\mathrm{M}=W\Omega V=\left[\begin{array}{cc}0&0\\ \overline{\omega}^{T}\overline{\upsilon}&-\rho\omega^{T}\upsilon\end{array}\right]\]  
It is easy to obtain that \(d\hbar/d\eta=0\) and \(d\lambda/d\eta=-\rho\overline{\omega}^{T}\upsilon<0\) are the two eigenvalues of \(\mathrm{M}\). That implies \(\hbar\) does not change; meanwhile, \(\lambda\)becomes smaller. Thus, there exist an upper bound \(\widetilde{\eta}\) for \(\eta\), such that \(|\lambda|<1\) if \(0<\eta<\widetilde{\eta}\). Therefore, except for \(\hbar=1\), the others are located strictly inside a unit disk. More importantly, it can be verified that \([\mathbf{1}^{T},\mathbf{0}^{T}]^{T}\) is the eigenvector of \(\Gamma+\eta\Omega\). That implies \(y_{i}^{qp}\) and \(y_{i}^{q}\) converge to zero; meanwhile, \(\lambda_{i}^{qp}\) and \(\lambda_{i}^{q}\) converge to a common value. Thus, the optimal conditions are satisfied, In addition, when the inequality constraints are taken into consideration, the projection operators Equation (11c) does not affect the convergence which has been verify in [20; 24]. Then, we can conclude that \(\lambda_{i}^{q}\rightarrow\lambda_{q}^{*}\), \(\lambda_{i}^{qp}\rightarrow\lambda_{q}^{*}\), \(x_{i}^{q}\rightarrow x_{i}^{q*}\) and \(x_{i}^{qp}\rightarrow x_{i}^{q*}\), as \(2k\rightarrow\infty\).  
By using the similar process of the proof above, we can obtain that \(y_{i}^{p}\) and \(y_{i}^{pq}\) converge to zero; meanwhile, \(\lambda_{i}^{p}\) and \(\lambda_{i}^{pq}\) converge to another common value. Then, we can get \(\lambda_{i}^{p}\rightarrow\lambda_{p}^{*}\), \(\lambda_{i}^{pq}\rightarrow\lambda_{p}^{*}\), \(x_{i}^{p}\rightarrow x_{i}^{p*}\) and \(x_{i}^{pq}\rightarrow x_{i}^{pq*}\) as \(2k+1\rightarrow\infty\).  
The proof is over.