Robust Optimization

Robust Optimization#

Real-world optimization problems rarely have perfectly known parameters. Demand forecasts carry error, material costs fluctuate, and physical measurements are imprecise. Robust optimization provides a principled framework for making decisions that remain feasible — and near-optimal — for every realization of uncertain parameters within a prescribed set [Ben-Tal et al., 2009, Ben-Tal and Nemirovski, 1999].

What is Robust Optimization?#

Given a nominal problem

\[ \min_x\; f(x, p) \quad \text{s.t.}\; g_i(x, p) \le 0,\; i=1,\ldots,m \]

where $p$ is an uncertain parameter taking values in an uncertainty set $\mathcal{U}$, the robust counterpart is:

\[ \min_x\; \max_{p \in \mathcal{U}}\; f(x, p) \quad \text{s.t.}\; g_i(x, p) \le 0\; \forall p \in \mathcal{U},\; i=1,\ldots,m \]

A key insight (Ben-Tal & Nemirovski Ben-Tal and Nemirovski [1999]) is that for standard uncertainty set shapes, this minimax problem can be reformulated as a deterministic program without nested optimization.

Uncertainty Sets#

discopt supports three standard uncertainty set families:

Set	Definition	Reformulation
Box	$\lvert p_j - \bar{p}_j \rvert \le \delta_j$	Substitute worst-case bound per component
Ellipsoidal	$\lVert \Sigma^{-1/2}(p-\bar{p}) \rVert_2 \le \rho$	Add SOC norm penalty [Ben-Tal and Nemirovski, 1999]
Polyhedral	$A\xi \le b,\; \xi = p - \bar{p}$	LP dual auxiliary variables [Bertsimas and Sim, 2004]

import discopt.modeling as dm
import numpy as np
from discopt.ro import BoxUncertaintySet, RobustCounterpart

Example 1 — Production Planning with Box Uncertainty#

A manufacturer produces three products with unit costs $c$ and must meet a total demand $d$:

\[ \min_x\; c^\top x \quad \text{s.t.}\; \mathbf{1}^\top x \ge d,\; 0 \le x \le 200 \]

Both costs $c$ and demand $d$ are uncertain: $c \in [\bar{c} - \delta_c, \bar{c} + \delta_c]$ and $d \in [\bar{d} - \delta_d, \bar{d} + \delta_d]$. We solve the nominal problem first, then compare with the robust counterpart.

c_bar = np.array([10.0, 15.0, 8.0])  # nominal unit costs
d_bar = 100.0  # nominal demand

# --- Nominal model (solve with nominal parameter values) ---
m_nom = dm.Model("nominal")
x_nom = m_nom.continuous("x", shape=(3,), lb=0, ub=200)
c_nom = m_nom.parameter("c", value=c_bar)
d_nom = m_nom.parameter("d", value=d_bar)
m_nom.minimize(dm.sum(c_nom * x_nom))
m_nom.subject_to(dm.sum(x_nom) >= d_nom, name="demand")

r_nom = m_nom.solve()
x_n = list(r_nom.x.values())[0]
print(f"Nominal:  obj = {r_nom.objective:.2f},  x = {np.round(x_n, 1)}")
print("  (All demand met by cheapest product: 8.0 * 100 = 800)")

# --- Robust model (costs +/- 10%, demand +/- 5%) ---
m_rob = dm.Model("robust")
x_rob = m_rob.continuous("x", shape=(3,), lb=0, ub=200)
c_rob = m_rob.parameter("c", value=c_bar)
d_rob = m_rob.parameter("d", value=d_bar)
m_rob.minimize(dm.sum(c_rob * x_rob))
m_rob.subject_to(dm.sum(x_rob) >= d_rob, name="demand")

unc_c = BoxUncertaintySet(c_rob, delta=0.10 * c_bar)
unc_d = BoxUncertaintySet(d_rob, delta=0.05 * d_bar)
rc = RobustCounterpart(m_rob, [unc_c, unc_d])
rc.formulate()

r_rob = m_rob.solve()
x_r = list(r_rob.x.values())[0]
print(f"\nRobust:   obj = {r_rob.objective:.2f},  x = {np.round(x_r, 1)}")
print("  (Worst-case cost 8.8, worst-case demand 105: 8.8 * 105 = 924)")
print(
    f"\nPrice of robustness: {r_rob.objective - r_nom.objective:.2f} "
    f"({(r_rob.objective / r_nom.objective - 1) * 100:.1f}% increase)"
)

Nominal:  obj = 800.00,  x = [  0.   0. 100.]
  (All demand met by cheapest product: 8.0 * 100 = 800)

Robust:   obj = 924.00,  x = [  0.   0. 105.]
  (Worst-case cost 8.8, worst-case demand 105: 8.8 * 105 = 924)

Price of robustness: 124.00 (15.5% increase)

Example 2 — Robust Demand with Capacity Constraints#

When capacity constraints force diversification, the robust counterpart changes the production plan. Here product $x$ is cheap but capacity-limited, and product $y$ is the expensive backup:

\[ \min_{x,y}\; 3x + 5y \quad \text{s.t.}\; x + y \ge d,\; x \le 40,\; 0 \le x,y \le 100 \]

With demand $d$ uncertain within $[40, 60]$ (nominal $\bar{d} = 50$, $\delta = 10$), the robust counterpart must ensure feasibility at $d = 60$.

# --- Nominal ---
m_nom = dm.Model("nominal")
x_n = m_nom.continuous("x", lb=0, ub=100)
y_n = m_nom.continuous("y", lb=0, ub=100)
d_n = m_nom.parameter("d", value=50.0)
m_nom.minimize(3.0 * x_n + 5.0 * y_n)
m_nom.subject_to(x_n + y_n >= d_n, name="demand")
m_nom.subject_to(x_n <= 40.0, name="x_cap")

r_nom = m_nom.solve()
print(f"Nominal (d=50):  obj = {r_nom.objective:.2f},  x = {r_nom.x}")

# --- Robust (demand uncertain +/- 10) ---
m_rob = dm.Model("robust")
x_r = m_rob.continuous("x", lb=0, ub=100)
y_r = m_rob.continuous("y", lb=0, ub=100)
d_r = m_rob.parameter("d", value=50.0)
m_rob.minimize(3.0 * x_r + 5.0 * y_r)
m_rob.subject_to(x_r + y_r >= d_r, name="demand")
m_rob.subject_to(x_r <= 40.0, name="x_cap")
RobustCounterpart(m_rob, BoxUncertaintySet(d_r, delta=10.0)).formulate()

r_rob = m_rob.solve()
print(f"Robust  (d<=60):  obj = {r_rob.objective:.2f},  x = {r_rob.x}")

print("\nThe robust solution orders 10 more units of the expensive product")
print("to guarantee feasibility if demand reaches 60.")

Nominal (d=50):  obj = 170.00,  x = {'x': array(40.), 'y': array(10.)}
Robust  (d<=60):  obj = 220.00,  x = {'x': array(40.), 'y': array(20.)}

The robust solution orders 10 more units of the expensive product
to guarantee feasibility if demand reaches 60.

Example 3 — Portfolio with Ellipsoidal Uncertainty#

Portfolio optimization with uncertain expected returns [Ben-Tal and Nemirovski, 1999]:

\[ \max_x\; \mu^\top x \quad \text{s.t.}\; \mathbf{1}^\top x = 1,\; x \ge 0 \]

where returns $\mu$ lie in an ellipsoidal set $\{\mu: \lVert \mu - \bar{\mu} \rVert_2 \le \rho\}$. The robust counterpart adds a norm penalty: $\max_x\; \bar{\mu}^\top x - \rho \lVert x \rVert_2$.

Robust portfolios are usually introduced with “the robust investor diversifies”, but that is only true past a threshold, and it is worth seeing where the threshold is. Concentrating everything in the best asset gives $\lVert x \rVert_2 = 1$; spreading equally over $k$ assets gives $\lVert x \rVert_2 = 1/\sqrt{k}$. So diversifying buys you at most $\rho\,(1 - 1/\sqrt{k})$ of penalty relief, and it costs you the return you give up by moving out of the best asset. Concentration stays optimal until $\rho$ is large enough to outweigh the spread in $\bar\mu$.

The sweep below shows exactly that transition rather than asserting it.

from discopt.ro import EllipsoidalUncertaintySet

rng = np.random.RandomState(42)
n_assets = 5
mu_bar = rng.uniform(0.05, 0.15, n_assets)

mu_sorted = np.sort(mu_bar)[::-1]
print(f"Expected returns: {np.round(mu_bar, 4)}")
print(f"Best asset {mu_sorted[0]:.4f}, runner-up {mu_sorted[1]:.4f}, "
      f"gap {mu_sorted[0] - mu_sorted[1]:.4f}")
print()


def solve_portfolio(rho=None):
    """Nominal portfolio when rho is None, else the robust counterpart."""
    m = dm.Model("portfolio")
    x = m.continuous("x", shape=(n_assets,), lb=0, ub=1)
    mu = m.parameter("mu", value=mu_bar)
    m.minimize(-(mu @ x))
    m.subject_to(dm.sum(x) == 1.0, name="budget")
    if rho is not None:
        RobustCounterpart(m, EllipsoidalUncertaintySet(mu, rho=rho)).formulate()
    r = m.solve()
    return -r.objective, np.asarray(list(r.x.values())[0])


ret_nom, x_nom = solve_portfolio()
n_held = int((x_nom > 1e-4).sum())
print(f"Nominal:  return = {ret_nom:.4f}  x = {np.round(x_nom, 3)}  "
      f"({n_held} asset{'s' if n_held > 1 else ''} held)")
print()
print(f"{'rho':>6}  {'return':>8}  {'assets':>6}  weights")
for rho in (0.02, 0.05, 0.10, 0.20, 0.50):
    ret, x = solve_portfolio(rho)
    held = int((x > 1e-4).sum())
    print(f"{rho:>6.2f}  {ret:>8.4f}  {held:>6d}  {np.round(x, 3)}")

Expected returns: [0.0875 0.1451 0.1232 0.1099 0.0656]
Best asset 0.1451, runner-up 0.1232, gap 0.0219

Nominal:  return = 0.1451  x = [0. 1. 0. 0. 0.]  (1 asset held)

   rho    return  assets  weights
  0.02    0.1251       1  [0. 1. 0. 0. 0.]
  0.05    0.1011       3  [0.    0.587 0.295 0.117 0.   ]
  0.10    0.0710       4  [0.091 0.408 0.287 0.214 0.   ]
  0.20    0.0212       5  [0.156 0.291 0.24  0.209 0.104]
  0.50   -0.1157       5  [0.183 0.235 0.215 0.203 0.163]

Reading the sweep. At $\rho = 0.02$ the robust portfolio is still fully concentrated — identical to the nominal one — because the uncertainty radius is smaller than the 0.0219 gap between the best asset and the runner-up. All robustness buys you there is a flat $\rho \lVert x \rVert_2 = 0.02$ haircut on the guaranteed return.

Diversification appears at $\rho = 0.05$ (three assets), broadens to all five by $\rho = 0.20$, and tends toward equal weights as $\rho$ grows — in the limit the penalty term dominates the return term entirely and the optimizer simply minimizes $\lVert x \rVert_2$.

Note also that the guaranteed return goes negative for large $\rho$. That is not a solver failure: if returns really could be off by 0.5 in 2-norm, no budget-feasible portfolio can guarantee a profit, and the model is correctly reporting that the worst case is a loss.

Multiple Uncertain Parameters#

Multiple parameters from the same uncertainty family can be passed as a list:

m_multi = dm.Model("multi_uncertain")
x = m_multi.continuous("x", shape=(3,), lb=0, ub=100)
c = m_multi.parameter("c", value=[10.0, 15.0, 8.0])
d = m_multi.parameter("d", value=50.0)

m_multi.minimize(dm.sum(c * x))
m_multi.subject_to(dm.sum(x) >= d, name="demand")
m_multi.subject_to(x[0] + 2 * x[1] <= 80.0, name="resource")

rc_multi = RobustCounterpart(
    m_multi,
    [
        BoxUncertaintySet(c, delta=[1.0, 1.5, 0.8]),
        BoxUncertaintySet(d, delta=5.0),
    ],
)
rc_multi.formulate()

r_multi = m_multi.solve()
xv = list(r_multi.x.values())[0]
print(f"Robust obj = {r_multi.objective:.2f},  x = {np.round(xv, 2)}")
print("\nRobust constraints (worst-case constants substituted):")
for con in m_multi._constraints:
    print(f"  [{con.name}]  {con.body} {con.sense} {con.rhs}")

Robust obj = 484.00,  x = [ 0.  0. 55.]

Robust constraints (worst-case constants substituted):
  [demand]  (55 - sum(x[3])) <= 0.0
  [resource]  ((x[3][0] + (2 * x[3][1])) - 80) <= 0.0

Example 4 — Drug Production (Ben-Tal et al. 2009, Example 1.1.1)#

This is the opening example from the textbook Robust Optimization [Ben-Tal et al., 2009]. A company produces two drugs from two raw materials. The active agent concentrations in the raw materials are uncertain, and the robust counterpart protects against worst-case yields.

\[ \max\; 5500 x_1 + 6100 x_2 - 100 r_1 - 199.9 r_2 \]

\[ \text{s.t.}\; a_1 r_1 + a_2 r_2 \ge 0.5 x_1 + 0.6 x_2, \quad r_1 + r_2 \le 1000, \quad 90 x_1 + 100 x_2 \le 2000 \]

\[ 40 x_1 + 50 x_2 \le 800, \quad 100 r_1 + 199.9 r_2 + 700 x_1 + 800 x_2 \le 100{,}000 \]

The uncertain concentrations are $a_1 = 0.01 \pm 0.5\%$ and $a_2 = 0.02 \pm 2\%$. The published solutions are: nominal profit = $8,819.66 (using only RawII), robust profit = $8,294.57 (switching to the more reliable RawI), a 5.95% price of robustness.

def drug_model(name):
    """Build the Ben-Tal Example 1.1.1 drug production model."""
    m = dm.Model(name)
    d1 = m.continuous("DrugI", lb=0, ub=1000)
    d2 = m.continuous("DrugII", lb=0, ub=1000)
    r1 = m.continuous("RawI", lb=0, ub=1000)
    r2 = m.continuous("RawII", lb=0, ub=1000)
    aI = m.parameter("aI", value=0.01)
    aII = m.parameter("aII", value=0.02)
    # Maximize profit (net revenue minus raw material cost)
    m.minimize(-(5500 * d1 + 6100 * d2 - 100 * r1 - 199.90 * r2))
    m.subject_to(aI * r1 + aII * r2 >= 0.5 * d1 + 0.6 * d2, name="balance")
    m.subject_to(r1 + r2 <= 1000, name="storage")
    m.subject_to(90 * d1 + 100 * d2 <= 2000, name="manpower")
    m.subject_to(40 * d1 + 50 * d2 <= 800, name="equipment")
    m.subject_to(100 * r1 + 199.90 * r2 + 700 * d1 + 800 * d2 <= 100000, name="budget")
    return m, aI, aII


# --- Nominal ---
m_nom, _, _ = drug_model("nominal")
r_nom = m_nom.solve()
print("Nominal solution:")
print(f"  Profit = ${-r_nom.objective:,.2f}")
for name, val in r_nom.x.items():
    if abs(val) > 0.01:
        print(f"  {name} = {val:.4f}")

# --- Robust (box uncertainty on active agent concentrations) ---
m_rob, aI, aII = drug_model("robust")
# aI = 0.01 +/- 0.5% = 0.01 +/- 0.00005
# aII = 0.02 +/- 2.0% = 0.02 +/- 0.0004
RobustCounterpart(
    m_rob, [BoxUncertaintySet(aI, delta=0.00005), BoxUncertaintySet(aII, delta=0.0004)]
).formulate()
r_rob = m_rob.solve()
print("\nRobust solution:")
print(f"  Profit = ${-r_rob.objective:,.2f}")
for name, val in r_rob.x.items():
    if abs(val) > 0.01:
        print(f"  {name} = {val:.4f}")

nom_profit = -r_nom.objective
rob_profit = -r_rob.objective
print(f"\nPrice of robustness: {(nom_profit - rob_profit) / nom_profit * 100:.2f}%")
print(
    "\nThe robust solution switches from RawII to RawI because RawI has"
    "\ntighter uncertainty (0.5% vs 2%), making it more reliable despite"
    "\nhigher cost per gram."
)

Nominal solution:
  Profit = $8,819.66
  DrugI = 17.5516
  RawII = 438.7889

Robust solution:
  Profit = $8,294.57
  DrugI = 17.4669
  RawI = 877.7319

Price of robustness: 5.95%

The robust solution switches from RawII to RawI because RawI has
tighter uncertainty (0.5% vs 2%), making it more reliable despite
higher cost per gram.

API Summary#

from discopt.ro import (
    AffineDecisionRule,        # adjustable robust (affine policy)
    BoxUncertaintySet,         # per-component interval
    EllipsoidalUncertaintySet, # ellipsoidal ball
    PolyhedralUncertaintySet,  # general polytope
    RobustCounterpart,         # static robust counterpart
    budget_uncertainty_set,    # Bertsimas-Sim convenience constructor
)

# -- Static robust optimization -----------------------------------------------
m = dm.Model()
p = m.parameter("p", value=nominal_value)
# ... build objective and constraints using p ...
rc = RobustCounterpart(m, BoxUncertaintySet(p, delta=0.1 * nominal_value))
rc.formulate()          # modifies m in-place; p replaced by worst-case constant
result = m.solve()

# -- Adjustable robust optimization (two-stage) -------------------------------
m = dm.Model()
x = m.continuous("x", ...)    # here-and-now decisions
y = m.continuous("y", ...)    # wait-and-see recourse
p = m.parameter("p", value=nominal_value)
# ... build objective and constraints using x, y, p ...

# Stage 1: affine rule  y(xi) = y0 + Y0*(p - p_bar)
adr = AffineDecisionRule(y, uncertain_params=p)
adr.apply()                    # retires y; adds y0, Y0 as decision variables

# Stage 2: handle remaining xi terms
rc = RobustCounterpart(m, BoxUncertaintySet(p, delta=delta))
rc.formulate()                 # fully deterministic model

result = m.solve()             # optimises over (x, y0, Y0)

Limitations#

Coefficient uncertainty (bilinear p * x) requires LP dual auxiliary variables; planned for a future release.
Recourse shape: AffineDecisionRule supports scalar and 1-D recourse variables; 2-D support is planned.
Non-affine recourse (polynomial, piecewise-linear) is not supported.
Distributionally robust optimization (moment-based or Wasserstein ambiguity sets) is not yet supported.

For a full adjustable RO tutorial including multi-parameter recourse and vector recourse variables, see Adjustable Robust Optimization.

Set	Definition	Reformulation
Box	\(\lvert p_j - \bar{p}_j \rvert \le \delta_j\)	Substitute worst-case bound per component
Ellipsoidal	\(\lVert \Sigma^{-1/2}(p-\bar{p}) \rVert_2 \le \rho\)	Add SOC norm penalty [Ben-Tal and Nemirovski, 1999]
Polyhedral	\(A\xi \le b,\; \xi = p - \bar{p}\)	LP dual auxiliary variables [Bertsimas and Sim, 2004]