"""GEMSEO problem formulations for the EX-link sizing problem.
The quasi-static formulation is
.. math::
\\begin{cases}
l_b \\le X \\le u_b \\\\
\\min_X f(X) = (-\\eta, H, B)^T \\\\
c(X) \\le 0 \\\\
c_{eq}(X) = 0
\\end{cases}
with ``c = (mra - 10, W - 0.985, g - 0.01, 10 - d, \\gamma - 0.02)^T`` and
``c_eq = (STE - 74, \\epsilon - 16)^T``.
Five entry points cover the strategies that apply to it, in the order a study
would normally use them:
:func:`maximise_efficiency`
Single objective, all constraints, global search.
:data:`DEFAULT_GLOBAL_ALGORITHM` is differential evolution, which -- unlike
NSGA-II -- handles the equality constraints directly.
:func:`refine`
Augmented Lagrangian polish from a starting design. Worth having because
the external penalty function is accurate only for small ``r`` and badly
conditioned when ``r`` is small.
:func:`pareto_front`
Multi-objective search over all three objectives.
:func:`local_pareto`
The MOEA in a box shrunk around a design already known to be good, rather
than over the whole space -- which is what makes it work at all here.
:func:`sweep_moving_limits`
Keep efficiency as the only objective and walk a moving upper limit on
``H`` and ``B`` down, re-solving each time.
"""
from __future__ import annotations
import logging
from collections.abc import Iterable, Sequence
from dataclasses import dataclass, field, replace
from typing import Any, ClassVar
import numpy as np
from gemseo import create_scenario
from gemseo.algos.design_space import DesignSpace
from gemseo.algos.opt.factory import OptimizationLibraryFactory
from gemseo.algos.pareto.pareto_front import ParetoFront
from gemseo.core.discipline import Discipline
from gemseo.scenarios.base_scenario import BaseScenario
from gemseo.typing import StrKeyMapping
from .constants import DEFAULT_SPEC, DEFAULT_TARGETS, DesignTargets, EngineSpec
from .design import GLOBAL_BOUNDS, VARIABLE_NAMES, Bounds, Design
from .disciplines import (
COUPLED_SAMPLES,
DynamicsDiscipline,
ExlinkDiscipline,
RangeDiscipline,
StructureDiscipline,
)
from .dynamics import DEFAULT_SPEED_RPM
from .kinematics import DEFAULT_SAMPLES
from .materials import DEFAULT_MATERIAL, DEFAULT_SAFETY, Material, SafetyFactors
from .model import Analysis, analyse
from .reference import PUBLISHED_DESIGN
LOGGER = logging.getLogger(__name__)
DEFAULT_GLOBAL_ALGORITHM = "DIFFERENTIAL_EVOLUTION"
"""Global search over the full design box."""
DEFAULT_LOCAL_ALGORITHM = "Augmented_Lagrangian_order_0"
"""Local refinement, accurate on the active constraints."""
DEFAULT_MOEA = "PYMOO_NSGA2"
"""Multi-objective evolutionary algorithm."""
OBJECTIVE_NAMES: tuple[str, str, str] = ("neg_efficiency", "height", "width")
"""``f(X) = (-eta, H, B)^T``."""
INEQUALITY_OUTPUTS: tuple[str, ...] = (
"rod_angle_margin",
"compatibility_margin",
"tdc_gap_margin",
"clearance_margin",
"side_load_margin",
)
"""Discipline outputs that must stay ``<= 0``."""
EQUALITY_OUTPUTS: tuple[str, ...] = ("stroke_error", "compression_ratio_error")
"""Discipline outputs that must equal 0."""
DEFAULT_EQUALITY_TOLERANCE: dict[str, float] = {
"stroke_error": 0.05,
"compression_ratio_error": 0.05,
}
"""Half-widths used when relaxing the equalities into pairs of inequalities.
An algorithm such as NSGA-II cannot take equality constraints. Relaxing
``STE = 74`` to ``|STE - 74| <= 0.05 mm`` and ``epsilon = 16`` to
``|epsilon - 16| <= 0.05`` keeps the design within manufacturing tolerance of
the specification -- the same pragmatism already applied to ``g``, which is
constrained to a tolerance rather than to zero.
"""
MOEA_EQUALITY_TOLERANCE: dict[str, float] = {
"stroke_error": 1.0,
"compression_ratio_error": 0.5,
}
"""Looser half-widths for the equalities during a multi-objective run.
See :func:`moea_targets` for why the multi-objective stage needs a relaxed
problem at all.
"""
MOEA_TDC_GAP: float = 1.0
"""Relaxed bound on ``g`` during a multi-objective run [mm]."""
[docs]
def moea_targets(
targets: DesignTargets = DEFAULT_TARGETS, tdc_gap: float = MOEA_TDC_GAP
) -> DesignTargets:
"""Relax ``g`` for a population-based run.
Why this is necessary is worth stating plainly, because it is the same wall
a population method hits on this problem: it returns solutions no better,
and often worse, than a local search from the same starting point.
Sampling 2000 designs uniformly from a box around a known-good solution and
testing each constraint on its own gives:
================== ============================
constraint satisfied by
================== ============================
``d >= 10`` 96 % of analysable designs
``W <= 0.985`` 76 %
``mra <= 10`` 28 %
``|eps - 16|`` 8 %
``gamma <= 0.02`` 6 %
``|STE - 74|`` 4 %
``g <= 0.01`` **0.1 %**
================== ============================
``g``, the gap between the two top dead centres, is nominally an inequality
but at 0.01 mm it is a *third equality in disguise*. With three equalities
in eleven variables the feasible set is a thin sheet, and the joint hit rate
for a random population is of order 1e-7. NSGA-II duly returns a "front" of
one point, because one point is all it ever found that was feasible.
So the multi-objective stage runs on a relaxed problem: it is a scouting
device that maps the *shape* of the efficiency-versus-size trade-off. The
design finally quoted comes from :func:`refine`, which drives ``g``, ``STE``
and ``epsilon`` back onto their true targets -- exactly the two-stage
two-stage workflow::
front = pareto_front()
final = refine(front.design)
On this problem :func:`sweep_moving_limits` is the more reliable way to trace
the trade-off, since every one of its steps is a warm-started local solve on
the *unrelaxed* problem.
Args:
targets: The targets to relax.
tdc_gap: The relaxed bound on ``g`` [mm].
Returns:
A copy of ``targets`` with ``max_tdc_gap`` widened.
"""
return replace(targets, max_tdc_gap=tdc_gap)
[docs]
@dataclass
class Outcome:
"""Result of one optimization run."""
design: Design
"""The best design found (for multi-objective runs, the highest-efficiency
point of the front)."""
analysis: Analysis
"""Its full analysis."""
algorithm: str
scenario: BaseScenario = field(repr=False)
front: list[Design] = field(default_factory=list, repr=False)
"""Pareto-optimal designs; empty for single-objective runs."""
@property
def feasible(self) -> bool:
"""Whether the returned design satisfies every constraint."""
return is_feasible(self.analysis)
@property
def n_evaluations(self) -> int:
"""Design points recorded in the top-level problem database.
For the augmented Lagrangian this counts *outer* iterations only: each
one solves a sub-problem with its own database, so the true number of
analyses is larger.
"""
return len(self.scenario.formulation.optimization_problem.database)
[docs]
def summary(self) -> str:
"""A short human-readable report of the outcome."""
return format_analysis(self.analysis, title=f"{self.algorithm} result")
INEQUALITY_TOLERANCE = 1.0e-4
"""Slack allowed on the inequality constraints when judging feasibility.
GEMSEO's own default (``ineq_tolerance``), and it matters here: an augmented
Lagrangian drives the active constraints *onto* their bounds, so a converged
design routinely lands a few parts in 1e5 outside one of them -- ``g``, in
practice, which the optimizer pushes hard against 0.01 mm. Demanding ``<= 0``
exactly would label designs infeasible that every solver in the chain considers
converged.
"""
EQUALITY_TOLERANCE = 0.05
"""Half-width allowed on ``STE - 74`` [mm] and ``epsilon - 16``.
Looser than GEMSEO's 0.01 default because both are physical dimensions: 0.05 mm
on a 74 mm stroke is well inside what the linkage could be manufactured to.
"""
[docs]
def is_feasible(
analysis: Analysis,
targets: DesignTargets = DEFAULT_TARGETS,
tolerance: float = EQUALITY_TOLERANCE,
inequality_tolerance: float = INEQUALITY_TOLERANCE,
) -> bool:
"""Whether an analysed design satisfies every constraint.
Args:
analysis: The analysis to check.
targets: Constraint right-hand sides.
tolerance: Half-width allowed on the two equality constraints.
inequality_tolerance: Slack allowed on the inequality constraints; see
:data:`INEQUALITY_TOLERANCE`.
Returns:
``True`` if the design was analysable and meets every constraint within
those tolerances.
"""
from .model import equality_constraints, inequality_constraints
if not analysis.valid:
return False
if np.any(inequality_constraints(analysis, targets) > inequality_tolerance):
return False
return bool(np.all(np.abs(equality_constraints(analysis, targets)) <= tolerance))
[docs]
def build_design_space(
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
fixed: Sequence[str] = (),
) -> DesignSpace:
"""Build the GEMSEO design space for the eleven variables.
Args:
bounds: The box ``l_b <= X <= u_b``.
initial: Starting point; defaults to the centre of the box.
fixed: Variables to leave *out* of the design space, because an outer
loop has already chosen them. The disciplines then read them from
their default input data, which
:func:`build_range_scenario` sets to the chosen value. This is how
the discrete gear variables pin ``I``: once a module and a tooth
count are chosen, ``I = 1.5 m z`` is not free.
Returns:
A design space with one scalar variable per entry of
:data:`exlink.design.VARIABLE_NAMES` that is not in ``fixed``.
"""
if initial is None:
start = 0.5 * (bounds.lower + bounds.upper)
else:
start = np.clip(initial.to_array(), bounds.lower, bounds.upper)
space = DesignSpace()
for index, name in enumerate(VARIABLE_NAMES):
if name in fixed:
continue
space.add_variable(
name,
lower_bound=float(bounds.lower[index]),
upper_bound=float(bounds.upper[index]),
value=float(start[index]),
)
return space
def _attach_constraints(
scenario: BaseScenario,
relax_equalities: bool,
equality_tolerance: dict[str, float] | None,
max_height: float,
max_width: float,
) -> None:
"""Attach the formulation's constraints to a scenario."""
for name in INEQUALITY_OUTPUTS:
scenario.add_constraint(name, constraint_type="ineq")
tolerance = dict(DEFAULT_EQUALITY_TOLERANCE)
if equality_tolerance:
tolerance.update(equality_tolerance)
for name in EQUALITY_OUTPUTS:
if relax_equalities:
half_width = tolerance[name]
# |residual| <= half_width, written as two one-sided inequalities.
scenario.add_constraint(
name,
constraint_type="ineq",
value=half_width,
constraint_name=f"{name}_upper",
)
scenario.add_constraint(
name,
constraint_type="ineq",
value=-half_width,
positive=True,
constraint_name=f"{name}_lower",
)
else:
scenario.add_constraint(name, constraint_type="eq")
if np.isfinite(max_height):
scenario.add_constraint(
"height", constraint_type="ineq", value=max_height, constraint_name="height_limit"
)
if np.isfinite(max_width):
scenario.add_constraint(
"width", constraint_type="ineq", value=max_width, constraint_name="width_limit"
)
[docs]
def build_scenario(
objective: str | Sequence[str] = "neg_efficiency",
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
samples: int = DEFAULT_SAMPLES,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets = DEFAULT_TARGETS,
relax_equalities: bool = False,
equality_tolerance: dict[str, float] | None = None,
max_height: float = float("inf"),
max_width: float = float("inf"),
) -> BaseScenario:
"""Assemble a constrained GEMSEO scenario.
Args:
objective: Output name, or several for a multi-objective run.
bounds: The design box.
initial: Starting design.
samples: Crank angles per revolution.
spec: Fixed engine data.
targets: Constraint right-hand sides.
relax_equalities: Express ``STE = 74`` and ``epsilon = 16`` as pairs of
inequalities. Required for algorithms that cannot take equality
constraints, such as NSGA-II.
equality_tolerance: Overrides for :data:`DEFAULT_EQUALITY_TOLERANCE`.
max_height: Moving upper limit on ``H`` [mm].
max_width: Moving upper limit on ``B`` [mm].
Returns:
A scenario ready to :meth:`~gemseo.scenarios.base_scenario.BaseScenario.execute`.
"""
discipline = ExlinkDiscipline(samples=samples, spec=spec, targets=targets)
space = build_design_space(bounds, initial)
scenario = create_scenario(
[discipline],
objective,
space,
formulation_name="DisciplinaryOpt",
)
_attach_constraints(scenario, relax_equalities, equality_tolerance, max_height, max_width)
return scenario
def _best_feasible_vector(problem: Any, objective: str) -> np.ndarray | None:
"""The database entry with the best objective among those satisfying every constraint.
A gradient solver can, and on this problem does, terminate marginally
outside a constraint: SLSQP's own "solution" is then the least-infeasible
point it saw, not a feasible one. Reporting that as the result publishes a
design that misses its brief -- here by 41 % on the top-dead-centre gap --
while the run very likely visited perfectly good feasible points on the way.
So the database is scanned directly. Every constraint is recorded there
under the name it was attached with, and GEMSEO's convention is uniform:
a constraint is satisfied when its recorded value is non-positive,
whichever sign convention it was declared under.
Args:
problem: The scenario's optimization problem.
objective: The objective that was minimised.
Returns:
The best feasible design vector, or ``None`` if the run never found
one.
"""
try:
values = problem.database.get_function_history(objective)
except (KeyError, ValueError):
return None
if values is None or len(values) == 0:
return None
feasible = np.ones(len(values), dtype=bool)
for constraint in problem.constraints:
try:
history = problem.database.get_function_history(constraint.name)
except (KeyError, ValueError):
continue
if history is None or len(history) != len(values):
return None
margin = np.atleast_2d(np.asarray(history, dtype=float).reshape(len(values), -1))
feasible &= np.all(margin <= INEQUALITY_TOLERANCE, axis=1)
if not np.any(feasible):
return None
scores = np.where(feasible, np.asarray(values, dtype=float).ravel(), np.inf)
return np.asarray(problem.database.get_x_vect(int(np.argmin(scores)) + 1))
def _best_design(
scenario: BaseScenario,
objective: str = "neg_efficiency",
fallback: Design | None = None,
) -> Design:
"""Extract the best design from a solved single-objective scenario.
The design space does not always carry all eleven variables: the range
problem pins ``I`` to the gear lattice and leaves it out, so the optimizer
returns ten. A missing entry is read back from the disciplines' own
default input data, which is where :func:`build_range_scenario` put the
pinned value and therefore what the disciplines were actually run with --
rather than from the caller, who would have to re-derive it and could get
it wrong.
Args:
scenario: A solved scenario.
objective: The objective that was minimised.
fallback: Last resort for a variable that is in neither the design
space nor any discipline's defaults; the reference design if
omitted.
Returns:
The best design found.
"""
problem = scenario.formulation.optimization_problem
vector = _best_feasible_vector(problem, objective)
if vector is None:
solution = problem.solution
if solution is not None and solution.x_opt is not None:
vector = solution.x_opt
else:
vector = problem.database.get_x_vect(-1)
names = list(problem.design_space.variable_names)
if len(names) == len(VARIABLE_NAMES):
return Design.from_array(vector)
base = (PUBLISHED_DESIGN if fallback is None else fallback).to_array()
values = dict(zip(names, np.asarray(vector, dtype=float).ravel(), strict=True))
for discipline in scenario.disciplines:
defaults = discipline.default_input_data
for index, name in enumerate(VARIABLE_NAMES):
if name not in values and name in defaults:
base[index] = float(np.ravel(defaults[name])[0])
return Design.from_array(
[values.get(name, base[index]) for index, name in enumerate(VARIABLE_NAMES)]
)
def _run(
scenario: BaseScenario,
algorithm: str,
settings: dict[str, Any],
samples: int,
spec: EngineSpec,
targets: DesignTargets,
) -> Outcome:
"""Execute a scenario and package the outcome."""
scenario.execute(algo_name=algorithm, **settings)
design = _best_design(scenario)
return Outcome(
design=design,
analysis=analyse(design, samples=samples, spec=spec, targets=targets),
algorithm=algorithm,
scenario=scenario,
)
[docs]
def maximise_efficiency(
algorithm: str = DEFAULT_GLOBAL_ALGORITHM,
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
max_iter: int = 4000,
samples: int = DEFAULT_SAMPLES,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets = DEFAULT_TARGETS,
max_height: float = float("inf"),
max_width: float = float("inf"),
relax_equalities: bool = False,
seed: int | None = 1,
**settings: Any,
) -> Outcome:
"""Maximise ``eta`` subject to every constraint.
The single-objective reading of the geometric problem, with the size
objectives either dropped or turned into the moving limits ``max_height`` /
``max_width``.
Args:
algorithm: Any GEMSEO optimizer; the default is differential evolution.
bounds: The design box.
initial: Starting design.
max_iter: Evaluation budget.
samples: Crank angles per revolution.
spec: Fixed engine data.
targets: Constraint right-hand sides.
max_height: Moving upper limit on ``H`` [mm].
max_width: Moving upper limit on ``B`` [mm].
relax_equalities: See :func:`build_scenario`.
seed: Random seed, for the stochastic algorithms that accept one.
**settings: Extra algorithm settings passed through to GEMSEO.
Returns:
The outcome, whose :attr:`Outcome.feasible` flag should be checked.
"""
scenario = build_scenario(
"neg_efficiency",
bounds=bounds,
initial=initial,
samples=samples,
spec=spec,
targets=targets,
relax_equalities=relax_equalities,
max_height=max_height,
max_width=max_width,
)
options: dict[str, Any] = {"max_iter": max_iter}
if seed is not None and algorithm in {"DIFFERENTIAL_EVOLUTION", "DUAL_ANNEALING"}:
options["seed"] = seed
if algorithm.startswith("Augmented_Lagrangian"):
options["sub_algorithm_name"] = "NELDER-MEAD"
options["sub_algorithm_settings"] = {"max_iter": 300}
options.update(settings)
return _run(scenario, algorithm, options, samples, spec, targets)
[docs]
def refine(
design: Design = PUBLISHED_DESIGN,
algorithm: str = DEFAULT_LOCAL_ALGORITHM,
bounds: Bounds | None = None,
relative: float = 0.25,
max_iter: int = 400,
samples: int = DEFAULT_SAMPLES,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets = DEFAULT_TARGETS,
sub_algorithm_name: str = "NELDER-MEAD",
sub_algorithm_settings: dict[str, Any] | None = None,
**settings: Any,
) -> Outcome:
"""Polish a design with the augmented Lagrangian.
The natural last step of a penalty-based workflow: take a point off the
Pareto front and run an augmented Lagrangian from it, precisely because the
external penalty function cannot be both accurate and well conditioned.
Args:
design: The starting design.
algorithm: Local algorithm; the default is the augmented Lagrangian.
bounds: The box to search in; defaults to a window around ``design``.
relative: Half-width of that default window, as a fraction of ``|X_0|``.
max_iter: Outer-iteration budget.
samples: Crank angles per revolution.
spec: Fixed engine data.
targets: Constraint right-hand sides.
sub_algorithm_name: Solver for the augmented Lagrangian sub-problems.
sub_algorithm_settings: Settings for that solver.
**settings: Extra algorithm settings passed through to GEMSEO.
Returns:
The refined outcome.
"""
box = Bounds.around(design, relative=relative) if bounds is None else bounds
scenario = build_scenario(
"neg_efficiency",
bounds=box,
initial=design,
samples=samples,
spec=spec,
targets=targets,
)
options: dict[str, Any] = {"max_iter": max_iter}
if algorithm.startswith("Augmented_Lagrangian"):
options["sub_algorithm_name"] = sub_algorithm_name
options["sub_algorithm_settings"] = sub_algorithm_settings or {"max_iter": 300}
options.update(settings)
return _run(scenario, algorithm, options, samples, spec, targets)
def _supported_settings(algorithm: str, candidates: dict[str, Any]) -> dict[str, Any]:
"""Keep only the settings the given GEMSEO algorithm actually declares.
Lets one set of "run the full budget" defaults be offered to every
multi-objective algorithm without breaking those that do not have, say, a
hypervolume criterion.
Args:
algorithm: The GEMSEO algorithm name.
candidates: Settings to filter.
Returns:
The subset of ``candidates`` accepted by that algorithm.
"""
factory = OptimizationLibraryFactory()
library = factory.get_class(factory.algo_names_to_libraries[algorithm])
fields = set(library.ALGORITHM_INFOS[algorithm].Settings.model_fields)
return {name: value for name, value in candidates.items() if name in fields}
def _extract_front(scenario: BaseScenario) -> list[Design]:
"""Collect the Pareto-optimal designs of a solved multi-objective scenario.
GEMSEO builds the front from the *feasible* points in the database, so an
empty result means the run never found a feasible design -- see
:data:`MOEA_EQUALITY_TOLERANCE`.
"""
problem = scenario.formulation.optimization_problem
pareto = ParetoFront.from_optimization_problem(problem)
inputs = np.atleast_2d(np.asarray(pareto.x_optima, dtype=float))
if inputs.size == 0:
return []
return [Design.from_array(row) for row in inputs]
[docs]
def pareto_front(
algorithm: str = DEFAULT_MOEA,
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
pop_size: int = 200,
max_gen: int = 60,
samples: int = 360,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets | None = None,
equality_tolerance: dict[str, float] | None = None,
seed: int = 1,
**settings: Any,
) -> Outcome:
"""Approximate the Pareto front of ``(-eta, H, B)``.
Args:
algorithm: A multi-objective algorithm; NSGA-II by default.
bounds: The design box.
initial: Starting design, used only by algorithms that take one.
pop_size: Population size. Over the full space this needs hundreds of
individuals to find anything; start smaller and grow.
max_gen: Generation budget. On this problem it matters more than
``pop_size``: seeded in a box around a good design, 20 generations
still return a single point whatever the population, because the
run has not yet worked its way into the thin feasible region, while
35 generations give fronts of several dozen. Spend extra budget
here before spending it on population.
samples: Crank angles per revolution. Lowered by default, because a
MOEA spends tens of thousands of evaluations here.
spec: Fixed engine data.
targets: Constraint right-hand sides; defaults to :func:`moea_targets`,
which relaxes ``g``. Read its docstring before overriding this.
equality_tolerance: Half-widths for the relaxed equalities; defaults to
:data:`MOEA_EQUALITY_TOLERANCE`.
seed: Random seed.
**settings: Extra algorithm settings passed through to GEMSEO.
Returns:
The outcome, with :attr:`Outcome.front` populated and
:attr:`Outcome.design` set to the most efficient point of the front.
Its designs satisfy the *relaxed* problem, so
:attr:`Outcome.feasible` will usually be ``False`` until the chosen
point has been through :func:`refine`.
"""
targets = moea_targets() if targets is None else targets
scenario = build_scenario(
list(OBJECTIVE_NAMES),
bounds=bounds,
initial=initial,
samples=samples,
spec=spec,
targets=targets,
relax_equalities=True,
equality_tolerance=equality_tolerance or MOEA_EQUALITY_TOLERANCE,
)
# GEMSEO's default convergence tests stop NSGA-II after about five
# generations on this problem -- long before a front has formed -- because
# the objectives barely move while the population is still hunting for
# feasible designs. Disable them so the generation budget is the budget
# that actually applies, keeping only the ones this algorithm declares.
options: dict[str, Any] = _supported_settings(
algorithm,
{
"ftol_rel": 0.0,
"ftol_abs": 0.0,
"xtol_rel": 0.0,
"xtol_abs": 0.0,
"hv_tol_rel": 0.0,
"hv_tol_abs": 0.0,
"stop_crit_n_hv": max_gen + 1,
"pop_size": pop_size,
"max_gen": max_gen,
"seed": seed,
},
)
options["max_iter"] = pop_size * (max_gen + 1)
options.update(settings)
scenario.execute(algo_name=algorithm, **options)
front = _extract_front(scenario)
if not front:
LOGGER.warning(
"%s found no feasible design in %d evaluations; widen "
"equality_tolerance, enlarge the population, or seed the run with "
"local_pareto()",
algorithm,
len(scenario.formulation.optimization_problem.database),
)
analyses = [analyse(d, samples=samples, spec=spec, targets=targets) for d in front]
if analyses:
best_index = int(np.argmax([a.metrics.efficiency for a in analyses]))
best_design, best_analysis = front[best_index], analyses[best_index]
else: # pragma: no cover - only if the algorithm returns nothing
best_design = _best_design(scenario)
best_analysis = analyse(best_design, samples=samples, spec=spec, targets=targets)
return Outcome(
design=best_design,
analysis=best_analysis,
algorithm=algorithm,
scenario=scenario,
front=front,
)
[docs]
def local_pareto(
design: Design,
relative: float = 0.1,
absolute_angle: float = 20.0,
**kwargs: Any,
) -> Outcome:
"""Run the MOEA in a box shrunk around a known-good design.
This is the move that makes the multi-objective step work at all: over the
full eleven-dimensional space the MOEA returned solutions worse than the
gradient-based ones, but seeded near an existing optimum it produced a
usable local front. Combining several such fronts into one starting
population, then re-running over the full box, is what eventually gave the
published front.
Args:
design: The centre of the shrunken box.
relative: Half-width as a fraction of ``|X_0|`` for the lengths.
absolute_angle: Half-width in degrees for ``theta_f`` and ``theta_r``.
**kwargs: Forwarded to :func:`pareto_front`.
Returns:
The outcome of the local multi-objective run.
"""
box = Bounds.around(design, relative=relative, absolute_angle=absolute_angle)
return pareto_front(bounds=box, initial=design, **kwargs)
[docs]
def sweep_moving_limits(
limits: Iterable[float],
algorithm: str = DEFAULT_LOCAL_ALGORITHM,
on_height: bool = True,
**kwargs: Any,
) -> list[Outcome]:
"""Trace a Pareto front by walking a moving size limit downwards.
The epsilon-constraint approach to the multi-objective problem: treat ``H`` (or
``B``) as a constraint rather than an objective, solve for maximum
efficiency, tighten the limit, and repeat. Each solve is an ordinary
single-objective problem, and the sequence of solutions traces the
efficiency-versus-size trade-off.
Args:
limits: The successive upper limits [mm], typically decreasing.
algorithm: The single-objective algorithm to use at each step. The
default is the augmented Lagrangian: each step has to satisfy the
two equalities *and* the relaxed-nothing inequalities exactly, and a
derivative-free method such as COBYLA does not get there reliably
within a sensible budget.
on_height: Apply the limit to ``H``; otherwise to ``B``.
**kwargs: Forwarded to :func:`maximise_efficiency`.
Returns:
One outcome per limit, in the order given.
"""
outcomes: list[Outcome] = []
previous: Design | None = kwargs.pop("initial", None)
for limit in limits:
limit_settings: dict[str, Any] = (
{"max_height": limit} if on_height else {"max_width": limit}
)
outcome = maximise_efficiency(
algorithm=algorithm, initial=previous, **limit_settings, **kwargs
)
outcomes.append(outcome)
if outcome.feasible:
# Warm-start the next, tighter problem from this solution.
previous = outcome.design
LOGGER.info(
"moving limit %s <= %.1f mm -> eta = %.4f (%s)",
"H" if on_height else "B",
limit,
outcome.analysis.metrics.efficiency,
"feasible" if outcome.feasible else "infeasible",
)
return outcomes
# =============================================================================
# The coupled problem
#
# Everything above optimizes the quasi-static formulation, in which geometry
# determines performance one way. Adding the sizing discipline and inertia
# closes a loop between them, so these scenarios use the MDF formulation: every
# objective evaluation runs an MDA to convergence before the optimizer sees a
# number.
# =============================================================================
COUPLED_OBJECTIVE_NAMES: tuple[str, str, str, str] = (
"neg_efficiency",
"height",
"width",
"total_mass",
)
"""``f(X) = (-eta, H, B, M)^T``: the three geometric objectives plus mass.
Structural mass cannot appear in the quasi-static formulation at all -- nothing
in it determines a cross-section. It only becomes an objective once the parts are
sized, and it is the objective the dynamics most affects.
"""
COUPLED_INEQUALITY_OUTPUTS: tuple[str, ...] = (
"saturation_margin",
"slenderness_margin",
"bearing_margin",
)
"""Constraints that exist only once the loads are dynamic.
``saturation_margin``
Positive when some member has been driven to the diameter ceiling, i.e. the
sizing loop ran away rather than settling: no section is thick enough to
survive the inertia its own mass creates.
``slenderness_margin``
Positive when a connecting link has grown thicker than a third of its
length, at which point calling it a rod -- and sizing it as a beam -- has
stopped being credible.
``bearing_margin``
Peak crankshaft bearing reaction against its limit.
"""
DEFAULT_MDA = "MDAGaussSeidel"
"""Fixed-point sweep, which is what this coupling wants.
The loop gain is sub-linear (see :mod:`exlink.coupled`), so plain Gauss-Seidel
converges without needing the coupled Jacobians a Newton MDA would demand -- and
the sizing step, a bisection, has no useful derivative to give it anyway.
"""
DEFAULT_MDA_SETTINGS: dict[str, Any] = {
"tolerance": 1.0e-6,
"max_mda_iter": 300,
"warm_start": True,
}
"""Settings for the inner MDA, chosen against measured cost.
Every objective evaluation runs this MDA to convergence, so its settings set
the price of the whole optimization. From a cold start at 2000 rev/min and 360
crank angles:
========== ======= ====== ==========================
tolerance sweeps time section error vs 1e-8
========== ======= ====== ==========================
``1e-8`` 56 9.3 s --
``1e-6`` 44 6.9 s 2e-6 mm
``1e-4`` 31 5.1 s 2e-4 mm
========== ======= ====== ==========================
``1e-6`` buys sub-micron sections for three quarters of the cost, which is far
inside anything a part would be made to.
``warm_start`` matters more than the tolerance: successive design points during
an optimization are close together, so starting each MDA from the previous
converged sections turns a ~50-sweep solve into a handful, and takes the cost
per evaluation from around 9 s to under 2.
"""
[docs]
class BearingMarginDiscipline(Discipline):
"""Turns the peak bearing reaction into a signed constraint.
A one-line discipline rather than a constraint offset, so that the limit
travels with the problem and shows up in the coupling graph.
Args:
limit: Allowable peak bearing reaction [N].
name: Discipline name.
"""
auto_detect_grammar_files: ClassVar[bool] = False
def __init__(self, limit: float = DEFAULT_TARGETS.max_bearing_load, name: str = "") -> None:
super().__init__(name=name)
self.limit = limit
self.input_grammar.update_from_names(["peak_bearing_load"])
self.output_grammar.update_from_names(["bearing_margin"])
self.default_input_data = {"peak_bearing_load": np.zeros(1)}
def _run(self, input_data: StrKeyMapping) -> StrKeyMapping:
load = float(np.ravel(dict(input_data)["peak_bearing_load"])[0])
return {"bearing_margin": np.array([load / self.limit - 1.0])}
def _compute_jacobian(
self,
input_names: Iterable[str] = (),
output_names: Iterable[str] = (),
) -> None:
"""A linear map, so its Jacobian is the reciprocal of the limit."""
self._init_jacobian(input_names, output_names)
self.jac["bearing_margin"]["peak_bearing_load"] = np.array([[1.0 / self.limit]])
[docs]
def build_coupled_scenario(
objective: str | Sequence[str] = "total_mass",
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
speed_rpm: float = DEFAULT_SPEED_RPM,
samples: int = COUPLED_SAMPLES,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets = DEFAULT_TARGETS,
material: Material = DEFAULT_MATERIAL,
safety: SafetyFactors = DEFAULT_SAFETY,
relax_equalities: bool = True,
equality_tolerance: dict[str, float] | None = None,
max_height: float = float("inf"),
max_width: float = float("inf"),
mda_name: str = DEFAULT_MDA,
mda_settings: dict[str, Any] | None = None,
formulation_name: str = "MDF",
) -> BaseScenario:
"""Assemble the coupled scenario, under a chosen formulation.
Three disciplines take part. :class:`~exlink.disciplines.ExlinkDiscipline`
supplies the geometric objectives and constraints and is feed-forward
from the design variables. The other two --
:class:`~exlink.disciplines.DynamicsDiscipline` and
:class:`~exlink.disciplines.StructureDiscipline` -- are strongly coupled to
each other through the section diameters and the internal load histories,
and MDF puts an MDA around them so that every point the optimizer evaluates
is a converged one.
Args:
objective: Output name, or several for a multi-objective run.
bounds: The design box.
initial: Starting design.
speed_rpm: Crankshaft speed [rev/min].
samples: Crank angles per revolution.
spec: Fixed engine data.
targets: Constraint right-hand sides.
material: The material.
safety: The design factors.
relax_equalities: Express the two equalities as pairs of inequalities.
On by default here: the coupled problem is harder, and few
algorithms that cope with it also take equality constraints.
equality_tolerance: Overrides for the relaxed half-widths.
max_height: Moving upper limit on ``H`` [mm].
max_width: Moving upper limit on ``B`` [mm].
mda_name: Inner MDA.
mda_settings: Overrides for :data:`DEFAULT_MDA_SETTINGS`.
formulation_name: ``"MDF"`` (the default: converge the MDA at every
optimizer iteration, so every point evaluated is a real engine) or
``"IDF"`` (hand the diameters to the optimizer as design variables
with consistency constraints, so no inner iteration runs but the
design space grows by seven variables and seven equalities). See
:mod:`exlink.formulations` for a measured comparison.
Returns:
A scenario ready to execute.
"""
disciplines = [
ExlinkDiscipline(samples=samples, spec=spec, targets=targets),
DynamicsDiscipline(
speed_rpm=speed_rpm,
samples=samples,
material=material,
safety=safety,
spec=spec,
),
StructureDiscipline(samples=samples, material=material, safety=safety, spec=spec),
BearingMarginDiscipline(limit=targets.max_bearing_load),
]
settings = dict(DEFAULT_MDA_SETTINGS)
settings.update(mda_settings or {})
# Only MDF takes an inner MDA; IDF replaces it with consistency constraints.
formulation_settings: dict[str, Any] = (
{"main_mda_name": mda_name, "main_mda_settings": settings}
if formulation_name == "MDF"
else {}
)
scenario = create_scenario(
disciplines,
objective,
build_design_space(bounds, initial),
formulation_name=formulation_name,
**formulation_settings,
)
_attach_constraints(scenario, relax_equalities, equality_tolerance, max_height, max_width)
for name in COUPLED_INEQUALITY_OUTPUTS:
scenario.add_constraint(name, constraint_type="ineq")
return scenario
DEFAULT_COUPLED_ALGORITHM = "SLSQP"
"""Single-objective solver for the coupled problem.
Gradient-based, and not by preference. The feasible set is a sliver -- two
equalities 0.1 wide inside an eleven-dimensional box with sides of tens of
millimetres -- and a derivative-free method has no direction to move in: COBYLA
returns its starting point unchanged after 120 evaluations and 313 s. With the
analytic Jacobians of :mod:`exlink.jacobian` and
:mod:`exlink.dynamics_jacobian`, SLSQP takes the same problem from 1.04 kg to
0.23 kg in 40 evaluations.
"""
[docs]
def minimise_mass(
algorithm: str = DEFAULT_COUPLED_ALGORITHM,
bounds: Bounds | None = None,
initial: Design | None = None,
speed_rpm: float = DEFAULT_SPEED_RPM,
max_iter: int = 60,
relative: float = 0.25,
min_efficiency: float = 0.25,
**kwargs: Any,
) -> Outcome:
"""Minimise total moving mass at a floor on efficiency.
The natural single-objective reading of the coupled problem: efficiency
becomes a constraint to hold, and the mass the dynamics creates becomes the
thing to reduce.
Budget realistically, and note that the budget depends entirely on whether
the solver has gradients. With them (the default) a few dozen evaluations
suffice: 40 take the published geometry from 1.04 kg to 0.23 kg at 2000 rev/min.
Without them a derivative-free search over eleven variables needs several
hundred, and on this feasible set may never move at all.
Args:
algorithm: A single-objective GEMSEO optimizer.
bounds: The design box; defaults to a window around ``initial``.
initial: Starting design; defaults to the refined reference.
speed_rpm: Crankshaft speed [rev/min].
max_iter: Evaluation budget.
relative: Half-width of the default box, as a fraction of ``|X_0|``.
min_efficiency: Floor on ``eta``.
**kwargs: Forwarded to :func:`build_coupled_scenario`.
Returns:
The outcome; its :attr:`Outcome.design` is the lightest design found.
"""
from .reference import REFINED_DESIGN
start = REFINED_DESIGN if initial is None else initial
box = Bounds.around(start, relative=relative) if bounds is None else bounds
scenario = build_coupled_scenario(
"total_mass", bounds=box, initial=start, speed_rpm=speed_rpm, **kwargs
)
scenario.add_constraint(
"efficiency",
constraint_type="ineq",
value=min_efficiency,
positive=True,
constraint_name="efficiency_floor",
)
settings: dict[str, Any] = {"max_iter": max_iter}
scenario.execute(algo_name=algorithm, **settings)
design = _best_design(scenario)
return Outcome(
design=design,
analysis=analyse(design, samples=COUPLED_SAMPLES),
algorithm=algorithm,
scenario=scenario,
)
EFFICIENCY_MASS_OBJECTIVES: tuple[str, str] = ("neg_efficiency", "total_mass")
"""``(-eta, M)``: the trade that only exists once the parts are sized.
The four-way front over :data:`COUPLED_OBJECTIVE_NAMES` is available, but ``H``,
``B`` and ``M`` are strongly correlated -- a bigger mechanism is a heavier one --
so most of it is redundant. Efficiency against mass is the trade with real
tension in it, because the geometry that maximises the quasi-static lever arm is
the geometry that makes the parts heavy.
"""
[docs]
def sweep_efficiency_floor(
floors: Iterable[float],
speed_rpm: float = DEFAULT_SPEED_RPM,
algorithm: str = DEFAULT_COUPLED_ALGORITHM,
initial: Design | None = None,
relative: float = 0.30,
max_iter: int = 60,
**kwargs: Any,
) -> list[Outcome]:
"""Trace the efficiency-versus-mass trade by epsilon-constraint.
Minimise mass subject to ``eta >= floor``, for a ladder of floors, warm
starting each solve from the last. This is the practical way to get the
coupled front: every step is an ordinary single-objective problem on the
*unrelaxed* constraints, and a handful of local solves replaces the tens of
thousands of MDA-backed evaluations a population method would need.
Args:
floors: Efficiency floors to impose, typically decreasing so that each
solve starts from a feasible point found by the last.
speed_rpm: Crankshaft speed [rev/min].
algorithm: The single-objective algorithm used at each step.
initial: Starting design; defaults to the refined reference.
relative: Half-width of the design box, as a fraction of ``|X_0|``.
max_iter: Evaluation budget per step.
**kwargs: Forwarded to :func:`minimise_mass`.
Returns:
One outcome per floor, in the order given.
"""
from .reference import REFINED_DESIGN
previous = REFINED_DESIGN if initial is None else initial
outcomes: list[Outcome] = []
for floor in floors:
outcome = minimise_mass(
algorithm=algorithm,
initial=previous,
speed_rpm=speed_rpm,
max_iter=max_iter,
relative=relative,
min_efficiency=floor,
**kwargs,
)
outcomes.append(outcome)
if outcome.analysis.valid:
previous = outcome.design
LOGGER.info(
"efficiency floor %.3f -> eta = %.4f",
floor,
outcome.analysis.metrics.efficiency,
)
return outcomes
[docs]
def coupled_pareto(
objective: Sequence[str] = EFFICIENCY_MASS_OBJECTIVES,
algorithm: str = DEFAULT_MOEA,
bounds: Bounds | None = None,
initial: Design | None = None,
speed_rpm: float = DEFAULT_SPEED_RPM,
pop_size: int = 40,
max_gen: int = 20,
relative: float = 0.25,
seed: int = 1,
**kwargs: Any,
) -> Outcome:
"""Approximate a Pareto front of the *coupled* problem with a MOEA.
Budget with care. Every individual of every generation runs an MDA to
convergence, so the cost is ``pop_size * max_gen`` MDAs -- the defaults here
are an order of magnitude smaller than the standalone
:func:`pareto_front` for that reason, and are a starting point rather than a
converged front. :func:`sweep_efficiency_floor` reaches the same trade for
a small fraction of the work and is the better tool unless a genuinely
non-convex front is suspected.
Args:
objective: Output names to trade off.
algorithm: A multi-objective algorithm.
bounds: The design box; defaults to a window around ``initial``.
initial: Starting design; defaults to the refined reference.
speed_rpm: Crankshaft speed [rev/min].
pop_size: Population size.
max_gen: Generation budget.
relative: Half-width of the default box.
seed: Random seed.
**kwargs: Forwarded to :func:`build_coupled_scenario`.
Returns:
The outcome, with :attr:`Outcome.front` populated.
"""
from .reference import REFINED_DESIGN
start = REFINED_DESIGN if initial is None else initial
box = Bounds.around(start, relative=relative) if bounds is None else bounds
scenario = build_coupled_scenario(
list(objective),
bounds=box,
initial=start,
speed_rpm=speed_rpm,
targets=moea_targets(),
**kwargs,
)
options: dict[str, Any] = _supported_settings(
algorithm,
{
"ftol_rel": 0.0,
"ftol_abs": 0.0,
"xtol_rel": 0.0,
"xtol_abs": 0.0,
"hv_tol_rel": 0.0,
"hv_tol_abs": 0.0,
"stop_crit_n_hv": max_gen + 1,
"pop_size": pop_size,
"max_gen": max_gen,
"seed": seed,
},
)
options["max_iter"] = pop_size * (max_gen + 1)
scenario.execute(algo_name=algorithm, **options)
front = _extract_front(scenario)
if not front:
LOGGER.warning(
"%s found no feasible design in %d evaluations",
algorithm,
len(scenario.formulation.optimization_problem.database),
)
analyses = [analyse(d, samples=COUPLED_SAMPLES) for d in front]
if analyses:
best = int(np.argmax([a.metrics.efficiency for a in analyses]))
design, analysis = front[best], analyses[best]
else: # pragma: no cover - only if the algorithm returns nothing
design = _best_design(scenario)
analysis = analyse(design, samples=COUPLED_SAMPLES)
return Outcome(
design=design,
analysis=analysis,
algorithm=algorithm,
scenario=scenario,
front=front,
)
# =================================================================================
# The vehicle-level problem
#
# Every formulation above stops at the engine. This one carries through to what
# the competition scores, and in doing so it turns a multi-objective problem into
# a single-objective one -- not by picking weights, but because the physics fixes
# the exchange rates.
# =================================================================================
RANGE_INEQUALITY_OUTPUTS: tuple[str, ...] = ("runs_margin", "gear_margin")
"""Constraints that exist only once the vehicle is in the problem.
``runs_margin``
Positive when the engine produces net work: friction must not exceed the
indicated work. A near-singular linkage at speed can fail this outright,
and nothing in the geometric problem would have noticed.
``gear_margin``
Positive when the gear pair is inside its face-width limit, i.e. the mesh
can be kept aligned across the tooth face.
"""
[docs]
def build_range_scenario(
bounds: Bounds = GLOBAL_BOUNDS,
initial: Design | None = None,
speed_rpm: float = DEFAULT_SPEED_RPM,
samples: int = COUPLED_SAMPLES,
vehicle: Any = None,
module: float | None = None,
teeth: int | None = None,
spec: EngineSpec = DEFAULT_SPEC,
targets: DesignTargets = DEFAULT_TARGETS,
material: Material = DEFAULT_MATERIAL,
safety: SafetyFactors = DEFAULT_SAFETY,
relax_equalities: bool = True,
beta_target: float | None = None,
grade: int = 8,
equality_tolerance: dict[str, float] | None = None,
mda_name: str = DEFAULT_MDA,
mda_settings: dict[str, Any] | None = None,
) -> BaseScenario:
"""Assemble the range problem: five disciplines, one objective, MDF.
The discipline graph is the argument of the whole study:
.. code-block:: text
ExlinkDiscipline feed-forward: geometry -> eta, H, B, constraints
|
| +--------------------------------------+
| | DynamicsDiscipline <==> Structure | <- MDA (strong)
| +--------------------------------------+
| | diameters
v v
RangeDiscipline --> neg_range
:class:`~exlink.disciplines.RangeDiscipline` sits downstream of the MDA and
consumes its converged coupling variable, so every point the optimizer
evaluates is a range computed on a self-consistent structure.
The gear pair is held fixed
------------------------------
The module and tooth count are **always** pinned here, never left to the
automatic choice inside :func:`exlink.gears.size_pair`. This is not a
convenience, and getting it wrong is instructive.
Letting the module float makes the objective a step function of ``I``: the
lightest workable module changes at a threshold, and the range jumps 40
km/L across it. A central difference straddling such a threshold returns a
gradient of order 1e5 against a true gradient of order 1e3, and SLSQP --
handed a quadratic subproblem built from two constraint gradients that are
both quantisation noise -- rejects it as *"inequality constraints
incompatible"* and stops at the starting point without evaluating anything.
Pinning the pair is also simply the correct treatment. With a 2:1 ratio,
choosing ``(m, z)`` fixes ``I = 1.5 m z``, so ``I`` leaves the design space
entirely: the continuous problem has ten variables, not eleven, and the
discrete choice is enumerated outside it by :func:`maximise_range`.
Args:
bounds: The design box.
initial: Starting design.
speed_rpm: Crankshaft speed [rev/min]. Held fixed; see
:func:`maximise_range`, which sweeps it.
samples: Crank angles per revolution.
vehicle: The car; a default Prototype-class entry if omitted.
module: Gear module [mm]. Derived from ``initial`` if omitted, then
held fixed for the whole run.
teeth: Teeth on the small gear; derived likewise.
spec: Fixed engine data.
targets: Constraint right-hand sides.
material: The material.
safety: The design factors.
relax_equalities: Express the two equalities as pairs of inequalities.
equality_tolerance: Overrides for the relaxed half-widths.
mda_name: Inner MDA.
mda_settings: Overrides for :data:`DEFAULT_MDA_SETTINGS`.
beta_target: When given, the search is reliability-constrained: a
:class:`~exlink.robustness.FailureProbabilityDiscipline` is added
and ``system_beta >= beta_target`` attached, so the optimizer must
hold a system reliability index rather than merely satisfy the
nominal constraints. ``None`` leaves the problem deterministic.
grade: ISO 286 IT grade assumed for the machined dimensions when
``beta_target`` is set. Ignored otherwise.
Returns:
A scenario ready to execute, over ten continuous variables.
"""
from .gears import lattice_inter_axle, size_pair, tooth_count
start = PUBLISHED_DESIGN if initial is None else initial
if module is None:
module = size_pair(start.I, 1000.0).module
if teeth is None:
teeth = tooth_count(start.I, module)
# Choosing the pair pins the centre distance, so the design starts there.
inter_axle = lattice_inter_axle(module, teeth)
start = start.replace(I=inter_axle)
disciplines: list[Discipline] = [
ExlinkDiscipline(samples=samples, spec=spec, targets=targets),
DynamicsDiscipline(
speed_rpm=speed_rpm,
samples=samples,
material=material,
safety=safety,
spec=spec,
),
StructureDiscipline(samples=samples, material=material, safety=safety, spec=spec),
BearingMarginDiscipline(limit=targets.max_bearing_load),
RangeDiscipline(
speed_rpm=speed_rpm,
samples=samples,
vehicle=vehicle,
module=module,
teeth=teeth,
material=material,
safety=safety,
spec=spec,
),
]
if beta_target is not None:
from .robustness import FailureProbabilityDiscipline
# Reliability enters as one more discipline, not as a post-processing
# step: its output is a constraint the optimizer must hold at every
# iteration. It depends on the design variables only, so it does not
# join the MDA -- the coupling is unchanged and ยง4.3's count still
# holds.
disciplines.append(
FailureProbabilityDiscipline(
grade=grade, samples=samples, targets=targets, spec=spec
)
)
# ``I`` is not a design variable any more, so every discipline has to be
# told its pinned value; they would otherwise fall back to the reference.
for discipline in disciplines:
if "I" in discipline.input_grammar:
discipline.default_input_data = {
**discipline.default_input_data,
"I": np.array([inter_axle]),
}
settings = dict(DEFAULT_MDA_SETTINGS)
settings.update(mda_settings or {})
scenario = create_scenario(
disciplines,
"neg_range",
build_design_space(bounds, start, fixed=("I",)),
formulation_name="MDF",
main_mda_name=mda_name,
main_mda_settings=settings,
)
_attach_constraints(
scenario, relax_equalities, equality_tolerance, float("inf"), float("inf")
)
# Two opposite sign conventions meet here, so they are attached separately.
# The coupled margins are violations: positive means saturated, too stubby,
# or over the bearing limit, so they are bounded above by zero. The range
# margins are margins: positive means the engine runs and the gears fit, so
# they are bounded below by zero.
for name in COUPLED_INEQUALITY_OUTPUTS:
scenario.add_constraint(name, constraint_type="ineq")
for name in RANGE_INEQUALITY_OUTPUTS:
scenario.add_constraint(name, constraint_type="ineq", positive=True)
if beta_target is not None:
# The reliability index, not the probability: P_f saturates towards 0
# and 1, where its gradient vanishes and the optimizer stalls, while
# beta = -Phi^-1(P_f) stays well scaled across the whole range.
scenario.add_constraint(
"system_beta",
constraint_type="ineq",
value=beta_target,
positive=True,
constraint_name="reliability",
)
return scenario
[docs]
def maximise_range(
algorithm: str = DEFAULT_COUPLED_ALGORITHM,
bounds: Bounds | None = None,
initial: Design | None = None,
speeds: Sequence[float] = (800.0, 1000.0, 1250.0),
candidates: int = 3,
max_iter: int = 40,
relative: float = 0.30,
**kwargs: Any,
) -> tuple[Design, float, dict[tuple[float, float, int], float]]:
"""Maximise kilometres per litre over a mixed discrete/continuous problem.
Two variables are enumerated outside the gradient-based solve, for two
different reasons.
**The gear pair, because it is genuinely discrete.** A module comes off the
ISO 54 list and a tooth count is an integer, and together they pin
``I = 1.5 m z``. There is no meaningful derivative with respect to either,
and pretending otherwise is what breaks the optimizer -- see
:func:`build_range_scenario`. Each candidate pair defines its own
ten-variable continuous problem, and the equality constraints ``STE = 74``
and ``epsilon = 16`` then have to be re-satisfied by the remaining
variables, since ``I`` was one of the variables that used to satisfy them.
That "choose the integers, repair the continuum" structure is the whole
mixed-integer content of the problem.
**The speed, because it enters as a construction parameter.** The dynamics
discipline is built *at* a speed. A sweep is also the honest treatment: the
trade it governs is sharp and non-convex -- flywheel mass falls as
``1 / omega^2`` while structural mass climbs with ``omega^2`` -- and seeing
the curve is worth more than seeing one converged number.
Args:
algorithm: A single-objective GEMSEO optimizer. Gradient-based by
default; see :data:`DEFAULT_COUPLED_ALGORITHM` for why.
bounds: The design box; defaults to a window around ``initial``.
initial: Starting design; defaults to the coupled reference.
speeds: Crankshaft speeds to optimize at [rev/min].
candidates: Gear pairs to try, nearest first on the ``I`` lattice.
max_iter: Evaluation budget per (speed, gear pair).
relative: Half-width of the default box, as a fraction of ``|X_0|``.
**kwargs: Forwarded to :func:`build_range_scenario`.
Returns:
``(best_design, best_km_per_litre, history)``, where the history maps
``(speed, module, teeth)`` to the range reached. Combinations where
nothing feasible was found map to zero.
"""
import math
from .constants import DEFAULT_SPEC
from .coupled import solve_for_design
from .gears import buildable_neighbours, lattice_inter_axle
from .performance import evaluate
from .reference import COUPLED_DESIGN
start = COUPLED_DESIGN if initial is None else initial
box = Bounds.around(start, relative=relative) if bounds is None else bounds
# Enumerate the lattice by what can carry the load, not by what is nearest.
# The nearest centre distances are reached with the smallest modules, and a
# small module needs a wide face -- so ranking by distance alone routinely
# offers only unbuildable pairs.
reference = solve_for_design(start, speed_rpm=speeds[0])
tangential = float(np.max(np.abs(reference.loads.gear_force))) * math.cos(
DEFAULT_SPEC.pressure_angle
)
ranked = buildable_neighbours(start.I, tangential, shaft_bore=12.0, count=candidates + 1)
pairs = [(module, teeth, value) for module, teeth, value, _pair in ranked[:candidates]]
history: dict[tuple[float, float, int], float] = {}
best_design, best_value = start, 0.0
for speed in speeds:
for module, teeth, _inter_axle in pairs:
key = (float(speed), float(module), int(teeth))
scenario = build_range_scenario(
bounds=box,
initial=start,
speed_rpm=speed,
module=module,
teeth=teeth,
**kwargs,
)
try:
scenario.execute(algo_name=algorithm, max_iter=max_iter)
except Exception: # a failed combination is data, not an error
history[key] = 0.0
continue
pinned = start.replace(I=lattice_inter_axle(module, teeth))
candidate = _best_design(scenario, objective="neg_range", fallback=pinned)
outcome = evaluate(candidate, speed_rpm=speed, module=module, teeth=teeth)
value = outcome.km_per_litre if outcome.feasible else 0.0
history[key] = value
if value > best_value:
best_design, best_value = candidate, value
return best_design, best_value, history
[docs]
def project_onto_equalities(
design: Design,
samples: int = COUPLED_SAMPLES,
targets: DesignTargets = DEFAULT_TARGETS,
spec: EngineSpec = DEFAULT_SPEC,
iterations: int = 8,
tolerance: float = 1.0e-9,
fixed: Sequence[str] = ("I",),
) -> Design:
"""Restore ``STE = 74`` and ``epsilon = 16`` by the smallest design change.
A gradient solver stops within *its* convergence tolerance of the
constraints it was given, not exactly on them. On this problem SLSQP
routinely finishes a few parts in ten thousand outside the relaxed equality
band -- physically irrelevant, since the tolerance study puts the machining
standard deviation of the stroke two orders of magnitude higher, but enough
for a strict feasibility check to reject the design.
Rather than loosen the check, the design is projected back onto the
equality manifold. With the exact Jacobian rows for both equalities from
:mod:`exlink.jacobian`, the minimum-norm Newton step is
.. math:: \\Delta X = -J^{+} r
which is the smallest change in the design that restores the equalities, so
everything else about the design -- and hence its range -- moves as little
as it can. Two or three iterations converge to machine precision.
Args:
design: A design near, but not on, the equality manifold.
samples: Crank angles per revolution.
targets: Constraint right-hand sides.
spec: Fixed engine data.
iterations: Newton step limit. The residual is measured at
``samples``, so a caller that later checks feasibility at a
different crank-angle resolution should project at that resolution:
``STE`` and ``g`` both shift slightly with it.
tolerance: Residual norm to stop at.
fixed: Variables the step may not touch. ``I`` by default, because it
is pinned to the gear lattice: a projection free to move it would
hand back a mechanism whose gears cannot be cut, which is exactly
what the pinning exists to prevent.
Returns:
The projected design, or the input unchanged if it cannot be analysed
or the Jacobian is singular there.
Warning:
Restoring the equalities can *break* the top-dead-centre gap. The
minimum-norm step on this problem is a few hundredths of a millimetre,
and that is enough to move ``g`` from 0.001 mm to 0.020 mm -- twice its
bound. This is not a defect in the projection; it is the same
hypersensitivity that :mod:`exlink.robustness` measures against
machining tolerance, seen once more. Always re-check the inequalities
afterwards rather than assuming a projected design is feasible.
"""
from .jacobian import kinematic_jacobian, metric_jacobian
from .model import analyse, equality_constraints
current = design
for _ in range(iterations):
analysis = analyse(current, samples=samples, spec=spec)
if not analysis.valid:
return design
residual = equality_constraints(analysis, targets)
if float(np.linalg.norm(residual)) <= tolerance:
return current
kinematic = kinematic_jacobian(current, analysis.require_solved().kinematics, spec)
rows = metric_jacobian(current, analysis, kinematic, spec)
jacobian = np.stack([rows["expansion_stroke"], rows["compression_ratio"]])
# Zeroing a column removes that variable from the step: the minimum-norm
# solution is then taken over the remaining ones only.
for index, name in enumerate(VARIABLE_NAMES):
if name in fixed:
jacobian[:, index] = 0.0
try:
step = -np.linalg.pinv(jacobian) @ residual
except np.linalg.LinAlgError: # pragma: no cover - defensive
return design
current = Design.from_array(current.to_array() + step)
return current
[docs]
@dataclass(frozen=True)
class MultiStartOutcome:
"""What a multi-start search found, and how varied its local optima were."""
design: Design
"""The best feasible design found."""
value: float
"""Its objective value."""
starts: int
feasible_starts: int
"""How many restarts ended somewhere feasible."""
values: list[float] = field(default_factory=list)
"""Objective at every feasible local solution found."""
seconds: float = 0.0
@property
def spread(self) -> float:
"""Range of the feasible local optima, as a fraction of the best.
Near zero means every restart reached effectively the same point, which
is evidence -- not proof -- that the single-start answer was not merely
the nearest local optimum.
"""
if len(self.values) < 2 or self.value == 0.0:
return 0.0
return (max(self.values) - min(self.values)) / abs(self.value)
[docs]
def multistart(
objective: str = "neg_efficiency",
initial: Design | None = None,
bounds: Bounds | None = None,
n_start: int = 10,
spread: float = 0.08,
max_iter: int = 60,
algorithm: str = DEFAULT_COUPLED_ALGORITHM,
samples: int = DEFAULT_SAMPLES,
seed: int = 0,
build: Any = None,
evaluate_value: Any = None,
**kwargs: Any,
) -> MultiStartOutcome:
"""Restart a gradient solve from perturbed points on the feasible manifold.
Why not GEMSEO's ``MultiStart``
--------------------------------
Because on this problem it cannot work, and the reason is structural rather
than a matter of budget. ``MultiStart`` draws its restarts from an LHS over
the design box, and this feasible set contains two *equality* constraints --
``STE = 74`` and ``epsilon = 16`` -- so it is a codimension-two manifold and
has measure zero. Sampling it uniformly succeeds with probability zero, and
that is what happens: **no feasible point in 12 000 uniform samples**, not
even within 10 % of a design known to be feasible. Only 5.8 % of the global
box is analysable at all. Run anyway, ``MultiStart`` returns its starting
point unchanged.
What works instead
------------------
Restart from points that are *on* the manifold. Each restart perturbs the
incumbent, projects the perturbation back onto the two equalities with
:func:`project_onto_equalities` -- a deterministic Newton step from the
analytic Jacobians, so the start is feasible in the equalities by
construction -- and lets the optimizer restore the inequalities from there.
That is the standard treatment for an equality-constrained problem, and it
is what makes the question answerable at all: whether the single-start
answers this package reports are local optima or something better.
Args:
objective: Output to minimise.
initial: Incumbent to perturb around; the refined reference if omitted.
bounds: The design box.
n_start: Restarts, including the incumbent itself.
spread: Perturbation size, as a fraction of each variable's magnitude.
max_iter: Budget per restart.
algorithm: A single-objective GEMSEO optimizer.
samples: Crank angles per revolution.
seed: Random seed, so the search is reproducible.
build: Scenario builder taking ``(bounds, initial)``; the geometric
:func:`build_scenario` if omitted.
evaluate_value: Maps a design to its objective; reads the analysis if
omitted.
**kwargs: Forwarded to the builder.
Returns:
The best feasible design found, and the spread of the local optima.
"""
import time
from .reference import REFINED_DESIGN
start = REFINED_DESIGN if initial is None else initial
box = Bounds.around(start, relative=0.35) if bounds is None else bounds
rng = np.random.default_rng(seed)
began = time.perf_counter()
def default_build(box_: Bounds, initial_: Design) -> BaseScenario:
return build_scenario(
objective,
bounds=box_,
initial=initial_,
samples=samples,
relax_equalities=True,
**kwargs,
)
def default_value(design: Design) -> float:
analysis = analyse(design, samples=samples)
return float(-analysis.metrics.efficiency)
builder = default_build if build is None else build
scorer = default_value if evaluate_value is None else evaluate_value
best_design, best_value = start, float("inf")
values: list[float] = []
feasible_starts = 0
for index in range(n_start):
if index == 0:
candidate = start
else:
base = start.to_array()
step = rng.normal(0.0, spread, size=base.size) * np.maximum(np.abs(base), 1.0)
moved = Design.from_array(np.clip(base + step, box.lower, box.upper))
# Put the restart back on the equality manifold; a perturbation off
# it starts the optimizer outside a set it may not re-enter.
candidate = project_onto_equalities(moved, samples=samples)
try:
scenario = builder(box, candidate)
scenario.execute(algo_name=algorithm, max_iter=max_iter)
found = _best_design(scenario, objective=objective, fallback=candidate)
except Exception: # a restart that fails is data, not an error
continue
if not is_feasible(analyse(found, samples=samples)):
continue
feasible_starts += 1
value = scorer(found)
values.append(value)
if value < best_value:
best_design, best_value = found, value
return MultiStartOutcome(
design=best_design,
value=best_value,
starts=n_start,
feasible_starts=feasible_starts,
values=values,
seconds=time.perf_counter() - began,
)