import itertools
from math import prod
import pandas as pd
import pyomo.environ as pyo
# ---------------------------------------------------------------- data
JOBS = [1, 2, 3, 4]
CENTERS = [1, 2]
LOAD = {1: 10.0, 2: 8.0, 3: 6.0, 4: 12.0} # CPU-hours per job
RATE = {1: 1.0, 2: 2.5} # $ / CPU-hour
PENALTY_FACTOR = 1.6 # reroute penalty
STAGE = {1: 1, 2: 2, 3: 3, 4: 3} # job -> stage
# Markov chain on availability
PI0 = 0.60 # P(a_1 = 1)
P11 = 0.80 # P(a_t = 1 | a_{t-1} = 1)
P01 = 0.20 # P(a_t = 1 | a_{t-1} = 0)
# ----------------------------------------------- scenario tree -------
SCEN = list(itertools.product([0, 1], repeat=3)) # (a1, a2, a3)
def scen_prob(s):
a1, a2, a3 = s
p1 = PI0 if a1 == 1 else 1 - PI0
p2 = (P11 if a2 == 1 else 1 - P11) if a1 == 1 else (P01 if a2 == 1 else 1 - P01)
p3 = (P11 if a3 == 1 else 1 - P11) if a2 == 1 else (P01 if a3 == 1 else 1 - P01)
return p1 * p2 * p3
PROB = {s: scen_prob(s) for s in SCEN}
assert abs(sum(PROB.values()) - 1.0) < 1e-9
# Marginal P(a_t = 1) for the EV problem
P_MARG = {1: PI0}
P_MARG[2] = PI0 * P11 + (1 - PI0) * P01
P_MARG[3] = P_MARG[2] * P11 + (1 - P_MARG[2]) * P01
# ------------------------------------------- realized cost helper ----
def realized_cost(j, c, a_t):
"""Cost when job j is assigned to center c and stage availability is a_t."""
if c == 2:
return RATE[2] * LOAD[j]
# c == 1, two cases on availability
return RATE[1] * LOAD[j] if a_t == 1 else PENALTY_FACTOR * RATE[2] * LOAD[j]
# ============================================================== RP ===
def build_RP():
"""Multistage SP with explicit non-anticipativity constraints."""
m = pyo.ConcreteModel("RP")
# declare a Pyomo Set with dimen=3 so the 3-tuple is a single index
m.S = pyo.Set(initialize=SCEN, dimen=3)
m.x = pyo.Var(JOBS, CENTERS, m.S, domain=pyo.Binary)
# one center per (job, scenario)
m.Assign = pyo.ConstraintList()
for j in JOBS:
for s in SCEN:
m.Assign.add(sum(m.x[j, c, s] for c in CENTERS) == 1)
# ---- non-anticipativity ----
# Stage 1 (job 1): identical across all scenarios
s_ref = SCEN[0]
m.NA1 = pyo.ConstraintList()
for c in CENTERS:
for s in SCEN:
if s != s_ref:
m.NA1.add(m.x[1, c, s] == m.x[1, c, s_ref])
# Stage 2 (job 2): depends only on a1
m.NA2 = pyo.ConstraintList()
for a1 in (0, 1):
group = [s for s in SCEN if s[0] == a1]
ref = group[0]
for c in CENTERS:
for s in group[1:]:
m.NA2.add(m.x[2, c, s] == m.x[2, c, ref])
# Stage 3 (jobs 3 and 4): depend only on (a1, a2)
m.NA3 = pyo.ConstraintList()
for a1, a2 in itertools.product([0, 1], repeat=2):
group = [s for s in SCEN if s[0] == a1 and s[1] == a2]
ref = group[0]
for j in (3, 4):
for c in CENTERS:
for s in group[1:]:
m.NA3.add(m.x[j, c, s] == m.x[j, c, ref])
# expected total cost
def obj_rule(m):
total = 0.0
for s in SCEN:
cost_s = 0.0
for j in JOBS:
a_t = s[STAGE[j] - 1]
for c in CENTERS:
cost_s += realized_cost(j, c, a_t) * m.x[j, c, s]
total += PROB[s] * cost_s
return total
m.Obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)
return m
# ============================================================== WS ===
def build_WS(s):
"""Wait-and-see (perfect information) for a single scenario."""
m = pyo.ConcreteModel(f"WS_{s}")
m.x = pyo.Var(JOBS, CENTERS, domain=pyo.Binary)
m.Assign = pyo.Constraint(
JOBS, rule=lambda m, j: sum(m.x[j, c] for c in CENTERS) == 1
)
m.Obj = pyo.Objective(
expr=sum(
realized_cost(j, c, s[STAGE[j] - 1]) * m.x[j, c]
for j in JOBS for c in CENTERS
),
sense=pyo.minimize,
)
return m
# ============================================================== EV ===
def build_EV():
"""Expected-value problem: replace random availability by its
marginal probability and minimise the resulting deterministic cost.
"""
m = pyo.ConcreteModel("EV")
m.x = pyo.Var(JOBS, CENTERS, domain=pyo.Binary)
m.Assign = pyo.Constraint(
JOBS, rule=lambda m, j: sum(m.x[j, c] for c in CENTERS) == 1
)
def exp_cost(j, c):
if c == 2:
return RATE[2] * LOAD[j]
p = P_MARG[STAGE[j]]
return p * RATE[1] * LOAD[j] + (1 - p) * PENALTY_FACTOR * RATE[2] * LOAD[j]
m.Obj = pyo.Objective(
expr=sum(exp_cost(j, c) * m.x[j, c] for j in JOBS for c in CENTERS),
sense=pyo.minimize,
)
return m
# =========================================================== driver ==
def main():
solver = pyo.SolverFactory("glpk")
# ---------- RP ----------
rp = build_RP()
solver.solve(rp)
rp_val = pyo.value(rp.Obj)
rp_rows = []
for s in SCEN:
cost_s = 0.0
assign = {}
for j in JOBS:
assign[j] = next(c for c in CENTERS if pyo.value(rp.x[j, c, s]) > 0.5)
cost_s += realized_cost(j, assign[j], s[STAGE[j] - 1])
rp_rows.append({
"scenario (a1,a2,a3)": s,
"prob": PROB[s],
**{f"job{j}": f"C{assign[j]}" for j in JOBS},
"cost": cost_s,
})
# ---------- WS per scenario ----------
ws_rows = []
ws_exp = 0.0
for s in SCEN:
m = build_WS(s)
solver.solve(m)
c_val = pyo.value(m.Obj)
ws_exp += PROB[s] * c_val
assign = {j: next(c for c in CENTERS if pyo.value(m.x[j, c]) > 0.5) for j in JOBS}
ws_rows.append({
"scenario (a1,a2,a3)": s,
"prob": PROB[s],
**{f"job{j}": f"C{assign[j]}" for j in JOBS},
"cost": c_val,
})
# ---------- EV / EEV ----------
ev = build_EV()
solver.solve(ev)
ev_val = pyo.value(ev.Obj)
ev_assign = {j: next(c for c in CENTERS if pyo.value(ev.x[j, c]) > 0.5) for j in JOBS}
eev_exp = 0.0
eev_rows = []
for s in SCEN:
cost_s = sum(realized_cost(j, ev_assign[j], s[STAGE[j] - 1]) for j in JOBS)
eev_exp += PROB[s] * cost_s
eev_rows.append({
"scenario (a1,a2,a3)": s,
"prob": PROB[s],
**{f"job{j}": f"C{ev_assign[j]}" for j in JOBS},
"cost": cost_s,
})
# ---------- pretty print ----------
df_rp = pd.DataFrame(rp_rows)
df_ws = pd.DataFrame(ws_rows)
df_eev = pd.DataFrame(eev_rows)
pd.set_option("display.width", 120)
pd.set_option("display.float_format", lambda v: f"{v:8.4f}")
bar = "=" * 78
sub = "-" * 78
print(bar)
print(" MULTISTAGE STOCHASTIC OPTIMIZATION -- GENOMICS WORKFLOW")
print(bar)
print(f" Jobs : {JOBS} loads (CPU-h) = {LOAD}")
print(f" Centers : C1 rate {RATE[1]}/h, C2 rate {RATE[2]}/h,"
f" reroute penalty x{PENALTY_FACTOR}")
print(f" Markov chain : pi0={PI0}, P(1|1)={P11}, P(1|0)={P01}")
print(f" Marginal P(a_t=1): "
+ ", ".join(f"t{t}={P_MARG[t]:.3f}" for t in (1, 2, 3)))
print(f" Scenarios : {len(SCEN)}\n")
print(sub)
print(" 1. WAIT-AND-SEE (perfect information per scenario)")
print(sub)
print(df_ws.to_string(index=False))
print(f"\n E[WS] = sum_s p_s * WS(s) = {ws_exp:.4f}\n")
print(sub)
print(" 2. RECOURSE PROBLEM RP (multistage stochastic optimum)")
print(sub)
print(df_rp.to_string(index=False))
print(f"\n RP (objective value) = {rp_val:.4f}\n")
print(sub)
print(" 3. MEAN-VALUE EV / EEV (deterministic, then evaluated)")
print(sub)
print(f" EV optimum (deterministic, marginal probs) = {ev_val:.4f}")
print(f" EV first-stage policy : "
+ ", ".join(f"job{j}->C{ev_assign[j]}" for j in JOBS))
print(df_eev.to_string(index=False))
print(f"\n EEV (EV policy evaluated under SP) = {eev_exp:.4f}\n")
print(sub)
print(" 4. STOCHASTIC PROGRAMMING METRICS")
print(sub)
print(f" E[WS] = {ws_exp:10.4f} (lower bound on RP)")
print(f" RP = {rp_val:10.4f} (true SP optimum)")
print(f" EEV = {eev_exp:10.4f} (using mean-value plan)")
print(f" VPI = RP - E[WS] = {rp_val - ws_exp:10.4f} (value of perfect information)")
print(f" VSS = EEV - RP = {eev_exp - rp_val:10.4f} (value of stochastic solution)")
print(bar)
# ---------- export ----------
with pd.ExcelWriter("./genomics_results.xlsx") as xw:
df_ws.to_excel(xw, sheet_name="WS", index=False)
df_rp.to_excel(xw, sheet_name="RP", index=False)
df_eev.to_excel(xw, sheet_name="EV", index=False)
pd.DataFrame([
{"metric": "E[WS]", "value": ws_exp},
{"metric": "RP", "value": rp_val},
{"metric": "EEV", "value": eev_exp},
{"metric": "VPI", "value": rp_val - ws_exp},
{"metric": "VSS", "value": eev_exp - rp_val},
]).to_excel(xw, sheet_name="metrics", index=False)Ex. S2: Multi-Stage Stochastic
The multi-stage stochastic scheduling problem is a straightfoward extension of the two-stage cases that we have studied previously. It provides the possibility to model a more complex decision-making process, where uncertain, exogenous inputs act on the system at each stage of its evolution.
We formulate the Multistage Stochastic Optimization for a Genomics Workflow.
Workflow DAG
job1 --> job2 +--> job3
+--> job4Each job j is assigned to one of two computing centers c in {1, 2}:
- Center 1 : low cost rate, but uncertain availability;
- Center 2 : higher cost rate, always available.
If a job is sent to center 1 and center 1 turns out to be unavailable when the job runs, the job is forcibly rerouted to center 2 and incurs a penalty on top of the nominal center‑2 rate. The penalty take sinto account the operational cost of last‑minute migration: queue waits, data movement, restart overhead. It is assumed that the cost of the penalty is known in advance (cost = penalty_factor * rate2 * load).
Stage / information structure
- Stage 1 – assign job 1 (HERE-AND-NOW, before any availability is known)
- Stage 2 – assign job 2, knowing the realized availability a1
- Stage 3 – assign job 3 and job 4, knowing (a1, a2)
The realized availability a3 of stage-3 is observed AFTER the stage-3 decisions, so it shows up only in the cost.
Availability is modelled as a 2-state Markov chain to make the multi- stage structure non-trivial (under iid, the stochastic and mean-value solutions coincide and VSS = 0).
P(a1 = 1) = pi0
P(a_t = 1 | a_{t-1} = 1) = p11 (sticky available)
P(a_t = 1 | a_{t-1} = 0) = p01 (sticky unavailable)Models
RP multistage Recourse Problem with non-anticipativity --> EVSS
PI Perfect Information per scenario --> EVPI
EV Expected-Value (deterministic) problem using marginal probs --> EVMReported metrics
VPI = EVSS - EVPI (value of Perfect Information)
VSS = EVM - EVSS (value of the Stochastic Solution)The decision maker assigns each job to a center before observing whether center 1 will be available when that job runs. The information available at the time of each decision is the realised availability of center 1 at all earlier stages.
Numerical values:
| Symbol | Meaning | Value |
|---|---|---|
| \(L_j\) | computational load of job \(j\) (CPU‑h) | \(L_1=10,\;L_2=8,\;L_3=6,\;L_4=12\) |
| \(r_1, r_2\) | unit rates of centers 1 and 2 ($/CPU‑h) | \(1.0,\;2.5\) |
| \(\alpha\) | reroute penalty multiplier | \(1.6\) |
| \(\pi_0\) | \(\mathbb{P}(a_1=1)\) | \(0.6\) |
| \(p_{11}\) | \(\mathbb{P}(a_{t+1}=1\mid a_t=1)\) | \(0.8\) |
| \(p_{01}\) | \(\mathbb{P}(a_{t+1}=1\mid a_t=0)\) | \(0.2\) |
Define the binary availability variables \(a_t\in\{0,1\}\) for \(t\in\mathcal{T}\) with \(a_t=1\) meaning “center 1 available during stage \(t.\)”
The stochastic process \((a_1,a_2,a_3)\) is a 2‑state Markov chain with parameters \((\pi_0,p_{11},p_{01})\). On the probability space \((\Omega,\mathcal{F},\mathbb{P})\) it generates the filtration
\[\mathcal{F}_0=\{\emptyset,\Omega\},\qquad \mathcal{F}_t=\sigma(a_1,\dots,a_t),\quad t=1,2,3.\]
The decision at stage \(t\) is taken under the information \(\mathcal{F}_{t-1}\):
| Stage \(t\) | Job(s) decided | Information set |
|---|---|---|
| 1 | job 1 | \(\mathcal{F}_0\) — none |
| 2 | job 2 | \(\mathcal{F}_1=\sigma(a_1)\) |
| 3 | jobs 3, 4 | \(\mathcal{F}_2=\sigma(a_1,a_2)\) |
The third realization \(a_3\) is observed only after the stage‑3 decisions and enters the model purely through the realized cost.
The scenario tree
Because the chain has two states and three stages, the scenario tree has a single root, branches with multiplicity 2 at each non‑leaf node, and \(|\mathcal{S}|=2^3=8\) leaves.
Let \(\mathcal{N}\) be the node set, \(\mathcal{N}_t\) its stage‑\(t\) slice, and write \(a(n)\) for the parent of \(n\) and \(\pi_n\) for its unconditional probability. The tree has
\[|\mathcal{N}_1|=1,\qquad |\mathcal{N}_2|=2,\qquad |\mathcal{N}_3|=4,\qquad |\mathcal{L}|=8.\]
For a leaf \(\ell\) corresponding to scenario \(s=(a_1^s,a_2^s,a_3^s)\), write \(n_t(\ell)\) for its stage‑\(t\) ancestor. The scenario probability is
\[p_s \;=\; \mathbb{P}(a_1=a_1^s)\,\mathbb{P}(a_2=a_2^s\mid a_1=a_1^s)\,\mathbb{P}(a_3=a_3^s\mid a_2=a_2^s),\]
which evaluates, for our parameters, to:
| Scenario \((a_1,a_2,a_3)\) | \(p_s\) | details |
|---|---|---|
| \((1,1,1)\) | \(0.384\) | \(p(a_1=1)\cdot p(a_2=1\mid a_1=1)\cdot p(a_3=1\mid a_2=1)\) |
| \((1,1,0)\) | \(0.096\) | \(p(a_1=1)\cdot p(a_2=1\mid a_1=1)\cdot p(a_3=0\mid a_2=1)\) |
| \((1,0,1)\) | \(0.024\) | \(\vdots\) |
| \((1,0,0)\) | \(0.096\) | \(\vdots\) |
| \((0,1,1)\) | \(0.064\) | \(\vdots\) |
| \((0,1,0)\) | \(0.016\) | \(\vdots\) |
| \((0,0,1)\) | \(0.064\) | \(\vdots\) |
| \((0,0,0)\) | \(0.256\) | \(p(a_1=0)\cdot p(a_2=0\mid a_1=0)\cdot p(a_3=0\mid a_2=0)\) |
where, for example, \[p(a_1=1)\cdot p(a_2=1\mid a_1=1)\cdot p(a_3=1\mid a_2=1)=0.6 \times 0.8\times 0.8 = 0.384,\] \[p(a_1=0)\cdot p(a_2=0\mid a_1=0)\cdot p(a_3=0\mid a_2=0)=0.4 \times (1-0.2)\times (1-0.2) = 0.256. \]
The eight probabilities sum to one and the strong persistence in the chain (\(p_{11}=0.8\), \(p_{00}=0.8\)) makes the corner scenarios \((1,1,1)\) and \((0,0,0)\) disproportionately likely.
Cost model
Define the realized cost of running job \(j\) on center \(c\) when stage‑\(\tau(j)\) availability of center 1 is \(a\)
\[\kappa(j,c,a)\;=\;\begin{cases} r_1 L_j & c=1,\;a=1 \\ \alpha\,r_2 L_j & c=1,\;a=0 \\ r_2 L_j & c=2. \end{cases}\]
For the numerical instance, the per‑job cost matrix becomes
| Job \(j\) | \(\kappa(j,1,1)\) | \(\kappa(j,1,0)\) | \(\kappa(j,2,\cdot)\) |
|---|---|---|---|
| 1 | 10 | 40 | 25 |
| 2 | 8 | 32 | 20 |
| 3 | 6 | 24 | 15 |
| 4 | 12 | 48 | 30 |
The crucial inequalities \(r_1 L_j < r_2 L_j < \alpha r_2 L_j\) hold for every \(j\): center 1 granted is strictly cheapest, center 2 is the safe middle, and center 1 rerouted is strictly worst. This monotonicity is what gives the optimal policy its threshold character.
Extensive form — scenario indexing with explicit non‑anticipativity
The formulation duplicates variables per scenario and adds equality constraints. Let \(x_{j,c,s}\in[0,1]\) for \(j\in\mathcal{J}\), \(c\in\mathcal{C}\), \(s\in\mathcal{S}.\) Define the equivalence relation on scenarios
\[s\;\sim_t\;s'\quad\iff\quad (a_1^s,\dots,a_{t-1}^s)=(a_1^{s'},\dots,a_{t-1}^{s'}),\]
i.e. \(s\) and \(s'\) are indistinguishable using only the information available at the start of stage \(t\). The deterministic equivalent is
\[ \begin{aligned} \min_{x}\quad & \sum_{s\in\mathcal{S}} p_s \sum_{j,c}\kappa\!\bigl(j,c,a^s_{\tau(j)}\bigr)\,x_{j,c,s} \\ \text{s.t.}\quad & \sum_{c}x_{j,c,s}=1, && \forall j,s, \\ & x_{j,c,s}=x_{j,c,s'}, && \forall j,c,\;\forall s\sim_{\tau(j)} s' \quad(\text{non-anticipativity})\\ & 0\le x_{j,c,s}\le 1. \end{aligned} \]
#if __name__ == "__main__":
main()==============================================================================
MULTISTAGE STOCHASTIC OPTIMIZATION -- GENOMICS WORKFLOW
==============================================================================
Jobs : [1, 2, 3, 4] loads (CPU-h) = {1: 10.0, 2: 8.0, 3: 6.0, 4: 12.0}
Centers : C1 rate 1.0/h, C2 rate 2.5/h, reroute penalty x1.6
Markov chain : pi0=0.6, P(1|1)=0.8, P(1|0)=0.2
Marginal P(a_t=1): t1=0.600, t2=0.560, t3=0.536
Scenarios : 8
------------------------------------------------------------------------------
1. WAIT-AND-SEE (perfect information per scenario)
------------------------------------------------------------------------------
scenario (a1,a2,a3) prob job1 job2 job3 job4 cost
(0, 0, 0) 0.2560 C2 C2 C2 C2 90.0000
(0, 0, 1) 0.0640 C2 C2 C1 C1 63.0000
(0, 1, 0) 0.0160 C2 C1 C2 C2 78.0000
(0, 1, 1) 0.0640 C2 C1 C1 C1 51.0000
(1, 0, 0) 0.0960 C1 C2 C2 C2 75.0000
(1, 0, 1) 0.0240 C1 C2 C1 C1 48.0000
(1, 1, 0) 0.0960 C1 C1 C2 C2 63.0000
(1, 1, 1) 0.3840 C1 C1 C1 C1 36.0000
E[WS] = sum_s p_s * WS(s) = 59.8080
------------------------------------------------------------------------------
2. RECOURSE PROBLEM RP (multistage stochastic optimum)
------------------------------------------------------------------------------
scenario (a1,a2,a3) prob job1 job2 job3 job4 cost
(0, 0, 0) 0.2560 C1 C2 C2 C2 105.0000
(0, 0, 1) 0.0640 C1 C2 C2 C2 105.0000
(0, 1, 0) 0.0160 C1 C2 C1 C1 132.0000
(0, 1, 1) 0.0640 C1 C2 C1 C1 78.0000
(1, 0, 0) 0.0960 C1 C1 C2 C2 87.0000
(1, 0, 1) 0.0240 C1 C1 C2 C2 87.0000
(1, 1, 0) 0.0960 C1 C1 C1 C1 90.0000
(1, 1, 1) 0.3840 C1 C1 C1 C1 36.0000
RP (objective value) = 73.6080
------------------------------------------------------------------------------
3. MEAN-VALUE EV / EEV (deterministic, then evaluated)
------------------------------------------------------------------------------
EV optimum (deterministic, marginal probs) = 83.6160
EV first-stage policy : job1->C1, job2->C1, job3->C1, job4->C1
scenario (a1,a2,a3) prob job1 job2 job3 job4 cost
(0, 0, 0) 0.2560 C1 C1 C1 C1 144.0000
(0, 0, 1) 0.0640 C1 C1 C1 C1 90.0000
(0, 1, 0) 0.0160 C1 C1 C1 C1 120.0000
(0, 1, 1) 0.0640 C1 C1 C1 C1 66.0000
(1, 0, 0) 0.0960 C1 C1 C1 C1 114.0000
(1, 0, 1) 0.0240 C1 C1 C1 C1 60.0000
(1, 1, 0) 0.0960 C1 C1 C1 C1 90.0000
(1, 1, 1) 0.3840 C1 C1 C1 C1 36.0000
EEV (EV policy evaluated under SP) = 83.6160
------------------------------------------------------------------------------
4. STOCHASTIC PROGRAMMING METRICS
------------------------------------------------------------------------------
E[WS] = 59.8080 (lower bound on RP)
RP = 73.6080 (true SP optimum)
EEV = 83.6160 (using mean-value plan)
VPI = RP - E[WS] = 13.8000 (value of perfect information)
VSS = EEV - RP = 10.0080 (value of stochastic solution)
==============================================================================Numerical results
Solving with GLPK on the instance defined above yields:
| Quantity | Value |
|---|---|
| \(\mathbb{E}[\mathrm{WS}]\) (perfect information) | \(59.81\) |
| \(\mathrm{RP}\) (multistage stochastic optimum) | \(73.61\) |
| \(\mathrm{EEV}\) (mean‑value plan, evaluated under the SP) | \(83.62\) |
| \(\mathrm{VPI}=\mathrm{RP}-\mathbb{E}[\mathrm{WS}]\) | \(13.80\) |
| \(\mathrm{VSS}=\mathrm{EEV}-\mathrm{RP}\) | \(10.01\) |
The optimal RP policy is the threshold policy:
- job 1 always goes to center 1;
- job 2 follows \(a_1\);
- jobs 3 and 4 follow \(a_2.\)
The mean‑value plan, in contrast, sends every job to center 1 — which happens to match RP at stage 1 but is wrong at stages 2 and 3, costing roughly 13% in expectation.
The perfect information (wait-and-see) solutions reveal an important contrast: with full information, even job 1 sometimes goes to center 2 (any scenario starting with \(a_1=0\)). When you know in advance that center 1 is unavailable, paying the rerouting penalty is strictly worse than choosing center 2 outright. RP cannot replicate this because the job‑1 decision must be made before \(a_1\) is observed, which is precisely what the VPI of \(13.80\) measures.