## 4 Solving Method
### The Augmented \(\epsilon\)-Constraint Method  
Given the limitations of the weighted sum method, especially in scenarios where objective importance is unclear or in a posteriori decision making, it becomes necessary to explore alternative approaches. As an a posteriori method, the \(\epsilon\)-constraint method has been extensively employed in the literature for approximating the Pareto frontier in multi-objective problems [52]. In this approach, one objective is selected as the primary goal, with the remaining objectives reformulated as constraints. However, identifying extreme value points typically requires sorting, expansion, and compression steps, which are time-consuming and make it difficult to estimate the maximum attainable value of the objective to be maximized. To address the bi-objective MIP model, we utilize the augmented \(\epsilon\)-constraint approach, which enhances the standard method by adding a penalty term to the objective function. This innovation brings multiple key advancements over classical \(\epsilon\)-constraint techniques. In contrast to the basic \(\epsilon\)-constraint method, which often generates weakly efficient solutions that allow one objective function to be improved while keeping the other objective functions unaffected, the augmented \(\epsilon\)-constraint method effectively eliminates such weak solutions by integrating both objectives into a single penalized objective function. Specifically, non-Pareto-optimal solutions require manual removal through pairwise comparisons. The augmented approach integrates both objectives into a single optimized function, inherently filtering out non-Pareto-optimal candidates during computation and eliminating the need for post-processing steps. Additionally, the penalty term's design allows dynamic adjustment of the tolerance level \(\epsilon\), enabling more flexible and precise control over solution density. For instance, when objective functions take integer values, setting \(\epsilon<1\) ensures the augmented method captures the full set of non-dominated objective vectors, making it particularly suitable for complex scenarios like transportation route optimization where fine-grained management of trade-offs between objectives (such as time and risk) is critical [57]. This improvement not only helps to better manage infeasible solutions but also accelerates convergence toward the Pareto frontier. Moreover, compared with traditional \(\epsilon\)-constraint variants (such as the revised or improved revised methods) that require solving two separate single-objective models for each iteration, the augmented method achieves equivalent results with only one model, significantly reducing redundant calculations and boosting computational efficiency. Overall, this approach effectively  
Figure 4: Graphical comparison of the Pareto frontier for an example: (**a**) weighted sum method with evenly distributed weight combinations. (**b**) \(\epsilon\)-constraint method with evenly distributed \(\epsilon\). balances solution quality and computational efficiency while maintaining the theoretical rigor of Pareto optimality. The augmented \(\epsilon\)-constraint method is elaborated as follows. Based on the bi-objective optimization model (2), let \(f_{1}^{\max}:=\underset{g\in G}{\max}f_{1}(g)\) and \(f_{2}^{\max}:=\underset{g\in G}{\max}f_{2}(g)\). Let \(f_{1}^{\min}:=\underset{g\in G}{\max}f_{1}(g)\) subject to \(f_{2}(g)\geq f_{2}^{\max}\), \(g\in G\). Let \(\epsilon\) be a small positive number such that \(f_{1}^{\max}-f_{1}^{\min}\) is an integer multiple of \(\epsilon\). Define a multiset \(\Omega:=\emptyset\), where \(\Omega\) represents a multiset composed of obtained solutions by the augmented \(\epsilon\)-constraint method. Define \(K\) as the parameter that determines the maximum number of iterations in the subsequent iterative process. It essentially quantifies how many discrete steps can be taken based on the chosen \(\epsilon\). \(K:=\big{(}f_{1}^{\max}-f_{1}^{\min}\big{)}/\epsilon\).