## 2 Materials and Methods
After the data are normalized, the gray correlation coefficient (\(\xi_{i}\)), which quantifies the relationship between the ideal and actual normalized results, is computed using Equation (9) [48; 49; 50; 51]. In this equation, the resolution coefficient \(\mathcal{O}\) was set to 0.5, following the standard practice in GRA applications. This value, originally proposed by [PERSON] [57], ensures a balanced level of sensitivity and has been widely adopted in similar engineering optimization studies [48; 49; 50]. \[\xi_{i}(\text{k})=\frac{\text{\Delta}_{\text{min}}+\text{\phi}\text{ \Delta}_{\text{max}}}{\text{\Delta}_{0i}(\text{k})+\text{\phi}\text{ \Delta}_{\text{max}}} \tag{9}\]  
\[\text{\Delta}_{0i}(\text{k})=\|\text{y}_{0}(\text{k})-\text{y}_{i}(\text{k})  
In Equation (10), \(\Delta_{00}\)(k) represents the deviation between \(y_{0}\)(k) and \(y_{i}\)(k). The values \(\Delta\)max and \(\Delta\)min are determined using Equations (11) and (12), respectively. The gray relational degree (\(\gamma_{i}\)) is calculated using Equation (13) with the normalized weighting factor.  
\[\gamma_{i}=\frac{1}{n}\ \sum_{k=1}^{n}w_{k}\xi_{i}(k) \tag{13}\]  
In this equation, a high \(\gamma_{i}\) indicates a strong correlation between \(y_{0}\)(k) and \(y_{i}\)(k). If the two compared series have identical values, then \(\gamma_{i}=1\). This implies that \(\gamma_{i}\) is used to measure the similarity between the compared series and the reference series. Additionally, in Equation (13), the total weighting factor satisfies the condition \(w_{k}=1\), where \(w_{k}\) represents the normalized weighting factor for each response. The weighting factor defines the influence of each output parameter on the objective function, allowing the transformation of the multi-objective optimization problem into a single equation. It is calculated using Equation (14) [48; 49; 50; 51], where t represents the number of responses, p is the number of parameters, and Delta denotes the S/N range. This weighting strategy was chosen to reflect the relative sensitivity of each output to the control factors based on empirical data, thus ensuring objectivity. While physical importance-based weights are also valid, the Delta-based approach avoids subjective bias and maintains consistency between the Taguchi and GRA stages. \[w_{i}=\frac{\sum_{j=1}^{p}Delta_{ij}}{\sum_{i=1}^{t}\sum_{j=1}^{p}Delta_{ij}} \tag{14}\]  
The Taguchi method is an efficient statistical tool that reduces the number of simulations by employing orthogonal arrays, enabling systematic evaluation of factor effects across multiple levels. However, it does not inherently support multi-objective optimization. To address this, GRA was employed as a decision-making method to consolidate multiple performance indicators into a single response metric. This combined Taguchi-GRA framework is known for its simplicity, low computational cost, and applicability in engineering systems involving multi-response problems. Nevertheless, GRA has certain limitations, such as sensitivity to weighting factors and normalization schemes, which are acknowledged in the literature. This integrated Taguchi-GRA optimization framework provides a novel and practical approach for enhancing the performance of CNG refueling stations. It enables a systematic evaluation of multi-parameter and multi-objective trade-offs, such as maximizing vehicle fuel delivery while minimizing compressor energy consumption. By transforming multiple performance outputs into a single comparable metric, this approach simplifies complex decision-making processes and allows for effective optimization under constrained conditions. The combined methodology is computationally efficient and highly interpretable, making it well suited for engineering systems like CNG refueling stations, where both operational efficiency and system level design considerations must be balanced.