Robust optimization

This section describes the theoretical framework behind the robust optimisation features of the platform.

General formulation

Let J(x, \theta) be a cost function that depends on a design variable x \in \mathbb{R}^{n_x} and a random parameter \theta \sim \mathcal{D} with known distribution. Let G(x, \theta) be a constraint function.

The robust optimization problem reads:

\begin{aligned} & \underset{x}{\text{minimize}} & & \rho_{\mathcal{D}}[J(x, \cdot)] \\ & \text{subject to} & & \lambda_{\mathcal{D}}[G(x, \cdot)] \geq 0 \end{aligned}

where \rho_{\mathcal{D}} is a robustness measure on the objective and \lambda_{\mathcal{D}} a reliability measure on the constraint. Both are statistical functionals evaluated over the random parameter \theta.

Robustness measures

These measures quantify the performance of a design x in the presence of uncertainty. They are used as the optimisation objective.

Mean

The MeanMeasure computes the expectation of the objective:

\rho_{\mathcal{D}}[J(x, \cdot)] = \mathbb{E}_{\mathcal{D}}[J(x, \theta)] = \int J(x, \theta) \, p(\theta) \, d\theta

where p is the probability density function of \theta. This is the standard risk-neutral approach.

For discrete distributions, the integral reduces to a weighted sum:

\rho_{\mathcal{D}}[J(x, \cdot)] = \sum_{i=1}^n J(x, \theta_i) \, p_i

Variance

The VarianceMeasure computes the variance of each output component:

\rho_{\mathcal{D}}[J(x, \cdot)] = \mathbb{V}_{\mathcal{D}}[J(x, \theta)] = \mathbb{E}_{\mathcal{D}}[J(x, \theta)^2] - \mathbb{E}_{\mathcal{D}}[J(x, \theta)]^2

This measure can be used for risk-averse optimisation.

Mean-Standard deviation tradeoff

The MeanStandardDeviationTradeoffMeasure computes a convex combination of the mean and the standard deviation for each output component:

\rho_{\mathcal{D}}[J(x, \cdot)]_k = (1 - \alpha_k) \, \mu_k(x) + \alpha_k \, \sigma_k(x)

where \mu_k = \mathbb{E}[J_k], \sigma_k = \sqrt{\mathbb{V}[J_k]} and \alpha_k \in [0, 1] is a trade-off coefficient.

Special cases:

  • \alpha_k = 0 gives the pure mean (risk-neutral),

  • \alpha_k = 1 gives the pure standard deviation,

  • \alpha_k = 0.5 gives equal weight to both.

Worst-case

The WorstCaseMeasure computes the extremal value of the function over the support of the distribution:

\rho_{\mathcal{D}}[J(x, \cdot)] = \inf_{\theta \in \operatorname{Supp}(\mathcal{D})} J(x, \theta)

for minimization, or the supremum for maximisation.

For continuous distributions, this is a nested optimisation problem solved with an inner solver (default: TNC). For discrete distributions, the support is scanned exhaustively.

Quantile

The QuantileMeasure computes the \alpha-quantile of the output distribution:

\rho_{\mathcal{D}}[J(x, \cdot)] = \inf \left\{ s \in \mathbb{R} \; : \; \mathbb{P}_{\mathcal{D}}(J(x, \theta) \leq s) \geq \alpha \right\}

This requires solving F(s) = \alpha where F(s) = \mathbb{P}(J(x, \theta) \leq s) is the cumulative distribution function of the output. A Brent root-finding algorithm is used.

This measure is restricted to scalar (single-output) functions.

Reliability measures

These measures quantify the feasibility of a design x under uncertainty. They are used as constraints in the optimisation problem.

Joint chance constraint

The JointChanceMeasure computes the probability that all constraint components are satisfied simultaneously, mapped through the probit transform \Phi^{-1} (standard normal quantile function) and shifted by a target level \alpha:

\lambda_{\mathcal{D}}[G(x, \cdot)] = \Phi^{-1}\left( \mathbb{P}_{\mathcal{D}} \left( \bigcap_{k=1}^{n_g} \{ G_k(x, \theta) \geq 0 \} \right) \right) - \Phi^{-1}(\alpha)

The probit transform stretches the probability scale near 0 and 1, improving the conditioning of the measure for gradient-based optimization. The integral is:

P(x) = \int \left( \prod_{k=1}^{n_g} \mathbf{1}_{\{ G_k(x, \theta) \geq 0 \}} \right) p(\theta) \, d\theta

\lambda = \Phi^{-1}(P(x)) - \Phi^{-1}(\alpha)

When a comparison operator other than Greater is used, the sign is adjusted so that the constraint \lambda \geq 0 corresponds to the desired reliability level.

Individual chance constraint

The IndividualChanceMeasure computes the probability that each constraint component is satisfied individually, mapped through the probit transform \Phi^{-1} and shifted by per-component target levels:

\lambda_{\mathcal{D}}[G(x, \cdot)]_k = \Phi^{-1}\left( \mathbb{P}_{\mathcal{D}}( G_k(x, \theta) \geq 0 ) \right) - \Phi^{-1}(\alpha_k)

for k = 1, \dots, n_g. The constraint \lambda_{\mathcal{D}} \geq 0 is interpreted component-wise.

Aggregated measure

The AggregatedMeasure concatenates the outputs of several measures into a single vector-valued measure:

\mathcal{M}(x) = \left[ \rho_1(x), \rho_2(x), \dots, \rho_m(x) \right]

This makes it possible to combine, for instance, a robustness measure on the objective with a reliability measure on the constraint into a single measure evaluation. It is used internally by the RobustOptimizationProblem.

Robust optimization problem

The RobustOptimizationProblem class assembles the robustness and reliability measures into a mathematical program:

\begin{aligned} & \underset{x}{\text{minimize}} & & \rho_{\mathcal{D}}[J(x, \cdot)] \\ & \text{subject to} & & \lambda_{\mathcal{D}}[G(x, \cdot)] \geq 0 \end{aligned}

Either the robustness measure or the reliability measure (or both) can be omitted, in which case the corresponding component remains deterministic. Concretely:

  • If \rho is omitted, the original deterministic objective is used.

  • If \lambda is omitted, the original deterministic constraint is used.

  • If both are provided, they must share the same uncertainty distribution.

The problem may also include deterministic bound constraints via setBounds().

Discretisation of measures

MeasureFactory converts a measure with a continuous distribution into a measure with a discrete (empirical) distribution:

\mathcal{D} \;\longrightarrow\; \sum_{i=1}^N \omega_i \, \delta_{\theta_i}

where \{\theta_i\} and \{\omega_i\} are generated by a WeightedExperiment (e.g. Monte Carlo, LHS, Sobol’ sequence). The discrete measure then replaces the integral by a finite sum, making the optimisation problem tractable via standard deterministic solvers.

Sequential Monte Carlo algorithm

The SequentialMonteCarloRobustAlgorithm is the principal solver. It iteratively refines the discretisation and re-optimises until convergence.

Algorithm

\begin{array}{l} \textbf{Given: } N_0 \text{ (initial sample size)}, \; \varepsilon \text{ (convergence tolerance)} \\ \textbf{Initialize: } N \leftarrow N_0 \\ \textbf{Repeat:} \\ \qquad \Delta N \leftarrow \text{samplingSizeIncrement}(N) \\ \qquad \text{Draw } \Delta N \text{ new samples } \{\theta_i\} \text{ from } \mathcal{D} \\ \qquad \text{Append to existing sample set; } N \leftarrow N + \Delta N \\ \qquad \text{Discretise the measures via MeasureFactory} \\ \qquad \text{Solve the resulting deterministic problem from } x_{k-1} \\ \qquad \text{Let } x_k \text{ be the new optimum} \\ \qquad \text{If } \|x_k - x_{k-1}\| < \varepsilon \text{, stop} \end{array}

Key aspects:

  • Sample size: At each iteration, new samples are added to improve the accuracy of the discretisation. The increment function can be configured (default: add N new samples each iteration).

  • Multi-start: On the first iteration, an optional multi-start search (using a Sobol’ low-discrepancy sequence over the bounds) helps find a promising starting point.

  • Convergence: The algorithm stops when the design point change falls below the absolute error tolerance.

  • Solver tolerance: The inner deterministic solver tolerance is tightened as \varepsilon_k = C / \sqrt{N} where C is a configurable factor.

This sequential approach avoids committing to a fixed discretisation level a priori: the approximation improves as the optimisation progresses, leading to efficient and robust convergence.