# Double-Consensus Based Distributed Optimal Energy Management for Multiple Energy Hubs
## Appendix A Constraints and Optimality Conditions  
1. The conversion relationship between inputs and outputs can be expressed as the following coupling matrix form, \[\left[\begin{array}{c}E_{out}^{\epsilon}\\ E_{out}^{h}\end{array}\right]\ =\ \left[\begin{array}{cc}\eta^{\epsilon \epsilon}&\eta^{\epsilon,chp}\rho\\ 0&\eta^{h,chp}\rho+\eta^{boil}(1-\rho)\end{array}\right]\left[\begin{array}{ c}E_{in}^{\epsilon}\\ E_{in}^{\delta}\end{array}\right]\] (11) where the dispatch factor \(\rho\in[0,1]\) is introduced to specify the share of natural gas by the CHP unit and boiler.
2. The total objective function of the EMP for multiple energy hubs is modeled as the following form, \[\min\ F\ =\ \sum\limits_{i\ =\ 1}^{n}\left(C_{i}(E_{i,in}^{\epsilon})+C_{i}(E_{ i,in}^{\delta})+W_{i}(\eta,E_{i,in}^{\delta})\right)\] (12) where \[\left\{\begin{array}{l}C_{i}(E_{i,in}^{\epsilon})\ =\ \sigma_{i}^{ \epsilon}\left(E_{i,in}^{\epsilon}\right)^{2}+b_{i}^{\epsilon}E_{i,in}^{ \epsilon}\\ C_{i}(E_{i,in}^{\delta})\ =\ \sigma_{i}^{\delta}\left(E_{i,in}^{\delta} \right)^{2}+b_{i}^{\delta}E_{i,in}^{\delta}\\ W_{i}(\eta,E_{i,in}^{\delta})\ =\ \omega_{i}^{\epsilon}\left(q_{i}^{ \epsilon,chp}\rho_{i}E_{i,in}^{\delta}\right)^{2}+\omega_{i}^{h}\left(\eta_{i }^{h,chp}\rho_{i}E_{i,in}^{\delta}+\eta_{i}^{boil}(1-\rho_{i})E_{i,in}^{\delta }\right)^{2}\end{array}\right.\] (13)
3. The local constraints arise from the limitations of the energy hubs' capability and the dispatch factors, which are expressed as \[E_{i,in}^{\epsilon,min}\leq E_{i,in}^{\epsilon}\leq E_{i,in}^{\epsilon,max}\] (14a) \[E_{i,in}^{\delta,min}\leq E_{i,in}^{\delta}\leq E_{i,in}^{\delta,max}\] (14b) \[0\leq\rho_{i}\leq 1\] (14c)
4. The Lagrange function for problem (3)-(6) is expressed as the following equation: \[L= F-\lambda_{\epsilon}\left[\sum\limits_{i\ =\ 1}^{n}\left(\eta_{i}^{ \epsilon\epsilon}E_{i}^{\epsilon}+\eta_{i}^{\epsilon,chp}E_{i}^{\epsilon \delta^{\prime}}\right)-\sum\limits_{i\ =\ 1}^{n}\left.t_{i}^{\epsilon}\right]-\lambda_{\delta}\left[\sum\limits_{i\ =\ 1}^{n}\left(\eta_{i}^{h,chp}E_{i}^{\epsilon \delta^{\prime}}+\eta_{i}^{boil}E_{i}^{\epsilon 1-\rho}\right)-\sum\limits_{i\ =\ 1}^{n}\left.t_{i}^{h}\right]\right.\] (14d) \[-\sum\limits_{i\ =\ 1}^{n}u_{i}^{\epsilon,min}(E_{i}^{\epsilon }-E_{i,in}^{\epsilon,min})-\sum\limits_{i\ =\ 1}^{n}u_{i}^{\epsilon,max}(E_{i,in}^{\epsilon,max}-E_{i}^{\epsilon})\] \[-\sum\limits_{i\ =\ 1}^{n}\sum\limits_{m\ =\ 1}^{4}u_{i,m}^{\delta}( \widetilde{d}_{i,m}E_{i}^{\epsilon\delta^{\prime}}+\widetilde{c}_{i,m}E_{i}^ {\delta,1-\rho}+\widetilde{f}_{i,m}),\]  
The optimal operating point is determined by the following KKT conditions:  
\[\left\{\begin{array}{l}\partial F/\partial E_{i}^{\epsilon}\ =\ \lambda_{ \epsilon}\eta_{i}^{\epsilon\epsilon}+u_{i}^{\epsilon,min}-u_{i}^{\epsilon, max}\\ \partial F/\partial E_{i}^{\epsilon\delta^{\prime}}\ =\ \lambda_{\epsilon}\eta_{i}^{ \epsilon,chp}+\lambda_{\delta}\eta_{i}^{h,chp}+\sum\limits_{m\ =\ 1}^{4}u_{i,m}^{\delta}\widetilde{d}_{i,m}\\ \\ \partial F/\partial E_{i}^{\epsilon 1-\rho}\ =\ \lambda_{\delta}\eta_{i}^{boil}+\sum \limits_{m\ =\ 1}^{4}u_{i,m}^{\delta}\widetilde{e}_{i,m}\\ \\ u_{i}^{\epsilon,min}(E_{i}^{\epsilon}-E_{i,in}^{\epsilon,min})\ =\ 0,u_{i}^{ \epsilon,max}(E_{i,in}^{\epsilon,max}-E_{i}^{\epsilon})\ =\ 0\\ \\ u_{i,m}^{\delta}(\widetilde{d}_{i,m}E_{i}^{\epsilon\delta^{\prime}}+\widetilde{c} _{i,m,i}E_{i}^{\delta 1-\rho}+\widetilde{f}_{i,m})\ =\ 0\\ \\ u_{i}^{\epsilon,min}\geq 0,u_{i}^{\epsilon,max}\geq 0,u_{i,m}^{\delta}\geq 0,\end{array}\right. \tag{14e}\]  
and Equation (4a) along with Equation (4b).