Adaptive Scenario Generation for KL-DRO Around a Learned Density Model

A DAPTIVE S CENARIO G ENERATION FOR KL-DRO A ROUND A L EARNED D ENSITY M ODEL

arXiv:2610.40644v1 [eess.SY] 29 Oct 2026

Aron Brenner & Saurabh Amin Laboratory for Information & Decision Systems Massachusetts Institute of Technology Cambridge, MA 02139, USA {abrenner,amins}@mit.edu

A BSTRACT Distributionally robust optimization (DRO) hedges data-driven two-stage stochastic programs against error in the distribution estimated from a finite record, and standard ambiguity sets differ in how they treat realizations absent from that record. A divergence ball around the empirical distribution contains only reweightings of the record, whereas a transport ball also contains distributions supported off the record, including on realizations the data-generating process cannot produce, against which the decision is then hedged as well. We instead propose a Kullback– Leibler (KL) ball centered at a density model fitted to the record, the learned nominal, which lets the adversary up-weight costly unobserved realizations in proportion to the probability the model assigns them. Solving the resulting problem poses a second difficulty, because it can accommodate only a few scenarios, and an independent sample of that size is likely to omit the rare realizations on which the worst case concentrates. We therefore build the scenario set adaptively, alternating between solving the finite problem and drawing new realizations from the learned nominal, steered toward the worst case at the current decision, with importance weights that correct the draws back to the model. On a multi-item newsvendor problem and a capacity expansion problem on realistic grid data, centering the ball at the learned nominal rather than at the record lowers cost under distribution shift, at comparable in-distribution cost, when the record is scarce and the model class captures structure that lets it extrapolate beyond the record. The advantage fades as the record grows or as the shift departs from the structure the model has learned.

1

I NTRODUCTION

Two-stage stochastic programs model decisions that must be committed before uncertainty is realized, with recourse actions taken after, in applications spanning finance, supply chains, telecommunications, and energy systems (Wallace & Ziemba, 2005; Birge & Louveaux, 2011). We consider their datadriven, distributionally robust form, n  o min c(x) + max RP Q(x, ξ) , Q(x, ξ) := min q(y), (1) x∈X

y∈Y(ξ,x)

P∈P(Q)

where x denotes here-and-now decisions with cost c(x), ξ ∈ RD the random vector, Q(x, ξ) the cost of the wait-and-see decisions y, P(Q) an ambiguity set of distributions considered plausible around a nominal distribution Q (in the standard setting, the empirical distribution Q̂N of a record of N past realizations), and RP a risk measure, such as the expectation or the conditional value at risk (CVaR). We are motivated by problems whose cost is dominated by tail events, in which several dimensions of uncertainty move against the plan at once, such as a demand surge that arrives while a supplier is down. Such problems expose two gaps in how data-driven distributionally robust optimization (DRO) treats realizations absent from the record. The first lies in the ambiguity set (Sec. 3). Divergence balls around the empirical distribution contain only reweightings of the record, while Wasserstein balls and moment-based sets admit unobserved realizations but judge them only by their distance from the 1

record or by a few moments. Each family carries a statistical guarantee tied to the record (Duchi et al., 2021; Mohajerin Esfahani & Kuhn, 2018), but none uses how the components of the uncertainty move together in the tail, the structure that distinguishes plausible unobserved realizations from implausible ones. The second gap lies in the finite scenario set on which the problem is actually solved, whose size is limited because each scenario adds a copy of the recourse problem. A divergence ball can only reweight this set. Scenario reduction and problem-driven selection choose it from the record or from a sample drawn once, so the set is likely to miss the rare realizations on which the worst case concentrates. Both gaps concern realizations that the record lacks, and both could be narrowed by a model that generates plausible ones. High-resolution operational records, often spanning only a few years, are available in many applications of stochastic programming, and generative models trained on such records can capture their dependence structure and respect known properties of the process, including the behavior of its extremes (Boulaguiem et al., 2022; Pandey et al., 2025; Chang & Sapsis, 2026). Such a model offers both a new center for the ambiguity set, addressing the first gap, and a new source of scenarios, addressing the second. This raises the question of how a generative model learned from a short record can be integrated into risk-averse stochastic programming, so that a limited set of scenarios hedges against tail risks the record underrepresents. We propose Kullback–Leibler DRO (KL-DRO) centered at a learned nominal: we take Q in equation 1 to be a density model Qθ trained on the record, such as a normalizing flow, a Gaussian copula, or a diffusion model, and the ambiguity set to be  Pρ (Qθ ) := P : KL(P ∥ Qθ ) ≤ ρ . (2) Over this ball, the adversary can up-weight costly realizations absent from the record in proportion to the probability Qθ assigns them, so what it can reach is governed by the learned nominal rather than by a distance from the record, which addresses the first gap. To close the second, we build the scenario set adaptively, alternating between solving the finite problem and drawing new realizations from the learned nominal, steered toward the worst case at the current decision, with importance weights that center the ball at an estimate of Qθ . Centering the ball at a learned nominal does, however, give up the finite-sample guarantees of balls centered at the record, so we choose the radius empirically (Sec. 5). Contributions.

Our contributions are as follows.

• Distributionally robust optimization with a learned nominal. We pose KL-DRO with the ambiguity set centered at a learned generative model, extending balls centered at parametric models or Bayesian posteriors for single-stage objectives (Watson & Holmes, 2016; Dellaporta et al., 2025) to high-dimensional nominals and two-stage objectives. • Finite scenario approximations for tail-dominated problems. We build the scenario set by adaptive importance sampling toward the ball’s worst case at each decision, so that the scenarios are drawn where the adversary places its weight. • Evaluations in controlled and realistic settings. We evaluate the framework on a multiitem newsvendor problem with a known distribution and controlled record size and shift, and on a capacity expansion problem over a 20-zone network with realistic demand and renewable data and an out-of-distribution test year.

2

BACKGROUND

Risk-averse stochastic programming. Given a limited record of past realizations, we seek decisions that are cost-effective both in distribution, even when the record is too small to learn the data-generating process, and out of distribution, when that process changes at evaluation time. To weigh both typical and tail costs, we take R to be the mean–CVaR objective, which can be defined as the minimum over a scalar ζ of an expectation (Rockafellar & Uryasev, 2000), RP [Q(x, ξ)] = (1 − β) EP [Q(x, ξ)] + β CVaRα,P [Q(x, ξ)]  = min ζ + EP [ψζ (Q(x, ξ))] , (3) ζ∈R

β ψζ (u) = (1 − β)(u − ζ) + 1−α (u − ζ)+ ,

2

where α is the tail level, β ∈ [0, 1] the weight on the tail, and the minimizing ζ the value at risk at level α. Throughout, we assume that the recourse problem is feasible (relatively complete recourse) and bounded for every x ∈ X and every ξ in the support of the nominal, so that Q(x, ξ) is finite. Kullback–Leibler DRO. For a nominal distribution Q, empirical or otherwise, the Kullback– Leibler (KL) ambiguity set {P : KL(P ∥ Q) ≤ ρ} is a standard choice for two-stage problems (Ben-Tal et al., 2013; Bayraksan & Love, 2015). Because every member assigns zero probability wherever Q does, the ball reweights the realizations Q supports and adds none, and when centered at the empirical record, it changes only the probabilities of the recorded scenarios. The ball is computationally attractive because dualizing the KL constraint reformulates the inner maximization of equation 1 as a minimization (Hu & Hong, 2013; App. B.5), n 
o max RP [Q(x, ξ)] = min ζ + λρ + λ log EQ exp ψζ (Q(x, ξ))/λ . (4) P : KL(P ∥ Q)≤ρ

ζ, λ>0

For ρ > 0, both sides are finite only if EQ [exp(Q(x, ξ)/λ)] < ∞ for some λ > 0, which holds for any finite nominal and any bounded cost. The dual is jointly convex in ζ and λ, and in x as well when the recourse cost is convex in x, as it is for linear recourse. In particular, for a finite nominal PK Q̂ = i=1 wi δξi , a set of K scenarios with weights wi > 0 summing to one, the expectation is a weighted sum, and equation 1 over the ball becomes, when X and c are convex, the single convex program min

x∈X , ζ, λ>0

c(x) + ζ + λρ + λ log

K X

wi exp ψζ (Q(x, ξi ))/λ ,

(5)

i=1

in which each scenario contributes one copy of the recourse problem and the risk measure adds one variable. The empirical measure Q̂N is the special case with K = N and wi = 1/N . Since the reformulation holds for any finite nominal, two modeling choices remain open: which realizations the adversary can reach, which is set by Q, and which finite Q̂ the problem is solved on. We discuss them in turn. Ambiguity centered at a learned nominal. The first choice determines which realizations the adversary can reach. Ideally, the nominal places mass on realizations the record lacks but not on those the data-generating process cannot produce, so that the decision is hedged against unrecorded events without paying for implausible ones. A density model Qθ fitted to the record places mass between and beyond the recorded realizations, distributed according to the dependence structure it has learned, so the ball around it contains distributions that weight events the record lacks. How well it withholds mass from implausible realizations depends on how faithfully it captures the dependence structure of the process. Because Qθ is itself an estimate, the ball around it also hedges against error in the model, as in Watson & Holmes (2016) and Dellaporta et al. (2025). We therefore adopt the ambiguity set equation 2, the KL ball above centered at Qθ rather than at the record. To make the nominal easy to sample and its draws easy to weight, we take Qθ to be the pushforward of the standard normal distribution φ := N (0, I), whose density we also write φ, by a deterministic decoder Gθ , so that sampling amounts to drawing z ∼ φ and applying Gθ . Normalizing flows, Gaussian copulas, and diffusion models with deterministic samplers are of this form, and the known latent density is what gives the sampler of Sec. 4 exact importance weights. Finite scenario approximation. The second choice is the finite nominal on which the problem is solved. The finite program equation 5, which we call the master problem, grows with each scenario, so the number of scenarios K must be small, and the decision is then largely determined by which K scenarios the master sees. When the realizations that drive the decision are rare, such as the few periods in which a demand surge coincides with a supplier outage, K independent draws are likely to miss them. Scenario reduction and problem-driven selection mitigate this by choosing the K scenarios from a larger pool (Sec. 3), but a pool fixed in advance cannot supply realizations it lacks. A learned nominal, by contrast, can be sampled wherever it has mass, so the K scenarios can be generated where the decision is sensitive. Because such scenarios are not drawn from Qθ itself, each carries a base weight wi , an importance weight back to Qθ , and we solve equation 5 with these weights normalized to sum to one, so that the ball is centered at an importance-weighted estimate of Qθ rather than at the distribution the scenarios were drawn from. 3

3

R ELATED W ORK

Ambiguity sets. Divergence-based sets (Ben-Tal et al., 2013; Bayraksan & Love, 2015) reweight the empirical record and so preclude unobserved realizations, or, with the arguments of the divergence reversed, admit them only on a fixed finite support (Van Parys et al., 2021). Transport-based sets (Mohajerin Esfahani & Kuhn, 2018; Gao & Kleywegt, 2023) admit them, but judge them by a cost that measures how far a recorded point moves and says nothing directly about whether the moved point is plausible. The Sinkhorn variant smooths this cost against a reference measure chosen in advance, such as a Gaussian, which sets where the worst case can place mass (Wang et al., 2021). Moment-based sets (Delage & Ye, 2010) constrain only low-order statistics. Closest in motivation are works that learn part of the set: the transport cost, from labels or the data manifold (Blanchet et al., 2019; Liu et al., 2022), or the nominal itself, a fitted model or Bayesian posterior at which a KL-based set is centered (Watson & Holmes, 2016; Dellaporta et al., 2025).

Scenario construction. Distribution-driven methods match specified moments (Høyland et al., 2003) or reduce a pool under a probability metric (Dupačová et al., 2003; Heitsch & Römisch, 2003) without reference to the decision. Problem-driven methods cluster by decision-based dissimilarity or minimize implementation error (Hewitt et al., 2022; Keutchayan et al., 2023), or, for tail risk measures, sample only the region of the support on which the risk measure depends (Fairbrother et al., 2022). All of these methods construct the set once, before the problem is solved, whereas the importance sampling methods of Sec. 4 adapt the proposal to the decision.

Generative models in optimization. A more recent line of work samples scenarios from a generative model fitted to the record (Chen et al., 2018; Jiang et al., 2018) or trains the generator on the downstream decision loss (Zhao et al., 2026; Islip et al., 2025), producing scenarios but no ambiguity set. A second line makes the generator the adversary, through a learned transport map to a continuous least favorable distribution (Xu et al., 2024; Cheng et al., 2025), a family of parametric generative models within a KL budget (Michel et al., 2021), a family of generators of bounded reconstruction loss (Wen & Yang, 2026), a Sinkhorn comparison between conditional generators (Zhang et al., 2026), or, for two-stage robust optimization, a search of a generator’s latent space for a single worst-case realization (Brenner et al., 2025). In each, the worst case is obtained by optimizing over that parameterization, and all but the last are posed for single-stage objectives. We instead center a divergence ball at the generator, which leaves the generator as trained, gives the worst case in closed form, and keeps the likelihood ratios exact.

4

A DAPTIVE I MPORTANCE S AMPLING A LGORITHM

Section 2 called for a finite scenario set generated where the decision is sensitive and weighted back to Qθ . We propose building it iteratively by adaptive importance sampling toward a target known up to a normalizing constant (Bugallo et al., 2017), re-targeting as the decision changes, as Dantzig & Glynn (1990), Infanger (1992), Parpas et al. (2015), and Aolaritei et al. (2026) do for a known density. Each round tilts a sample of the model toward the worst case at the current decision and forms the batch either by resampling that sample, which needs only draws of the model, or by drawing from a Gaussian fitted to the tilt by weighted moment matching, a cross-entropy step (Rubinstein & Kroese, 2004) that can reach realizations beyond the initial sample. Our approach differs from these methods in its target and its nominal. First, our target is the worst case of the KL ball, so the draws supply the realizations the adversary up-weights (Sec. 4.2), whereas Aolaritei et al. (2026) target the variance of a stochastic gradient. Second, we draw from a learned density model (Sec. 4.3), whereas the methods above assume a known density and Parpas et al. (2015) fit a kernel density estimate to Markov chain samples, which degrades in high dimension. Kim et al. (2025) also steer a generative model toward costly inputs with importance weights, but to train a predictor rather than to make a decision. Figure 1 summarizes the resulting loop, whose components the rest of this section develops: the master solved each round (Sec. 4.1), the target toward which scenarios are drawn (Sec. 4.2), and the two ways a batch is drawn (Sec. 4.3). 4

(a) Adaptive loop

(b) Evaluate sample

(d) Draw

Gθ

Gθ

Latent space z

(c) Tilt

Sampling step at xt Sample Sample nominal, evaluate cost Tilt Reweight by exp(Q/λ) Draw Resample, or sample rt ; weight by w

xt+1

Q Gθ

St ∪ {(ξj , wj )}

Q

Scenario space ξ

Master problem Solve equation 6 on St+1 ⇒ xt+1

Figure 1: (a) The adaptive loop. (b–d) One sampling step in a 2-D slice of the latent space (contours of φ; cost rises to the upper right): (b) an initial sample is evaluated; (c) it is tilted by exp(Q/λ), and a Gaussian proposal gt (orange) is fitted; (d) new scenarios are resampled (rings) or drawn from the proposal (filled), with weights back to the nominal (marker area). Bottom row. The same steps after decoding by Gθ : new scenarios reach beyond the record (gray), and the target moves with the decision (dotted).

4.1

A LTERNATING M ASTER AND S AMPLER

The master on the scenario set. Each round draws a batch of n scenarios at the incumbent xt with base weights back to the nominal (Sec. 4.3), adds it to the set held, and solves equation 5 on that set for xt+1 (Algorithm 1). Scenarios drawn at earlier decisions can be kept at later ones, because a batch’s weights depend only on how it was drawn. Hence, for any decision x fixed in advance, the weighted sum of the batch’s costs, before normalization, estimates n times the nominal expected cost without bias. Since equation P 5 requires weights that sum to one, the weights of all rounds are normalized together, w̄i = wi / j wj over the Kt+1 scenarios held. The master’s objective is therefore a biased estimate of its counterpart under Qθ , through this division by a random total and through the logarithm and the minimizations over ζ and λ below. With these weights, the master on the held scenarios is Kt+1

min

x∈X , ζ, λ>0

c(x) + ζ + λρ + λ log

X

w̄i exp ψζ (Q(x, ξi ))/λ .

(6)

i=1

At a decision x, the adversary’s worst case is available in closed form, pi ∝ w̄i exp(ψζ ∗ (Q(x, ξi ))/λ∗ ), where ζ ∗ and λ∗ attain the minimum at x, so that KL(p ∥ w̄) = ρ. When equation 6 is too large to solve jointly, as in the capacity expansion study, this closed form offers a cheaper alternative, in which the master is solved with the weights fixed at the worst case Algorithm 1 Adaptive scenario acquisition. Require: Nominal Q that can be sampled; radii ρ and ρref ; initial sample size M ; batch size n; rounds T ; re-targeting rounds T ∋ 0; mixture weight a 1: Choose an initial decision x0 ; set S0 ← ∅ 2: for t = 0, . . . , T − 1 do 3: if t ∈ T then ▷ Re-target; otherwise keep the last 4: Draw ξ1 , . . . , ξM ∼ Q and evaluate Qi ← Q(x ▷ Initial sample P t , ξi ) 5: Bisect λt so that τi ∝ exp(Qi /λt ) satisfies i τi log(M τi ) = ρref ▷ Tilt 6: end if 7: if resampling then 8: Resample n indices ∝ τ ; Bt ← {(ξi , ci /(M τi ))} over distinct i ▷ Resampling 9: else 10: Fit gt to the tilted latent points of the initial sample ▷ Fitted proposal 11: Draw from rt = (1 − a) φ + a gt with weights φ/rt ; keep n draws as Bt , or resample n by the tilt as above 12: end if 13: St+1 ← St ∪ Bt 14: Solve equation 6 on St+1 for xt+1 ▷ Master 15: end for 16: return xT

5

computed at xt , which requires Kt recourse evaluations. Weak duality then bounds the resulting optimality gap on the held scenarios at the cost of one more master solve (App. B.5). Sampler and radii. The sampler, which chooses the scenarios that enter the set, tilts an initial sample of M draws from P the nominal toward the worst case at a reference radius ρref , τi ∝ exp(Q(xt , ξi )/λt ) with i τi log(M τi ) = ρref , and forms the batch from this tilt (Sec. 4.3). The two radii thus play different roles, ρref determining which scenarios enter the set and ρ how the master weights them, so that ρ = 0 with ρref > 0 solves the nominal problem with a decision-aware sampler. We take ρref ≥ ρ as a heuristic, because a sampler that reaches beyond the adversary only costs efficiency, whereas one that reaches less leaves tail scenarios out of the set, and no reweighting can restore them. The tilt is the worst case under the expectation rather than under the mean–CVaR objective, but because ψζ is nondecreasing, both worst cases rank scenarios by cost in the same way. Since the weights correct any proposal back to the nominal, this choice affects efficiency but not the quantity being estimated. Cost of a round. Each initial sample costs M recourse evaluations, which run in parallel, whereas only the batch of n enters the master, whose size grows with each scenario. To save evaluations, the initial sample need only be redrawn and re-tilted in the rounds t ∈ T , because the tilt depends on the decision only through the costs at xt and so changes little once the incumbent settles. Later batches reuse the last initial sample and proposal, and their weights remain exact, although resampling the same initial sample again does not reduce its own sampling error. 4.2

T HE S AMPLING TARGET

Latent space. The sampler works in the latent space, where the nominal is the standard Gaussian φ. A simple parametric family would match Qθ poorly in the ambient space, but because φ is known, a proposal from such a family has likelihood ratios in closed form in the latent space, without evaluating the density of Qθ . These ratios correct the expectation of any function of z, including the cost Q(x, Gθ (z)) of the decoded scenario. The whole problem can also be posed in the latent space, n  o min c(x) + max Rz∼p Q(x, Gθ (z)) , (7) x∈X

p : KL(p ∥ φ)≤ρ

in which the adversary reweights φ. This latent problem is equivalent to the original one, because decoding maps the latent ball onto the ball around Qθ for any measurable Gθ . Pushing two distributions through the same map cannot increase the divergence between them (Cover & Thomas, 2006), so the latent ball decodes into the ambient ball. Conversely, any P in the ambient ball is the image dP of φ reweighted by dQ ◦ Gθ , whose divergence from φ equals KL(P ∥ Qθ ). For a law-invariant R, θ such as the expectation and CVaR (Kuhn et al., 2025), the two problems therefore have the same value and minimizers, whether or not Gθ is invertible. The tilt as a proposal. exponential tilt

Under the expectation, the worst case of equation 7 at a fixed decision is the

π(z | x) ∝ φ(z) exp Q(x, Gθ (z))/λ , (8) with λ fixed by the radius. Used as a sampling target, the tilt weights a realization by the level of its cost, whereas the decision responds only to the part of the cost that x can change. The level is therefore a faithful proxy when the costly realizations are those the decision can mitigate, such as unserved load that added capacity would avert, but it fails when the costliest realizations are costly whatever the decision, since the tilt then concentrates the sample where the master cannot respond. The tilt also links the two roles the scenario set plays. Read as robustness, π(· | x) is the least favorable distribution, and read as estimation, it is the exponential tilting of rare-event sampling (Blanchet & Lam, 2011). For large λ, which corresponds to a small radius, it is close to a mixture of φ and the minimum-variance proposal for the expected cost, which reweights φ by the cost itself when the cost is nonnegative, and for small λ, it concentrates on the costliest scenarios, where the master’s worst case places its weight. The sampler draws from a proposal that approximates the tilt, and importance weights correct for the approximation, however rough it is. For any proposal g positive wherever φ is, draws zj ∼ g with the weight  w(z) = φ(z) g(z) (9) P 1 give an unbiased estimate n j w(zj ) Q(x, Gθ (zj )) of Eφ [Q(x, Gθ (z))]. 6

4.3

D RAWING THE BATCH

Sampling the target equation 8 directly is expensive, because it is known only through evaluations of Q, each a second-stage solve. Of the existing alternatives, guided sampling of the generator yields draws without likelihood ratios, and Markov chain Monte Carlo (Parpas et al., 2015) needs many sequential evaluations of Q per retained sample. We instead form the batch from the tilt τ of the initial sample, in one of two ways. Resampling. The simpler variant resamples n indices with probabilities τ , and a scenario drawn ci times enters with base weight wi = ci /(M τi ), which corrects the selection back to the nominal. The batch then holds only realizations in the initial sample, so its reach is governed by M . In return, resampling needs only draws, so it also applies to models available only as samplers and to the empirical measure. Fitted proposal. To reach beyond the initial sample, the second variant fits a Gaussian gt = P N (µ̂, Σ̂) to the latent points zi of the initial sample, ξi = Gθ (zi ), weighted by the tilt, µ̂ = i τi zi P and Σ̂ = i τi (zi − µ̂)(zi − µ̂)⊤ . When the M evaluations are few relative to the latent dimension, Σ̂ can be regularized, for example by shrinkage toward its diagonal or, at the extreme, by fixing Σ̂ = I so that only the mean shifts. The sampler then draws from the defensive mixture rt = (1 − a) φ + a gt with a ∈ (0, 1), decodes ξ = Gθ (z), and attaches the weight h φ(z) gt (z) i−1 1 w(z) = = (1 − a) + a ≤ , (10) rt (z) φ(z) 1−a with gt /φ in closed form. The weight is exact however gt was fitted, and since the component φ bounds it, a poor fit lowers efficiency but cannot produce an unbounded weight. The draws either enter the master directly with the weights equation 10 or form a pool from which n are resampled by the tilt as above. In that case, the uniform weights 1/M are replaced by the weights of equation 10 normalized within the pool, which makes the batch’s weighted sum only approximately unbiased (App. B.4). Choosing a variant. Resampling needs only draws of the nominal, but its weights vary with the tilt, so a concentrated tilt and a small n give a noisy correction. The fitted proposal, in contrast, reaches beyond the initial sample, with every weight below 1/(1 − a) when its draws enter the master directly, but it needs latent coordinates and a proposal family that fits the tilt, a fit that degrades with dimension. It should therefore help most when the events that drive the worst case are rare in an initial sample of size M and the tilted target is close to a single Gaussian in the latent space, and least when decoded draws cannot leave the record’s range or the map from z to a scenario is strongly nonlinear (Sec. 5.1).

5

E XPERIMENTS

We evaluate the framework in two settings. The first, a multi-item newsvendor problem and a standard testbed for DRO (Scarf, 1958; Hanasusanto et al., 2015) and data-driven optimization (Ban & Rudin, 2019), has a known demand distribution, so optimality gaps can be computed accurately, and a record size and shift under our control, which lets us isolate the regimes in which a learned nominal improves on the record. The second, a capacity expansion problem for the Texas grid on realistic data, tests whether these findings carry over to a planning problem of practical scale. Methods. All methods solve the master equation 5 on a finite scenario set and differ in where the scenarios come from, how they are chosen, and how the master hedges. Our method draws scenarios from a learned nominal, either by resampling its draws (KL-Gen) or from a proposal fitted to the tilt (KL-Gen-P), chooses them by the tilt at the current decision, and hedges with a KL ball. We group the comparators by which of these components they change. Source: KL-DRO centered at the record (KL-Rec) removes generation; KL-DRO over record points moved toward higher cost by optimal transport (KL-OT) extrapolates by distance from the record instead of by a model; alternative nominals isolate the model class; and in the newsvendor an oracle applies our sampler to the true distribution. Selection: cost-space scenario clustering (CSSC; Keutchayan et al., 2023) chooses 7

Table 1: Newsvendor: cost above the optimum (%, median over records) under the level (L), tail (T), support (S), and coincidence (C) shifts. (a) Out-of-distribution cost at each method’s selection rule. (b) In-distribution cost needed to halve SAA’s out-of-distribution cost; ·: never reached; † : lowest in under 80% of bootstrap resamples. Bold: lowest among the rows shown. (a) OOD cost at SAA’s ID budget N = 16 Method SAA KL-Rec W-DRO KL-OT

L

T

S

C

L

643 154 434 244 129 647 155 435 227 141 80 39 92 154 24 88 43 99 159 45

KL-Gen-P Log-normal 52 52 47 Flow 52 52 43 Gaussian 229 104 158 KL-Gen (flow) 229 103 152

(b) ID cost to halve SAA’s OOD cost

N = 128 T

S

N = 16 C

L

T

S

88 98 39 · · · 89 111 39 · · · 45 30 18 55 68 55 57 44 25 57 64 57

31 35 56 35 15 33 41 58 35 18 49 250 107 174 51 49 206 101 138 46

N = 128 C

S

C

· · · · · · · · 97 12 42 12 97 12 29 12

L

· · 20 ·

21 23 21 24 9† 19 19 19 23† 11 40 · 40 42 · 36 · 36 33 ·

T

· 9† · 9 · · · ·

15 19 · ·

scenarios once, by their effect on the decision. Hedge: sample average approximation (SAA) applies none, and Wasserstein DRO (W-DRO; Mohajerin Esfahani & Kuhn, 2018) hedges with a transport ball around the record. Both studies use the mean–CVaR objective with β = 0.25 and α = 0.9, and Tables 2 and 9 give the implementations. Hedging can lower out-of-distribution cost at the expense of in-distribution cost, so a comparison at a single setting would confound how much a method hedges with how well it hedges. We therefore compare the frontiers of out-of-distribution (OOD) against in-distribution (ID) cost that each method’s settings trace. 5.1

M ULTI -I TEM N EWSVENDOR WITH S UBSTITUTION

Problem. A retailer orders five products at a cost of 1 per unit before demand is observed. As recourse, unmet demand for a product can be served from the next product’s stock at a penalty of 0.5 per unit (substitution) or lost at a cost of 100 per unit, and leftover stock is salvaged at 0.25. For a single product, the order that minimizes expected cost is the quantile of demand at the critical ratio (100 − 1)/(100 − 0.25) ≈ 0.99, and the mean–CVaR objective used here raises it further, so the orders lie in the upper tail of demand, which a short record samples poorly. Log demand follows a known Gaussian mixture with a rare shock that raises demand for every product at once (App. A). A decision-maker sees a record of N ∈ {16, 32, 64, 128} periods and is scored under one of four shifted distributions: a level shift in every product’s demand, a tail shift to a heavier-tailed shock, a support shift that adds a shock on different products, and a coincidence shift, in which the record’s shocks strike one product at a time and the test distribution’s strike all products together, with unchanged marginals. Setup. For each record, we fit five learned nominals: a Gaussian copula; a log-normal; a two-component log-normal mixture, the family of the base distribution; a RealNVP flow; and a Gaussian fitted to raw demand, which captures dependence only through its covariance. KL-Gen-P, KL-Gen, and the oracle run a single round (T = 1) from the risk-neutral decision on the record. Because both the distribution the records are drawn from and the shifted distributions are known, we report each cost as a percentage above the optimum under the distribution it is scored on, and Table 1 reads the resulting frontiers in two ways: (a) the OOD cost within an ID budget of 5 points over SAA, chosen on held-out records without the shifted distribution, and (b) the ID cost needed to halve SAA’s OOD cost, which prices protection but is not a deployable rule. Record size and shift. The optimal order depends on each product’s top 1% of demand, which lies above the record’s maximum in 85% of records at N = 16. Reweighting the record cannot reach this region, so KL-Rec costs as much as SAA at every N . The methods that extrapolate beyond the record can reach this region. At N = 16, where cross-validation keeps the flow at its log-normal initialization (App. A.3), KL-Gen-P with the log-normal or the flow does so most cheaply, halving SAA’s OOD cost at about a third of the ID cost that W-DRO and KL-OT need (Table 1b). Which method leads at a fixed ID budget, however, depends on the shift (Table 1a). W-DRO leads or ties under the tail shift at every N and overtakes KL-Gen-P under the level and support shifts as N grows. 8

OOD (real) KL-OT

35

36 34

32

Hour of the worst day

Wind

Solar

0.6 0.8

0.4 0.2

KL-OT δ = 1

0.9

1

1.1

Peak net load / OOD KL-OT δ = 3

W-DRO

KL-Gen

KL-Gen-P

25 20 15 10 5

OOD OT Gen

In-distribution cost (bn$/yr) KL-Rec

1

0.6

0.6 OOD OT Gen

30

1.5

System wind (GW)

0

Lull

38

(d) Most severe scenarios

KL-Gen

Shape distance to real days

Wind (±2 SD)

40

Spatial coherence

Out-of-distribution cost (bn$/yr)

Windy

42

31

(c) Reach vs. realism

(b) Worst admitted days

(a) Cost trade-off

KL-Gen (copula)

Record days

0 0.6

1

1.5

2

Demand / record max 20 worst OOD days

1 (panels defined in App. B.6). 2

Figure 2: Capacity expansion at k = (a) OOD against ID cost (bn$/yr) of each setting, with Pareto frontiers; KL-Rec, CSSC, and SAA lie outside the range. (b) Zonal wind on the most severe OOD day and admitted scenarios (top), and their spatial coherence (bottom). (c) Severity against demand-rescaled distance to the nearest real day, for the most severe admitted scenarios. (d) Demand relative to its 2001–2010 maximum against system wind; stars: the 20 most severe OOD days.

KL-Gen-P stays lowest under the coincidence shift because a fitted nominal places mass on joint shocks, which transport reaches only by moving every product of a record point at once. The proposal and the nominal. At the same 356 cost evaluations, KL-Gen-P cuts the out-ofdistribution cost of KL-Gen by a factor of two to six under the level and support shifts for every nominal fitted to log demand other than the copula (Tables 3–6). This is because only about four of 356 draws from a nominal exceed a product’s 99th percentile, whereas the proposal draws there directly (Sec. 4.3). The copula is the exception, because its draws stay within the record’s range whatever the latent proposal, and for it the proposal gains at most a factor of two. Nor does a larger generated pool substitute for sampling toward the worst case, since CSSC, which reduces a pool of 1,024 generated draws at 95 times the evaluations, costs as much as KL-Gen. Finally, the nominal’s tail matters more than its exact family. With the proposal, the log-normal costs less out of distribution than the Gaussian in every cell, and less than the log-normal mixture, the family of the base distribution, under the level and support shifts up to N = 64 and under the coincidence shift at every N . These results support training generative models to reproduce the extremes of the process accurately (Jaini et al., 2020; Laszkiewicz et al., 2022; Pandey et al., 2025). 5.2

C APACITY E XPANSION P LANNING

Problem. A planner sizes wind, solar, gas, and battery capacity in 20 zones of the Texas system, and transmission between them, within a capital budget, then dispatches hourly, with a linearized unit commitment of gas plants, over a day of uncertain zonal demand and renewable availability (D = 1,440). Unserved load costs $10,000/MWh, so cost is dominated by the few days on which high demand coincides with low renewable output across zones, and a record of a year or less rarely holds the worst of them. Setup. A record covers k years of 2001–2010: every second day of one year (k = 21 ), a full year (k = 1), or, in App. B.7, five years (k = 5). The learned nominal is a seasonally conditioned normalizing flow, fitted to the record alone and decoded through the record’s marginal distributions (App. B.3). Every method minimizes annualized capital cost plus the mean–CVaR of operating cost under the same capital budget (B = 21.5 bn$/yr), holds at most K = 48 scenarios, and, if adaptive, re-targets every round. For the transport methods, KL-OT and W-DRO, the tables report the transport budget δ = 1 (App. B.4). Each plan is scored by the same objective, in bn$/yr, on every day of 2011–2014 (ID), the four years after the record decade, and of 2020 (OOD), the most extreme year of 2001–2020 in the tail of demand (App. B.2). Frontiers across record lengths. Adaptive sampling alone cuts cost well below that of a static set of record days. At ρ = 0, KL-Rec is far cheaper than SAA at every k (Table 10). At k = 12 and 9

the five radius pairs of Table 10, every setting of KL-Rec, and every setting of KL-OT except its cheapest in distribution, is dominated in median cost by a setting of KL-Gen (Figure 2(a)). Over the three k, KL-Gen is cheaper than KL-OT out of distribution at all 15 radius pairs, at a small in-distribution premium for k ≥ 1 (Table 11). The advantage narrows as the record grows: KL-Gen’s out-of-distribution margin over KL-OT falls from 2.3–6.5 bn$/yr at k = 12 to 0.07–0.76 at k = 5, and KL-Rec joins the Pareto front at k = 5. KL-OT reaches the front at every k only through its cheapest in-distribution setting, and SAA, CSSC, and W-DRO never reach it. KL-Gen-P and KL-Gen with a seasonal Gaussian copula stay close to KL-Gen. The copula appears to capture what these decisions need, although the flow fits held-out days better (App. B.3), and at this dimension the fitted proposal adds little over resampling (Sec. 4.3). What the generator extrapolates. The most severe days of 2020 lie within the demand range of 2001–2010 and are severe because high demand coincides with low renewable output across zones (panel d), a pattern present even in a k = 12 record. The flow reproduces it: its most severe scenarios reach the OOD days’ severity with spatial shapes as close to real days as the record’s own severe days (panel c) and wind as spatially coherent as on real days (panel b). What transport extrapolates. KL-OT makes scenarios costly by inflating demand zone by zone, up to 2.4 times the highest system demand of 2001–2010 at δ = 3 (panel d), whereas on real days each zone’s demand moves with the system total. At δ = 1 its most severe scenarios fall just short of the OOD days’ severity, and at δ = 3 they reach it only with shapes farther from real days than the record’s own severe days (panel c). Widening the transport budget brings its OOD cost toward that of KL-Gen only at a higher ID cost (panel a). This premium is paid every year to hedge against days unlike those in the data.

6

C ONCLUSION

For two-stage stochastic programs with severe tail risk, we proposed centering a KL ball at a generative model of a short record, so that the adversary judges unobserved realizations by learned dependence rather than distance from the record, and solving the resulting problem on scenarios drawn adaptively toward its worst case. Across both case studies, we found that the learned nominal helps most when the record is short and the model class extrapolates along structure the record shows, such as the co-movement of demand across products or of renewable shortfalls across zones. The advantage fades as the record grows, and transport around the record hedges more cheaply against shifts the model class does not anticipate. As for the sampler, a fitted proposal adds most when the realizations driving the worst case are rare among the model’s draws. Limitations. Both case studies have linear recourse, so integer and nonlinear convex recourse remain untested. On the theoretical side, the ball around Qθ contains a distribution P, shifted or not, only if ρ ≥ KL(P ∥ Qθ ). This divergence is set by the model’s misspecification, which the record does not reveal, and is infinite if the model assigns no mass where P has some, as nominals decoded through the record’s empirical marginals do beyond its range. We therefore choose the radius empirically, without a finite-sample guarantee. Statistically, the master’s objective is a biased estimate of its counterpart under Qθ , and we do not show that the returned plan converges: under the alternation, it is a best response whose optimality gap on the held scenarios is bounded only after the fact (App. B.5).

AI U SE S TATEMENT The authors received assistance from Claude Code in implementing the methods and baselines, producing the tables and figures from the data of the computational experiments, reviewing the literature, and editing and revising the text, and take full responsibility for the content of the paper.

R EFERENCES Liviu Aolaritei, Bart P. G. Van Parys, Henry Lam, and Michael I. Jordan. Stochastic Optimization with Optimal Importance Sampling, February 2026. URL http://arxiv.org/abs/2504. 03560. arXiv:2504.03560 [math]. 10

Gah-Yi Ban and Cynthia Rudin. The Big Data Newsvendor: Practical Insights from Machine Learning. Operations Research, 67(1):90–108, 2019. doi: 10.1287/opre.2018.1757. Güzin Bayraksan and David K. Love. Data-Driven Stochastic Programming Using Phi-Divergences. In Dionne Aleman, Aurélie Thiele, J. Cole Smith, and Harvey J. Greenberg (eds.), The Operations Research Revolution, pp. 1–19. INFORMS, September 2015. ISBN 978-0-9843378-80. doi: 10.1287/educ.2015.0134. URL http://pubsonline.informs.org/doi/10. 1287/educ.2015.0134. Aharon Ben-Tal, Dick den Hertog, Anja De Waegenaere, Bertrand Melenberg, and Gijs Rennen. Robust Solutions of Optimization Problems Affected by Uncertain Probabilities. Management Science, 59(2):341–357, February 2013. doi: 10.1287/mnsc.1120.1641. URL https: //pubsonline.informs.org/doi/10.1287/mnsc.1120.1641. John R. Birge and François Louveaux. Introduction to Stochastic Programming. Springer, second edition, 2011. doi: 10.1007/978-1-4614-0237-4. Jose Blanchet and Henry Lam. Rare event simulation techniques. In S. Jain, R. R. Creasey, J. Himmelspach, K. P. White, and M. Fu (eds.), Proceedings of the 2011 Winter Simulation Conference, pp. 146–160, 2011. doi: 10.1109/WSC.2011.6147747. Jose Blanchet, Yang Kang, Karthyek Murthy, and Fan Zhang. Data-driven optimal transport cost selection for distributionally robust optimization. In Proceedings of the Winter Simulation Conference, pp. 3740–3751, 2019. doi: 10.1109/WSC40007.2019.9004785. Younes Boulaguiem, Jakob Zscheischler, Edoardo Vignotto, Karin van der Wiel, and Sebastian Engelke. Modeling and simulating spatial extremes by combining extreme value theory with generative adversarial networks. Environmental Data Science, 1:e5, 2022. doi: 10.1017/eds.2022.4. Cameron Bracken, Travis Thurber, and Nathalie Voisin. Hourly wind and solar generation profiles at 1/8th degree resolution, 2023. URL https://zenodo.org/records/10214348. Dataset. Aron Brenner, Rahman Khorramfar, Jennifer Z. Sun, and Saurabh Amin. A deep generative learning approach for two-stage adaptive robust optimization. In International Conference on Learning Representations, 2025. Mónica F. Bugallo, Víctor Elvira, Luca Martino, David Luengo, Joaquín Míguez, and Petar M. Djurić. Adaptive importance sampling: The past, the present, and the future. IEEE Signal Processing Magazine, 34(4):60–79, 2017. doi: 10.1109/MSP.2017.2699226. Kai Chang and Themistoklis P. Sapsis. Extreme event aware (η-) learning. Nature Communications, 17:9919, 2026. doi: 10.1038/s41467-026-76811-x. arXiv:2510.19161. Yize Chen, Yishen Wang, Daniel Kirschen, and Baosen Zhang. Model-free renewable scenario generation using generative adversarial networks. IEEE Transactions on Power Systems, 33(3): 3265–3275, 2018. Xiuyuan Cheng, Yao Xie, Linglingzhi Zhu, and Yunqin Zhu. Worst-case generation via minimax optimization in Wasserstein space, December 2025. URL http://arxiv.org/abs/2512. 08176. arXiv:2512.08176 [stat]. Paul G. Constantine. Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies. Society for Industrial and Applied Mathematics, 2015. doi: 10.1137/1.9781611973860. Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. Wiley, 2nd edition, 2006. George B. Dantzig and Peter W. Glynn. Parallel processors for planning under uncertainty. Annals of Operations Research, 22:1–21, 1990. doi: 10.1007/BF02023045. Erick Delage and Yinyu Ye. Distributionally Robust Optimization Under Moment Uncertainty with Application to Data-Driven Problems. Operations Research, 58(3):595–612, 2010. doi: 10.1287/opre.1090.0741. 11

Charita Dellaporta, Patrick O’Hara, and Theodoros Damoulas. Decision making under the exponential family: Distributionally robust optimisation with Bayesian ambiguity sets. In Proceedings of the 42nd International Conference on Machine Learning, 2025. Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In International Conference on Learning Representations, 2017. URL https://openreview. net/forum?id=HkpbnH9lx. John C. Duchi, Peter W. Glynn, and Hongseok Namkoong. Statistics of robust optimization: A generalized empirical likelihood approach. Mathematics of Operations Research, 46(3):946–969, 2021. doi: 10.1287/moor.2020.1085. J. Dupačová, N. Gröwe-Kuska, and W. Römisch. Scenario reduction in stochastic programming. Mathematical Programming, 95(3):493–511, March 2003. doi: 10.1007/s10107-002-0331-0. URL http://link.springer.com/10.1007/s10107-002-0331-0. Jamie Fairbrother, Amanda Turner, and Stein W. Wallace. Problem-driven scenario generation: an analytical approach for stochastic programs with tail risk measure. Mathematical Programming, 191(1):141–182, January 2022. doi: 10.1007/s10107-019-01451-7. URL https://link. springer.com/10.1007/s10107-019-01451-7. Rui Gao and Anton Kleywegt. Distributionally Robust Stochastic Optimization with Wasserstein Distance. Mathematics of Operations Research, 48(2):603–655, 2023. doi: 10.1287/moor.2022. 1275. Paul J. Goulart and Yuwen Chen. Clarabel: An interior-point solver for conic programs with quadratic objectives, 2024. URL https://arxiv.org/abs/2405.12762. Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2026. URL https://www. gurobi.com. Grani A. Hanasusanto, Daniel Kuhn, Stein W. Wallace, and Steve Zymler. Distributionally robust multi-item newsvendor problems with multimodal demand distributions. Mathematical Programming, 152(1-2):1–32, 2015. doi: 10.1007/s10107-014-0776-y. Holger Heitsch and Werner Römisch. Scenario Reduction Algorithms in Stochastic Programming. Computational Optimization and Applications, 24(2–3):187–206, February 2003. doi: 10.1023/A: 1021805924152. URL https://doi.org/10.1023/A:1021805924152. Mike Hewitt, Janosch Ortmann, and Walter Rei. Decision-based scenario clustering for decisionmaking under uncertainty. Annals of Operations Research, 315(2):747–771, August 2022. doi: 10.1007/s10479-020-03843-x. URL https://link.springer.com/10.1007/ s10479-020-03843-x. Zhaolin Hu and L. Jeff Hong. Kullback-Leibler divergence constrained distributionally robust optimization. Optimization Online, 2013. URL https://optimization-online.org/ ?p=12225. Kjetil Høyland, Michal Kaut, and Stein W. Wallace. A Heuristic for Moment-matching Scenario Generation. Computational Optimization and Applications, 24(2–3):169–185, 2003. doi: 10.1023/ A:1021853807313. Gerd Infanger. Monte Carlo (importance) sampling within a Benders decomposition algorithm for stochastic linear programs. Annals of Operations Research, 39(1):69–95, December 1992. doi: 10.1007/BF02060936. URL http://link.springer.com/10.1007/BF02060936. David Islip, Roy H. Kwon, Sanghyeon Bae, and Woo Chang Kim. Contextual scenario generation for two-stage stochastic programming, 2025. URL https://arxiv.org/abs/2502.05349. Priyank Jaini, Ivan Kobyzev, Yaoliang Yu, and Marcus Brubaker. Tails of Lipschitz triangular flows. In Proceedings of the 37th International Conference on Machine Learning, 2020. 12

Congmei Jiang, Yongfang Mao, Yi Chai, Mingbiao Yu, and Songbing Tao. Scenario generation for wind power using improved generative adversarial networks. IEEE Access, 6:62193–62203, 2018. doi: 10.1109/ACCESS.2018.2875936. Julien Keutchayan, Janosch Ortmann, and Walter Rei. Problem-driven scenario clustering in stochastic optimization. Computational Management Science, 20(1):13, December 2023. doi: 10.1007/s10287-023-00446-2. URL https://link.springer.com/10.1007/ s10287-023-00446-2. Heasung Kim, Taekyun Lee, Hyeji Kim, and Gustavo de Veciana. Generating Informative Samples for Risk-Averse Fine-Tuning of Downstream Tasks. In Advances in Neural Information Processing Systems, 2025. Daniel Kuhn, Soroosh Shafiee, and Wolfram Wiesemann. Distributionally robust optimization. Acta Numerica, 34:579–804, 2025. doi: 10.1017/S0962492924000084. Mike Laszkiewicz, Johannes Lederer, and Asja Fischer. Marginal tail-adaptive normalizing flows. In Proceedings of the 39th International Conference on Machine Learning, 2022. Jiashuo Liu, Jiayun Wu, Bo Li, and Peng Cui. Distributionally robust optimization with data geometry. In Advances in Neural Information Processing Systems, volume 35, 2022. Casey R. McGrath, Casey D. Burleyson, Zarrar Khan, Aowabin Rahman, Travis Thurber, Chris R. Vernon, Nathalie Voisin, and Jennie S. Rice. tell: a Python package to model future total electricity loads in the United States. Journal of Open Source Software, 7(79):4472, 2022. doi: 10.21105/ joss.04472. Paul Michel, Tatsunori Hashimoto, and Graham Neubig. Modeling the second player in distributionally robust optimization. In International Conference on Learning Representations, 2021. Brian Mirletz, Laura Vimmerstedt, Tyler Stehly, Dana Stright, Stuart Cohen, Wesley Cole, Patrick Duffy, David Feldman, Parthiv Kurup, Vignesh Ramasamy, Jarett Zuboy, Gbadebo Oladosu, Jeffrey Hoffmann, Annika Eberle, Owen Roberts, Daniel Mulas Hernando, Greg Avery, Evan Rosenlieb, Anna Schleifer, Dayo Akindipe, Eric Witter, Becca Fuchs, Gabriel Zuckerman, Alex Zolan, Abdalla Abou Jaoude, Levi Larson, Chris Lohse, Nahuel Guaita, Ishita Trivedi, Fred Joseck, Ty Hedalen, Ben Rakov, and Ashok Sekar. 2024 Annual Technology Baseline (ATB) Cost and Performance Data for Electricity Generation Technologies. Technical Report 6006, DOE Open Energy Data Initiative (OEDI); National Renewable Energy Laboratory (NREL), June 2024. URL https://www.osti.gov/biblio/2377191. Peyman Mohajerin Esfahani and Daniel Kuhn. Data-driven distributionally robust optimization using the Wasserstein metric: performance guarantees and tractable reformulations. Mathematical Programming, 171(1-2):115–166, September 2018. doi: 10.1007/s10107-017-1172-1. URL http://link.springer.com/10.1007/s10107-017-1172-1. Kushagra Pandey, Jaideep Pathak, Yilun Xu, Stephan Mandt, Michael Pritchard, Arash Vahdat, and Morteza Mardani. Heavy-tailed diffusion models. In International Conference on Learning Representations, 2025. Panos Parpas, Berk Ustun, Mort Webster, and Quang Kha Tran. Importance Sampling in Stochastic Programming: A Markov Chain Monte Carlo Approach. INFORMS Journal on Computing, 27(2): 358–377, 2015. doi: 10.1287/ijoc.2014.0630. R. Tyrrell Rockafellar and Stanislav Uryasev. Optimization of conditional value-at-risk. Journal of Risk, 2(3):21–41, 2000. Reuven Y. Rubinstein and Dirk P. Kroese. The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte-Carlo Simulation, and Machine Learning. Springer, New York, 2004. Herbert Scarf. A min-max solution of an inventory problem. In Kenneth J. Arrow, Samuel Karlin, and Herbert Scarf (eds.), Studies in the Mathematical Theory of Inventory and Production, pp. 201–209. Stanford University Press, 1958. 13

Maurice Sion. On general minimax theorems. Pacific Journal of Mathematics, 8(1):171–176, 1958. Bart P. G. Van Parys, Peyman Mohajerin Esfahani, and Daniel Kuhn. From Data to Decisions: Distributionally Robust Optimization Is Optimal. Management Science, 67(6):3387–3402, 2021. doi: 10.1287/mnsc.2020.3678. Stein W. Wallace and William T. Ziemba (eds.). Applications of Stochastic Programming. MPS-SIAM Series on Optimization. SIAM, 2005. Jie Wang, Rui Gao, and Yao Xie. Sinkhorn distributionally robust optimization. arXiv preprint arXiv:2109.11926, 2021. James Watson and Chris Holmes. Approximate models and robust decisions. Statistical Science, 31 (4):465–489, 2016. Jiaqi Wen and Jianyi Yang. Distributionally Robust Optimization via Generative Ambiguity Modeling, 2026. URL https://arxiv.org/abs/2602.08976. arXiv:2602.08976. Chen Xu, Jonghyeok Lee, Xiuyuan Cheng, and Yao Xie. Flow-based Distributionally Robust Optimization, February 2024. URL http://arxiv.org/abs/2310.19253. arXiv:2310.19253 [cs]. Ziwei Zhang, Jonathan Yu-Meng Li, and Zhihao Jin. Generative distributionally robust optimization. arXiv preprint arXiv:2607.24983, 2026. Zihao Zhao, Christopher Yeh, Lingkai Kong, and Kai Wang. Diffusion-DFL: Decision-focused Diffusion Models for Stochastic Optimization. In International Conference on Learning Representations, 2026. URL https://openreview.net/forum?id=uhv3f80jmG.

A

N EWSVENDOR

A.1

P ROBLEM

A retailer chooses order quantities x ∈ R5+ at unit cost, c(x) = 1⊤ x, before demand ξ ∈ R5+ is observed. Stock of product i ships yij units to demand for product j along the arcs E = {(j, j)}5j=1 ∪ {(j + 1, j)}4j=1 , and the recourse cost is Q(x, ξ) = min 0.5 y,u,v≥0

s.t.

4 X

yj+1,j + 100

j=1

X

5 X

uj − 0.25

j=1

5 X

vi

(11)

i=1

X

yij + vi = xi ,

j:(i,j)∈E

yij + uj = ξj ,

i, j = 1, . . . , 5,

i:(i,j)∈E

where uj is unmet demand and vi leftover stock. The critical ratio is (100−1)/(100−0.25) = 0.9925, and recourse is relatively complete, since y = 0 is always feasible. The objective is c(x) plus the mean–CVaR equation 3 of Q(x, ξ) with β = 0.25 and α = 0.9. A.2

D ISTRIBUTIONS AND R ECORDS

Under the base distribution P0 , log demand is 0.95 N (µ, Σ) + 0.05 N (µ + ∆, Σ), with µj = log 100, standard deviations 0.35, and correlations 0.6. The shock ∆ is 0.9 on products 1, 3, and 5 and 0.225 on products 2 and 4. Each shift family pairs the distribution the records are drawn from with a shifted distribution (Figure 3). Level. Records from P0 ; the shifted distribution adds 0.35 to every component of µ. Tail. Records from P0 ; the shifted distribution multiplies the covariance of the shock component by 8, keeping expected demand unchanged. Support. Records from P0 ; the shifted distribution adds a third component of weight 0.05 whose shock is 1.35 on products 2 and 4 and 0.34 on the others. 14

Coincidence. Records from a mixture of a calm component of weight 0.75 and five components of weight 0.05, the jth adding ∆j to product j only; the shifted distribution is P0 , with the same marginals. There are 300 records at every N , nested across N , and the level, tail, and support families share them.

(b) Tail

(d) Coincidence

500

500

500

500

200

200

100

100

100

200

500 1000

200 100

50 50

Product 3

1000

Product 2

1000

50

100

200

500 1000

200 100

50 50

50 50

100

200

500 1000

1000

1000

500

500

500

500

200

200

200

100

100

100

50

50

50

50

100

200

500 1000

Product 1

50

100

200

500 1000

Product 1

Product 3

1000

Product 2

1000

Product 3

Shifted distribution Product 3

(c) Support

1000

Product 3

Training distribution Product 3

(a) Level 1000

50

100

50

100

200

500 1000

200

500 1000

200 100 50

50

100

200

Product 1

500 1000

Product 1

Figure 3: The four newsvendor shift families: joint densities of log demand for products 1 and 3 (1 and 2 in (c)) under the records’ distribution (top) and the shifted distribution (bottom), with contours spanning 3.5 orders of magnitude so that the shocks are visible. (a) Level: the density translates. (b) Tail: the shock spreads into a wide, faint band. (c) Support: a new component appears where the records’ distribution has almost no mass. (d) Coincidence: shocks to one product at a time become joint shocks, with unchanged marginals.

A.3

L EARNED N OMINALS

Each nominal, and a heavy-tailed variant of the flow, is fitted to all five products of a record and sampled as a map of z ∼ N (0, I5 ), so the weights φ/rt are exact whether or not the map is invertible (Sec. 4.2). Figure 4 compares them on single records. Gaussian. Sample mean and covariance of raw demand, with draws clipped at zero. Copula. The correlation of the normal scores Φ−1 ((m + 12 )/N ) of each coordinate’s rank m ∈ {0, . . . , N − 1} in the record, where Φ is the standard normal distribution function. Draws are decoded through the record’s empirical quantiles, so they stay within its range. Log-normal. A maximum-likelihood Gaussian on log demand. Log-normal mixture. Two components fitted by expectation maximization, the family of the base distribution, sampled through the inverse of its Rosenblatt transform. Flow. A RealNVP flow (Dinh et al., 2017) on whitened log demand, with six affine coupling layers of width 16, initialized at the log-normal and fitted by maximum likelihood with Adam. Its number of steps (0 to 400) and weight decay are chosen by five-fold cross-validation, which keeps zero steps on every record at N = 16. There the flow remains a log-normal, but with the unbiased covariance (divisor N − 1) rather than the maximum-likelihood one (divisor N ), so its rows differ from the log-normal’s by this 16/15 wider covariance and by independent sampler draws. Heavy-tailed flow. A tail-adaptive flow after Jaini et al. (2020): the same flow on a Student-t base with a learned number of degrees of freedom per coordinate, obtained by mapping each latent coordinate through Φ and then the Student-t quantile function. Cross-validation chooses among the log-normal, a log-t, and the t-base flow.

15

True distribution

Gaussian

Copula

Log-normal

Log-normal mixture

Flow

Heavy-tailed flow

N = 16 Product 3

1000 200 50

50

200

1000 50

200

1000 50

200

1000 50

200

1000 50

200

1000 50

200

1000 50

200

1000

N = 128 Product 3

1000 200 50

50

200 1000 50 200 1000 50 200 1000 50 200 1000 50 200 1000 50 200 1000 50 200 1000 Product 1 Product 1 Product 1 Product 1 Product 1 Product 1 Product 1

Figure 4: What each nominal can represent: 4,000 draws from each nominal fitted to one record of N = 16 (top) and one of N = 128 (bottom), shown for products 1 and 3 (black points: the record; dashed box: its range). The copula’s draws stay within the record’s range, while the other nominals extrapolate beyond it. Here, cross-validation keeps the flow at its log-normal initialization, and the heavy-tailed flow selects a log-t with 30 degrees of freedom.

A.4

M ETHODS

Table 2 summarizes the methods. Methods with a KL ball are run at every radius of the grid ρ ∈ {0, 0.003, 0.005, 0.007, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1, 0.2, 0.4, 0.8}, and methods that depend on a decision use the risk-neutral decision x̂ on the record and tilt at ρref = 0.8. KL-Rec and SAA. KL-DRO on the empirical distribution (Hu & Hong, 2013). SAA is KL-Rec at ρ = 0. KL-Gen-P and the oracle. One round of the sampler at x̂: M = 256 latent draws are evaluated, a Gaussian is fitted to their tilt, with its covariance shrunk toward its diagonal with weight 0.25, and n = 100 scenarios are drawn from the defensive mixture with a = 0.9, once per record for all radii. The oracle applies the same sampler to the true distribution and serves only as a reference. KL-Gen. Resamples 100 of 356 draws of the nominal by the tilt, with replacement, matching the 256 + 100 recourse evaluations of the sampler. CSSC. Reduces 1,024 draws of the log-normal or the log-normal mixture to K = 100 medoids weighted by cluster mass, clustering each draw’s costs under x̂ and 32 single-scenario decisions with k-medoids, a scalable variant of the program of Keutchayan et al. (2023). Transport step. W-DRO and KL-OT move each record point ξ 0 to an approximate maximizer of Q(x, ξ) − γ∥ω(ξ) − ω(ξ 0 )∥22 , the Lagrangian of a type-2 Wasserstein ball, with ω the standardized log demand, by damped fixed-point steps on the recourse duals. The multiplier γ is set at x̂ so that the mean squared displacement is ε times the median squared nearest-neighbor spacing of the record, and we call ε the displacement. Both methods alternate three rounds of moving the points and solving the master. W-DRO. Moves all N record points, with ε ∈ {0.01, 0.03, 0.1, 0.3, 1, 3, 10}, and has no KL ball. KL-OT. Moves min(100, N ) record points selected by the tilt, at the displacements of W-DRO, within the ball at every radius of the grid. Parity. All methods share the records, the evaluation distributions, the selection rule, and the master. KL-Gen-P and KL-Gen evaluate the recourse 356 times per record, and KL-Gen’s master holds a median of 69 distinct scenarios against the 100 of KL-Gen-P. CSSC evaluates it 33 times per pool draw, and W-DRO and KL-OT about 1,300 times per setting at N = 16 and 8,400 to 9,700 at N = 128. KL-OT also chooses among 91 settings, against 13 for the other KL methods and 7 for W-DRO. The differences in evaluations and in KL-OT’s grid favor the baselines. 16

Table 2: Newsvendor methods. KL-Gen-P and KL-Gen are run with each nominal of Figure 4, and CSSC with the log-normal and the log-normal mixture. Weights are normalized to sum to one. Method

Scenarios in the master

Weights

Radius

SAA KL-Rec KL-Gen-P, oracle KL-Gen

the N record points the N record points n = 100 draws of the sampler 356 draws of the nominal, resampled 100 times by the tilt K = 100 medoids of 1,024 draws of the nominal the N record points, moved min(100, N ) record points selected by the tilt, moved

uniform uniform φ/r ci /(356 τi )

0 grid grid grid

cluster mass

grid

uniform uniform

0 grid

CSSC W-DRO KL-OT

A.5

S OLVING THE M ASTER

The master is a linear program at ρ = 0, solved by Gurobi (Gurobi Optimization, LLC, 2026), and otherwise an exponential-cone program solved by Clarabel (Goulart & Chen, 2024). A Clarabel solution is accepted if its value agrees with the objective recomputed at its decision to 10−5 relative and is no worse than the solution for ρ = 0. Otherwise, a sequential least squares programming solve starts from the solution for ρ = 0. This fallback served 10% of the master solves with ρ > 0, and 22% of those of KL-Gen-P. A.6

S ELECTION AND E VALUATION

Selection rule. Records are split into odd- and even-numbered halves. On each half, each method’s setting is chosen as the largest whose mean in-distribution cost is within 5 percentage points of that of SAA, or the cheapest if none is, and it is scored on the other half, so the shifted distribution is never consulted. Settings are ordered by radius, by displacement for W-DRO, and by displacement and then radius for KL-OT. Table 1(b) instead reads the frontiers of Figure 5 at a fixed target: the lowest median in-distribution cost among a method’s settings whose median out-of-distribution cost is at most half that of SAA, with the stability of the lowest method assessed over 200 bootstrap resamples of the records. Under relative rules that admit in-distribution costs up to 1.10 or 1.25 times that of SAA, the lowest method is unchanged in 7 of the 8 columns of Table 1(a) at 1.10, where the log-normal replaces W-DRO at N = 128 under the tail shift, and in all 8 at 1.25. Scoring. In- and out-of-distribution costs are the objective of App. A.1 under the records’ distribution and the shifted distribution, as a percentage above the optimum under that distribution, computed on a scrambled Sobol quadrature with 4,096 points √ per mixture component. Tables report the median over records with the standard error 1.2533 σ̂/ n, where σ̂ is the interquartile range divided by 1.349. A.7

F ULL R ESULTS

Figure 5 and Tables 3–6 give the full results.

B

C APACITY E XPANSION

B.1

M ODEL

Variables. Let Z be the 20 zones, indexed by n, L the candidate lines with zone-line incidence matrix A, G the technologies, of which Gth are gas and Gcont = G \ Gth the others, S the scenarios, and t = 1, . . . , 24 the hours of a day. The first stage builds generation capacity xgen ng and transmission capacity xtr from zero. Gas capacity is P̄ u for a plant count u , which is relaxed to be continuous g ng ng ℓ throughout. The recourse variables of scenario s, collected in ys , are dispatch pstng , flows fstℓ , load ch dis shedding ℓshed stn , and battery charge pstn , discharge pstn , and state of charge estn . 17

(a) Level

(b) Tail

200

(c) Support

(d) Coincidence

500 200

N = 16

500 100

200 100

50 50

100

200

100

100

50

50

50 20

200

20

50

100

200

20

50

100

200

100

200 100

50

50

30

30

20

15 20 30

50

100

100

50

15 20 30

50

100

30

50

100

200

15 20 30

50

100

50

100

20 15

30 10

20

50

100

100

50

20

30

50

15

20

30

20

50

100

10

20

30

50

100 30

50

50 30

30

10 20

30

50

ID cost above optimum (%) KL-Rec KL-OT W-DRO

20

50

20

10

15 50

100 100

20 15

30 10

20 15 15 20 30

200

50

50

N = 128

100

30

50

100

200

50

100

200

N = 64

OOD cost above optimum (%)

N = 32

500

20

20 200

10

20

50

100

ID cost above optimum (%)

KL-Gen-P (copula) KL-Gen-P (log-normal)

10 10

15

20

ID cost above optimum (%)

KL-Gen-P (log-normal mixture) KL-Gen-P (flow)

30

ID cost above optimum (%)

SAA half SAA's OOD

Figure 5: Newsvendor frontiers: out-of-distribution against in-distribution cost above the optimum (%, median over records) as each method’s settings vary. Each line is a method’s Pareto frontier (lower left is better). Rows are record sizes and columns shifts. The open square is SAA, and the dashed line is half of SAA’s out-of-distribution cost, the target of Table 1(b).

18

Table 3: Newsvendor, every method, N = 16 (300 records): cost above the optimum (%) under the level (L), tail (T), support (S), and coincidence (C) shifts. (a) Median out-of-distribution cost (± 2 SE) at each method’s selection rule. ID cost at ρ = 0: median in-distribution cost of the unhedged fit (none for W-DRO and KL-OT). (b) Lowest median in-distribution cost at which a setting halves SAA’s median out-of-distribution cost; ·: never. Bold: lowest per column, excluding the oracle; † : lowest in fewer than 80% of bootstrap resamples. (a) OOD cost at SAA’s ID budget Method

L

SAA 643 ±59 KL-Rec 647 ±55 W-DRO 80 ±10 KL-OT 88 ±11 Gaussian KL-Gen-P 229 ±31 KL-Gen 480 ±43 Copula KL-Gen-P 490 ±46 KL-Gen 503 ±48 Log-normal KL-Gen-P 52 ±11 KL-Gen 238 ±37 CSSC 240 ±35 Log-normal mixture KL-Gen-P 65 ±13 KL-Gen 216 ±34 CSSC 193 ±35 Flow KL-Gen-P 52 ±9 KL-Gen 229 ±34 Heavy-tailed flow KL-Gen-P 55 ±11 KL-Gen 178 ±29 Oracle (reference)

18 ±3

ID cost

(b) ID cost to halve SAA’s OOD

C

ρ=0

L

154 ±8 434 ±39 244 ±21 155 ±7 435 ±36 227 ±20 39 ±4 92 ±10 154 ±20 43 ±4 99 ±11 159 ±20

184 184 – –

104 ±4 158 ±19 49 ±7 135 ±5 316 ±27 128 ±12

115 40 116 ·

136 ±6 338 ±30 126 ±12 137 ±6 344 ±30 133 ±13

125 127

T

S

S

C

· · · · · · 55 68 55 57 64 57

T

· · 97 97

· 40 · ·

42 ·

· ·

· ·

· ·

· ·

52 ±6 47 ±8 105 ±5 162 ±23 104 ±5 163 ±22

31 ±5 56 ±9 66 ±9

48 21 23 21 48 43 · 43 57 44 · 44

24 40 45

52 ±6 59 ±10 102 ±5 151 ±20 97 ±6 136 ±24

44 ±8 59 ±8 63 ±9

41 26 27 26 41 40 · 40 49 39 · 39

39 49 52

52 ±5 43 ±7 103 ±5 152 ±21

33 ±6 49 ±7

41 19 19 19 38 36 · 36

23† 33

39 ±4 54 ±11 96 ±5 126 ±20

49 ±8 39 ±6

28 22 22 22 36 31 · 31

25 31

21 ±2

22 ±3

3

32 ±3

Table 4: Newsvendor, every method, N = 32 (300 records); as Table 3. (a) OOD cost at SAA’s ID budget Method

L

SAA 412 ±41 KL-Rec 412 ±41 W-DRO 46 ±5 KL-OT 46 ±5 Gaussian KL-Gen-P 222 ±27 KL-Gen 414 ±30 Copula KL-Gen-P 285 ±29 KL-Gen 313 ±29 Log-normal KL-Gen-P 36 ±7 KL-Gen 190 ±24 CSSC 223 ±26 Log-normal mixture KL-Gen-P 49 ±10 KL-Gen 168 ±26 CSSC 160 ±28 Flow KL-Gen-P 41 ±9 KL-Gen 195 ±28 Heavy-tailed flow KL-Gen-P 44 ±8 KL-Gen 170 ±24 Oracle (reference)

19 ±3

ID cost

(b) ID cost to halve SAA’s OOD

C

ρ=0

L

T

S

C

126 ±5 273 ±26 149 ±14 126 ±5 284 ±27 143 ±13 23 ±3 57 ±6 85 ±11 24 ±3 61 ±6 98 ±12

87 87 – –

· · 26 24

· · 43 37

· · 26 24

· · 56 56

104 ±4 159 ±17 127 ±3 278 ±19

50 ±7 119 ±9

93 90

· ·

· ·

· ·

37 ·

111 ±4 195 ±18 114 ±4 214 ±19

58 ±6 66 ±8

61 60

· ·

· ·

· ·

44 53

49 ±5 35 ±5 98 ±4 138 ±16 103 ±4 151 ±18

23 ±3 44 ±7 51 ±7

37 33 44

13† 28 ·

17 13† · 28 · ·

19 29 35

43 ±5 51 ±9 94 ±5 124 ±17 90 ±5 119 ±18

39 ±5 48 ±7 54 ±6

30 32 32

20 28 30

26 · ·

20 28 30

28 35 38

55 ±5 33 ±6 98 ±4 135 ±18

24 ±3 45 ±7

34 30

13 27

17 ·

13 27

17 29

37 ±5 46 ±7 94 ±4 120 ±16

43 ±6 34 ±5

24 28

14 17† 25 ·

14 25

20 25

21 ±2

23 ±3

3

T

S

32 ±3

19

Table 5: Newsvendor, every method, N = 64 (300 records); as Table 3. (a) OOD cost at SAA’s ID budget Method

L

SAA 245 ±30 KL-Rec 264 ±30 W-DRO 31 ±4 KL-OT 39 ±8 Gaussian KL-Gen-P 245 ±20 KL-Gen 406 ±21 Copula KL-Gen-P 157 ±18 KL-Gen 205 ±21 Log-normal KL-Gen-P 39 ±7 KL-Gen 204 ±21 CSSC 209 ±21 Log-normal mixture KL-Gen-P 41 ±9 KL-Gen 143 ±21 CSSC 127 ±21 Flow KL-Gen-P 34 ±6 KL-Gen 204 ±22 Heavy-tailed flow KL-Gen-P 35 ±7 KL-Gen 177 ±20 Oracle (reference)

19 ±3

ID cost

(b) ID cost to halve SAA’s OOD L

S

C

ρ=0

T

S

C

107 ±4 169 ±18 108 ±4 182 ±18 34 ±3 44 ±5 45 ±4 48 ±6

76 ±10 78 ±8 37 ±4 51 ±7

38 38 – –

· · · · 16 32 16 29

· · 16 16

· · 36 36

107 ±3 172 ±13 55 ±6 126 ±2 267 ±13 113 ±8

90 87

· ·

· ·

· ·

· ·

92 ±3 115 ±11 101 ±3 151 ±14

24 ±4 39 ±5

30 33

· ·

· ·

· ·

22 26

52 ±4 31 ±5 101 ±3 145 ±13 101 ±3 151 ±14

21 ±3 43 ±6 45 ±6

35 34 42

9† 20 9† · · · · · ·

17 · ·

45 ±5 42 ±6 89 ±4 109 ±13 85 ±5 99 ±13

35 ±5 42 ±6 44 ±5

22 25 26

15 22 · · 24 ·

15 · ·

23 · ·

51 ±4 31 ±4 100 ±3 144 ±14

22 ±3 47 ±6

34 35

10 17 · ·

10 ·

15 ·

43 ±5 35 ±5 96 ±3 128 ±13

28 ±4 35 ±5

25 30

12 28 · ·

12 ·

18 23

22 ±2

25 ±3

3

T

30 ±3

Table 6: Newsvendor, every method, N = 128 (300 records); as Table 3. (a) OOD cost at SAA’s ID budget Method

L

SAA 129 ±17 KL-Rec 141 ±19 W-DRO 24 ±5 KL-OT 45 ±8 Gaussian KL-Gen-P 250 ±18 KL-Gen 389 ±22 Copula KL-Gen-P 77 ±9 KL-Gen 151 ±16 Log-normal KL-Gen-P 35 ±7 KL-Gen 189 ±17 CSSC 202 ±18 Log-normal mixture KL-Gen-P 34 ±8 KL-Gen 106 ±17 CSSC 96 ±16 Flow KL-Gen-P 41 ±7 KL-Gen 206 ±20 Heavy-tailed flow KL-Gen-P 56 ±11 KL-Gen 184 ±19 Oracle (reference)

9 ±2

ID cost

(b) ID cost to halve SAA’s OOD L

S

C

ρ=0

T

S

C

88 ±3 98 ±11 89 ±3 111 ±13 45 ±3 30 ±3 57 ±4 44 ±6

39 ±5 39 ±4 18 ±3 25 ±3

15 15 – –

· · · · 12 42 12 29

· · 12 12

· · 20 ·

107 ±2 174 ±11 51 ±4 124 ±2 261 ±13 103 ±6

81 82

· ·

· ·

· ·

· ·

75 ±3 60 ±7 91 ±3 108 ±11

10 ±1 21 ±3

15 18

· ·

· ·

· ·

13 ·

56 ±4 35 ±5 99 ±3 134 ±12 101 ±3 140 ±12

15 ±2 40 ±5 41 ±5

29 34 40

9† · ·

· 9† · · · ·

15 · ·

54 ±5 82 ±4 79 ±4

35 ±5 85 ±12 81 ±11

23 ±4 39 ±5 36 ±4

13 19 17

13 26 · · · ·

13 · ·

· · ·

58 ±4 35 ±5 101 ±3 138 ±13

18 ±3 46 ±5

33 31

11 ·

· ·

9 ·

19 ·

68 ±4 45 ±7 98 ±3 128 ±12

49 ±6 51 ±6

24 30

13 32 · ·

13 ·

· ·

37 ±2

14 ±2

3

T

18 ±1

20

Objective.

Given scenario weights ps summing to one (App. B.5), the master is   X 1 X min C inv (x) + (1 − β) ps ν Csop (x, ys ) + β ζ + p s vs 1−α s∈S

s.t.

vs ≥ ν Csop (x, ys ) − ζ,

(12)

s∈S

vs ≥ 0,

∀s,

C inv (x) ≤ B,

together with the constraints below, where ν = 8760/24 annualizes a day’s operating cost, B is the capital budget below, β = 0.25, α = 0.9, and ζ and vs are the auxiliary variables of the Rockafellar– Uryasev representation of CVaR (Rockafellar & Uryasev, 2000). The investment and operating costs are X X X gen cinv cinv ctr |ℓ| xtr C inv (x) = g xng + g P̄g ung + ℓ , n, g∈Gcont

Csop (x, ys ) =

X

n, g∈Gth

shed cvar g pstng + c

t,n,g

X

ℓ∈L

ℓshed stn ,

t,n

tr var where |ℓ| is the length of line ℓ, cinv = cfuel hg + cvom + g and c are annualized, and for gas cg g CO2 CO2 vom c r hg , with variable operation and maintenance cost cg , heat rate hg , and emissions rate rCO2 .

Constraints.

For all s, t, and n, zonal balance is X ch shed pstng + pdis stn − pstn + ℓstn − dstn = (Afst )n , g

tr where dstn is the scenario’s demand in zone n at hour t, and flows are limited by −xtr ℓ ≤ fstℓ ≤ xℓ , a network-flow model without power-flow constraints. Wind and solar dispatch is bounded by the scenario’s capacity factors, pstng ≤ astng xgen ng , and gas dispatch by installed capacity, pstng ≤ P̄g ung . For batteries, gen dis 0 ≤ pch stn , pstn ≤ xn,bat ,

0 ≤ estn ≤ Hxgen n,bat ,

dis dis estn = es,t−1,n + η ch pch , stn − pstn /η

where xgen n,bat is the battery power capacity of zone n, with the last hour linked cyclically to the first. Load shedding satisfies 0 ≤ ℓshed stn ≤ max{dstn , 0}. Gas plants follow a linearized clustered unit commitment: the committed count cstng ∈ [0, ung ] of each gas technology in each zone changes by start-ups and shut-downs, output satisfies pg P̄g cstng ≤ pstng ≤ P̄g cstng , started and stopped units respect minimum up- and down-times (cyclic over the day), and each start-up adds csu g P̄g to Csop . The count cstng is continuous, so the recourse remains a linear program. There are no ramping constraints. Capital budget. B = 21.5 bn$/yr for every method, radius, record, seed, and k. It was fixed before any run, inside the range (19.5 to 23.6 bn$/yr) of the investment cost of unbudgeted plans from a preliminary configuration of KL-Rec at k = 1 (ρ = 0, ρref = 0.4). At B = 19.35 and 23.65 bn$/yr (±10 %), KL-Gen remains cheaper out of distribution than KL-Rec and KL-OT at (ρref , ρ) = (0.2, 0.05) at every k. Network and parameters. The 214 counties of the region of the Electric Reliability Council of Texas (ERCOT) are grouped into 20 zones by k-means on county locations. Dallas and Rockwall counties, which have no coordinates in the source table, take their 2020 Census internal points. Zonal demand is the sum of county demand, and zonal capacity factors are county averages. Candidate lines join zone centroids along the edges of their Gabriel graph (40 lines, 116 to 279 km, median 178 km). Figure 6 shows the zones and lines, and Table 7 gives the cost and technology parameters. B.2

DATA AND R ECORDS

Data. The source is hourly county-level ERCOT data: wind and solar capacity factors from Bracken et al. (2023) and demand simulated by the Total ELectricity Load (TELL) model (McGrath et al., 2022), aggregated to the 20 zones. A scenario is one day of 20 zonal demands and 40 zonal capacity factors, so D = 1,440, cut at local midnight. Short gaps are interpolated, and days with longer gaps are dropped. Table 8 gives the periods. 21

Zone Centroid 2% of Demand 5% of Demand 10% of Demand Figure 6: The 20 zones of the capacity expansion model. Counties of the ERCOT region are shaded by zone (other Texas counties in gray) and zone boundaries are drawn in dark gray. Lines join zone centroids along the 40 candidate transmission corridors, and centroid markers are sized by the zone’s share of 2010 demand.

Table 7: Cost and technology parameters, in real 2024 dollars. Generation costs are from the 2024 Annual Technology Baseline of the National Renewable Energy Laboratory (Mirletz et al., 2024) (Moderate scenario, 2030, converted from 2022 dollars by ×1.061), and heat rates are assumed. The gas price and emissions rate are from the Energy Information Administration (Annual Energy Outlook 2024, 2030). Capital costs are annualized with the capital recovery factor CRF(r, L) = r(1 + r)L /((1 + r)L − 1) over the lifetime L, plus fixed operation and maintenance (O&M). VOM is variable O&M. Technology Wind Solar Gas, combined cycle Gas, combustion turbine Battery 4h Parameter shed

Lost load c Carbon price cCO2 Gas price cfuel Emissions rate rCO2 Transmission ctr Minimum output pg (CC, CT) Start-up cost csu g (CC, CT)

Capital ($/kW)

Fixed O&M ($/kW-yr)

VOM ($/MWh)

Heat rate (MMBtu/MWh)

Lifetime (yr)

1494 1266 1489 1348 1540

31.0 19.1 33.7 26.5 34.5

0.0 0.0 2.2 7.4 0.0

– – 6.4 9.5 –

30 30 30 30 15

Value

Parameter

Value

$10,000/MWh $150/tCO2 $3.43/MMBtu 0.05306 t/MMBtu $1,500/MW-km 40%, 30% $79, $179/MW

Discount rate r Gas units (combined cycle, turbine) Battery duration H Efficiency η ch = η dis Line lifetime, fixed O&M Minimum up/down time (CC, CT)

7% 500, 150 MW 4√h 0.85 40 yr, 2%/yr 6 h, 1 h

Table 8: Data periods. Records A to J are the years 2010, 2009, 2007, 2008, 2005, 2001, 2002, 2003, 2004, and 2006. The years 2015–2019 are used only to place 2020 among the other years. Period

Years

Days

Record In distribution (ID) Out of distribution (OOD)

2001–2010 2011–2014 2020

22

10 × 364 1,458 364

Role records and fitted nominals evaluation evaluation

Records. At k = 1 a record is one year, at k = 12 every second day of it (182 days), and at k = 5 one of five sets of five years. Every record is run with five seeds. Out-of-distribution year. 2020 was designated out of distribution before any method was scored on it, although its days later also served for development (App. B.6). It is the most extreme year of 2001–2020 in the tail of demand: its annual peak (77.4 GW) lies 2.5 standard deviations of year-to-year variation above the 2001–2019 trend. Net load is demand less the output of a fixed reference fleet of 62.6 GW of wind and 24.3 GW of solar, the mean zonal capacities of 25 KL-OT plans at k = 12 . The 99th percentile of daily peak net load is 50.4 GW in 2020, against at most 46.7 GW in any other year. B.3

L EARNED N OMINALS

KL-Gen and KL-Gen-P use a conditional normalizing flow, and KL-Gen (copula) and KL-Gen-P (copula) a Gaussian copula with a separate mean and covariance for each season. Both are fitted to the record alone and share their marginal transform and smoothing. Copula. Each coordinate of a day, in log demand and logit capacity factors, is mapped to its normal score among the record’s days, and the scores of each season define a Gaussian. A draw uses one extra latent coordinate z0 to choose the season in proportion to the record’s season counts, forms the seasonal Gaussian draw from the other coordinates of z, maps each coordinate back through the record’s empirical quantiles, and smooths the day over adjacent hours with weights (0.2, 0.6, 0.2). Flow. The flow replaces the seasonal Gaussian with a model of the 16 leading whitened principal components of the scores, a conditional RealNVP (Dinh et al., 2017) with an affine layer and two coupling layers of width 16 (2,720 parameters), conditioned on two harmonics of the day of the year and on the season. The remaining components follow a seasonal Gaussian. The flow is fitted by maximum likelihood with AdamW, with Gaussian noise of standard deviation 0.3 added to the components, for a number of steps chosen on a season-stratified 20% of the record’s days, and then refitted on all of them. All seeds of a record share one flow. In a draw, z0 selects a record day, whose calendar day and season are the covariates, and the other coordinates drive the flow and the Gaussian remainder. The day is then decoded and smoothed as for the copula. Properties and fidelity. Before smoothing, each coordinate of a draw of either nominal lies within the record’s range, although a draw can combine coordinates into a system-wide peak that no record day reaches. In a separate study on 2001–2010, the flow’s held-out log-likelihood exceeded that of a seasonal Gaussian of the copula’s form, with tuned shrinkage, in each of the ten years, by a median of 0.97 nats per day with 182 training days and 1.85 with 320 (three training seeds). Its configuration was fixed from that study before any plan was scored. The most severe scenarios KL-Gen admits are less spatially coherent in solar than real days (mean correlation of zonal solar profiles 0.82, against 0.87 for the 20 most severe days of 2020). Both nominals overstate the frequency of the decade’s tail: at k = 21 , the share of their draws whose peak net load exceeds the 2001–2010 99th percentile is about 1.7 times the decade’s own share for the copula and 2.1 times for the flow. B.4

M ETHODS

Table 9 summarizes the methods. KL-Rec, KL-Gen, KL-Gen-P, their copula variants, and KL-OT share the adaptive loop of Algorithm 1: Rounds. The initial plan solves the unbudgeted, risk-neutral problem on one scenario fixed by the seed. Each of twelve rounds evaluates M = 256 candidates at the current plan (for the proposal methods, an initial sample of M = 64, to which the proposal is fitted, and 192 draws from it), tilts them at ρref , and admits n = 4, so the master holds at most 48 scenarios. Base weights. A new scenario’s base weight, ci /(M τi ), or ci bi /τi for the proposal methods, where bi ∝ φ(zi )/rt (zi ) and τi ∝ bi exp(Qi /λt ), is fixed at admission, and the weights of all rounds are normalized together. Master. Each round recomputes the worst case at the current plan (App. B.5) and solves for the next plan. The last plan is reported. 23

Radii. KL-Rec, KL-OT, and the copula variants are run at k ≤ 1 at the eleven pairs with ρref ∈ {0.2, 0.4}, ρ ∈ {0, 0.02, 0.05, 0.1, 0.2, 0.4}, and ρ ≤ ρref . KL-Gen, KL-Gen-P, and every method at k = 5 are run at the five pairs (0.2, 0), (0.2, 0.05), (0.2, 0.2), (0.4, 0), and (0.4, 0.05). The methods differ as follows. KL-Rec. Candidates are record days drawn with replacement. KL-Gen. Candidates are draws of the flow. KL-Gen-P. Candidates are 192 draws from a proposal that shifts the mean of φ (the mean-shift case of Sec. 4.3) in a 10-dimensional active subspace (Constantine, 2015) of the cost gradients of 64 fitting draws, mixed with φ at a = 0.5. Four are resampled by the tilt. KL-Gen (copula) and KL-Gen-P (copula). KL-Gen and KL-Gen-P with the copula in place of the flow, with all else identical. The tables report KL-Gen (copula). KL-OT (δ). Candidates are 256 record days. From the second round, each is moved before evaluation to an approximate maximizer of Q(xt , ξ) − γ∥ω(ξ) − ω(ξ 0 )∥22 , the Lagrangian of a type-2 Wasserstein ball, where ξ 0 is the record day and ω its standardized log demand and logit capacity factors, by three damped fixed-point steps on the dispatch subgradient. A moved day keeps its record day’s base weight. The multiplier γ is set once, at the second round, so that the mean cost increase of moved days is δ times the mean cost of the unmoved candidates. The transport budget δ ∈ {0.3, 1, 3} is thus a cost-increase target at one plan rather than a radius. W-DRO (δ). Holds 48 record days with uniform weights and no KL ball. Each round moves every day by the step of KL-OT, with γ set as for KL-OT, to the target δ, at the first plan, and re-solves the master on each day’s largest penalized cost, for at most four rounds. At δ = 1, 68 of 100 runs reach the four-round cap, with a median relative gap of 2 to 4 × 10−3 across k. CSSC. Selects 48 days from 256 record days (all 182 at k = 12 ) by the clustering program of Keutchayan et al. (2023), with opportunity costs C inv (x∗i ) + ν Q(x∗i , ξj ) for the plan x∗i that solves the problem on day i alone, and weights them by cluster mass. It has no KL ball. SAA. Solves the master on 48 record days with uniform weights. Table 9: Capacity expansion methods. All use the objective in equation 12 with β = 0.25, α = 0.9, and the budget B. Method

Candidates per round

Scenarios in the master

Weights

KL-Rec

256 record days

KL worst case at ρ

KL-Gen KL-Gen-P

256 draws of the flow 192 draws from a fitted proposal 256 copula draws 192 draws from a proposal fitted to the copula 256 record days, moved from the second round a moved copy of each day, from the second round

at most 48, four per round same same

same same

same same

same same

same

same

48 record days, each valued at its largest penalized cost over the day and its moved copies 48 selected

1/48

48 record days

1/48

KL-Gen (copula) KL-Gen-P (copula) KL-OT (δ) W-DRO (δ)

CSSC SAA

B.5

256 record days (all 182 at k = 1 ) 2 none

cluster masses

S OLVING THE M ASTER

The worst case. For the mean–CVaR risk measure, with Qi := Q(x, ξi ) and normalized base P weightsP w̄, the adversary maximizes minζ {ζ + i pi ψζ (Qi )} over the ball Bρ of probability vectors p with i pi log(pi /w̄i ) ≤ ρ. The objective is linear in p and convex in ζ and the feasible set is convex and compact, so the maximum and minimum may be exchanged (Sion, 1958). For fixed ζ the 24

maximizer is the tilt pi ∝ w̄i exp(ψζ (Qi )/λ), with λ bisected so that the KL constraint binds, and ζ is found by golden-section search over [mini Qi , maxi Qi ]. The alternation. Each round solves the linear program equation 12 with the weights fixed at the worst case computed at the current plan, so the reported plan xT is a best response to the worst case at xT −1 rather than a solution of the KL-DRO problem on the final scenario set. Weak duality bounds its optimality gap. With J(x, p) the value of equation 12 at plan x and weights p, minimized over the recourse and auxiliary variables, pT the final weights, and p′ the worst case at xT , the gap is at most maxp∈Bρ J(xT , p) − max{minx J(x, pT ), minx J(x, p′ )}, which costs one more master solve. Over 200 runs of KL-Rec, KL-OT, KL-Gen, and KL-Gen (copula) at (ρref , ρ) = (0.2, 0.05) and k = 12 , the bound is small for a typical run (median 0.02 to 0.21 bn$/yr by method, under 0.7% of the objective) but not for every run (90th percentile 0.4 to 1.1 bn$/yr, and above 0.5 bn$/yr in 36 runs). Individual plans are therefore not settled, and one further step moves a run’s OOD cost by more than 0.5 bn$/yr in 100 of the 200 runs. These moves are small against the variation over seeds (median 0.5 against a seed standard deviation of 4.1 bn$/yr), and they leave the comparisons intact. Re-scored at the further step, every comparison of KL-Gen with KL-Rec, KL-OT, and KL-Gen (copula) at this radius pair keeps its sign, in and out of distribution, and KL-Gen’s out-of-distribution advantage is 15.4 bn$/yr over KL-Rec and 2.5 over KL-OT, against 18.3 and 2.4 at the reported plans. Solver and computation. Masters are built directly in the Gurobi Python interface and solved by the Gurobi barrier method (Gurobi Optimization, LLC, 2026) on 8 threads, and the candidates’ recourse costs are evaluated in parallel on 8 cores. B.6

E VALUATION AND S TATISTICS

Scoring. The final plan of every run is scored by the objective in equation 12 with uniform weights over all ID or OOD days. Statistics. Paired differences between methods at the same seed are summarized by their median over the five seeds, so the record is the unit of analysis: ten records at k = 12 , five (A to E) at k = 1, and five five-year blocks at k = 5. Table 11 gives the median over records with its standard error, the standard deviation over records divided by the square root of their number, and a two-sided Wilcoxon signed-rank test. Development. The radii, the risk weight β, the marginals and smoothing of the nominals, and the construction of the base weights were chosen after comparing configurations on the ID and OOD days, which therefore served for development as well as for evaluation. The flow’s architecture and training were fixed before any plan was scored (App. B.3), and it shares the copula’s marginals and smoothing. Severity. A day’s severity is its system peak net load, with net load as in App. B.2, and its peak hour is the hour of that peak. The OOD days of Figure 2 are the 20 most severe days of 2020 (42.2 to 54.3 GW), and a method’s most severe scenarios are the top 10% by severity of the scenarios it admits, 48 per run, pooled over record years and seeds at (ρref , ρ) = (0.2, 0.05). Figure 2(a). Each point is a setting’s median over records of the median over seeds, and a method’s frontier joins its non-dominated settings. KL-OT at δ = 3 was run at three of the five radius pairs, and W-DRO is a single point. Further radius pairs of KL-Rec and KL-OT are dominated by the KL-Gen points shown, except KL-OT at (0.2, 0.02), which is cheaper in distribution than every KL-Gen setting. Figure 2(b). The maps show zonal wind capacity factors, standardized per zone over the day and clipped at ±2, on the most severe day of 2020 and on each method’s most severe admitted scenario over all 50 runs. Spatial coherence is the mean correlation of the 24-hour profiles over the 190 pairs of zones, over daylight hours for solar. Boxes span the interquartile range, whiskers the 10th to 90th percentiles, and the gray band is the interquartile range over the 20 OOD days. 25

Figure 2(c). Each day is represented by its 20 zonal net loads at its peak hour, standardized per zone over 2001–2010. The horizontal axis is peak net load relative to the median OOD day. The vertical axis is the distance to the nearest real day of 2001–2010 outside the record’s year, after the scenario’s demand is scaled so that its system net load equals that day’s, so only the spatial shape counts. Small markers are medians per record, large markers medians over records, the star the median OOD day, and the dotted lines the OOD days’ median severity and the median distance of the decade’s own most severe days. A small distance marks a shape common among real days but does not prove realism. Over the full day (24 hours of demand, wind, and solar), the most severe scenarios of KL-Gen lie 0.424 from the nearest real day of 2011–2014 and 2020, against 0.367 for real days of the same severity and 0.372 for KL-OT (δ = 1). Compared with each OOD day, the nearest admitted scenarios of KL-Gen lie within 0.2 of the OOD days’ demand level, against 0.5 to 1.0 below for KL-OT and KL-Rec. With demand rescaled, the three methods’ shape distances are nearly the same. Figure 2(d). Demand is system demand at the scenario’s peak hour, relative to the largest such value over the 2001–2010 days (72.6 GW), and wind is system wind output at that hour under the reference fleet. Gray points are the 2001–2010 days; some of them, and seven scenarios of KL-OT (δ = 3), lie outside the axes. Reported settings. KL-OT and W-DRO were also run at δ ∈ {0.3, 3}, KL-OT at four radius pairs for k ≤ 1 and two for k = 5. These settings never reach the Pareto front, so the tables report δ = 1, and Figure 2 also shows KL-OT at δ = 3. B.7

F ULL R ESULTS

20 15 10 5 0

0.2

0.4

(b) Worst-case effective sample size vs. ρ

(c) Capacity added vs. ρ Capacity added vs. ρ = 0 (GW)

Admitted above record q99 (%)

(a) Admitted tail share vs. ρref

Worst-case ESS (% of scenarios)

Tables 10 and 11 give the full results, in which the copula ablation is the two KL-Gen (copula) columns of Table 10 and the two KL-Gen − KL-Gen (copula) columns of Table 11. Figure 7 shows that the two radii act as designed, and Figure 8 summarizes KL-Gen’s advantage across record lengths.

ρref (reference radius; ρ = 0)

100

80

60

40 0

.05

.2

ρ (objective radius)

KL-Rec

KL-OT δ = 1

KL-Gen (copula)

10

5

0

0

.05

.2

ρ (objective radius)

KL-Gen

Figure 7: The two radii act as designed (capacity expansion, k = 12 ). Boxes span the interquartile range over the

ten records (each the median over five seeds; in (a), seeds are pooled), and whiskers the full range. (a) Share of admitted scenarios above the record’s 99th percentile of peak net load, at ρ = 0 (dotted: the record’s own 1 %). Raising ρref admits more severe scenarios, in 37 of 40 (method, record) pairs. (b) Effective sample size of the worst-case weights, as a share of the admitted scenarios, at ρref = 0.2. It falls as ρ grows, and the realized divergence equals ρ in every run. (c) Capacity added relative to ρ = 0 (GW) at ρref = 0.2, paired within record and seed. Raising ρ moves capital toward the budget and changes the mix.

Paired comparisons. Table 11 compares methods at the same radius pair (ρref , ρ). With five records at k = 1 and five blocks at k = 5, no difference can reach p < 0.05, so we report the number of records in which the first method is cheaper. KL-Gen vs. KL-Rec. KL-Gen’s out-of-distribution cost is lower at every radius pair: • at k = 12 , by 12.7–22.2 bn$/yr, in 10 of 10 records (p < 0.01 at all five pairs); • at k = 1, by 8.4–14.5 bn$/yr, in 5 of 5 records; • at k = 5, by 3.8–13.8 bn$/yr, in 5 of 5 blocks. 26

(a) KL-Gen − KL-Rec

(b) KL-Gen − KL-OT (δ = 1)

−10 −20 −30 −40 −50 −60

KL-Gen − KL-OT (bn$/yr)

5

KL-Gen − KL-OT (bn$/yr)

KL-Gen − KL-Rec (bn$/yr)

0

(c) KL-Gen − KL-OT (δ = 3)

0 −5 −10 −15 −20

p=.03 p=.002

p=.62 p=.06

p=.12 p=.06

p=.85 p=.01

p=.06 p=.62

p=.31 p=.19

½

1

5

½

1

5

Record length k (years)

5

0

−5

−10

p=.01 p=.28

p=.06 p=.06

p=.06 p=.06

½

1

5

Record length k (years)

Out of distribution (2020)

Record length k (years)

In distribution (2011–2014)

Figure 8: KL-Gen’s advantage across record lengths (capacity expansion). Paired cost difference (bn$/yr; negative: KL-Gen is cheaper) against (a) KL-Rec, (b) KL-OT (δ = 1), and (c) KL-OT (δ = 3), out of distribution (filled) and in distribution (hatched). Each record contributes the median, over the radius pairs of Table 10 at which both methods were run (all five; for δ = 3, three at k ≤ 1 and two at k = 5), of its median over seeds. Boxes span the interquartile range over records and whiskers the full range; the numbers are two-sided Wilcoxon signed-rank p-values over records (with five records, the smallest attainable is 0.0625). The out-of-distribution gap to KL-Rec and to KL-OT (δ = 1) shrinks as the record grows; KL-OT (δ = 3) narrows it at k = 21 at a higher in-distribution cost.

The gap shrinks as the record grows in the median over radius pairs (Figure 8), though not at every pair. In distribution, KL-Gen is cheaper at k = 12 (by 0.46–2.92 bn$/yr), mixed at k = 1, and up to 0.73 bn$/yr more expensive at k = 5. KL-Gen vs. KL-OT. KL-Gen’s out-of-distribution cost is lower at all 15 radius pairs: • at k = 12 , by 2.3–6.5 bn$/yr, in 7–9 of 10 records (p < 0.05 at four pairs); • at k = 1, by 0.4–1.9 bn$/yr, in 3–5 of 5 records; • at k = 5, by 0.07–0.76 bn$/yr, in 3–4 of 5 blocks. In distribution the two are indistinguishable at k = 21 . At k ≥ 1, KL-Gen costs up to 0.45 bn$/yr more, and is cheaper in at most 3 of 5 records. At k ≥ 1 and the same radii, KL-Gen therefore trades a small in-distribution premium for lower out-of-distribution cost. Ablations. The proposal variant stays within 0.65 bn$/yr of KL-Gen on both axes, and none of its ten differences at k = 12 is significant. Out of distribution, KL-Gen and KL-Gen (copula) are within two standard errors of each other at every radius pair and k except two pairs at k = 5, where KL-Gen is cheaper (by 0.31 and 0.84 bn$/yr). In distribution, the copula is up to 0.36 bn$/yr cheaper at k ≥ 1, beyond two standard errors at one pair at k = 1 and two at k = 5. Replacing the copula with the flow thus changes out-of-distribution cost little and adds a small in-distribution premium when the record is longer.

27

Table 10: Capacity expansion: in-distribution (2011–2014) and out-of-distribution (2020) cost (bn$/yr) at the five radius pairs (ρref , ρ) run at every k, and of SAA, CSSC, and W-DRO, which have no KL radii (median over records of the median over five seeds). Standard errors over records are 0.02–1.6 in distribution and 0.05–6.5 out of distribution for the methods with radii, and up to 6.7 and 15.6 for the others. Table 11 gives paired differences. Bold: on the Pareto front at that k. KL-OT and W-DRO use δ = 1. In-distribution ρref

ρ

Out-of-distribution

KL-Rec KL-OT KL-Gen KL-Gen-P KL-Gen (copula) KL-Rec KL-OT KL-Gen KL-Gen-P KL-Gen (copula)

k = 12 (10 record years) 0.2 0 33.44 30.35 0.05 32.67 30.43 0.2 32.44 30.73

30.34 30.89 31.14

30.71 30.65 31.00

30.39 31.09 31.08

56.96 53.62 50.19

38.42 35.97 38.86

33.91 33.07 33.19

33.35 32.90 33.17

38.79 33.51 33.66

30.29 30.73

30.53 30.88

30.48 31.00

30.34 30.66

55.71 47.46

40.57 37.35

33.58 33.08

33.60 32.97

34.30 32.88

0.4

0 0.05

32.31 31.70

SAA CSSC W-DRO

42.79 34.12 31.42

97.06 49.10 41.42

k = 1 (5 record years) 0.2 0 30.83 30.27 0.05 30.56 30.49 0.2 30.82 30.96

30.34 30.75 31.08

30.41 31.14 31.08

30.45 31.07 31.32

49.67 43.43 42.76

35.98 34.46 34.16

32.87 33.02 33.21

33.11 32.77 33.18

34.83 33.11 33.29

30.02 30.41

30.20 30.91

30.44 31.13

30.26 30.69

46.48 46.60

35.70 33.67

32.96 32.85

32.88 33.02

33.10 32.79

0.4

0 0.05

30.36 30.58

SAA CSSC W-DRO

40.99 33.65 30.69

k = 5 (5 record blocks) 0.2 0 30.11 30.18 0.05 30.73 30.53 0.2 30.83 30.71 0.4

0 0.05

29.78 29.88

29.89 30.29

83.64 45.32 36.43

30.05 30.80 31.02

30.24 30.60 30.68

29.98 30.52 31.02

42.59 46.28 44.18

34.77 33.26 34.69

33.54 32.72 33.27

32.65 32.81 33.13

33.97 33.16 33.58

30.37 30.87

30.19 30.92

30.03 30.54

39.75 37.09

35.24 33.12

32.76 32.61

33.05 32.62

33.34 33.04

SAA CSSC W-DRO

42.97 32.73 34.78

88.00 42.12 38.51

Table 11: Capacity expansion: paired cost differences (bn$/yr; negative: the first method is cheaper) at each radius pair of Table 10: the median over records (standard error) of each record’s median over seeds. ∗ /∗∗ : two-sided Wilcoxon signed-rank p < 0.05/0.01 over records. With five records (k = 1) or blocks (k = 5), the smallest attainable p is 0.0625. In-distribution

ρref

ρ

KL-Gen − KL-Rec

KL-Gen − KL-OT

KL-Gen-P − KL-Gen

Out-of-distribution KL-Gen − KL-Gen (copula)

KL-Gen − KL-Rec

KL-Gen − KL-OT

KL-Gen-P − KL-Gen

KL-Gen − KL-Gen (copula)

k = 12 (10 record years) 0.2 0 −2.92 (1.02)∗ +0.14 (0.30) −0.07 (0.12) +0.30 (0.15) −22.21 (4.64)∗∗ −2.65 (2.22)∗ −0.27 (0.75) −1.44 (0.87) 0.05 −1.72 (1.06) +0.28 (0.30) −0.05 (0.14) +0.19 (0.16) −18.29 (5.13)∗∗ −2.41 (2.85) −0.07 (0.30) −0.55 (0.90) 0.2 −1.19 (0.78)∗ +0.31 (0.15) +0.20 (0.15) −0.06 (0.25) −14.95 (3.86)∗∗ −2.28 (1.48)∗∗ −0.10 (0.58) +0.05 (0.58) 0.4

0 −1.93 (1.14)∗ +0.12 (0.45) +0.04 (0.18) +0.16 (0.15) −17.66 (4.45)∗∗ −6.51 (2.96)∗ −0.13 (1.36) −0.03 (1.10) 0.05 −0.46 (1.18) +0.16 (0.24) +0.10 (0.14) −0.02 (0.15) −12.71 (5.57)∗∗ −4.09 (1.46)∗∗ +0.01 (0.66) +0.30 (1.33)

k = 1 (5 record years) 0.2 0 −0.47 (0.34) +0.30 (0.50) −0.16 (0.05) +0.19 (0.04) −14.51 (1.34) 0.05 +0.45 (0.26) +0.45 (0.23) −0.19 (0.29) +0.24 (0.20) −8.42 (2.78) 0.2 +0.16 (0.33) +0.17 (0.12) −0.00 (0.28) +0.06 (0.13) −9.23 (1.25)

−1.93 (1.93) −0.03 (1.47) −0.16 (1.33) −1.46 (3.47) −0.65 (3.39) +0.16 (2.22) −0.46 (1.87) −0.02 (1.59) −0.03 (1.25)

0 −0.43 (0.33) +0.14 (0.25) +0.33 (0.24) −0.07 (0.21) −12.63 (2.71) 0.05 −0.23 (0.59) +0.42 (0.18) +0.20 (0.08) +0.22 (0.18) −14.31 (1.87)

−1.33 (1.18) −0.29 (0.66) +0.00 (0.57) −0.39 (0.85) −0.11 (0.30) −0.00 (0.43)

k = 5 (5 record blocks) 0.2 0 +0.25 (0.25) −0.11 (0.09) +0.18 (0.04) +0.07 (0.09) −8.90 (2.50) 0.05 −0.07 (0.31) +0.09 (0.08) +0.01 (0.16) +0.31 (0.26) −13.80 (2.06) 0.2 +0.18 (0.17) +0.38 (0.36) −0.08 (0.22) −0.04 (0.19) −10.19 (1.77)

−0.07 (0.34) −0.14 (0.89) −0.79 (0.82) −0.39 (0.81) −0.28 (0.13) −0.84 (0.35) −0.33 (0.58) −0.35 (0.17) +0.17 (0.85)

0.4

0.4

0 +0.57 (0.34) +0.33 (0.10) −0.13 (0.08) +0.26 (0.07) 0.05 +0.73 (0.15) +0.24 (0.19) +0.31 (0.16) +0.36 (0.15)

28

−6.50 (1.87) −3.79 (1.97)

−0.76 (0.50) +0.37 (0.31) −0.15 (0.67) −0.18 (0.74) −0.08 (0.11) −0.31 (0.15)

Related documents

• Learning, fast and slow: a two-fold data-based model adaptation algorithm for processes under varying operating conditions#1014530
• Bayesian Optimization under Uncertainty for Training a Scale Parameter in Stochastic Models#1343443
• Non-asymptotic Convergence of Stochastic Gradient Descent in Score-based Generative Models#1233913
• tidyHEBO: Improving the Robustness of Heteroskedastic Evolutionary Bayesian Optimization#1473706
• BoLT: A Benchmark to Democratize Black-box Optimization Research for Expensive LLM Tasks#1467153