Propellant Optimizer

Prometheus includes a gradient-based propellant formulation optimizer that searches composition space for the mixture maximising a user-defined figure of merit — typically a weighted combination of specific impulse and bulk density.

Algorithm overview 

The optimizer wraps scipy.optimize.minimize (SLSQP) with a multi-start strategy: n_starts independent SLSQP runs are launched from random feasible starting points, and the best result across all starts is returned. Gradients are estimated by forward finite differences with step size fd_step (default 1 × 10⁻⁴ in mass-fraction units).

Objective 

The optimizer maximises the log-scale figure of merit

\[\ln \text{FoM} = \ln I_{sp} + n \cdot \ln \rho\]

where \(I_{sp}\) is the selected Isp variant (actual, vacuum, or sea-level), \(\rho\) is the bulk mixture density, and \(n\) is the user-supplied density exponent (typically 0 for pure Isp, 0.25 for volumetric propulsion index, or 1.0 for density-weighted Isp).

Constraints 

All constraints are linear and are enforced directly by SLSQP:

Mass balance (equality): \(\sum x_i = 1\).
Ratio locks (equality, fixed-proportion groups): internal ingredient ratios are held constant while their combined mass fraction varies.
Sum bounds (inequality, sum-to-total groups): a sub-set of ingredients is constrained to lie within a specified total range.
Box bounds: each ingredient has independent lower and upper bounds.

Convergence tolerance 

SLSQP stops when the absolute change in the objective falls below ftol (default 1 × 10⁻⁴). Tighter values (e.g. 1 × 10⁻⁶) improve accuracy at the cost of more iterations; looser values (e.g. 1 × 10⁻³ or 1 × 10⁻²) exit earlier and are useful for broad exploratory sweeps or when the objective landscape is flat near the optimum.

For high-fidelity final runs after an exploratory staged pass, tightening to 1 × 10⁻⁵ or 1 × 10⁻⁶ is recommended.

Expansion modes 

Three nozzle expansion modes are available, controlling both the Isp evaluation and (in staged mode) the overall optimization strategy.

Frozen 

The nozzle expansion uses frozen composition: the product mixture is fixed at the chamber equilibrium state throughout expansion. This is a fast approximation — each performance evaluation requires only a single equilibrium solve.

Suitable for quick screening runs or propellant families where the frozen/ shifting difference is small.

Shifting 

The nozzle expansion uses shifting composition: the product mixture re-equilibrates at each pressure step along the nozzle. This is physically accurate (the optimizer targets the quantity that actually matters) but roughly 10× slower than frozen because each performance evaluation requires many sequential equilibrium solves along the expansion path.

Shifting (Staged)

A two-stage strategy that combines the speed of frozen exploration with the accuracy of shifting refinement:

Stage 1 — Frozen exploration: n_starts random starting compositions are optimised with frozen expansion using the full iteration budget (max_iter / start). Because frozen evaluation is ~10× cheaper, this stage surveys the composition landscape at low cost and produces n_starts locally optimal frozen compositions.
Stage 2 — Shifting refinement: The top n_refine stage-1 results (ranked by frozen objective) are warm-started into SLSQP with full shifting expansion and a separate iteration budget (iter/start (stage 2)). Because the warm starts are already close to the optimum, far fewer shifting iterations are needed to converge.

Benchmark result on a 5-ingredient Bi₂O₃/APCP formulation:

The staged approach delivered 2.5× speedup with identical solution quality and a higher convergence rate (4/4 vs 8/12) because the frozen stage pre-identifies the feasible basins.

Recommended defaults for staged mode: n_starts = 12, max_iter / start = 10 (stage 1), n_refine = 4, iter/start (stage 2) = 10. Increasing n_starts beyond ~16 with the same n_refine = 4 rarely improves quality and increases stage-1 cost proportionally.

Progress display 

In staged mode the GUI shows two separate graphs: one for the frozen stage-1 traces and one for the shifting stage-2 traces. The progress bar spans the total start count across both stages, and the status label shows [stage 1] / [stage 2] tags to indicate which phase is active.

API reference 

class prometheus_equilibrium.optimization.MultiStartGradientOptimizer(problem: OptimizationProblem, objective: ObjectiveSpec, operating_point: OperatingPoint, prop_db, spec_db, solver, enabled_databases: list[str], max_atoms: int)

Bases: object

Multi-start SLSQP optimizer for constrained propellant formulations.

Launches n_starts independent SLSQP runs from random feasible starting compositions (generated via FormulationConstraintCompiler) and returns the best result found across all starts.

When n_workers != 1 the starts are distributed across a ProcessPoolExecutor. Each worker process receives the database state once (via the pool initializer) and handles many FD evaluations without repeated pickle overhead. Set n_workers=1 to force sequential execution (useful for debugging or when the pool startup cost outweighs the parallelism gain for very small problems).

Parameters:

problem – Constraint definition for composition generation.
objective – Objective settings.
operating_point – Chamber / expansion conditions.
prop_db – Loaded propellant database.
spec_db – Loaded species database.
solver – Equilibrium solver instance.
enabled_databases – Enabled thermo database labels for species selection.
max_atoms – Product-species max-atoms filter.

optimize(n_starts: int = 8, max_iter_per_start: int = 100, fd_step: float = 0.0001, ftol: float = 0.0001, seed: int | None = None, n_workers: int = 0, progress_callback: Callable[[dict], None] | None = None, should_stop: Callable[[], bool] | None = None, warm_starts: list[ndarray] | None = None) → OptimizationResult

Run multi-start SLSQP and return the best result.

Parameters:

n_starts – Number of independent random starting points.
max_iter_per_start – Maximum SLSQP iterations per start.
fd_step – Finite-difference step size for gradient estimation (in mass-fraction units, i.e. 0–1 scale).
seed – Optional random seed for reproducible starts.
n_workers – Number of parallel worker processes. 0 (default) means automatic: min(n_starts, os.cpu_count()). 1 forces sequential execution.
progress_callback – Optional callback receiving per-start progress dicts with keys start, converged, objective_value, start_trace, best_value, isp, n_starts.
should_stop – Optional callback; if it returns True the run stops after the current start completes (sequential) or cancels pending futures (parallel).
warm_starts – Optional list of pre-computed starting vectors (0–1 mass fractions, aligned with problem.variables). When provided, n_starts is ignored and these are used directly instead of random sampling.

Returns:

OptimizationResult summarising the best start.

Raises:

RuntimeError – If no start produces a feasible result.

class prometheus_equilibrium.optimization.StagedGradientOptimizer(problem: OptimizationProblem, objective: ObjectiveSpec, operating_point: OperatingPoint, prop_db, spec_db, solver, enabled_databases: list[str], max_atoms: int, n_refine: int = 4, max_iter_stage2: int = 20)

Bases: object

Two-stage optimizer: frozen exploration followed by shifting refinement.

Stage 1 runs n_starts random starts using frozen nozzle expansion (fast, ~10× cheaper per evaluation). The top n_refine final compositions from stage 1 are then warm-started into stage 2, which re-optimises with full shifting expansion to reach the true optimum without paying the shifting cost for every exploratory evaluation.

The optimize() method signature is intentionally identical to MultiStartGradientOptimizer.optimize() so that the two classes are drop-in substitutes in the GUI worker and headless runner. The staged-specific parameters (n_refine, max_iter_stage2) are captured at construction time.

When operating_point.shifting is False this class behaves identically to MultiStartGradientOptimizer (staged mode provides no benefit over a single frozen pass).

Parameters:

problem – Constraint definition for composition generation.
objective – Objective settings.
operating_point – Target chamber / expansion conditions (shifting=True for the stage-2 evaluation).
prop_db – Loaded propellant database.
spec_db – Loaded species database.
solver – Equilibrium solver instance.
enabled_databases – Enabled thermo database labels for species selection.
max_atoms – Product-species max-atoms filter.
n_refine – Number of stage-1 optima to carry forward into stage 2. Clamped to the number of converged stage-1 starts.
max_iter_stage2 – Maximum SLSQP iterations for each stage-2 refinement start.

optimize(n_starts: int = 8, max_iter_per_start: int = 30, fd_step: float = 0.0001, ftol: float = 0.0001, seed: int | None = None, n_workers: int = 0, progress_callback: Callable[[dict], None] | None = None, should_stop: Callable[[], bool] | None = None, warm_starts: list[ndarray] | None = None) → OptimizationResult

Run staged optimization: frozen stage 1 then shifting stage 2.

max_iter_per_start is used as the iteration budget for stage 1. Stage 2 uses max_iter_stage2 supplied at construction. When operating_point.shifting is False, stage 2 is skipped and the stage-1 result is returned directly (frozen-only fallback).

Parameters:

n_starts – Number of random starting points for stage 1.
max_iter_per_start – SLSQP iteration budget per start in stage 1.
fd_step – Finite-difference step size (0–1 mass-fraction scale).
seed – Optional random seed for reproducible starts.
n_workers – Number of parallel worker processes (0 = auto, 1 = seq).
progress_callback – Optional progress callback. Payloads include a stage key (1 or 2) and n_starts reflects the combined total across both stages.
should_stop – Optional cancellation callback checked between stages and after each parallel future.
warm_starts – Ignored (accepted for API compatibility).

Returns:

OptimizationResult from stage 2 (shifting), or stage 1 (frozen) if shifting is disabled or stage 1 is cancelled before completion.

Raises:

RuntimeError – If stage 1 produces no converged starts.

class prometheus_equilibrium.optimization.OptimizationProblem(variables: list[~prometheus_equilibrium.optimization.problem.VariableBound], fixed_proportion_groups: list[~prometheus_equilibrium.optimization.problem.FixedProportionGroup] = <factory>, sum_to_total_groups: list[~prometheus_equilibrium.optimization.problem.SumToTotalGroup] = <factory>, total_mass_fraction: float = 1.0)

Bases: object

Constraint and normalization definition for a formulation search.

The compiler uses a two-level hierarchical sampling strategy:

Level 1 – allocate the total mass fraction among compositional units (fixed-proportion groups, sum-to-total groups, and standalone ingredients) using budget-aware sequential sampling. Mass balance is satisfied by construction; no closure ingredient is required.
Level 2 – distribute each unit’s allocated total among its members. Fixed groups use the ratio lock; sum groups use sequential sampling with the last member absorbing the group residual.

Each ingredient must belong to at most one group (fixed or sum).

Parameters:

variables – Per-ingredient bounds (0–1 fractions).
fixed_proportion_groups – Ratio-locked ingredient groups.
sum_to_total_groups – Groups constrained to a bounded total mass fraction.
total_mass_fraction – Required total mass-fraction sum (default 1.0).

variables: list[VariableBound]

fixed_proportion_groups: list[FixedProportionGroup]

sum_to_total_groups: list[SumToTotalGroup]

total_mass_fraction: float = 1.0

validate() → None

Validate problem consistency.

Raises:: ValueError – If bounds or group definitions are inconsistent.

class prometheus_equilibrium.optimization.ObjectiveSpec(isp_variant: str = 'isp_actual', rho_exponent: float = 0.0)

Bases: object

Objective settings for scalar optimization.

Parameters:

isp_variant – One of "isp_actual", "isp_vac", or "isp_sl".
rho_exponent – Exponent n in Isp * rho**n constrained to [0, 1].

isp_variant: str = 'isp_actual'

rho_exponent: float = 0.0

class prometheus_equilibrium.optimization.OperatingPoint(chamber_pressure_pa: float, expansion_type: str, expansion_value: float, ambient_pressure_pa: float, shifting: bool = True)

Bases: object

Performance evaluation operating point.

Parameters:

chamber_pressure_pa – Chamber pressure in Pa.
expansion_type – Either "pressure" or "area_ratio".
expansion_value – Pe in Pa for pressure mode, or Ae/At for area-ratio mode.
ambient_pressure_pa – Ambient pressure in Pa.
shifting – Whether to score shifting (True) or frozen (False) expansion.

chamber_pressure_pa: float

expansion_type: str

expansion_value: float

ambient_pressure_pa: float

shifting: bool = True

class prometheus_equilibrium.optimization.OptimizationResult(best_objective: float, best_log_objective: float, best_composition: dict[str, float], best_isp: float, best_density: float, trial_history: list[tuple[int, float]], start_history: dict[int, list[tuple[int, float]]], completed_trials: int, pruned_trials: int, start_history_meta: dict[int, list[dict[str, object]]] = <factory>, start_compositions: dict[int, dict[str, float]] = <factory>)

Bases: object

Summary payload returned by an optimizer run.

best_objective

Best raw objective value Isp * rho**n.

Type:: float

best_log_objective

Best log-objective value used for optimization.

Type:: float

best_composition

Best mass-fraction composition.

Type:: dict[str, float]

best_isp

Best-start Isp corresponding to selected Isp variant.

Type:: float

best_density

Best-start bulk density in kg/m^3.

Type:: float

trial_history

List of (start_index, best_log_value_so_far) points.

Type:: list[tuple[int, float]]

start_history

Mapping start_index -> [(iteration_index, log_value), ...] containing objective traces for each converged start.

Type:: dict[int, list[tuple[int, float]]]

start_history_meta

Optional mapping start_index -> [point_meta, ...] with per-iteration feasibility diagnostics for each converged start.

Type:: dict[int, list[dict[str, object]]]

completed_trials

Number of starts that produced a valid result.

Type:: int

pruned_trials

Number of starts that failed (solver error or infeasible).

Type:: int

best_objective: float

best_log_objective: float

best_composition: dict[str, float]

best_isp: float

best_density: float

trial_history: list[tuple[int, float]]

start_history: dict[int, list[tuple[int, float]]]

completed_trials: int

pruned_trials: int

start_history_meta: dict[int, list[dict[str, object]]]

start_compositions: dict[int, dict[str, float]]

Configuration file format 

Optimization runs can be saved and reloaded as .prop-opt.json files. All mass-fraction values are stored as percentages (0–100).

{
  "schema_version": 1,
  "problem": {
    "variables": [
      {"ingredient_id": "AP", "minimum": 50.0, "maximum": 75.0, "pinned": false},
      {"ingredient_id": "AL", "minimum": 5.0,  "maximum": 20.0, "pinned": false},
      {"ingredient_id": "HTPB", "minimum": 8.0, "maximum": 16.0, "pinned": false},
      {"ingredient_id": "IPDI", "minimum": 1.0, "maximum": 4.0, "pinned": false}
    ],
    "fixed_proportion_groups": [
      {"group_id": "binder", "members": ["HTPB", "IPDI"], "ratios": [0.12, 0.04]}
    ],
    "sum_to_total_groups": [
      {"group_id": "solids", "members": ["AP", "AL"],
       "minimum_total": 80.0, "maximum_total": 88.0}
    ],
    "total_mass_fraction": 100.0
  },
  "objective": {"isp_variant": "isp_actual", "rho_exponent": 0.25},
  "operating_point": {
    "chamber_pressure_pa": 6894757.0,
    "expansion_type": "area_ratio",
    "expansion_value": 40.0,
    "ambient_pressure_pa": 101325.0,
    "shifting": true
  },
  "run": {
    "n_starts": 12,
    "max_iter_per_start": 10,
    "fd_step": 0.0001,
    "ftol": 0.0001,
    "n_workers": 0,
    "seed": 42
  },
  "staged": {
    "enabled": true,
    "n_refine": 4,
    "max_iter_stage2": 10
  },
  "solver": {
    "type": "gmcb",
    "enabled_databases": ["NASA-7", "NASA-9", "TERRA"],
    "max_atoms": 6
  }
}

Key fields:

run.ftol: SLSQP convergence tolerance. Omit to use the default (1 × 10⁻⁴).
staged.enabled: Set to true to activate the two-stage frozen→shifting strategy. Requires operating_point.shifting = true; ignored otherwise.
staged.n_refine: Number of stage-1 optima to carry into stage-2 refinement.
staged.max_iter_stage2: SLSQP iteration budget per start in stage 2.

Headless runner 

Run an optimization from a saved config without launching the GUI:

prometheus-optimize my_propellant.prop-opt.json
prometheus-optimize my_propellant.prop-opt.json --output results.json
prometheus-optimize my_propellant.prop-opt.json --n-starts 20 --seed 0

Results are written as JSON to stdout (or --output) and a human-readable summary is printed to stderr.