## 3 Methodology
### Feature Importance and Interpretability Using SHAP  
While tree-based ensemble learning algorithms such as GBR offer superior predictive accuracy in high-dimensional and nonlinear domains, their internal decision structures are inherently complex and often criticized for lacking transparency. This opacity can hinder the interpretability of the model's behavior, particularly in applied engineering domains where traceable decision logic is essential for adoption [40]. To address this issue and improve the explainability of the developed predictive framework, the SHapley Additive explanations (SHAP) methodology was implemented as a post hoc interpretability tool. SHAP is grounded in cooperative game theory and quantifies the marginal contribution of each input feature to the model's output by evaluating the change in prediction when a feature is included or excluded across all possible permutations of feature subsets [41]. This is mathematically formalized by the Shapley value formulation given in Equation (15) [42]. \[\phi_{i}(f)=\sum_{S\subseteq\left\{x_{i},\ldots,x_{p}\right\}\smallsetminus \left\{y_{i}\right\}}\frac{|S|!(p-|S|-1)!}{p!}\left[f(S\cup\{x_{i}\})-f(S)\right] \tag{15}\]  
where x\({}_{i}\) is the feature variable; p is the number of the features; S denotes a subset of the features; and f(x\({}_{i}\)) is the output of the model. In SHAP, a feature's importance is determined by altering its value and then quantifying the prediction error. The SHAP analysis was carried out at two levels: global and local. Global SHAP values were aggregated over the entire dataset to assess the average importance of each feature across all prediction instances. This allowed the identification of dominant input variables that systematically influenced the predictions of total cost and CO\({}_{2}\) emissions. For example, features such as transport distance and energy intensity were consistently observed to have the highest SHAP values in both targets, confirming their central role in the logistics-emissions-cost relationship [43; 44]. This global interpretability analysis not only validated the relevance of the selected input features but also informed the prioritization of decision variables in the subsequent optimization phase. In parallel, local SHAP values were calculated for individual samples to understand how specific feature values contributed to the model's prediction in particular scenarios. These instance-level explanations are especially valuable for scenario-specific decision support, where engineers or policymakers may need to justify why certain configurations result in higher emissions or costs [45]. For example, in cases where grid CO\({}_{2}\) intensity was low but the predicted emission was still high, SHAP force plots revealed that elevated moisture content and long transport distances outweighed the benefit of clean electricity, thus offering actionable insight for process adjustment. SHAP is also used to describe the working process of a trained model by using an additional feature attribution method that treats the output model as an interpretable model of the linear addition of the input variables. The prediction f(x) of the original model can be represented, as shown in Equation (16) [42]. \[\text{f(x)}=\text{h(x_{\text{s}})}=\varphi_{0}+\sum_{i=1}^{\text{k}}\varphi_{i} \text{x}_{i} \tag{16}\]  
where k is the number of input features and \(\varphi_{0}\) denotes the constant value when no inputs are used. Inputs x and xs are related by a mapping function, x = mx(xs). To further explore the interaction effects between variables, SHAP dependence plots were generated. These plots highlighted key nonlinear relationships, such as the amplification of energy intensity effects in regions with high grid emission factors, or the compounded cost impact of combined long transport distance and high gradation adjustment levels [46].