Multi-Objective Mixed-Integer LP

• MOMILP is an optimization framework that minimizes several linear objectives under mixed integer-continuous constraints, yielding a Pareto front of nondominated solutions.
• Exact algorithms like facet-driven dichotomic search and objective space methods rigorously compute supported extreme points for clear trade-off identification.
• Scalarization techniques, combined with outer approximation and dual methods, offer practical approaches to tackle complex engineering and operations research problems.

A multi-objective mixed-integer linear program (MOMILP) is an optimization problem in which multiple linear objectives are optimized simultaneously over a feasible region defined by linear constraints, with a subset of variables required to be integer-valued and the rest continuous. The field encompasses a rigorous theoretical framework, exact and approximate algorithms for Pareto front computation, specialized scalarization methods, outer-approximation and duality theory, and a growing range of applications in engineering and operations research.

1. Mathematical Formulation and Pareto Optimality

Given $k \geq 2$ objectives, a typical MOMILP is expressed as: $$\min_{x \in X} \; f(x) = (f_1(x), \dots, f_k(x)), \quad X = \{ x \in \mathbb{Z}^n \times \mathbb{R}^p : A x \leq b \}$$ where each $f_i(x) = c_i^\top x$ is linear, $A \in \mathbb{R}^{m \times (n+p)}$, $b \in \mathbb{R}^m$. The feasible set $X$ imposes mixed-integer and continuous decision variables. The outcome set in objective space is $Z = f(X) \subset \mathbb{R}^k$.

Dominance is defined via component-wise order: $y \leq y'$ iff $y_i \leq y'_i$ for all $i$, and $y < y'$ means strictly better in at least one component. A point $y \in Z$ is nondominated if no $y' \in Z$ satisfies $y' \leq y$ and $y' \neq y$; the set of all such outcomes is the Pareto front $Z_N$.

Supported nondominated points lie on the boundary of $\operatorname{conv}{Z}$, the convex hull, and are obtained by minimizing strictly positive weighted sums, i.e., for $\lambda \in \mathbb{R}^k_{>0}$: $\min_{x \in X} \lambda^\top f(x)$ Nondominated extreme points are the supported nondominated solutions which are vertices of $\operatorname{conv}{Z}$ (Przybylski et al., 2019).

2. Exact Algorithms for MOMILP Pareto Fronts

Facet-Driven Dichotomic Search

A key approach for enumerating all nondominated extreme points is the dichotomic search, generalized to arbitrary $k$ objectives (Przybylski et al., 2019). The algorithm iteratively identifies new extreme points by:

1. Maintaining a set $S$ of current extreme points and computing the convex hull polytope $P = \operatorname{conv}{S}$.
2. For each unexplored facet $F$ of $P$ with inward normal $\lambda$, the weighted-sum problem $\min_{x \in X} \lambda^\top f(x)$ is solved.
3. If a new extreme point $y^*$ is found (i.e., $\lambda^\top y^* < \lambda^\top y'$ for $y' \in F$), it is added to $S$ and the hull is updated; otherwise, the facet is validated.
4. The process iterates until all facets are validated.

Complexity is incremental polynomial in $|Z_{NE}|$ for fixed $k$, provided the lexicographic MILP subroutine is polynomial time (Przybylski et al., 2019).

Objective Space Algorithms

Pettersson & Özlen introduce a three-phase objective-space algorithm to enumerate the complete Pareto front, including unsupported points (Pettersson et al., 2019):

1. Phase 1: Slice the feasible set by fixing integer variables, compute convex hulls via Benson’s algorithm in each slice, and assemble a covering collection $\mathcal{C}$ of simplices.
2. Phase 2: Eliminate polytope overlaps in $\mathcal{C}$ by set subtraction and convex decomposition, yielding a collection with disjoint interiors.
3. Phase 3: For each pair of polytopes, check and excise dominated regions by projecting dominance regions, leaving an exact polyhedral decomposition of the Pareto front.

Correctness is guaranteed, with covered complexity exponential in the number of objectives and integer variables (Pettersson et al., 2019).

Outer Approximation and Dual Benson

Outer-approximation algorithms construct the upper image (Edgeworth–Pareto hull) as a polyhedron, incrementally refining an outer superset via supporting hyperplanes obtained through calls to a weighted-sum MILP "oracle" (Bökler et al., 2021). The number of required oracle calls is polynomial in the number of facets for fixed $k$, and the procedure runs in incremental-polynomial delay if the subproblems can be solved in polynomial time.

PaMILO adapts the Dual-Benson algorithm: it builds and refines an outer approximation in dual space by shooting rays (i.e., weighted-sum directions) and updating with supporting inequalities, thereby discovering all supported extreme points even in mixed-integer and (convex) quadratic settings (Bökler et al., 2022).

3. Scalarization and Pareto Front Approximation

Several scalarization strategies are effective for MOMILPs:

• Weighted-sum scalarization: Recovers supported Pareto points. In non-convex/disconnected feasible regions, this can miss unsupported points (Burachik et al., 2019).
• $\epsilon$-constraint method: Fixes all but one objective via constraints, optimizing the remaining one. Ensures weak efficiency but may require careful refinement to obtain only nondominated points. Efficient for problems where the secondary objectives are discrete-valued (e.g., cardinality in bi-objective subset selection (Charkhgard et al., 2018)).
• Weighted-constraint scalarization: For each objective $k$, minimize $f_k(x)$ subject to $w_i f_i(x) \leq w_k f_k(x)$ for $i \neq k$. This surjects onto the set of weakly efficient points and can recover Pareto points even in disconnected fronts (Burachik et al., 2019).
• Pascoletti–Serafini scalarization: Minimizes the maximum excess over a utopia point under a weight direction, practical mainly for convex or connected fronts (Burachik et al., 2019).

For large-scale problems or applications requiring representative subsets of the Pareto front, advanced $\epsilon$-constraint strategies provide coverage-oriented, uniformity-oriented, or cardinality-controlled representations with acceleration mechanisms such as redundancy checking and LP relaxation early exit (Mesquita-Cunha et al., 2021).

4. Relaxations, Duality, and Bound Sets

Theoretical foundations encompass a hierarchy of relaxations (Dunbar et al., 2023):

• Continuous Relaxation (MOLP): Dropping integer constraints yields an upper bound on the nondominated set.
• Convex Hull Relaxation: Tight at supported points but can be slack for unsupported Pareto points.
• Lagrangian Relaxation: By dualizing complicating constraints, one can obtain tighter bounds for unsupported points. The corresponding Lagrangian dual set offers weak duality in general and strong duality at supported points under additional convexity conditions.

Superadditive duality theory extends from the single-objective value function to the multiobjective case, with set-valued and vector-valued superadditive duals offering structural insights into nondominated outcomes and relaxation tightness (Dunbar et al., 2023).

5. Practical Implementations and Applications

Implementation relies on high-performance MILP solvers (CPLEX, Gurobi, SCIP) and convex hull libraries (CGAL, Qhull). Static or incremental convex hull algorithms may be used, with practical limits as $k$ grows (typically tractable for $k \leq 5$) (Przybylski et al., 2019). Data structures must efficiently maintain the sets of candidate extreme points, hull facets, and facilitate rapid duplicate detection (Przybylski et al., 2019).

Empirical studies demonstrate scalability to hundreds of variables in classical assignment and knapsack MOMILPs (Przybylski et al., 2019), and industry-relevant models include supply chain planning for 3D printing (Talebi, 2024), feature subset selection (Charkhgard et al., 2018), and engineering design (Burachik et al., 2019). For large-scale or time-limited scenarios, interval Pareto front representations via Chebyshev scalarization with incumbent and relaxation-based bounds can guide decision support (Filcek et al., 2023).

6. Heuristics and Approximation Schemes

Approximation methods such as LP relaxation–based matheuristics provide fast, high-quality trade-off sets for tri-objective MOIPs by combining outer-approximation lower bound sets (via Bensolve) with rounding and path-relinking in the integer domain (An et al., 2021). Such schemes can be extended to mixed-integer settings and enriched with local-search or feasibility-pump heuristics.

Approximate dichotomic schemes for extreme point generation can be terminated early to yield $\epsilon$-approximate convex hulls (Przybylski et al., 2019). Selective Pareto front representations for multi-criteria decision support can be tailored by coverage, uniformity, or cardinality goals using modified $\epsilon$-constraint loops (Mesquita-Cunha et al., 2021).

7. Complexity, Limitations, and Future Directions

The computational complexity of exact Pareto front enumeration in MOMILPs grows exponentially with the number of objectives and integer variables, a reflection of the combinatorial nature of the feasible set and the potential for non-convexity in the Pareto front (Przybylski et al., 2019, Pettersson et al., 2019). Output-sensitivity and incremental polynomial-time behavior can be achieved for supported points and fixed $k$ (Przybylski et al., 2019, Bökler et al., 2021).

Outstanding challenges include handling higher-dimensional ($k > 5$) objective spaces, scaling to industrial problem sizes, and developing robust numerical techniques for convex hull and relaxation computations in high dimensions. The integration of branch-and-bound or branch-and-cut methods, as well as hybrid schemes combining dual-Benson enumeration with facet-driven approaches, offer promising directions for improved scalability and flexibility (Przybylski et al., 2019). Enhanced approximation frameworks, interval-bound navigation on the Pareto front, and robust relaxations for mixed-integer convex settings represent expanding frontiers in the MOMILP research landscape.