Skip to content

Cost, optimisation and surrogates

model

Cost and revenue modelling on top of simulation metrics.

A :class:CostModel is a list of terms, each turning one metric into money:

  • fixed costs fixed("rent", 5000)
  • variable costs / revenue variable("energy", "machine.busy_time", 0.4)
  • resource costs resource("teller", per_capacity_time=30)
  • waiting costs waiting("teller", per_time=2.5)
  • penalties penalty("sla", "resource.teller.wait.p95", above=10, amount=1000)
  • revenue revenue("sales", "entity.order.disposed", 12)

Costs are computed per replication from that replication's metrics, so experiments get a full distribution of cost/profit (with CIs) for free.

CostModel dataclass 

CostModel(terms: list[CostTerm] = list(), prefix: str = 'cost')

resource 

resource(resource: str, *, per_capacity_time: float = 0.0, per_busy_time: float = 0.0) -> CostModel

Staffing/ownership cost (per unit of capacity per time) and usage cost (per busy unit-time).

waiting 

waiting(resource: str, *, per_time: float) -> CostModel

Cost of customer waiting: total wait time at resource times per_time.

penalty 

penalty(name: str, metric: str, *, amount: float = 0.0, above: float | None = None, below: float | None = None, per_unit: float | None = None) -> CostModel

A charge when metric breaches a threshold: flat amount plus per_unit x excess.

calculate 

calculate(result: Any) -> CostBreakdown

Cost breakdown for a :class:SimulationResult or a metrics mapping.

metrics 

metrics(metrics: Mapping[str, float]) -> dict[str, float]

Derived metrics (cost.total, cost.revenue, cost.profit, cost.<term>).

from_dict classmethod 

from_dict(data: Mapping[str, Any] | Sequence[Mapping[str, Any]]) -> CostModel

Build from configuration: {"terms": [{"name":..., "kind":..., "metric":..., ...}]}.

objective

A thin bridge from simulation models to optimisation libraries.

:class:Objective turns a model + metric into a plain function that scipy.optimize, OR-Tools callbacks, evolutionary or Bayesian optimisation libraries can call. SimulSI's own :func:~simulsi.optimization.optimize is built on it.

from simulsi.models import mmc obj = Objective(mmc, "resource.server.wait.mean", ["servers"], replications=2) obj.bounds [(1.0, None)]

Objective 

Objective(model: Model, metric: str | None, parameters: Sequence[str], *, fixed: Mapping[str, Any] | None = None, replications: int = 5, seed: int = 0, minimize: bool = True, transform: Callable[[Mapping[str, float]], float] | None = None, invalid_value: float = math.inf)

f(x) -> float for a vector (or mapping) of decision variables.

  • x is ordered like parameters; integer parameters are rounded.
  • Every evaluation uses the same replication seeds (common random numbers), so differences between candidates are not swamped by noise. The objective is still stochastic across seeds: validate the optimum with a fresh :class:~simulsi.Experiment using a different seed.
  • transform maps the per-run metrics to a scalar, e.g. a cost model: transform=lambda m: costs.calculate(m).total_cost.
  • Invalid parameter vectors (outside bounds) return invalid_value instead of raising, which suits derivative-free optimisers.

evaluate 

evaluate(x: Sequence[float] | Mapping[str, Any]) -> Evaluation

Run the replications for x and return the (unsigned) objective value.

search

Simulation optimisation: search a parameter box, then confirm the winner.

:func:optimize searches by grid, Latin hypercube random search or Bayesian optimisation (a Gaussian-process surrogate with expected improvement), all evaluated with common random numbers through :class:Objective. Because every estimate is noisy, the best few candidates are then re-run with fresh seeds and compared with ranking and selection (:func:~simulsi.analysis.selection.select_best), so the reported winner is not just the luckiest draw.

optimize 

optimize(model: Model, metric: str | None, ranges: Mapping[str, tuple[float, float]], *, method: Method = 'bayes', budget: int = 30, replications: int = 5, minimize: bool = True, transform: Callable[[Mapping[str, float]], float] | None = None, fixed: Mapping[str, Any] | None = None, seed: int = 0, grid_levels: int = 5, n_initial: int | None = None, confirm_top: int = 3, indifference: float | None = None) -> OptimizationResult

Search ranges ({parameter: (low, high)}) for the best metric.

  • method="grid" evaluates every combination (integers: every value; floats: grid_levels points); it fails if that exceeds budget.
  • "random" evaluates a Latin hypercube of budget points.
  • "bayes" starts from a Latin hypercube of n_initial points (default max(5, 2d + 1)) and then picks each next point by expected improvement on a Gaussian-process surrogate.
  • "cem" is the cross-entropy method: sample a population from a Gaussian over the box, keep the best fifth, refit the Gaussian to them and repeat. It suits policy parameters (thresholds, trigger levels) and copes with noisy, non-smooth responses.

Each candidate is the mean of replications runs with shared seeds. With indifference set (and a metric), the best confirm_top distinct candidates are compared afresh by :func:select_best (new seeds), and its choice becomes best_parameters.

pareto

Multi-objective optimisation: the Pareto front of trade-offs.

Real decisions trade objectives against each other - cost against service level, on-time performance against passenger delay. :func:pareto_search evaluates candidate designs (a grid, a Latin hypercube over ranges, or named scenarios) with common random numbers and returns the non-dominated ones: designs for which no other design is at least as good on every objective and strictly better on one.

from simulsi.models import disruption_recovery res = pareto_search( ... disruption_recovery, ... {"cost": "min", "otp": "max"}, ... grid={"policy": ["delay", "cancel", "spares"], "spares": [1, 3]}, ... replications=2, ... ) len(res.front) >= 1 True

ParetoResult dataclass 

ParetoResult(objectives: dict[str, Sense], designs: list[Design], replications: int, notes: list[str] = list())

front property 

front: list[Design]

Non-dominated designs, sorted by the first objective.

plot 

plot(x: str | None = None, y: str | None = None, ax: Any = None) -> Any

Scatter of two objectives with the front highlighted (needs matplotlib).

dominates 

dominates(a: Sequence[float], b: Sequence[float], senses: Sequence[Sense]) -> bool

True if a is at least as good as b everywhere and strictly better somewhere.

pareto_front 

pareto_front(points: Sequence[Sequence[float]], senses: Sequence[Sense]) -> list[int]

Indices of the non-dominated points (NaN points are never on the front).

pareto_search(model: Model, objectives: Mapping[str, Sense], *, grid: Mapping[str, Sequence[Any]] | None = None, ranges: Mapping[str, tuple[float, float]] | None = None, scenarios: Mapping[str, Mapping[str, Any]] | None = None, budget: int = 30, replications: int = 5, seed: int = 0, fixed: Mapping[str, Any] | None = None, workers: int = 1) -> ParetoResult

Evaluate candidate designs and mark the Pareto-optimal ones.

Candidates come from exactly one of grid (full factorial), ranges (a Latin hypercube of budget points) or scenarios (named parameter sets). objectives maps metric names to "min" or "max". Means over replications (common random numbers) decide dominance; designs whose intervals overlap may swap places with more replications, which the notes point out.

surrogate

Gaussian-process surrogates (metamodels) of simulation outputs.

A surrogate is a cheap statistical approximation of an expensive model, fitted to a few dozen simulation runs. Use it to explore a response surface, to find promising regions before confirming them with full experiments, or inside :func:~simulsi.optimization.optimize (Bayesian optimisation).

import numpy as np X = np.linspace(0, 1, 12)[:, None] gp = GaussianProcess(seed=0).fit(X, np.sin(6 * X[:, 0])) round(float(gp.predict(np.array([[0.5]]))[0]), 2) 0.14

The GP uses a squared-exponential kernel with one length scale per input (automatic relevance determination) and a fitted noise term, so it smooths simulation noise instead of interpolating it. Always check loo_r2() (leave-one-out accuracy) before trusting predictions, and never extrapolate outside the sampled ranges.

GaussianProcess 

GaussianProcess(*, noise: bool = True, restarts: int = 4, seed: int = 0)

GP regression with an ARD squared-exponential kernel and Gaussian noise.

Inputs are rescaled to [0, 1] and outputs standardised internally; hyperparameters maximise the log marginal likelihood (L-BFGS-B with restarts random starts).

noise_std property 

noise_std: float

Estimated noise standard deviation, in output units.

predict 

predict(X: ArrayLike, *, return_std: bool = False) -> Any

Mean prediction (and the standard deviation of the mean if return_std).

loo_residuals 

loo_residuals() -> Array

Leave-one-out prediction errors, in closed form (no refitting).

loo_r2 

loo_r2() -> float

Leave-one-out R^2: share of output variance the surrogate predicts out of sample.

Surrogate dataclass 

Surrogate(parameters: list[str], ranges: dict[str, tuple[float, float]], outputs: dict[str, GaussianProcess], design: list[dict[str, Any]], data: dict[str, Array])

GP surrogates of one or more simulation outputs over a parameter box.

predict 

predict(points: Sequence[Mapping[str, Any]] | Mapping[str, Any], output: str | None = None, *, return_std: bool = False) -> Any

Predict output at one parameter mapping or a list of them.

accuracy 

accuracy() -> dict[str, float]

Leave-one-out R^2 per output (1 = perfect; below ~0.8, add runs).

latin_hypercube 

latin_hypercube(n: int, d: int, seed: int) -> Array

n points in [0, 1]^d, one per stratum in every dimension.

fit_surrogate 

fit_surrogate(model: Model | Callable[..., Any], ranges: Mapping[str, tuple[float, float]], outputs: Sequence[str], n: int = 40, *, replications: int = 1, seed: int = 0, fixed: Mapping[str, Any] | None = None) -> Surrogate

Run a Latin hypercube design of n points and fit a GP per output.

Each point averages replications runs; all points share the same seeds (common random numbers), which makes the surface smoother.