## 2 Matrix-Based Model Formulation for MES
### District Gas Network
#### 2.3.1 Matrix-Based Universal DGN Model
Once the compressor is driven by the gas turbine, the consumed gas flow \(\tau_{\rm C}\) that can be extracted from any gas node is given by  
\[\mbox{$\tau_{\rm C}$}=\mbox{$\tau_{\rm C,st}$}\mbox{$\cdot$}\mbox{$\eta_{\rm C }$}^{-1}=\left(\mbox{$\alpha_{\rm C}$}+\mbox{$\beta_{\rm C}$}\mbox{$\cdot$} \mbox{$H_{\rm C}$}+\mbox{$\gamma_{\rm C}$}\mbox{$\cdot$}\mbox{$H_{\rm C}$}^{2 }\right)\mbox{$\cdot$}\mbox{$\eta_{\rm C}$}^{-1} \tag{21}\]  
where \(\tau_{\rm C,st}\) is the equivalent energy consumption vector of the extracted gas under standard gross calorific value (GCV); \(q_{\rm C}\) is the GCV vector at the extraction node; and \(\mbox{$\alpha_{\rm C}$}\), \(\beta_{\rm C}\), and \(\gamma_{\rm C}\) are the efficiency coefficient vectors. Since the sum of the gas volume flow at each node is zero, the nodal volume flow balance equation of the DGN is expressed as  
\[-A_{\rm g}\mbox{$f_{\rm P}$}+B_{\rm g}\mbox{$f_{\rm E}$}+\mbox{$f_{\rm L}$}+ \mbox{$T_{\rm C}$}\mbox{$\tau_{\rm C}$}=0 \tag{22}\]  
where \(B_{\rm g}=[B_{\rm g}^{\rm S}\,B_{\rm g}^{\rm C}]\) is the node-equipment incidence matrix; \(B_{\rm g}^{\rm S}\) and \(B_{\rm g}^{\rm C}\) are the node-slack source incidence matrix and the node-compressor incidence matrix, respectively; \(\mbox{$f_{\rm E}$}=[\mbox{$f_{\rm C}$0};\mbox{$f_{\rm C}$1}]\) is the equipment gas flow vector; \(\mbox{$f_{\rm C}$2}\) is the slack source flow vector; \(\mbox{$f_{\rm C}$3}=\mbox{$f_{\rm C}$4}-\mbox{$f_{\rm C}$5}\) is the net load flow vector; \(\mbox{$f_{\rm C}$6}=\mbox{$f_{\rm C}$7}\mbox{$f_{\rm C}$8}\mbox{$f_{\rm C}$9}^{-1}\) is the load volume demand vector; \(\mbox{$\theta_{\rm C}$0}\) and \(q\) are the load energy demand vector and the nodal GCV vector, respectively; \(\mbox{$\theta_{\rm C}$1}\) is the non-slack source flow vector; and \(T_{\rm C}\) is the node-extraction incidence matrix. In this paper, five compressor modes are considered: (1) Mode I: the squared boost ratio \(R_{\rm C}^{0}\) is specified; (2) Mode II: the boost squared difference \(D_{\rm C}^{0}\) is specified; (3) Mode III: the compressor flow rate \(f_{\rm C}^{0}\) is specified; (4) Mode IV: the squared inlet pressure \(\Pi_{\rm I}^{0}\) is specified; (5) Mode V: the squared outlet pressure \(\Pi_{\rm O}^{0}\) is specified. To accommodate the diverse operating conditions of the compressor, the equipment control equation is formulated. For the \(k\)-th equipment, its control expression can be written as  
$\theta_{\rm C}$2  
where \(k/s/i/o\) is the index of the equipment/slack source node/compressor inlet node/compressor outlet node; \(w_{\rm N,k_{\rm S}}/w_{\rm N,k_{\rm I}}/w_{\rm N,k_{\rm O}}/w_{\rm E,k}/d_{k}\) is the control parameter of slack source node pressure/compressor inlet node pressure/compressor outlet node pressure/equipment flow/compressor, whose value for different equipment is shown in Table 1; \(\Pi_{s}/\Pi_{i}/\Pi_{o}\) is the squared pressure of slack source node \(s/\)compressor inlet node \(i\)/compressor outlet node \(o\); \(f_{\text{E,k}}\) is the gas flow of the \(k\)-th equipment; and \(\mathbb{G}^{\text{equ}}\) is the equipment set of the DGN. Equation (23) can be rewritten into matrix form as  
$\theta_{\rm C}$3  
where the nonzero element in the \(k\)-th row and \(s\)-th/\(i\)-th/\(o\)-th column of \(\mathbf{W}_{\text{N}}\) is \(w_{\text{N,k}s}/w_{\text{N,k}}/w_{\text{N,k}o}\); \(\mathbf{W}_{\text{E}}\) is the diagonal matrix of \(w_{\text{E,k}}\); and \(\mathbf{d}\) is the vector of \(d_{k}\). According to the aforementioned assumptions of the steady-state gas flow, the injection of heterogeneous gas sources will significantly alter the gas composition within the DGN, primarily affecting gas density and GCV. Hence, referring to the existing DGN models [27; 31], the deviation equations that follow the nodal mass and energy conservation laws are constructed to track the gas quality across the network:  
\[\left\{