Multi Objective Optimization

Multi-Objective Optimization

You are an expert in multi-objective combinatorial optimization. This skill covers the full a-posteriori workflow — define Pareto optimality, generate the nondominated front exactly (weighted sum, epsilon-constraint with Gurobi) or approximately (NSGA-II, from scratch and via pymoo), measure approximation quality (hypervolume, IGD/IGD+), and support the final single-solution decision. Use the protocols below to pick the right method for the objective count, problem class, and budget, and to report results that survive peer review.

Initial Assessment

Establish these facts before writing any model or algorithm:

• How many objectives M, and are they genuinely conflicting? Sample feasible solutions and check pairwise objective correlation. Two strongly positively correlated objectives collapse to (nearly) one; optimizing both separately wastes the entire multi-objective apparatus.
• What does the decision maker actually need? One solution under known priorities (a-priori articulation), the whole trade-off curve to choose from afterwards (a-posteriori), or an interactive loop? This single question selects most of the method.
• Is the model an explicit MILP or a black-box evaluation? Epsilon-constraint needs a solvable scalarized model; if each scalarized MILP solves in seconds-to-minutes, an exact bi-objective front is usually affordable. Black-box or very large models push you to evolutionary methods.
• Are the objective functions integer-valued? Integer objectives let the epsilon-constraint sweep step by exactly 1 and terminate with the provably complete front. Continuous objectives need a grid density decision and yield a representation, not the full front.
• Estimate the nondominated set size. Bi-objective integer programs can have very many nondominated points (growing exponentially in the worst case — see Ehrgott (2005), Multicriteria Optimization). The cost of an exact sweep is one MILP per front point; budget accordingly.
• Objective scales and units. Cost in euros vs CO2 in tons differ by orders of magnitude. Every distance-based mechanism (crowding, hypervolume, IGD) silently breaks without normalization. Fix ideal/nadir estimates and record them.
• Hard constraints vs objectives. A "constraint" with a negotiable bound is often better treated as an objective (and vice versa). Epsilon-constraint makes this conversion explicit; confirm with the stakeholder which quantities are negotiable.
• Solver availability and license. Gurobi for the exact scalarization loops here; for license-free settings the same loops run on HiGHS/SCIP with longer runtimes.
• Evaluation budget for evolutionary methods. Population size N and generations G give N×G evaluations. The front cannot hold more points than the population; size N at 2-5× the front cardinality you intend to report.
• Comparison protocol. If two or more algorithms (or parameter settings) will be compared, fix now: instances, seeds per algorithm (10-30), identical normalization bounds, one shared hypervolume reference point, and a reference front (exact if available, else pooled best-known).
• Reproducibility. Every stochastic component takes an explicit seed (np.random.default_rng(seed)); every exact solve logs status, gap, and runtime per epsilon grid point.
• Decision-support endgame. Plan how the front turns into a decision: knee points, pseudo-weights, or a stakeholder workshop. A 200-point front without a selection protocol is not a deliverable.

Pareto Optimality and Method Selection

Definitions

For minimization of $M$ objectives over feasible set $X$:

$$ \min{x \in X} \; F(x) = \big(f1(x), \dots, f_M(x)\big). $$

Solution $x$ dominates $y$ (written $x \prec y$) iff $fi(x) \le fi(y)$ for all $i$ and $f_j(x) 0) | one MILP per weight | quick supported-front sketch; convex problems; warm-starting other methods | | Augmented epsilon-constraint | yes | complete front for M=2 with integer objectives | one MILP per front point (+1) | the default exact method for bi-objective MILPs | | Lexicographic (hierarchical) | n/a — one point | efficient point under priority order | M chained solves | priorities known and strict; also builds payoff tables | | Achievement scalarizing | yes | one efficient point per reference point | one solve | decision maker supplies aspiration levels | | NSGA-II (Deb et al. 2002) | yes | approximation, no guarantee | N×G evaluations | M = 2-3, black-box, large instances, nonlinearities | | NSGA-III (Deb & Jain 2014) / MOEA/D (Zhang & Li 2007) | yes | approximation | N×G evaluations | M ≥ 4, or when a structured spread is wanted |

Complexity notes

• Fast non-dominated sorting: $O(M N^2)$ per population of size $N$; crowding distance $O(M N \log N)$. NSGA-II's per-generation cost is dominated by these plus evaluation.
• Exact hypervolume is polynomial for $M \le 3$ ($O(n \log n)$ for $M = 2$) and #P-hard in general; use Monte Carlo approximation or moocore for $M \ge 5$.
• An epsilon-constraint sweep for $M = 2$ solves exactly $|Y_N| + 1$ MILPs (the last one proves infeasibility). For $M = 3$ the AUGMECON2 grid multiplies solves by the grid resolution of the second constrained objective.

Exact Pareto Fronts with Gurobi

Scalarization math

Weighted sum. For weights $w \in \mathbb{R}^M{> 0}$, any optimum of $\minx \summ wm f_m(x)$ is efficient. The converse fails for non-convex feasible images: unsupported points are optimal for no weight vector. Sweeping weights therefore yields only the supported front, and evenly spaced weights yield unevenly spaced points (clustered where the hull is flat).

Epsilon-constraint (Haimes, Lasdon & Wismer 1971). Keep one objective, bound the others:

$$ \minx \; f1(x) \quad \text{s.t.} \quad fj(x) \le \varepsilonj \;\; (j = 2, \dots, M), \;\; x \in X. $$

Every efficient solution is optimal for some $\varepsilon$, including unsupported ones. Two practical refinements:

1. Augmentation (Mavrotas 2009, "Effective implementation of the ε-constraint method in multi-objective mathematical programming problems" — AUGMECON): a plain epsilon-constraint optimum can be only weakly efficient. Add the constrained objectives' slacks to the objective with a small premium $\delta$, or equivalently optimize $f1 + \delta f2$ directly. With integer-valued $f_1$ and $\delta list[tuple[float, float]]:

"""Supported nondominated points of max(p1·x, p2·x) s.t. w·x bestf2: front.append((z1, z2)) bestf2 = z2 return front

p1, p2, w = [10.0, 1.0, 6.0], [1.0, 10.0, 6.0], [4.0, 4.0, 4.0] print(weightedsumfront(p1, p2, w, capacity=8.0))

Expected: [(16, 7), (7, 16)] — only the two supported points. The nondominated

point (11, 11) lies inside the convex hull and is optimal for NO weight vector.

### Augmented epsilon-constraint (complete front)

```python
import gurobipy as gp
from gurobipy import GRB

def pareto_front_epsilon(
    p1: list[float], p2: list[float], w: list[float], capacity: float
) -> list[tuple[int, int]]:
    """Complete nondominated front of max(p1·x, p2·x) s.t. w·x = 0.0, name="epsilon")

    front: list[tuple[int, int]] = []
    eps = 0.0
    while True:
        eps_con.RHS = eps
        m.optimize()
        if m.Status != GRB.OPTIMAL:          # INFEASIBLE terminates the sweep
            break
        z1, z2 = round(f1.getValue()), round(f2.getValue())
        front.append((z1, z2))
        eps = z2 + 1                         # demand strictly better f2 next
    return front

p1, p2, w = [10.0, 1.0, 6.0], [1.0, 10.0, 6.0], [4.0, 4.0, 4.0]
print(pareto_front_epsilon(p1, p2, w, capacity=8.0))
# Expected: [(16, 7), (11, 11), (7, 16)] — the complete front in 4 MILP solves,
# including the unsupported point (11, 11) that the weighted-sum sweep misses.

For larger models, keep the same persistent-model pattern: build once, update only eps_con.RHS between solves so Gurobi reuses presolve work, set a per-solve TimeLimit, and record (Status, MIPGap, Runtime) per grid point. Any point returned at GRB.TIME_LIMIT with SolCount > 0 is feasible but possibly dominated — flag it in the output table rather than silently mixing it into an "exact" front.

Gurobi's native multi-objective API (a priori, one point)

Gurobi's setObjectiveN solves blended (weighted) or hierarchical (lexicographic with allowed degradation) multi-objective models. It returns one preferred solution, not a front — use it when priorities are known, and the sweep above when the trade-off curve is the deliverable.

import gurobipy as gp
from gurobipy import GRB

def lexicographic_knapsack(
    p1: list[float], p2: list[float], w: list[float], capacity: float, reltol: float
) -> tuple[float, float]:
    """Hierarchical solve: maximize p1·x first, then p2·x, allowing reltol
    relative degradation of the first objective."""
    n = len(w)
    m = gp.Model("bi_knapsack_lex")
    m.Params.OutputFlag = 0
    x = m.addVars(n, vtype=GRB.BINARY, name="x")
    m.addConstr(gp.quicksum(w[i] * x[i] for i in range(n))  9.6) lets
# the second level pick the balanced solution. With reltol=0.0 it returns (16, 7).

The same payoff-table pattern (lexicographically optimize each objective in turn) produces the ideal point and, for $M = 2$, the exact nadir — the bounds the epsilon sweep and all normalization need.

NSGA-II: Mechanics and a Complete Implementation

NSGA-II (Deb, Pratap, Agarwal & Meyarivan 2002, "A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II") replaces a GA's scalar fitness with a two-level criterion: nondominated rank first, crowding distance second. Everything else is a standard elitist (mu+lambda) GA — see genetic-algorithms for operator depth and selection-and-replacement-strategies for the pressure analysis that carries over unchanged.

NSGA-II generation (population P, |P| = N, minimization)
  1. rank   np.ndarray:
    """Front index per row of objective matrix F (minimization); 0 = best. O(M N^2)."""
    weak = (F[:, None, :]  np.ndarray:
    """Crowding distance of each row within ONE front (minimization)."""
    n, m = F.shape
    if n  0:
            d[order[1:-1]] += (fk[2:] - fk[:-2]) / span
    return d

F = np.array([[1.0, 5.0], [2.0, 4.0], [3.0, 3.0], [2.5, 4.5], [4.0, 4.5]])
ranks = fast_nondominated_sort(F)
print(ranks, crowding_distance(F[ranks == 0]))
# Expected: ranks [0 0 0 1 2] — (2.5,4.5) is dominated by (2,4); (4,4.5) also by
# (2.5,4.5). Crowding on front 0: [inf 2.0 inf] (boundaries infinite, middle 1+1).

End-to-end NSGA-II for a bi-objective permutation flow shop

The application minimizes (makespan, total flowtime) over job permutations — the standard bi-objective extension of the permutation flow shop (see flow-shop-scheduling for the single-objective toolkit this builds on). The block is self-contained: machinery, permutation operators, main loop, and a seeded demo.

import numpy as np

def fast_nondominated_sort(F: np.ndarray) -> np.ndarray:
    """Front index per row of objective matrix F (minimization); 0 = best."""
    weak = (F[:, None, :]  np.ndarray:
    """Crowding distance of each row within one front (minimization)."""
    n, m = F.shape
    if n  0:
            d[order[1:-1]] += (fk[2:] - fk[:-2]) / span
    return d

def nsga2_survival(F: np.ndarray, n_keep: int) -> np.ndarray:
    """Indices of the n_keep rows kept by rank-then-crowding truncation."""
    ranks = fast_nondominated_sort(F)
    keep: list[int] = []
    for r in range(ranks.max() + 1):
        idx = np.where(ranks == r)[0]
        if len(keep) + idx.size  np.ndarray:
    """Binary tournament under the crowded-comparison operator; returns indices."""
    a = rng.integers(0, len(ranks), n_draws)
    b = rng.integers(0, len(ranks), n_draws)
    a_wins = (ranks[a]  crowd[b]))
    return np.where(a_wins, a, b)

def evaluate_flowshop(perm: np.ndarray, proc: np.ndarray) -> tuple[float, float]:
    """(makespan, total flowtime) of permutation perm on proc[(n_jobs, n_machines)]."""
    p = proc[perm].astype(float)
    comp = np.zeros_like(p)
    comp[0] = np.cumsum(p[0])
    for j in range(1, p.shape[0]):
        comp[j, 0] = comp[j - 1, 0] + p[j, 0]
        for k in range(1, p.shape[1]):
            comp[j, k] = max(comp[j - 1, k], comp[j, k - 1]) + p[j, k]
    return comp[-1, -1], comp[:, -1].sum()

def order_crossover(p1: np.ndarray, p2: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """OX variant: copy a slice from p1, fill remaining positions in p2's order."""
    n = p1.size
    a, b = np.sort(rng.choice(n, size=2, replace=False))
    child = np.full(n, -1)
    child[a : b + 1] = p1[a : b + 1]
    child[child  np.ndarray:
    """Swap two random positions."""
    q = perm.copy()
    i, j = rng.choice(perm.size, size=2, replace=False)
    q[i], q[j] = q[j], q[i]
    return q

def nsga2_flowshop(proc: np.ndarray, pop_size: int = 60, n_gen: int = 150,
                   p_cx: float = 0.9, p_mut: float = 0.3,
                   seed: int = 0) -> tuple[np.ndarray, np.ndarray]:
    """NSGA-II minimizing (makespan, total flowtime). Returns (front F, permutations)."""
    rng = np.random.default_rng(seed)
    n = proc.shape[0]
    pop = np.array([rng.permutation(n) for _ in range(pop_size)])
    F = np.array([evaluate_flowshop(p, proc) for p in pop])
    for _ in range(n_gen):
        ranks = fast_nondominated_sort(F)
        crowd = np.zeros(pop_size)
        for r in np.unique(ranks):
            idx = np.where(ranks == r)[0]
            crowd[idx] = crowding_distance(F[idx])
        parents = crowded_tournament(ranks, crowd, rng, 2 * pop_size)
        kids = []
        for i in range(0, 2 * pop_size, 2):
            child = (order_crossover(pop[parents[i]], pop[parents[i + 1]], rng)
                     if rng.random()  None:
        p = self.proc[x.astype(int)].astype(float)
        comp = np.zeros_like(p)
        comp[0] = np.cumsum(p[0])
        for j in range(1, p.shape[0]):
            comp[j, 0] = comp[j - 1, 0] + p[j, 0]
            for k in range(1, p.shape[1]):
                comp[j, k] = max(comp[j - 1, k], comp[j, k - 1]) + p[j, k]
        out["F"] = [comp[-1, -1], comp[:, -1].sum()]

rng = np.random.default_rng(7)
proc = rng.integers(1, 20, size=(12, 4))
algorithm = NSGA2(
    pop_size=48,
    sampling=PermutationRandomSampling(),
    crossover=OrderCrossover(),
    mutation=InversionMutation(),
    eliminate_duplicates=True,
)
res = minimize(BiObjectiveFlowShop(proc), algorithm, ("n_gen", 120),
               seed=7, verbose=False)

ideal, nadir = res.F.min(axis=0), res.F.max(axis=0)
span = np.where(nadir > ideal, nadir - ideal, 1.0)
hv = HV(ref_point=np.array([1.1, 1.1]))((res.F - ideal) / span)
print(len(res.F), round(float(hv), 3))
# Expected: a front of a few nondominated permutations and a normalized
# hypervolume around 1.0-1.2 (upper bound 1.21 for ref point (1.1, 1.1)).

pymoo notes worth knowing: ElementwiseProblem evaluates one solution per call (simple, slower) while Problem receives the whole population matrix (vectorize the evaluation there); constraints go in out["G"] with the convention `G HV(B). It rewards both convergence and spread, needs no known true front, but depends on $r$ and on normalization.

• IGD: mean distance from each point of a reference front $Z$ to its nearest approximation point. Cheap, intuitive, but not Pareto-compliant — a dominated front can score better. Prefer IGD+ (Ishibuchi, Masuda, Tanigaki & Nojima 2015), which measures only the dominated-direction component and is weakly Pareto-compliant.
• Report HV (primary) plus IGD+ (when a reference front exists), with the normalization bounds and reference point stated explicitly. Standard protocol: normalize all fronts by the ideal/nadir of the pooled reference front, then use $r = (1.1, \dots, 1.1)$ (Ishibuchi et al. 2018 recommend a reference p

…

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: hajibabaie
• Source: hajibabaie/combinatorial-optimization-skills
• License: MIT

Install and usage instructions live in the source repository linked above.