1 Introduction↩︎

A co-design problem is the joint design of components whose specifications mutually depend on each other. A canonical example: a battery must store enough energy to power an actuator, but the actuator must lift the battery itself, so the battery’s mass appears in the energy it has to supply. There is no way to fix one without fixing the other.

[1] showed that under modest monotonicity assumptions, such problems have a clean algebraic theory: the relations between resources and functionalities form a category closed under series, parallel, and feedback composition, and solutions are computed by a Kleene fixed-point iteration in the lattice of antichains. The original publication included a software tool, MCDPL [2], distributed through co-design.science, which exposes the framework via a domain-specific language.

codesign-mcdp is a from-scratch Python implementation of the algorithmic core, designed around composable Python objects rather than a separate DSL. It targets users who want to formulate and solve MCDPs from inside their existing Python codebase, without context switching to a DSL or an external solver. The algorithmic core (posets, antichains, primitive design problems, composition, solver, and both builders) has no required runtime dependencies beyond the standard library; the uncertainty, online-learning, and visualisation layers pull in numpy, scipy, and matplotlib as optional extras. The whole package is roughly 8000 lines of Python and is licensed under MIT.

1.0.0.1 Attribution.

The theory implemented here is not the author’s. Monotone co-design was introduced by Andrea Censi [1], [3], who also established the treatment of cyclic constraints [4] and of uncertainty [5], and who wrote the reference implementation MCDPL [2]. The framework has been developed since by Gioele Zardini, Dejan Milojevic, and colleagues, across future mobility systems [6]–[8], autonomous systems and embodied intelligence [9]–[11], and the co-design of perception and compute [12], with more recent work on composable and distributional uncertainty [13]–[15], categorical and scalability extensions [16], [17], and online learning [18]. codesign-mcdp is an independent Python implementation of that framework. It contributes software, not theory. Every definition, theorem, and algorithm in Section 2, and the uncertainty and online-learning layers, belong to the authors cited above. Users of this library should cite their papers, not only this manual. Citing the software is not a substitute for citing the theory. Please cite [1] for the monotone co-design framework itself, and the paper behind whichever layer you use: [5] for uncertainty, [18] for the online-learning layer, and the relevant application paper from the co-design group [7], [12].

The temporal, sequential, vector-state, and receding-horizon layers of Section 15 are the exception. Those are the author’s own work, built on top of the monotone co-design framework rather than taken from it. They are being written up separately and are not yet published, so they carry no citation here. Section 18 maps the literature and marks the boundary again.

1.0.0.2 Scope of this document.

This manual covers the public API and the mathematical structure it implements. Section 2 recalls the formal background. Sections 4 through 7 describe the data types, primitive design problems, composition operators, and the solver in detail. Section 8 covers set-based and stochastic uncertainty, and Section 9 the online-learning layer; Section 10 documents the visualisation helpers. Sections 11 and 12 document the two high-level builders. Section 13 is a complete, per-module reference for every public symbol in the library. Section 14 walks through every example shipped with the library. Section 15 introduces the temporal, vector-state, and online co-design layers, and Section 16 surveys realistic problem domains for them. Section 17 collects modelling guidelines and common pitfalls, and Section 19 discusses limitations and future directions.

1.0.0.3 Conventions.

Code listings are typeset in monospaced font. Type signatures use Python’s standard notation (\(\texttt{f: F} \rightarrow \texttt{R}\), for example). Mathematical quantities use \(F\), \(R\), \(h\), \(\top\), \(\bot\), \(\mathsf{Min}\), etc.

2 Mathematical Background↩︎

This section recalls the formal definitions from [1] needed to specify the API. The definitions, the theorems, and the fixed-point algorithm are all his; nothing in this section is original to this library. Proofs are omitted; the reader is referred to the paper for the full development. Standard order-theoretic background is in [19]. The treatment of design problems whose constraint graph contains cycles is developed further in [4].

2.1 Posets and Lattices↩︎

2.2 Antichains as a Lattice↩︎

The set of antichains of a poset carries a natural order: \(A \leq A'\) in \(\mathsf{A}[P]\) iff every point in \(A'\) is dominated by some point in \(A\). Intuitively, \(A \leq A'\) means \(A\) is "lower" or "easier", having weaker Pareto-minimal demands.

This lattice is the state space of the Kleene iteration (Section 6.3).

2.3 Design Problems↩︎

Intuitively, monotonicity says that demanding more functionality cannot reduce the minimal resources needed. This corresponds to the common-sense engineering principle that you don’t get more for less.

2.4 Composition↩︎

The class of MCDPs is closed under three operators.

The recursion is solved by Kleene iteration in \(\mathsf{A}[R \times X]\), seeded at \(\{\bot\}\) and ascending until a fixed point is reached.

2.5 The Kleene Fixed-Point Algorithm↩︎

In implementation, the iteration terminates when one of three conditions holds: \(A_{k+1} = A_k\) (fixed point), \(A_{k+1} = \emptyset\) (infeasibility), or every point in \(A_{k+1}\) has \(r_X = \top_X\) (loop axis saturates).

3 Installation and Project Layout↩︎

codesign-mcdp requires Python 3.9 or newer. The algorithmic core has no required runtime dependencies; the optional layers listed below pull in numpy, scipy, matplotlib, or graphviz as needed.

3.0.0.1 Installation.

Clone the repository and install in editable mode:

git clone https://github.com/cbriat/codesign-mcdp.git
cd codesign-mcdp
pip install -e .

3.0.0.2 Optional extras.

• pip install -e ".[viz]" adds matplotlib for the visualisation helpers (Section 10) and the plotting examples.

• pip install -e ".[online]" adds numpy and scipy for the uncertainty (Section 8) and online-learning (Section 9) layers.

• pip install -e ".[diagram]" adds graphviz for the block-diagram renderer (Section 10.1).

• pip install -e ".[nb]" adds the Jupyter tooling needed to (re)generate the notebooks.

3.0.0.3 Repository layout.

codesign-mcdp/
  codesign/          the Python package (21 modules)
  examples/          25 runnable example scripts
  notebooks/         the same examples as executed .ipynb notebooks
  tests/             smoke tests
  docs/              this manual plus rendered trace images
  build_notebooks.py script to regenerate the notebooks
  pyproject.toml     packaging metadata, PEP 621
  .github/workflows/ CI configuration

4 Core Data Types↩︎

This section documents the data types provided by codesign.posets and codesign.antichains.

4.1 Posets↩︎

The Poset abstract base class has the following interface:

class Poset(ABC):
    name: str

    def leq(self, x, y) -> bool: ...      # is x <= y?
    def eq(self, x, y) -> bool: ...
    def lt(self, x, y) -> bool: ...
    def bottom(self): ...                  # least element, if any
    def top(self): ...                     # greatest element, if any
    def is_bottom(self, x) -> bool: ...
    def is_top(self, x) -> bool: ...
    def join(self, x, y): ...              # least upper bound
    def format(self, x) -> str: ...

Four concrete poset classes are provided.

4.1.0.1 Reals(unit="").

The non-negative reals augmented with \(+\infty\). Bottom is \(0.0\); top is math.inf. The optional unit string is used for pretty-printing only.

4.1.0.2 Naturals(unit="").

The non-negative integers augmented with \(+\infty\). Bottom is \(0\); top is math.inf.

4.1.0.3 Ports(components: dict[str, Poset]).

The componentwise-ordered product of named factor posets. Elements are Python dicts mapping component names to values. The bottom is the dict of all bottoms; the top is the dict of all tops. Ports can be nested freely: a component poset may itself be a Ports. This is the type used for every design problem’s \(F\) and \(R\).

The name reflects the library’s vocabulary: each named component is a port of the design problem, addressable by name in the constraint DSL. For backward compatibility the older name NamedProduct remains exported as an alias, so existing code that imports NamedProduct continues to work; new code should prefer Ports.

4.1.0.4 Discrete(elements, leq=None).

A finite poset with explicit element set and a caller-supplied leq predicate. Useful for catalog-like enumerations where the order is not the default discrete order.

4.2 Antichains↩︎

The Antichain class wraps a set of poset elements together with the poset they live in. The constructor normalises the input by removing dominated points and deduplicating.

class Antichain:
    poset: Poset
    points: list

    @classmethod
    def of_bottom(cls, poset) -> "Antichain": ...
    @classmethod
    def empty(cls, poset) -> "Antichain": ...
    @classmethod
    def singleton(cls, poset, point) -> "Antichain": ...
    @classmethod
    def from_set(cls, poset, points) -> "Antichain": ...

    @classmethod
    def union_min(cls, poset, antichains) -> "Antichain": ...

    def filter_above(self, floor) -> "Antichain": ...
    def leq(self, other) -> bool: ...
    def eq(self, other) -> bool: ...
    def has_any_top(self) -> bool: ...
    def is_empty(self) -> bool: ...

The two key operations on antichains are union_min, a class method that takes the union of a collection of antichains and re-applies \(\mathsf{Min}\), and filter_above, which keeps only points dominating a given floor. Together they implement the body of the Kleene iteration in Theorem 1; the full method surface is documented in the reference card (Section 13.1.2.1).

5 Primitive Design Problem Types↩︎

Six concrete subclasses of DesignProblem are provided, each appropriate for a different style of relation between functionality and resources.

5.1 AlgebraicDP↩︎

For relations where each resource is a closed-form monotone function of the functionality.

AlgebraicDP(
    F: Ports,
    R: Ports,
    equations: dict[str, Callable[[dict], float]],
    name: str = "algebraic",
)

The equations dictionary maps each resource name to a callable taking a functionality dict and returning a scalar. h(f) returns a singleton antichain.

battery = AlgebraicDP(
    F=Ports({"capacity": Reals(unit="J")}),
    R=Ports({"mass": Reals(unit="kg")}),
    equations={"mass": lambda f: f["capacity"] / 1.8e6},
)

5.2 FunctionDP↩︎

For relations where the user supplies h directly as a function from functionality dicts to antichains. Use this when the relation is multi-valued (a true Pareto front) or has branching logic.

FunctionDP(
    F: Ports,
    R: Ports,
    h_fn: Callable[[dict], Antichain],
    name: str = "function",
)

F = Ports({"watts":    Reals(unit="W")})
R = Ports({"power_in": Reals(unit="W"),
           "thd":      Reals()})

def h_amp(f):
    watts = f["watts"]
    return Antichain.from_set(R, [
        {"power_in": watts / 0.25, "thd": 0.001},   # class-A
        {"power_in": watts / 0.90, "thd": 0.020},   # class-D
    ])

amp = FunctionDP(F, R, h_amp)

5.3 CatalogDP↩︎

For selecting from a finite catalog of implementations (motors, batteries, sensors).

CatalogDP(
    F: Ports,
    R: Ports,
    catalog: list[CatalogEntry],
    name: str = "catalog",
)

CatalogEntry(
    provides: dict[str, Any],
    costs: dict[str, Any],
    name: str = "",
)

For a request \(f\), \(h(f)\) is the antichain of entry.costs for every entry whose provides dominates \(f\) in \(F\). Entries that are dominated by another feasible entry are pruned by \(\mathsf{Min}\) in \(R\).

motor = CatalogDP(
    F=Ports({"torque": Reals(unit="N*m")}),
    R=Ports({"mass":   Reals(unit="kg"),
                    "cost":   Reals(unit="USD")}),
    catalog=[
        CatalogEntry(name="Tiny",   provides={"torque": 2.0},
                     costs={"mass": 0.2, "cost": 30.0}),
        CatalogEntry(name="Light",  provides={"torque": 8.0},
                     costs={"mass": 0.5, "cost": 200.0}),
        CatalogEntry(name="Heavy",  provides={"torque": 8.0},
                     costs={"mass": 0.8, "cost": 120.0}),
    ],
)
# motor.h({"torque": 7.0}) yields the antichain
# Antichain[(mass=0.5, cost=200), (mass=0.8, cost=120)]

5.4 ConstraintDP↩︎

A feasibility predicate plus a scalar cost, lifted to \(\mathsf{Min}\). Useful when the resource-functionality relation is most natural to express as is this choice good enough, and how much does it cost?.

ConstraintDP(
    F: Ports,
    R: Ports,
    sampler: Callable[[dict], Iterable[dict]],
    feasible: Callable[[dict, dict], bool],
    cost: Callable[[dict], dict],
    name: str = "constraint",
)

The sampler yields candidate implementations for the given functionality; feasible(impl, f) filters them; cost(impl) maps each feasible implementation into the resource poset, and the result is reduced by \(\mathsf{Min}\).

5.5 ODE_DP↩︎

Derives a monotone resource relation from a differential equation \(\dot{x} = \text{rhs}(x, t, f)\). Two modes are supported.

5.5.0.1 Final value

integrates the ODE from \(t = 0\) to \(t = t_{\text{end}}\) by explicit Euler and reports the extracted resources of the final state.

5.5.0.2 Steady state

solves \(\text{rhs}(x, \cdot, f) = 0\) for \(x\) by Newton iteration and reports the extracted resources of the steady-state value.

ODE_DP(
    F: Ports,
    R: Ports,
    rhs: Callable[[state, float, dict], state],
    extract: Callable[[state], dict],
    mode: str = "final_value",   # or "steady_state"
    t_end: float = 10.0,
    n_steps: int = 200,
    x0_fn: Callable[[dict], state] = None,   # default: lambda f: 0.0
    name: str = "ode",
)

F = Ports({"delta_T": Reals(unit="K")})
R = Ports({"power":   Reals(unit="W")})
k = 0.5   # heat-loss coefficient, W/K

heater = ODE_DP(
    F, R,
    rhs=lambda x, t, f: k * f["delta_T"] - x,
    extract=lambda x: {"power": float(x)},
    mode="steady_state",
    x0_fn=lambda f: 0.0,
)
# solve(heater, {"delta_T": 30.0}) -> Antichain[(power=15 W)]

5.6 UncertainDP↩︎

Brackets an unknown \(h\) with a known lower and upper bound. The optimistic mode (lower) yields a lower bound on the minimal resources; the pessimistic mode (upper) yields an upper bound. Section VII of [1] develops the theory, and [5] treats uncertainty in monotone co-design problems in full.

UncertainDP(
    F: Ports,
    R: Ports,
    lower: DesignProblem,
    upper: DesignProblem,
    mode: str = "upper",
    name: str = "uncertain",
)

with_mode("lower" | "upper") returns the bracket selected for solving.

F = Ports({"capacity": Reals(unit="J")})
R = Ports({"mass":     Reals(unit="kg")})

lo = AlgebraicDP(F, R, {"mass": lambda f: f["capacity"] / 2.0e6})
hi = AlgebraicDP(F, R, {"mass": lambda f: f["capacity"] / 1.6e6})

battery = UncertainDP(F, R, lower=lo, upper=hi)

# Solve both brackets and report the interval.
m_opt  = solve(battery.with_mode("lower"), {"capacity": 3.6e6})
m_pess = solve(battery.with_mode("upper"), {"capacity": 3.6e6})
# m_opt  -> Antichain[(mass=1.80 kg)]
# m_pess -> Antichain[(mass=2.25 kg)]

6 Composition Operators↩︎

The three operators of Section 2 are exposed as classes (Series, Parallel, Loop) with lowercase function aliases (series, par, loop) matching the paper’s notation. Each takes existing design problems and returns a new DesignProblem, so composites nest freely and are solved by the same solve entry point. Because the class of MCDPs is closed under all three, a composite is itself a monotone design problem. The compositional structure is the point of the framework, and is what the later literature builds on: modular co-design of vehicle control systems [11], structured co-design of embodied intelligence [10], and the categorical accounts of [14] and [16] all rest on it. [20] study what happens when monotonicity is dropped.

6.1 Series↩︎

series(dp1: DesignProblem, dp2: DesignProblem) -> Series
Series(dp1, dp2)

Requires dp1.R == dp2.F structurally (same component names and posets). Implements the series definition: for each point of \(h_1(f)\) run \(h_2\) and take the union \(\mathsf{Min}\).

battery = AlgebraicDP(
    F=Ports({"capacity": Reals(unit="J")}),
    R=Ports({"mass":     Reals(unit="kg")}),
    equations={"mass": lambda f: f["capacity"] / 1.8e6},
)
actuator = AlgebraicDP(
    F=Ports({"mass":     Reals(unit="kg")}),
    R=Ports({"power":    Reals(unit="W")}),
    equations={"power": lambda f: 10.0 * f["mass"] ** 2},
)
chain = series(battery, actuator)
# solve(chain, {"capacity": 3.6e6}) maps energy -> power directly.

6.2 Parallel↩︎

par(dp1: DesignProblem, dp2: DesignProblem) -> Parallel
Parallel(dp1, dp2)

Requires non-overlapping component names in \(F_1, F_2\) and in \(R_1, R_2\). The output \(F\) and \(R\) are the concatenated Portss. The antichain is the Cartesian product followed by \(\mathsf{Min}\).

propulsion = AlgebraicDP(
    F=Ports({"energy":      Reals(unit="J")}),
    R=Ports({"prop_mass":   Reals(unit="kg")}),
    equations={"prop_mass": lambda f: f["energy"] / 1.8e6},
)
sensor = AlgebraicDP(
    F=Ports({"sample_rate": Reals(unit="Hz")}),
    R=Ports({"sensor_pwr":  Reals(unit="W")}),
    equations={"sensor_pwr": lambda f: 0.05 * f["sample_rate"]},
)
combined = par(propulsion, sensor)
# combined.F has both energy and sample_rate;
# combined.R has both prop_mass and sensor_pwr.

6.3 Loop↩︎

loop(inner: DesignProblem, axis: str) -> Loop
Loop(inner, axis="...")

Closes the feedback loop on a single named axis that must appear in both inner.F and inner.R. The resulting DP has F = inner.F minus the axis and R = inner.R minus the axis. Calling h(f) triggers the Kleene fixed-point iteration described in Theorem 1.

inner = AlgebraicDP(
    F=Ports({"payload":       Reals(unit="kg"),
                    "battery_mass":  Reals(unit="kg")}),
    R=Ports({"battery_mass":  Reals(unit="kg")}),
    equations={
        "battery_mass":
            lambda f: 0.5 * (f["payload"] + f["battery_mass"]),
    },
)
drone = loop(inner, axis="battery_mass")
# drone.F = {"payload"}, drone.R = {}
# solve(drone, {"payload": 1.0}) iterates to the fixed point.

7 The Solver↩︎

All design problems, primitive or composite, are solved through a single entry point, solve. It evaluates the problem at a given functionality, running the Kleene fixed-point iteration of Theorem 1 whenever a feedback loop is present, and returns a SolveResult carrying the resulting antichain together with convergence and feasibility metadata. This section documents that entry point, the underlying iteration, scalar minimisation over the returned front, and the solver’s observability features.

7.1 Top-level entry point↩︎

solve(
    dp: DesignProblem,
    functionality: dict | None = None,
    max_iter: int = 200,
    *,
    trace: bool = False,
    verbose: int = 0,
    on_iteration: Callable[[TraceEntry], None] | None = None,
    start_from: SolveResult | Antichain | None = None,
    uncertainty: list[str] | None = None,
    n_samples: int = 1000,
    rng_seed: int | None = None,
) -> SolveResult

The keyword-only arguments trace, verbose, on_iteration, and start_from are described in Section 7.4; uncertainty, n_samples, and rng_seed in Section 8. The result is a dataclass with the following fields:

@dataclass
class SolveResult:
    antichain: Antichain      # the Pareto front (in dp.R)
    iterations: int           # Kleene steps taken (0 if no loop)
    status: str               # "converged", "max_iter", or "diverged"
    feasible: bool            # True if antichain has finite, nonempty points
    trace: list[TraceEntry]   # the per-step trace if trace=True, else None

    # converged is a read-only alias for status == "converged"

result = solve(drone, {"endurance":      300.0,
                       "extra_payload":  0.5,
                       "extra_power":    5.0},
               max_iter=200, trace=True)

print(result.status)         # "converged"
print(result.iterations)     # e.g. 17
print(result.feasible)       # True
print(result.antichain)      # Antichain[(total_mass=0.55 kg)]

# Scalarise: cheapest design under a custom cost function
best = minimize_cost(result, cost_fn=lambda r: r["total_mass"])
print(best)                  # {"total_mass": 0.5492}

7.2 The Kleene iteration↩︎

kleene_loop(
    loop_dp: Loop,
    f_outer: dict | None,
    max_iter: int = 200,
    *,
    trace: bool = False,
    verbose: int = 0,
    on_iteration: Callable[[TraceEntry], None] | None = None,
    start_from: Antichain | None = None,
    info_out: dict | None = None,
) -> Antichain

The body of the iteration is direct from Theorem 1. The implementation adds:

• a divergence cap of \(10^{30}\) that converts an iterate exceeding the cap to \(\top\), so floating-point overflow becomes infeasibility rather than inf/nan;

• try/except around the inner \(h\) to catch OverflowError, ValueError, and ZeroDivisionError, again mapping to \(\top\);

• three termination conditions: fixed point reached, antichain empty (infeasible), or every point’s loop axis equal to \(\top\) (provably infeasible).

7.3 Scalar minimisation over an antichain↩︎

minimize_cost(
    result: SolveResult,
    cost_fn: Callable[[dict], float],
) -> dict | None

Returns the single resource bundle in result.antichain that minimises cost_fn, or None if the result is infeasible. This is the scalarisation step that turns a Pareto front into a single engineering choice.

7.4 Observability: trace, verbose, callback, status↩︎

The solver can be made to report its progress in three independent ways, and its termination reason is exposed through a structured status field that is orthogonal to feasible.

7.4.0.1 The status field.

SolveResult.status is one of three strings:

"converged"

The iteration reached a fixed point (or provably hit \(\top\)) within max_iter steps. The answer in result.antichain is the algorithm’s verdict; feasible tells you whether that verdict is finite (a real design) or \(\top\) (no design exists).

"max_iter"

The iteration was cut off at max_iter. The current antichain may or may not be near a fixed point. Usually a sign to increase max_iter, or to look for an unintentional non-monotonicity in the model.

"diverged"

At least one numeric component crossed the divergence cap of \(10^{30}\) before the iteration could settle, and the solver stopped. Distinguishes pure numerical blow-up from clean \(\top\)-infeasibility.

The previous Boolean converged field is preserved as a backward-compatibility alias for status == "converged".

7.4.0.2 Structured trace.

Passing trace=True to solve populates result.trace with a list of TraceEntry dataclasses, one per Kleene step (plus the seed at iteration 0):

@dataclass
class TraceEntry:
    iteration: int           # 0 = seed, then 1, 2, ...
    antichain: Antichain     # snapshot at this step
    n_points:  int           # convenience: len(antichain)
    delta:     float | None  # max absolute change (numeric)
                             # or 0/1 (discrete); None at iter 0
    elapsed_ms: float        # wall time for this step alone

With trace=False (the default) the field is None, so the cost of unused tracing is zero.

7.4.0.3 Verbose printing.

The integer verbose controls live output:

verbose=0

(default) silent.

verbose=1

one summary line at the end, e.g. [solve] converged: 17 iters, |A|=1, total=2.0ms, feasible=True.

verbose=2

one line per iteration, showing the delta, the antichain size, and the step’s wall time. Useful for watching oscillating fixed points.

7.4.0.4 Iteration callback.

Passing on_iteration=callable registers a callback that receives each TraceEntry as soon as it is produced. Useful for live plots, custom logging, or interrupting an iteration that looks pathological:

def log_every_5(entry):
    if entry.iteration         print(f"iter {entry.iteration}: delta={entry.delta:.3e}")

solve(dp, f, on_iteration=log_every_5)

prev = None
results = []
for L in np.linspace(5.0, 30.0, 50):
    r = solve(dp, {"daily_load_kwh": float(L), ...},
              max_iter=400, start_from=prev)
    results.append(r)
    prev = r            # reuse the converged inner antichain

from codesign import viz

r = solve(drone, f, trace=True, max_iter=200)
ax = viz.plot_convergence(r)
ax.set_title(f"converged in {r.iterations} steps")

8 Uncertainty↩︎

Uncertainty in monotone co-design was introduced by [5], and the composable and distributional treatments of [13], [14], and [15] extend it. This section documents what the library implements, which is the set-based and stochastic part.

Two kinds of parameter uncertainty are supported, declared on Module instances and consumed by the same solve entry point. Both are independent of the model’s structure: a module can have set-based uncertainty, stochastic uncertainty, both, or neither.

8.1 Declaration↩︎

A Module carries optional attributes:

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}

    def __init__(self, specific_energy=2.0e6, efficiency=0.9):
        self.specific_energy = specific_energy
        self.efficiency = efficiency
        super().__init__()

    def h(self, f):
        # nominal product 2.0e6 * 0.9 = 1.8 MJ/kg delivered,
        # matching the canonical drone (examples 1/6/7)
        return {"mass": f["capacity"]
                       / (self.specific_energy * self.efficiency)}

b = Battery()
b.uncertain_set  = Box(...)          # deterministic, set-based
b.uncertain_dist = Stochastic(...)   # statistical

The uncertainty solver discovers every such module by walking the _codesign_modules attribute that System.build attaches to the returned DP. Before each solve, the module’s named parameters are reassigned; afterwards, their nominal values are restored. The user’s original Battery() instance is unchanged.

8.2 Set-based uncertainty (deterministic)↩︎

8.2.0.1 Box.

Axis-aligned interval product. Each parameter is declared as (lo, hi) or (lo, hi, direction), where direction is one of the four tokens "more_is_better", "more_is_worse", "less_is_better", "less_is_worse". With all directions declared, the worst case is a single corner read off in constant time:

b.uncertain_set = Box(
    specific_energy=(1.7e6, 2.3e6, "more_is_better"),
    efficiency=(0.83, 0.97, "more_is_better"),
)
# worst case = (specific_energy=1.7e6, efficiency=0.83)

When some directions are undeclared, all \(2^n\) endpoint combinations are evaluated and the worst is kept; cheap when \(n\) is modest.

8.2.0.2 Ellipsoid.

An \(n\)-D ellipsoid \((p - c)^\top \Sigma^{-1} (p - c) \leq 1\) in parameter space. The covariance \(\Sigma\) encodes both scale and correlation between parameters, so an ellipsoid is a more honest model than a box when parameters are believed to vary together:

b.uncertain_set = Ellipsoid(
    center={"specific_energy": 2.0e6, "efficiency": 0.9},
    cov=[
        [1.0e10, -2.0e3],
        [-2.0e3,  2.5e-3],
    ],
    params=["specific_energy", "efficiency"],
    directions={
        "specific_energy": "more_is_better",
        "efficiency":      "more_is_better",
    },
)

With all directions declared, the worst case is the boundary point \(c + L\,u^\star\), where \(L\) is the Cholesky factor of \(\Sigma\) and \(u^\star\) is the unit vector in the direction of badness. Otherwise, the boundary is sampled (controlled by boundary_samples) and the worst sample is returned.

8.2.0.3 2D conveniences.

Disk(center, radius, ...) and Circle(center, radius, ...) reduce to the isotropic ellipsoid \(\Sigma = r^2 I\). For monotone systems with declared directions, the two are equivalent: the worst case lies on the boundary either way.

8.3 Stochastic uncertainty↩︎

8.3.0.1 Stochastic.

A joint distribution built from named scipy-stats frozen marginals plus a copula:

from scipy import stats

b.uncertain_dist = Stochastic(
    marginals={
        "specific_energy": stats.uniform(loc=1.7e6, scale=0.6e6),
        "efficiency":      stats.uniform(loc=0.83, scale=0.14),
    },
    copula=GaussianCopula(correlation=[[1.0, 0.4],
                                        [0.4, 1.0]]),
)

Sampling is by the standard copula–marginal trick: the copula generates correlated uniforms in \([0,1]^d\); the marginal ppfs map them to the named parameter axes.

8.3.0.2 Copulas.

Independence() is the default (coordinatewise uniform). GaussianCopula(correlation) uses the Gaussian copula with the given correlation matrix; samples are drawn by Cholesky factorisation followed by \(\Phi\) on each coordinate. Subclassing Copula requires only one method, sample_uniform(n, d, rng), so other copulas (e.g.Clayton, Gumbel) drop in without further plumbing.

8.4 Querying↩︎

The same solve entry point accepts an uncertainty list and returns a richer object:

result = solve(
    dp, f,
    uncertainty=["worst_case", "mean", "p95", "cvar95", "samples"],
    n_samples=1000,
    rng_seed=42,
)

result.worst_case        # SolveResult at the worst point of the set
result.mean              # dict[r_port -> mean across MC samples]
result.p95               # dict[r_port -> 95th percentile]
result.cvar95            # dict[r_port -> CVaR at 95result.samples           # list[Antichain], one per MC sample
result.feasibility_rate  # fraction of feasible MC samples

The allowed labels are:

"worst_case"

The deterministic worst over every module’s uncertain_set. One canonical solve at the worst-case parameter point.

"mean", "p95", "cvar95"

Marginal statistics over Monte Carlo samples, drawn from every module’s uncertain_dist. The CVaR (conditional value at risk) at the 95% level is the mean of the worst 5% of samples.

"samples"

The raw antichain per MC sample, for the user to aggregate however they like.

For a single solve call you typically observe the ordering \[\text{nominal} < \text{mean} < \text{p95} < \text{CVaR95} < \text{worst\_case},\] which separates the four standard ways of reporting an answer "under uncertainty." Notebooks 11 and 12 visualise this for the example-7 drone.

bat = Battery()  # the Module subclass from example~7
bat.uncertain_set = Box(
    specific_energy=(1.7e6, 2.3e6, "more_is_better"),
    efficiency=(0.83, 0.97, "more_is_better"),
)
# ... wire `bat` into the drone system ...

result = solve(drone, f, uncertainty=["worst_case"])
worst_mass = list(result.worst_case.antichain.points)[0]["total_mass"]
# worst_mass = 0.5668 kg, versus 0.5492 kg nominal

from scipy import stats

bat.uncertain_dist = Stochastic(
    marginals={
        "specific_energy": stats.uniform(loc=1.7e6, scale=0.6e6),
        "efficiency":      stats.uniform(loc=0.83, scale=0.14),
    },
    copula=GaussianCopula(correlation=[[1.0, 0.4],
                                        [0.4, 1.0]]),
)

res = solve(drone, f,
            uncertainty=["mean", "p95", "cvar95", "samples"],
            n_samples=1000, rng_seed=42)
# res.mean["total_mass"]    -> 0.5506 kg
# res.p95["total_mass"]     -> 0.5632 kg
# res.cvar95["total_mass"]  -> 0.5647 kg
# len(res.samples)          -> 1000

9 Online Learning↩︎

When a co-design problem has many discrete candidates (catalog entries, robot types, component families) and each candidate’s inner solve is non-trivial, the naive “evaluate every entry” strategy is wasteful. The codesign.online module implements the compositional online learner of [18]: maintain history-dependent bounds on each candidate’s inner-solve output, evaluate the most promising one under the lcb (lower-confidence-bound) rule by default, then prune any candidate whose lower bound is already dominated by the incumbent.

9.1 Optimistic evaluators↩︎

The bound machinery is encapsulated in a small class hierarchy. OptimisticEvaluator is the abstract base; subclasses implement bound(candidate) returning a pair of dicts (lower, upper) mapping each numeric R component to its current lower and upper bound at the queried feature point. Tighter bounds prune more candidates; looser bounds are always safe.

Three concrete evaluators are provided:

9.1.0.1 Monotonicity.

MonotonicityEvaluator(features, r_components) assumes the inner-solve output is component-wise monotone in the named features. Given an observation at feature point \(f_0\) with antichain-min \(R_0\), any candidate whose features dominate \(f_0\) has resources \(\geq R_0\) (lower bound); any candidate dominated by \(f_0\) has resources \(\leq R_0\) (upper bound). Aggressive when applicable. Only correct if the monotonicity assumption genuinely holds: choosing a derived feature (e.g.cost_per_capacity = unit_cost / (speed * payload)) is often what makes this safe in practice.

9.1.0.2 Lipschitz.

LipschitzEvaluator(features, r_components, L) assumes \(|h(c_1) - h(c_2)| \leq L \cdot \|features(c_1) - features(c_2)\|\) (Euclidean distance, per R component). Each observation tightens the bound by a cone of slope \(L\). The constant can be a scalar (same for every R component) or a dict per component. Safe default: with a sensible \(L\) the bound never prunes a Pareto-optimal candidate.

9.1.0.3 Linear-parametric.

LinearParametricEvaluator(features, r_components, min_obs=3, prior_box=None, noise_bound=0.0, solver="highs") implements the certified optimistic bound of [18] (Section V-C3, eqs.–28 and Lemma V.5). It assumes each resource coordinate is an exact affine function of the features, \(\text{req}_k(c) = \phi(c)^\top \theta_k^*\) with \(\phi(c) = [1, features(c)]\) (the paper’s noiseless linear model). Rather than fit a single least-squares line, it maintains the whole confidence polytope \(\Theta(H)\) — the subset of the prior box \(\Theta_0\) consistent with every observation taken as a linear equality — and returns, per resource coordinate, the coordinatewise minimum of the predicted resource over that polytope, \([\text{req}_\text{opt}(i, H)]_k = \min_{\theta \in \Theta(H)} \phi(i)^\top \theta\). Each such minimum is one linear program, solved with scipy.optimize.linprog (the solver kwarg selects the method). Because the true parameter \(\theta_k^*\) is feasible for the LP, its optimum is always \(\leq \phi(i)^\top \theta_k^* = \text{req}_k(i)\): the returned value is a guaranteed lower bound (the optimism guarantee of Lemma V.5), so — unlike the former OLS \(\pm\) confidence-band heuristic — it can never wrongly eliminate a Pareto-optimal candidate. No upper bound is certified, so the returned upper bound stays at \(+\infty\). prior_box sets \(\Theta_0\) (a per-parameter box) to keep the LP bounded while the fit is under-determined: None leaves every parameter unbounded, which is always safe but yields the trivial bound until enough observations pin \(\theta\) down; noise_bound\({} > 0\) relaxes each observation equality into a two-sided band of that half-width, a documented extension for bounded observation noise. The old confidence kwarg is deprecated and ignored (passing it emits a DeprecationWarning). Net effect: this evaluator is both safe and effective when the affine assumption genuinely holds, and degrades safely — to the no-information bound, never to a wrong elimination — when it does not.

9.1.0.4 Gaussian process.

GaussianProcessEvaluator(features, r_components, length_scale=0.3, sigma_f=1.0, noise=1e-3, confidence=2.0, min_obs=3) fits a zero-mean GP with an RBF kernel and bounds new queries by a confidence band on the predictive standard deviation. Implemented in pure numpy. The right tool when the response surface has feature interactions or local nonlinearity: the GP still produces informative (if uncertified) bounds there, whereas the certified linear-parametric evaluator correctly degrades to the no-information bound once the affine assumption fails. When the map is affine, prefer the linear-parametric evaluator — its bound is then both certified-safe and tight.

9.1.0.5 Subclassing.

Custom evaluators only need to override bound(candidate); the base class handles observation recording and the default fallback bound of \((0, +\infty)\).

9.2 The online solver↩︎

result = solve_online(
    candidate_fn,         # candidate_fn(candidate) -> DP
    functionality,        # outer F vector
    *,
    candidates,           # list of feature dicts
    evaluator,            # OptimisticEvaluator instance
    budget=None,          # max inner solves; None = unbounded
    max_iter=200,         # forwarded to each inner solve
    verbose=0,
    warm_start=None,      # seed indices or n for farthest-point
    picker="lcb",         # "lcb", "ucb", "random", or callable
)

Algorithm (a simplified Algorithm 2 from Alharbi et al.):

1. Bound every remaining candidate via the evaluator.

2. Eliminate every candidate whose lower bound is already dominated by the current incumbent antichain.

3. Pick the most promising survivor with the configured picker (lcb by default) on its lower-bound sum.

4. Run the inner solve at that candidate; observe the result in the evaluator; merge into the incumbent via union_min.

5. Repeat until the candidates are exhausted or the budget is hit.

The returned OnlineResult exposes:

antichain

Min over the evaluated, surviving candidates.

n_evaluated, n_eliminated, n_candidates

Bookkeeping for the elimination cascade.

evaluated_ids, eliminated_ids, incumbent_ids

Lists of indices into the original candidate list.

history

Per-iteration log: {pick, antichain, remaining, evaluated, eliminated}.

def candidate_fn(robot):
    capacity = robot["speed"] * robot["payload"]
    return AlgebraicDP(F, R, {
        "total_cost":   lambda f, cap=capacity, uc=robot["unit_cost"]:
                        (f["target_throughput"] / cap) * uc,
        "total_energy": lambda f, ek=robot["energy_per_km"]:
                        f["target_range"] * ek * 24.0,
    })

ev = LipschitzEvaluator(
    features=["speed", "payload", "unit_cost", "energy_per_km"],
    r_components=["total_cost", "total_energy"],
    L={"total_cost": 300.0, "total_energy": 30.0},
)

result = solve_online(
    candidate_fn, mission,
    candidates=robot_catalog,        # list of 200 feature dicts
    evaluator=ev,
)
# result.n_evaluated   = 192 / 200
# result.n_eliminated  = 8
# result.antichain     = 5 Pareto-optimal robot types

9.3 Warm-start, picker strategies, and the GP evaluator↩︎

The basic online solver has three small extensions that significantly expand the tuning space without complicating the core algorithm.

9.3.0.1 Warm-start.

A new warm_start argument to solve_online pre-populates the evaluator with seed observations before the picker takes over. It accepts either a list of candidate indices (manually picked corner runs) or an integer \(n\) that triggers a greedy farthest-point heuristic over the feature space. The motivation is concrete: the MonotonicityEvaluator’s lower bounds only tighten for candidates whose features dominate every observation in the partial order. Without observations at the low-feature corner where the Pareto front typically lives, the picker explores indiscriminately and the evaluator remains uninformative throughout the budget. In example 16, evaluating the example with four corner warm-start runs lifts Monotonicity’s Pareto recovery from zero to one of four classes; a hybrid with one of the other evaluators is the natural next step for cases where this is not enough.

9.3.0.2 Pluggable picker strategies.

A new picker argument selects the candidate-scoring rule. Built-in options: "lcb" (the previous default, minimises the sum of lower-bound components, pure exploitation of the optimistic estimate); "ucb" (lower bound minus a \(\kappa \cdot\) (upper \(-\) lower) exploration bonus, tunable via the tuple form ("ucb", {"kappa": 1.0})); and "random" (uniform baseline for comparing the value of structural priors). Custom callables are accepted and receive (lower, upper, r_components, **kwargs). The strategies that matter in practice are LCB for confident priors and UCB with \(\kappa \approx 0.3\) to \(0.5\) when the prior is uncertain enough that pure exploitation gets stuck in a local region.

9.3.0.3 Gaussian-process evaluator.

A new GaussianProcessEvaluator class uses a zero-mean GP with an RBF kernel, implemented in pure numpy. Hyperparameters are length_scale (default \(0.3\), suitable for features in \([0, 1]\)), sigma_f (default \(1.0\), rescaled per-output by the empirical observation standard deviation), noise (default \(10^{-3}\) jitter), confidence (default \(2.0\sigma\) for roughly \(95\%\) coverage), and min_obs (default \(3\)). The GP is more expressive than the linear-parametric evaluator and captures local nonlinearity that a global linear fit misses. Since the linear-parametric evaluator became the certified confidence-polytope bound the two occupy clearly separated niches. The certified LinearParametricEvaluator is the tool of choice when the resource map is genuinely affine in the features: on the linear fleet catalogue of example 14 it recovers all five Pareto types (where the former OLS-band version wrongly dropped one). But the example 16 bioprocess effect model is markedly nonlinear (temperature U-shapes, a power-law failure rate), and there the certified bound correctly refuses to over-claim — it degrades to the no-information bound and recovers zero of the four Pareto classes rather than eliminate a candidate it cannot rule out. The GP, willing to extrapolate a smooth surrogate, stays informative on that same grid (recovering one of the four classes). Rule of thumb: certified linear-parametric when you trust the affine assumption (safe and effective), GP when you do not.

10 Visualisation↩︎

The codesign.viz module exposes four helpers built on matplotlib and a small GraphViz emitter:

plot_antichain(result, axes)

Scatter the antichain on the chosen 2D or 3D axes. Accepts a SolveResult, an UncertaintyResult (uses its worst case), or a bare Antichain. The 2D version optionally shades each point’s dominated region for visual clarity.

plot_convergence(result)

Semilog plot of the Kleene delta against iteration index. Accepts a SolveResult with a trace or a trace list directly. Useful for diagnosing oscillating fixed points and tuning max_iter.

plot_uncertainty(unc_result, port, nominal=None)

Histogram of the Monte Carlo samples for the named R port, with the nominal, mean, p95, and CVaR95 marked as vertical lines.

to_dot(dp, name="codesign")

GraphViz dot string for the system structure: modules as boxes, outer ports as labelled rectangles, constraint connections as edges. Render with dot -Tpng or paste into an online viewer.

All four accept an optional ax=... argument so they compose into larger figures. The module is imported as from codesign import viz.

10.1 Block diagrams of Systems↩︎

The to_dot function above emits a plain string with module-level wiring. For richer visualisation, the library ships a dedicated diagram renderer in codesign.diagram that produces Simulink-style block diagrams:

One box per subsystem

with the F (functionality) ports listed on the left and the R (resource) ports on the right. Each port is individually addressable, so edges attach to specific ports rather than to the module as a whole.

Outer F and outer R as separate nodes

on the diagram’s left and right margins. The system inputs and outputs are visible, not buried inside the constraint list.

Port-level edges

whenever the constraint was written in the operator-overloaded form (module1.r_port >= module2.r_port * ...). Constraints defined by callables that cannot be introspected are rendered with a dashed edge from a small “\(\lambda\)” marker, so they remain visible.

Cycle detection

via Tarjan’s strongly-connected-components algorithm on the module-level constraint graph. Edges within any cycle of size \(\geq 2\) are coloured amber, making the Kleene-iteration loop visible at a glance.

The renderer returns a graphviz.Digraph object, which can be displayed inline in a Jupyter notebook or written to SVG / PDF / PNG via dot.render(filename, format="svg"):

from codesign import draw_system

dot = system.draw_diagram()        # convenience method on System
# or equivalently:
dot = draw_system(system, rankdir="LR", highlight_cycles=True)

dot.render("bioprocess", format="svg", cleanup=True)

Optional dependency: install via pip install codesign-mcdp[diagram] (which pulls the graphviz Python package), and ensure the dot binary is on PATH (apt-get install graphviz on Debian / Ubuntu, brew install graphviz on macOS).

from codesign import draw_system
dp = make_bioprocess(cell_line=CHO_K1, glucose_setpoint_mm=8.0,
                     annual_demand_kg=100.0)
draw_system(dp, name="bioprocess").render("bioprocess",
                                          format="svg",
                                          cleanup=True)

from codesign import viz
import matplotlib.pyplot as plt

result = solve(vehicle, {"payload": 2.0, "mission_energy": 2.0e5})
viz.plot_antichain(result, axes=["total_mass", "total_cost"])
plt.show()

res = solve(drone, f,
            uncertainty=["mean", "p95", "cvar95", "samples"],
            n_samples=1000, rng_seed=42)
ax = viz.plot_uncertainty(res, port="total_mass",
                          nominal=0.5492, bins=30)

dot = viz.to_dot(microgrid, name="microgrid")
# Render externally:
# echo "<dot output>" | dot -Tpng > microgrid.png

11 The MCDP Builder↩︎

The MCDP class in codesign.mcdpl provides an MCDPL-style declarative syntax for a single self-contained design problem.

with MCDP("name") as m:
    m.provides("functionality_name", unit="...")    # outer F port
    m.requires("resource_name",     unit="...")    # outer R port

    m.constraint("resource_name", lambda f: <expr>)   # closed-form
    # or
    m.rule(lambda f: Antichain.from_set(...))         # multi-valued

    m.loop_on("axis_name")    # optional, may be repeated

dp = m.build()

The builder emits an AlgebraicDP or FunctionDP internally, then wraps it in Loops as requested by loop_on. Multiple constraint calls on the same target are joined (max). Mostly useful when the model has a small number of variables and reads naturally as a list of \(\geq\) inequalities, mirroring the paper’s mcdp { ... } blocks.

with MCDP("drone") as m:
    m.provides("endurance",     unit="s")
    m.provides("extra_payload", unit="kg")
    m.provides("extra_power",   unit="W")
    m.provides("battery_mass",  unit="kg")   # loop axis: in F and R
    m.requires("battery_mass",  unit="kg")
    m.requires("report_mass",   unit="kg")   # mirror for visibility

    def battery_mass_eq(f):
        return (f["extra_power"]
                + 10.0 * (9.81 * (f["battery_mass"]
                                  + f["extra_payload"])) ** 2
                ) * f["endurance"] / 1.8e6

    m.constraint("battery_mass", battery_mass_eq)
    m.constraint("report_mass",  battery_mass_eq)
    m.loop_on("battery_mass")

drone = m.build()

12 The System Builder↩︎

For larger designs the System class in codesign.system provides modular composition with named subsystems and algebraic connection constraints. This is the recommended way to build non-trivial co-design models. Two equivalent surface syntaxes are available: the operator-overloaded form (recommended for new code) and the legacy lambda form (still supported).

12.1 Operator-overloaded syntax (recommended)↩︎

provides, requires, and add each return a port handle. Arithmetic operators on handles build expression trees lazily; the >= operator registers a constraint with the parent system.

from codesign import Module, Reals, System, solve

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def h(self, f):
        return {"mass": f["capacity"] / 1.8e6}

class Actuator(Module):
    F = {"lift_force": Reals(unit="N")}
    R = {"power":      Reals(unit="W")}
    def h(self, f):
        return {"power": 10.0 * f["lift_force"] ** 2}

sys = System("drone")
endurance     = sys.provides("endurance",     unit="s")
extra_payload = sys.provides("extra_payload", unit="kg")
extra_power   = sys.provides("extra_power",   unit="W")
total_mass    = sys.requires("total_mass",    unit="kg")

battery  = sys.add("battery",  Battery())
actuator = sys.add("actuator", Actuator())

battery.capacity    >= (actuator.power + extra_power) * endurance
actuator.lift_force >= 9.81 * (battery.mass + extra_payload)
total_mass          >= battery.mass + extra_payload

drone = sys.build()

The four port kinds are:

• Outer F (from provides): can appear on the RHS of constraint expressions. Cannot be a constraint target.

• Outer R (from requires): can be a constraint target. Cannot appear in expression RHS.

• Module F (handle.f_port): can be a constraint target. Cannot appear in expression RHS (its value is determined by the constraints on it, not by the current iteration state).

• Module R (handle.r_port): can appear in expression RHS. Cannot be a constraint target (its value is determined by the module’s own h, not externally).

Misusing a port (e.g.trying to put a module F port on the RHS, or constraining an outer F) raises immediately with a message naming the port and explaining the rule. F and R port names within a single subsystem must be disjoint.

The supported operators are +, -, *, /, **, unary -, along with the helper functions codesign.sqrt, exp, and log for expressions that need transcendental terms. For demands that cannot be expressed algebraically, the legacy lambda form remains available and can be mixed freely with the operator form (see below).

12.2 Legacy lambda syntax↩︎

sys.constrain("battery.capacity",
              lambda x: (x["actuator.power"] + x["extra_power"]) * x["endurance"])
sys.constrain("actuator.lift_force",
              lambda x: 9.81 * (x["battery.mass"] + x["extra_payload"]))
sys.constrain("total_mass",
              lambda x: x["battery.mass"] + x["extra_payload"])

The argument x is a dict carrying outer F values under their bare names (x["endurance"]) and subsystem R ports under their dotted names (x["battery.mass"]). Both syntaxes compile to the same internal representation; the operator form is purely sugar.

12.3 Class-based modules↩︎

The Module base class lets a design problem be defined declaratively:

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}

    def __init__(self, specific_energy=1.8e6):
        self.specific_energy = specific_energy
        super().__init__()

    def h(self, f):
        return {"mass": f["capacity"] / self.specific_energy}

A Module subclass declares its \(F\) and \(R\) ports as class-level dicts and overrides h(self, f). The constructor wires everything into the underlying DesignProblem machinery. Parameterised modules are written by overriding __init__ and calling super().__init__() after setting any instance attributes. h may return a dict (treated as a singleton antichain), a list of dicts (multi-valued antichain), or an Antichain directly.

sys = System("drone")
endurance     = sys.provides("endurance",     unit="s")
extra_payload = sys.provides("extra_payload", unit="kg")
extra_power   = sys.provides("extra_power",   unit="W")
total_mass    = sys.requires("total_mass",    unit="kg")

b = sys.add("battery",  Battery())
a = sys.add("actuator", Actuator())

b.capacity   >= (a.power + extra_power) * endurance
a.lift_force >= 9.81 * (b.mass + extra_payload)
total_mass   >= b.mass + extra_payload

drone = sys.build()
r = solve(drone, {"endurance": 300.0,
                  "extra_payload": 0.5,
                  "extra_power": 5.0})
# r.antichain -> Antichain[(total_mass=0.55 kg)]

12.4 Semantics↩︎

The build method emits a Loop whose axis bundles every subsystem’s resource poset: \[\texttt{\_\_modules\_\_} = \prod_{m \in \mathrm{modules}} R_m.\] The inner DP performs one Kleene step:

1. Read the current bundle estimate to populate \(\mathtt{ctx[m.r]}\) for every subsystem \(R\) port.

2. For each subsystem \(m\), evaluate the constraint demands on its \(F\) ports to obtain \(f_m\), then call \(m.h(f_m)\) to obtain an antichain.

3. Take the Cartesian product across subsystems; for each combination, evaluate the outer-\(R\) constraint demands to produce a candidate point.

4. Return the resulting antichain.

The outer Kleene iteration converges over all subsystem \(R\) ports simultaneously.

12.5 Validation↩︎

build performs the following checks:

• every subsystem \(F\) port has at least one constrain target;

• every outer \(R\) has at least one constrain target;

• every constraint target refers to a known module \(F\) port or a declared outer \(R\);

• module names do not clash with the reserved axis name __modules__ and contain no ..

If any check fails, build raises ValueError with a message identifying the missing or invalid item.

13 API Reference↩︎

This section is the complete, per-module reference for the public API exported by import codesign. The six subsections below follow the architecture of the library: order structures; design problems and their solver; the two builders and the rendering helpers; the uncertainty and online-learning layers; the temporal planners; and the vector-state and closed-loop layers. Within each group one subsubsection documents each source module, and each public symbol is a reference card giving its purpose, signature, a minimal runnable snippet, its parameters, and the exceptions it raises. Cards are cross-linked to the tutorial chapters (Sections 2–16) and the worked examples (Section 14); the condensed export list in Section 20 remains a one-page index. All reference labels use the ref:<module>:<symbol> scheme.

13.1 Order structures↩︎

The two foundational modules: posets supplies the functionality and resource spaces, and antichains the Pareto-front value type returned by every design problem. See Section 2 for the underlying order theory.

13.1.1 posets↩︎

The codesign.posets module supplies the partially ordered sets that serve as functionality spaces \(F\) and resource spaces \(R\). Every concrete poset derives from the abstract Poset base and honours its contract: a least element bottom() that seeds the Kleene iteration, a greatest element top() that marks infeasibility, and a join() that returns the least upper bound used to combine constraints. The mathematical background is given in Section 2.

13.1.1.1 Poset

Abstract base class defining the poset contract. It is not instantiable directly; a concrete subclass must override leq, bottom and top, and inherits working defaults for the derived operators. The four shipped subclasses (Reals, Naturals, Ports, Discrete) cover essentially every modelling need.

Signature: class Poset(ABC) with class attribute name: str. The derived operators are exercised here through a concrete subclass, since Poset itself cannot be instantiated:

from codesign import Reals
P = Reals(unit="kg")
# Derived operators are provided by the Poset base class:
print(P.eq(1.0, 1.0))          # True
print(P.lt(1.0, 2.0))          # True
print(P.comparable(1.0, 2.0))  # True (a chain: all comparable)
print(P.is_bottom(0.0))        # True

The contract methods (the abstract three plus the inherited defaults):

Parameters: the constructor takes no arguments; name is a display label overridden by subclasses.

Errors: instantiating Poset directly raises TypeError (abstract methods unimplemented). The default join raises NotImplementedError when a and b are incomparable (neither leq holds); chain and product subclasses override join so it never fails on them.

13.1.1.2 Reals

The non-negative reals \(\mathbb{R}_{\geq 0}\) augmented with \(+\infty\), forming a CPO. This is the scalar chain for any real-valued port (mass, cost, energy). Bottom is 0.0, top is math.inf.

Signature: Reals(unit: str = "", name: str = "R+").

from codesign import Reals
R = Reals(unit="kg")
print(R.bottom(), R.top())   # 0.0 inf
print(R.leq(0.5, 1.2))       # True
print(R.join(0.5, 1.2))      # 1.2
print(R.format(0.56))        # 0.56 kg
print(R.format(R.top()))     # the top marker
print(R.name)                # R+[kg]

Parameters:

• unit — display-only unit string; travels through composition so format can print 0.56 kg. Never used in any algebraic check.

• name — display name; if left at the default and a unit is given, becomes R+[unit] automatically.

Deviations from the Poset contract: overrides is_top (direct \(+\infty\) test), join (returns max) and format (prints \(\top\) for \(+\infty\), else the value with its unit).

Errors: none raised directly.

13.1.1.3 Naturals

The non-negative integers \(\mathbb{N}\) augmented with a top (\(+\infty\)), forming a CPO. Mirrors Reals for integer-valued quantities such as part counts or bit budgets. Bottom is 0, top is math.inf.

Signature: Naturals(unit: str = "", name: str = "N+").

from codesign import Naturals
N = Naturals(unit="parts")
print(N.leq(3, N.top()))   # True
print(N.join(2, 5))        # 5
print(N.is_top(N.top()))   # True
print(N.bottom())          # 0
print(N.format(3))         # 3 parts
print(N.format(N.top()))   # ⊤

Parameters: unit and name as for Reals. format mirrors Reals: it appends the unit (so N.format(3) prints 3 parts) and renders the top as \(\top\); finite values keep their integer form (3.0 prints as 3). When unit is set and the default name is kept, the name is decorated as N+[unit].

Deviations from the Poset contract: overrides leq to treat \(+\infty\) as strictly above every finite natural (and equal only to itself), is_top, join, and format. Finite inputs are coerced with int().

Errors: none raised directly.

13.1.1.4 Ports

A named product poset: a typed bundle of named component posets, ordered component-wise. This is the workhorse type; every design problem’s \(F\) and \(R\) is a Ports. Elements are Python dicts mapping each component name to a value in its factor poset. The order is \(x \leq y \iff x[k] \leq y[k]\) for every component \(k\); bottom and top are the dicts of all component bottoms and tops.

Signature: Ports(components: Mapping[str, Poset]).

from codesign import Ports, Reals
F = Ports({"capacity": Reals(unit="J")})
x = F.make(capacity=3.6e6)
print(F.leq(x, F.top()))    # True
print(F.is_top(F.top()))    # True
print(F.any_top({"capacity": F.top()["capacity"]}))  # True
print(F.format(x))          # (capacity=3.6e+06 J)
print(list(F.keys()))       # ['capacity']

Parameters:

• components — non-empty {name: poset} mapping. Copied into a fresh dict, so later mutation of the argument cannot corrupt the instance. A component poset may itself be a Ports (nesting is allowed).

Beyond the contract, Ports adds three public methods:

Note the distinction between is_top (every component at top — total infeasibility) and any_top (at least one component at top — enough to disqualify a candidate in the antichain machinery).

Errors: the constructor raises ValueError if components is empty. make raises ValueError if any declared component is missing from the keyword arguments, and ValueError if any unknown component name is supplied. leq, join and format may raise KeyError if an element dict lacks a declared component.

13.1.1.5 Discrete

A finite poset over an explicit element list with a caller-supplied order predicate. Useful for enumerated design spaces (part numbers, operating modes). The default order is equality: nothing is below anything but itself.

Signature: Discrete(elements, leq_fn=None, name: str = "D").

from codesign import Discrete
modes = Discrete(["idle", "cruise", "boost"], name="mode")
print(modes.leq("idle", "idle"))    # True
print(modes.leq("idle", "cruise"))  # False (equality order)
# bottom()/top() raise ValueError: no canonical extremum.

Parameters:

• elements — iterable of elements; stored as a list.

• leq_fn — order predicate leq_fn(a, b) -> bool; defaults to a == b (the discrete order).

• name — display name.

Errors: bottom() and top() each raise ValueError, because an arbitrary discrete order has no canonical least or greatest element; subclass and override them, or wrap in a poset that adds explicit extrema.

13.1.1.6 NamedProduct

Backward-compatible alias: NamedProduct = Ports. Existing code importing NamedProduct continues to work unchanged; new code should prefer Ports (see 13.1.1.4). The two names are the same object.

from codesign import NamedProduct, Ports
print(NamedProduct is Ports)   # True

Errors: as for Ports.

13.1.2 antichains↩︎

The codesign.antichains module provides the Antichain type, the return value of every design problem’s resource map \(h : F \to \mathsf{A}[R]\). Antichains generalise the Pareto front to an arbitrary partial order; their lattice structure (order, join by union followed by \(\mathsf{Min}\)) is developed in Section 2.

13.1.2.1 Antichain

A finite set of mutually incomparable elements of a poset, paired with the poset they live in. Construction is self-normalising: dominated points and duplicates are discarded, so the constructor is idempotent under \(\mathsf{Min}\). Instances are iterable and support len() and truthiness (empty is falsey).

Signature: Antichain(poset: Poset, points: Iterable = ()). The class-method constructors are the primary construction API; the bare constructor is rarely called directly.

from codesign import Antichain, Ports, Reals
R = Ports({"mass": Reals(unit="kg"), "cost": Reals(unit="$")})
pts = [{"mass": 1.0, "cost": 5.0},
       {"mass": 2.0, "cost": 3.0},
       {"mass": 1.5, "cost": 6.0}]   # dominated by the first
A = Antichain.from_set(R, pts)
print(len(A))    # 2 (dominated point dropped)
B = Antichain.singleton(R, {"mass": 0.5, "cost": 7.0})
U = Antichain.union_min(R, [A, B])
print(len(U))    # 3

Parameters:

• poset — the underlying Poset (typically a Ports) supplying leq, eq and, where present, any_top.

• points — iterable of poset elements; normalised to an antichain on insertion.

Public methods (class-method constructors marked cm):

The antichain order is \(A \leq A'\) iff \(\mathord{\uparrow} A \supseteq \mathord{\uparrow} A'\): every point of other is dominated by some point of self. union_min and filter_above together implement the body of the Kleene iteration (Section 2); has_any_top uses the poset’s any_top when available (as on Ports), otherwise falls back to is_top.

Errors: none raised directly. Operations delegate to the underlying poset, so a malformed point may surface a KeyError from Ports.leq rather than an antichain-level exception.

13.2 Design problems and solving↩︎

The six primitive design-problem types (dp), the scalar primitives built on them, the three composition operators, and the Kleene solver that computes fixed points. The tutorial treatment is in Sections 5–7.

13.2.1 dp↩︎

The dp module defines the central abstraction, the DesignProblem, together with the concrete primitives that realise it in closed form, from a catalogue, from a constraint, or from an ODE. Every design problem exposes a single method h(f) returning the antichain of minimal resources able to deliver functionality f.

13.2.1.1 DesignProblem

The abstract base class for all design problems: a monotone relation from a functionality poset F to an antichain in a resource poset R. Subclasses must implement h; the class is never instantiated directly.

class DesignProblem(ABC):
    F: Poset            # functionality space
    R: Poset            # resource space
    name: str = "dp"
    def h(self, f) -> Antichain: ...   # abstract

Concrete instances expose h, F, R and a readable __repr__:

from codesign import AlgebraicDP, Ports, Reals
dp = AlgebraicDP(
    F=Ports({"speed": Reals()}),
    R=Ports({"power": Reals(unit="W")}),
    equations={"power": lambda f: 2.0 * f["speed"]},
)
print(repr(dp))              # DP(algebraic: speed:R+ -> A[...])
print(dp.h({"speed": 3.0}))  # Antichain[(power=6 W)]

• F — Poset — the functionality space.

• R — Poset — the resource space.

• h(f) — returns the minimal-resource Antichain for f.

Errors: being an ABC with an abstract h, direct instantiation raises TypeError; subclasses raise their own (below).

13.2.1.2 AlgebraicDP

A functional design problem whose every resource is a closed-form (monotone) expression in the functionality. The relation is a singleton antichain. See Section 5 for the tutorial treatment. Monotonicity is the caller’s responsibility and is not verified.

AlgebraicDP(F: Poset, R: Ports,
            equations: Mapping[str, Callable | Any],
            name: str = "algebraic")

from codesign import AlgebraicDP, Ports, Reals
motor = AlgebraicDP(
    F=Ports({"torque": Reals(unit="Nm")}),
    R=Ports({"mass": Reals(unit="kg")}),
    equations={"mass": lambda f: 0.5 * f["torque"] + 1.0},
)
print(motor.h({"torque": 4.0}))  # Antichain[(mass=3 kg)]

• R — Ports — named resource bundle (required).

• equations — one entry per resource port, a callable f_dict -> value or a constant.

Errors: raises TypeError if R is not a Ports; raises ValueError if any resource port declared in R has no entry in equations.

13.2.1.3 FunctionDP

Wraps an arbitrary user callable h_fn. The return value is normalised to an antichain: an Antichain is passed through, a dict or a bare point becomes a singleton, and any other iterable is treated as a set of points.

FunctionDP(F: Poset, R: Poset,
           h_fn: Callable[[Any], Any],
           name: str = "function")

from codesign import FunctionDP, Ports, Reals
fdp = FunctionDP(
    F=Ports({"load": Reals()}),
    R=Ports({"cost": Reals()}),
    h_fn=lambda f: {"cost": f["load"] ** 2},
)
print(fdp.h({"load": 3.0}))  # Antichain[(cost=9)]

• h_fn — callable returning an Antichain, a dict, a point, or an iterable of points.

Errors: none raised directly; the caller guarantees monotonicity, and any exception raised inside h_fn propagates unchanged.

13.2.1.4 CatalogEntry

A dataclass describing one implementation in a catalogue: the functionality it can deliver and the resources it costs.

from codesign import CatalogEntry
e = CatalogEntry(provides={"speed": 5.0},
                 costs={"cost": 10.0}, name="A")

Errors: none raised directly.

13.2.1.5 CatalogDP

Selects the cheapest catalogue entries able to deliver f. An entry is feasible when its provides dominates f port-wise; the result is the Min of the feasible costs. If no entry works, h returns the singleton at \(\top\) (infeasible), rather than raising.

CatalogDP(F: Ports, R: Ports, catalog: Sequence,
          name: str = "catalog")

from codesign import CatalogDP, CatalogEntry, Ports, Reals
cat = CatalogDP(
    Ports({"speed": Reals()}), Ports({"cost": Reals()}),
    catalog=[
        CatalogEntry({"speed": 5.0}, {"cost": 10.0}, "A"),
        CatalogEntry({"speed": 9.0}, {"cost": 25.0}, "B"),
    ])
print(cat.h({"speed": 6.0}))  # only B works: (cost=25)

• catalog — a sequence of CatalogEntry objects or plain dicts with provides/costs/name keys.

Errors: raises TypeError if F or R is not a Ports; ValueError (naming the DP) if the catalogue is empty, or if two entries share the same non-empty name (the clashing name(s) are listed). Multiple entries left unnamed (the default name="") are allowed.

13.2.1.6 ConstraintDP

Defines a problem by a sampler over an implementation space, a feasible predicate, and a cost function; h returns the Min over all feasible samples (or the singleton at \(\top\) when none is feasible).

ConstraintDP(F: Poset, R: Ports,
             sampler: Callable[[Any], Iterable],
             feasible: Callable[[Any, Any], bool],
             cost: Callable[[Any], Mapping],
             name: str = "constraint")

from codesign import ConstraintDP, Ports, Reals
cdp = ConstraintDP(
    F=Ports({"span": Reals()}), R=Ports({"area": Reals()}),
    sampler=lambda f: [0.5, 1.0, 2.0],
    feasible=lambda impl, f: impl >= f["span"],
    cost=lambda impl: {"area": impl},
)
print(cdp.h({"span": 0.8}))  # Antichain[(area=1)]

• sampler(f) — yields candidate implementations.

• feasible(impl, f) — whether impl delivers f.

• cost(impl) — resource vector for a feasible sample.

Errors: none raised directly (no upfront interface check); exceptions inside the callables propagate.

13.2.1.7 ODE_DP

Derives a monotone relation from an ODE \(dx/dt = \texttt{rhs}(x, t, f)\), reading resources from the trajectory with extract. In final_value mode it integrates (explicit Euler) to t_end; in steady_state mode it Newton-solves \(\texttt{rhs}(x)=0\). See Section 7 for the wider solve pipeline.

ODE_DP(F: Poset, R: Ports,
       rhs: Callable[[Any, float, Any], Any],
       extract: Callable[[Any], Mapping],
       mode: str = "final_value",
       t_end: float = 10.0, n_steps: int = 200,
       x0_fn: Callable[[Any], Any] | None = None,
       name: str = "ode")

from codesign import ODE_DP, Ports, Reals
tank = ODE_DP(
    F=Ports({"inflow": Reals()}), R=Ports({"level": Reals()}),
    rhs=lambda x, t, f: f["inflow"] - 0.5 * x,
    extract=lambda x: {"level": x},
    mode="steady_state", x0_fn=lambda f: 0.0,
)
print(tank.h({"inflow": 2.0}))  # Antichain[(level=4)]

• mode — "final_value" or "steady_state".

• x0_fn(f) — initial state; a scalar or a positional sequence of floats (defaults to 0.0).

• t_end, n_steps — integration horizon and step count (final_value only).

Errors: the constructor raises ValueError if mode is not "final_value" or "steady_state". On the first h call the state is validated: a dict-valued (named) state raises TypeError (name the keys onto a list and index positionally), as does a str/bytes state.

13.2.1.8 UncertainDP

Brackets a nominal problem between an optimistic lower DP lower and a pessimistic upper DP upper, returning whichever the mode selects (paper Sec. VII). See Section 8 for the uncertainty solver.

UncertainDP(F: Poset, R: Poset,
            lower: DesignProblem, upper: DesignProblem,
            mode: str = "upper", name: str = "uncertain")

from codesign import AlgebraicDP, UncertainDP, Ports, Reals
F = Ports({"x": Reals()}); R = Ports({"y": Reals()})
lo = AlgebraicDP(F, R, {"y": lambda f: f["x"]})
hi = AlgebraicDP(F, R, {"y": lambda f: 1.5 * f["x"]})
u = UncertainDP(F, R, lower=lo, upper=hi, mode="upper")
print(u.h({"x": 2.0}))                     # (y=3)
print(u.with_mode("lower").h({"x": 2.0}))  # (y=2)

• lower, upper — the two bounding DPs.

• mode — "lower" (optimistic) or "upper" (pessimistic).

13.2.1.8.1 with_mode(mode)

Returns a fresh UncertainDP sharing the same bounds but with the given mode; the convenient way to obtain both fronts.
Errors: the constructor (and hence with_mode) raises ValueError if mode is not "lower" or "upper".

13.2.2 primitives↩︎

The primitives module provides small factory functions for the everyday plumbing DPs used to glue diagrams together (paper Fig. 26, Fig. 35). Each returns a ready-made AlgebraicDP (13.2.1.2) over scalar Reals ports, so its result is always a singleton antichain. Every factory accepts an optional poset argument; when omitted a plain Reals() is used for all ports. As thin wrappers over AlgebraicDP, they raise nothing directly; the ports they build are always Ports, so the wrapped constructor never faults.

13.2.2.1 adder

Sums several scalar inputs into one scalar output.

adder(in_names: list[str], out_name: str, poset=None)

from codesign import adder
a = adder(["p1", "p2"], "total")
print(a.h({"p1": 2.0, "p2": 3.0}))  # Antichain[(total=5)]

13.2.2.2 multiplier

Multiplies two scalar inputs, e.g. current times voltage gives power.

multiplier(in_a: str, in_b: str, out_name: str, poset=None)

from codesign import multiplier
m = multiplier("current", "voltage", "power")
print(m.h({"current": 2.0, "voltage": 5.0}))  # (power=10)

13.2.2.3 scale

Multiplies one input by a constant factor.

scale(in_name: str, out_name: str, factor: float, poset=None)

from codesign import scale
s = scale("capacity", "mass", factor=0.02)
print(s.h({"capacity": 1000.0}))  # Antichain[(mass=20)]

13.2.2.4 constant

Ignores its (trivial) functionality port "_" and emits a fixed value.

constant(out_name: str, value: float, poset=None)

from codesign import constant
c = constant("base_cost", value=42.0)
print(c.h({"_": 0.0}))  # Antichain[(base_cost=42)]

13.2.2.5 identity

Passes a single named scalar through unchanged (same port name on both sides).

identity(name: str, poset=None)

from codesign import identity
i = identity("speed")
print(i.h({"speed": 7.0}))  # Antichain[(speed=7)]

13.2.3 composition↩︎

The composition module supplies the three operators that close design problems under composition (Section 6): Series, Parallel and Loop, each with a lowercase functional alias (series, par, loop) matching the paper’s notation. All three preserve monotonicity, so a composite is itself a DesignProblem and solves through the same solve entry point (13.2.4.1).

13.2.3.1 Series

Series composition: dp1 then dp2, with each resource of dp1 fed in as a functionality of dp2. The composite has F = dp1.F and R = dp2.R. When both interfaces are Ports, the connectable check set(dp2.F.keys()) <= set(dp1.R.keys()) is enforced upfront; extra resource ports on dp1 are permitted and ignored.

Series(dp1: DesignProblem, dp2: DesignProblem,
       name: str | None = None)
# alias: series(dp1, dp2, name=None) -> Series

from codesign import series, scale
s1 = scale("torque", "current", factor=2.0)
s2 = scale("current", "mass", factor=0.5)
chain = series(s1, s2)
print(chain.h({"torque": 4.0}))  # Antichain[(mass=4)]

• dp1, dp2 — the two stages; dp1’s resources must cover dp2’s functionality ports.

Errors: raises ValueError (naming the missing ports) when both interfaces are Ports and some port dp2 requires is not produced by dp1.

13.2.3.2 Parallel

Parallel composition: dp1 and dp2 stacked side by side. The composite has F = F1 x F2 and R = R1 x R2 (Ports concatenation), and its antichain is the Cartesian product of the two.

Parallel(dp1: DesignProblem, dp2: DesignProblem,
         name: str | None = None)
# alias: par(dp1, dp2, name=None) -> Parallel

from codesign import par, scale
p = par(scale("x", "cost_x", 1.0),
        scale("y", "cost_y", 2.0))
print(p.h({"x": 3.0, "y": 4.0}))
# Antichain[(cost_x=3, cost_y=8)]

• dp1, dp2 — must have disjoint F port names and disjoint R port names.

Errors: raises TypeError if either F or either R is not a Ports; raises ValueError if the F port names overlap, or if the R port names overlap.

13.2.3.3 Loop

Feedback composition: closes a named axis present on both inner.F and inner.R, feeding the produced resource back as the consumed functionality (paper Def. 16, Section 6.3). The outer F and R drop the looped axis. Evaluating h defers to the solver’s Kleene iteration (13.2.4.2).

Loop(inner: DesignProblem, axis: str,
     name: str | None = None)
# alias: loop(inner, axis, name=None) -> Loop

from codesign import loop, AlgebraicDP, solve, Ports, Reals
inner = AlgebraicDP(
    F=Ports({"payload": Reals(), "mass": Reals()}),
    R=Ports({"mass": Reals(), "power": Reals(unit="W")}),
    equations={
        "mass": lambda f: 0.5 * f["payload"] + 0.3 * f["mass"],
        "power": lambda f: 2.0 * f["mass"],
    })
lp = loop(inner, axis="mass")     # feed R.mass back to F.mass
print(solve(lp, {"payload": 2.0}).antichain)
# Antichain[(power=2.857 W)]

• inner — the DP to close; needs Ports on both sides.

• axis — the port name (on both F and R) to feed back.

Errors: raises TypeError if inner.F or inner.R is not a Ports; raises ValueError if axis is not declared on both F and R of inner.

13.2.4 solver↩︎

The solver module drives every design problem through one entry point, solve, which runs the Kleene fixed-point iteration whenever a loop is present and returns a SolveResult. See Section 7 for the tutorial treatment and Section 7.4 for the observability features.

13.2.4.1 solve

Solves a (possibly composite) design problem at a given functionality. For a Loop it iterates to the least fixed point; for any other DP it is a single evaluation of h. This is the single most-used function in the library.

solve(
    dp: DesignProblem,
    functionality: Mapping | None = None,
    max_iter: int = 200,
    *,
    trace: bool = False,
    verbose: int = 0,
    on_iteration: Callable[[TraceEntry], None] | None = None,
    start_from: SolveResult | Antichain | None = None,
    record_trace: bool = False,
    uncertainty: list[str] | None = None,
    n_samples: int = 1000,
    rng_seed: int | None = None,
) -> SolveResult

from codesign import AlgebraicDP, solve, minimize_cost, Ports, Reals
dp = AlgebraicDP(
    F=Ports({"speed": Reals()}),
    R=Ports({"power": Reals(unit="W")}),
    equations={"power": lambda f: 2.0 * f["speed"]},
)
res = solve(dp, {"speed": 10.0}, trace=True)
print(res.status, res.iterations, res.feasible)
# converged 0 True
print(res.antichain)  # Antichain[(power=20 W)]

• functionality — outer F values; when None the solver uses dp.F.bottom().

• max_iter — cap on Kleene steps; reaching it sets status="max_iter".

• trace — collect a per-step trace on result.trace (13.2.4.5).

• verbose — 0 silent, 1 final summary, 2 per-iteration feed.

• on_iteration — callback invoked with each TraceEntry as it is produced.

• start_from — warm-start seed: a prior SolveResult (its inner antichain is reused) or an Antichain; ignored for non-loop DPs.

• record_trace — legacy alias for trace.

• uncertainty — if given (labels "worst_case", "mean", "p95", "cvar95", "samples"), dispatches to the uncertainty solver (Section 8) and returns an UncertaintyResult.

• n_samples, rng_seed — Monte Carlo sample count and optional seed for stochastic uncertainty.

Errors: raises TypeError if start_from is not a SolveResult, an Antichain, or None. Non-convergence is not an error: it is reported through status ("max_iter" or "diverged") and feasible, so callers must inspect those fields.

13.2.4.2 kleene_loop

The raw fixed-point iterator behind solve: it returns the loop’s antichain directly (without the SolveResult wrapper), writing metadata into an optional info_out dict. Most callers should prefer solve.

kleene_loop(
    loop_dp: Loop,
    f_outer: Mapping | None,
    max_iter: int = 200,
    *,
    trace: bool = False, verbose: int = 0,
    on_iteration: Callable | None = None,
    info_out: dict | None = None,
    start_from: Antichain | None = None,
) -> Antichain

from codesign import loop, kleene_loop, AlgebraicDP, Ports, Reals
inner = AlgebraicDP(
    F=Ports({"payload": Reals(), "mass": Reals()}),
    R=Ports({"mass": Reals(), "power": Reals(unit="W")}),
    equations={
        "mass": lambda f: 0.5 * f["payload"] + 0.3 * f["mass"],
        "power": lambda f: 2.0 * f["mass"],
    })
lp = loop(inner, axis="mass")
print(kleene_loop(lp, {"payload": 2.0}))
# Antichain[(power=2.857 W)]

• info_out — receives iterations, status, inner_antichain, and (if trace) trace.

• start_from — seed antichain in the inner R poset; a cold start (bottom singleton) is used when None.

Errors: raises TypeError if start_from is not None and lacks a points attribute (i.e.is not an Antichain). Runaway numeric values are capped at DIVERGENCE_CAP (\(10^{30}\)) and reported as status="diverged" rather than raised.

13.2.4.3 minimize_cost

Scalarises a solved antichain: returns the single resource bundle minimising a scalar cost_fn, or None when the result is infeasible or empty. This is the step that turns a Pareto front into one engineering choice.

minimize_cost(result: SolveResult,
              cost_fn: Callable[[Mapping], float]) -> Mapping | None

from codesign import AlgebraicDP, solve, minimize_cost, Ports, Reals
dp = AlgebraicDP(Ports({"speed": Reals()}),
                 Ports({"power": Reals(unit="W")}),
                 {"power": lambda f: 2.0 * f["speed"]})
res = solve(dp, {"speed": 10.0})
print(minimize_cost(res, cost_fn=lambda r: r["power"]))
# {'power': 20.0}

• result — a SolveResult, or any object with feasible and antichain attributes.

• cost_fn(r) — scalar cost of a resource bundle.

Errors: none raised directly; returns None on an infeasible or empty result.

13.2.4.4 SolveResult

The dataclass returned by solve, carrying the front plus convergence and feasibility metadata.

from codesign import SolveResult, AlgebraicDP, solve, Ports, Reals
dp = AlgebraicDP(Ports({"x": Reals()}), Ports({"y": Reals()}),
                 {"y": lambda f: f["x"]})
ac = solve(dp, {"x": 1.0}).antichain
sr = SolveResult(antichain=ac, iterations=3, status="max_iter")
print(sr.converged)  # False

The read-only property converged is a backward-compatibility alias for status == "converged". Errors: none raised directly.

13.2.4.5 TraceEntry

A dataclass recording one step of a Kleene iteration; a list of these is placed on result.trace when trace=True.

from codesign import solve, AlgebraicDP, Ports, Reals
dp = AlgebraicDP(Ports({"speed": Reals()}),
                 Ports({"power": Reals()}),
                 {"power": lambda f: 2.0 * f["speed"]})
te = solve(dp, {"speed": 10.0}, trace=True).trace[-1]
print(te.iteration, te.n_points, te.delta)  # 1 1 None

Errors: none raised directly.

13.3 Builders and rendering↩︎

The two high-level builders — the declarative mcdpl front end and the modular system builder, with its module base class and sugar constraint DSL — together with the diagram and viz rendering helpers. See Sections 11–12 and 10.

13.3.1 mcdpl↩︎

The codesign.mcdpl module supplies a single builder class, MCDP, a thin declarative front-end that mirrors the paper’s mcdp { ... } concrete syntax in pure Python. The tutorial in Section 11 walks through the surface syntax; the card below is the reference summary.

13.3.1.1 MCDP

Fluent builder for a single self-contained monotone design problem. One declares functionalities with provides and resources with requires, attaches one equation per resource with constraint (or a full antichain-valued rule), optionally closes feedback axes with loop_on, and emits a plain DesignProblem with build. Every declaration method returns self, so calls chain. It also works as a context manager.

Signature: MCDP(name="mcdp").

from codesign import MCDP, solve

with MCDP("battery") as m:
    m.provides("capacity", unit="J")
    m.requires("mass", unit="kg")
    m.constraint("mass", lambda f: f["capacity"] / 1.8e6)
battery = m.build()

res = solve(battery, {"capacity": 3.6e6})
print(res.antichain)          # Antichain[(mass=2 kg)]

Constructor parameters:

• name — str, optional (default "mcdp"). Names the emitted DesignProblem and appears in its errors.

Public methods:

Errors: loop_on raises ValueError naming the missing side and listing the current provides/requires declarations if axis is absent from either side. build raises ValueError if no functionality was declared, if no resource was declared, or (in the closed-form path) if any required resource has neither a constraint nor a rule. The other methods raise nothing directly.

13.3.2 system↩︎

The codesign.system module supplies the System builder, the recommended surface for non-trivial models. Section 12 is the full tutorial (both the operator-overloaded and legacy-lambda syntaxes); the card below summarises the API.

13.3.2.1 System

Modular MCDP composition with named subsystems. One declares outer functionalities (provides) and outer resources (requires), adds subsystems by name (add), wires them with connection constraints (constrain), and emits a single DesignProblem with build. The result solves like any other DP and can itself be nested as a subsystem. Internally build closes one feedback loop over the bundle of every subsystem’s \(R\) ports (the hidden __modules__ axis).

Signature: System(name="system").

from codesign import Module, Reals, System, solve

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def h(self, f):
        return {"mass": f["capacity"] / 1.8e6}

sys = System("drone")
endurance  = sys.provides("endurance", unit="s")
total_mass = sys.requires("total_mass", unit="kg")
battery    = sys.add("battery", Battery())
battery.capacity >= endurance * 5.0     # module F demand
total_mass       >= battery.mass        # outer R demand
res = solve(sys.build(), {"endurance": 300.0})
print(res.antichain)          # total_mass = 0.0008333 kg

Constructor parameters:

• name — str, optional (default "system"). Names the emitted DP.

Public methods:

The operator-overloaded >= form (as in the snippet) is sugar for constrain: Port.__ge__ registers the same internal constraint, so both syntaxes yield identical solves.

Errors: provides/requires raise ValueError on a duplicate name. add raises ValueError if the name is reused, contains a dot, equals the reserved __modules__, or if the DP’s \(F\)/\(R\) are not Ports. constrain raises TypeError if demand is neither callable nor an Expr, or if target is a non-constrainable port kind. build (via _validate) raises ValueError when there are no subsystems and no outer resources, when no outer resource is declared, when a constraint targets an unknown module or port, or when a subsystem \(F\) port or outer \(R\) has no constraint. draw_diagram raises ImportError if graphviz is not installed.

13.3.3 module↩︎

The codesign.module module supplies Module, a class-based way to declare a design problem by subclassing rather than by wiring up posets and an h_fn by hand. It is the natural unit to hand to System.add (see 13.3.2.1).

13.3.3.1 Module

Base class for declarative design problems. A subclass provides two class-level dict[str, Poset] attributes, F and R, and an instance method h(self, f); the constructor wraps these into the underlying DesignProblem machinery, so an instance is itself a DesignProblem. The subclass h may return a dict (a singleton antichain), a list of dicts (multi-valued), or an Antichain directly — the return value is normalised automatically. An optional class attribute module_name overrides the default name (the lower-cased class name).

Signature: Module() (subclasses may override __init__ to take parameters, calling super().__init__() last).

from codesign import Module, Reals, solve

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def h(self, f):
        return {"mass": f["capacity"] / 1.8e6}

battery = Battery()                       # a DesignProblem
print(battery.h({"capacity": 3.6e6}))     # Antichain[(mass=2 kg)]
print(solve(battery, {"capacity": 1.8e6}).antichain)  # mass=1 kg

Class-level declarations (set by the subclass, not passed to __init__):

• F — dict[str, Poset], functionality ports.

• R — dict[str, Poset], resource ports. Units are carried by the individual posets (e.g. Reals(unit="kg")).

• module_name — str, optional name override.

13.3.3.1.1 h(self, f)

The relation to override. Takes a functionality dict f keyed by the F port names, returns the resource requirement as a dict, a list of dicts, or an Antichain. The framework wraps non-antichain returns via Antichain.singleton / Antichain.from_set.

Errors: __init__ raises ValueError if neither F nor R is declared. The default h raises NotImplementedError if a subclass fails to override it. The return normaliser raises TypeError if h returns something other than a dict, a list/tuple of dicts, or an Antichain.

13.3.4 sugar↩︎

The codesign.sugar module implements the operator-overloaded constraint DSL used by the System builder (Section 12). A Port is a typed handle onto one of a system’s ports; arithmetic on ports builds an Expr tree lazily; >= on a constrainable port registers a constraint; ModuleHandle (returned by System.add) yields ports by attribute access; and sqrt/exp/log lift the corresponding math functions to operate on Expr trees.

The following listing exercises all six symbols; individual cards refer back to it.

from codesign import Module, Reals, System, solve, sqrt, exp, log

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def h(self, f):
        return {"mass": f["capacity"] / 1.8e6}

sys = System("s")
endurance = sys.provides("endurance", unit="s")   # outer_f Port
total     = sys.requires("total_mass", unit="kg")  # outer_r Port
bat       = sys.add("battery", Battery())          # ModuleHandle
e = (bat.mass + 0.5) * sqrt(endurance)             # an Expr
print(e.pretty())   # ((battery.mass + 0.5) * sqrt(endurance))
bat.capacity >= endurance * 5.0    # >= registers a constraint
total        >= bat.mass
print(solve(sys.build(), {"endurance": 300.0}).antichain)

13.3.4.1 Port

A named handle onto a port of a subsystem or an outer \(F\)/\(R\) name. Its kind is one of "module_f", "module_r", "outer_f", "outer_r", which governs where it may appear: module_f and outer_r ports may be the LHS of a >= constraint; module_r and outer_f ports may appear on the demand (RHS) side. Ports are not built directly — they are returned by System.provides/requires and by ModuleHandle attribute access. Arithmetic operators (inherited from Expr) build expression trees; >= registers a constraint with the parent system.

• full_name — property; "module.port" for module ports, the bare name otherwise.

• kind — property; the port-kind string above.

• pretty() — the full_name, used in expression rendering.

Errors: Port.__ge__ raises TypeError when the LHS is a module_r port (set by the module’s h) or an outer_f port (a system input); both name the offending port.

13.3.4.2 Expr

Base class for constraint-expression AST nodes (Add, Sub, Mul, Div, Pow, Neg, Func, Const, and Port). Users rarely name Expr directly; instances arise from arithmetic on ports and numbers. The arithmetic operators (+ - * / ** and unary -) return new Expr nodes rather than evaluating numerically, so the whole demand is captured as a tree that System later compiles to a callable.

• pretty() — a readable, fully parenthesised string of the tree (e.g. ((battery.mass + 0.5) * sqrt(endurance))); also the __repr__.

Errors: Expr.__ge__ raises TypeError (only ports, not compound expressions, may be a constraint LHS); __le__ raises TypeError (constraints are written target >= demand, not reversed); __bool__ raises TypeError to catch an expression used in a conditional. Coercing a non-numeric, non-Expr operand raises TypeError.

13.3.4.3 ModuleHandle

The handle returned by System.add. Attribute access (handle.port_name) returns a Port of kind module_f or module_r depending on which side of the subsystem’s DP the port lives on. Within one subsystem the \(F\) and \(R\) port names must be disjoint (enforced at construction).

• name — property; the subsystem’s registered name.

• dp — property; the underlying DesignProblem.

Its __repr__ lists the \(F\) and \(R\) ports, e.g. <ModuleHandle ’battery’: F=[’capacity’] R=[’mass’]>.

Errors: the constructor raises ValueError if the subsystem’s \(F\) and \(R\) port names overlap; attribute access raises AttributeError (listing the available \(F\)/\(R\) ports) for an unknown port name.

13.3.4.4 sqrt

Lift math.sqrt onto an expression: sqrt(x) -> Expr. The argument may be a Port, an Expr, or a number (coerced to a constant); the result is a Func node preserving the tree, as in the listing above (sqrt(endurance)). Errors: raises TypeError if the argument is neither a number nor an Expr.

13.3.4.5 exp

Lift math.exp onto an expression: exp(x) -> Expr. Same argument rules as sqrt (see 13.3.4.4); e.g. exp(endurance).pretty() is "exp(endurance)". Errors: raises TypeError on a non-number, non-Expr argument.

13.3.4.6 log

Lift math.log (natural log) onto an expression: log(x) -> Expr. Same argument rules as sqrt (see 13.3.4.4); e.g. log(endurance).pretty() is "log(endurance)". Errors: raises TypeError on a non-number, non-Expr argument.

13.3.5 diagram↩︎

The codesign.diagram module renders a System as a GraphViz block diagram — one box per subsystem with its \(F\) ports on the left and \(R\) ports on the right, outer \(F\)/\(R\) as margin nodes, constraint wiring as port-to-port edges, and any Kleene-iteration cycle coloured. See Section 10.1 for the rendered figures.

13.3.5.1 draw_system

Build a graphviz.Digraph for a system. It accepts either a live System or the DesignProblem returned by System.build() (which carries the modules and constraints on its _codesign_* attributes). Constraints written in the operator-overloaded form resolve to individual ports; lambda-based constraints are drawn as a dashed edge from a \(\lambda\) marker so none are silently dropped. Building the Digraph object needs only the graphviz Python package; the dot binary is required only to render it to a file.

Signature:

draw_system(system, *, name=None, rankdir="LR", show_ports=True,
            highlight_cycles=True, graph_attrs=None)

from codesign import Module, Reals, System, draw_system

class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def h(self, f):
        return {"mass": f["capacity"] / 1.8e6}

sys = System("drone")
endurance  = sys.provides("endurance", unit="s")
total_mass = sys.requires("total_mass", unit="kg")
battery    = sys.add("battery", Battery())
battery.capacity >= endurance * 5.0
total_mass       >= battery.mass

dot = draw_system(sys, rankdir="LR")   # a graphviz.Digraph
# dot.render("drone", format="svg")    # writes drone.svg (needs dot)

Parameters:

• system — a live System or a built System DP.

• name — str or None; diagram title, defaults to system.name.

• rankdir — "LR" (default, horizontal) or "TB" (vertical, for tall systems).

• show_ports — bool; when False, only module names are shown and edges attach to the box body.

• highlight_cycles — bool; colour edges inside strongly-connected components of size \(>1\).

• graph_attrs — dict of extra GraphViz graph attributes (e.g. {"ranksep": "0.6"}).

Errors: raises ImportError if the graphviz package is not installed; raises TypeError if system is neither a System nor a built System DP (i.e. lacks the _codesign_modules attribute). The convenience method System.draw_diagram forwards to this function.

13.3.6 viz↩︎

The codesign.viz module collects ready-made diagnostics. It is reached via the namespace import from codesign import viz (its functions are deliberately not in the top-level __all__) and called as viz.plot_antichain(...), and so on. The three plot_* functions require matplotlib (imported lazily, so the rest of the package works without it) and return a matplotlib.axes.Axes; to_dot has no plotting dependency and returns a string. Section 10 shows the figures in context.

The snippet below builds a small drone system used by every card that follows; the plot_* calls run headless under MPLBACKEND=Agg.

from codesign import Module, Reals, System, solve, viz

# ... define Battery (F: capacity -> R: mass) and Actuator
# ... (F: lift_force -> R: power) Modules as in earlier cards
sys = System("drone")
endurance = sys.provides("endurance", unit="s")
power_out = sys.requires("power", unit="W")
mass_out  = sys.requires("total_mass", unit="kg")
b = sys.add("battery", Battery())
a = sys.add("actuator", Actuator())
b.capacity   >= a.power * endurance
a.lift_force >= 9.81 * b.mass
power_out >= a.power
mass_out  >= b.mass
dp = sys.build()
res = solve(dp, {"endurance": 300.0}, trace=True)

13.3.6.1 viz.plot_antichain

Scatter an antichain’s Pareto front on two or three chosen \(R\)-port axes, optionally shading the dominated region (2D only). Accepts a SolveResult, an UncertaintyResult (uses its worst-case antichain), or an Antichain.

Signature:

plot_antichain(result, axes, *, ax=None, title=None,
               shade_dominated=True, point_size=60.0, label=None)

ax = viz.plot_antichain(res, ["power", "total_mass"])

Parameters: result the object to plot; axes a sequence of two or three \(R\)-port names; ax an existing axes (a new figure is made if omitted); title; shade_dominated; point_size; label (a legend entry is drawn only when named).

Errors: raises TypeError if result is not one of the three accepted types; ValueError if axes has other than 2 or 3 names, or if no point is plottable (empty, infeasible, or non-numeric); ImportError if matplotlib is absent.

13.3.6.2 viz.plot_convergence

Semilog plot of the per-iteration Kleene delta from a solve trace. Accepts a SolveResult whose trace was populated (call solve(..., trace=True)) or the trace list directly; zero deltas are clamped to floor for the log scale.

Signature:

plot_convergence(result_or_trace, *, ax=None, floor=1e-18,
                 title="Kleene-iteration convergence", label=None)

ax = viz.plot_convergence(res)   # res from solve(..., trace=True)

Parameters: result_or_trace; ax; title; floor (lower clamp for zero deltas); label.

Errors: raises ValueError if the trace has no numeric deltas; ImportError if matplotlib is absent.

13.3.6.3 viz.plot_uncertainty

Histogram of Monte-Carlo samples on one \(R\) port, with vertical lines for each requested summary (mean, p95, CVaR95, worst case) and an optional nominal marker. The UncertaintyResult must have been produced with "samples" in the uncertainty list.

Signature:

plot_uncertainty(result, port, *, ax=None, bins=30, title=None,
                 nominal=None, show_summaries=True, xlabel=None,
                 legend_loc="best")

ures = solve(dp, {"endurance": 300.0},
             uncertainty=["mean", "p95", "samples"],
             n_samples=200, rng_seed=0)
ax = viz.plot_uncertainty(ures, "total_mass",
                          xlabel="total mass (kg)")

Parameters: result an UncertaintyResult with samples; port the \(R\) port to histogram; bins; nominal; show_summaries; xlabel (defaults to the raw port name); legend_loc.

Errors: raises ValueError if the result carries no samples, or if the chosen port has no finite values; ImportError if matplotlib is absent.

13.3.6.4 viz.to_dot

Emit a GraphViz dot string describing a DP’s structure: a System-built DP renders its subsystems and constraint edges, a composition tree (Series/Parallel/Loop) renders as nested boxes, and a leaf DP renders as a single box. Returning a string keeps the module free of the graphviz dependency; save it and run dot -Tsvg.

Signature: to_dot(dp, *, name="codesign") -> str.

print(viz.to_dot(dp, name="drone").splitlines()[0])
# digraph drone {
print(viz.to_dot(dp, name="my drone").splitlines()[0])
# digraph "my drone" {   (non-identifier name is quoted)
}}

Parameters: dp the design problem to describe; name the graph id — used bare when it is a valid GraphViz identifier, otherwise double-quoted with embedded quotes and backslashes escaped.

Constraint edges track the modules referenced on the right-hand side even when a module port is wrapped in a function or negation — the reference walker descends through .inner (sqrt(...), exp, a leading minus) as well as .left/.right/.arg, so a constraint like bat.capacity >= sqrt(act.power) still draws the act\(\rightarrow\)bat edge.

Errors: none raised directly; the walker swallows introspection failures and returns a best-effort graph.

13.4 Uncertainty and online learning↩︎

The set-based and stochastic uncertainty layer and the compositional online-learning layer of optimistic evaluators. The modelling discussion is in Sections 8 and 9.

13.4.1 uncertainty↩︎

Reference cards for the uncertainty layer of module codesign.uncertainty. Both regimes are declared on a Module through its uncertain_set (set-based) and uncertain_dist (stochastic) attributes, and both are consumed by the ordinary solve entry point; see Section 8 for the modelling discussion. The set-based classes (UncertaintySet, Box, Ellipsoid, Disk, Circle) are grouped first, the stochastic ones (Independence, GaussianCopula, Stochastic) after, followed by the result dataclass and the solver entry point. Every listing requires numpy; the stochastic cards additionally require scipy.

13.4.1.1 UncertaintySet

Abstract base for deterministic, set-based parameter uncertainty. A concrete subclass declares which parameters vary (param_names) and where the set’s worst point lies (worst_case_values); the worst-case solve then reassigns those named module attributes before re-solving. Signature: UncertaintySet() (no constructor of its own).

from codesign import Box
from codesign.uncertainty import UncertaintySet
s = Box(mass=(1.0, 2.0, "more_is_worse"))
print(isinstance(s, UncertaintySet))   # True
print(s.param_names())                 # ['mass']

Public methods (both abstract, overridden by subclasses):

• param_names() -> list[str] — names of the uncertain parameters carried by the set.

• worst_case_values(h_callable, base_values, f_inner) -> dict[str, float] — the parameter dict at the worst point of the set. h_callable(overrides) -> Antichain is evaluated only for undeclared directions (boundary search); base_values supplies nominal values for out-of-set parameters.

Errors: the base methods raise NotImplementedError; instantiate a concrete subclass instead.

13.4.1.2 Box

Axis-aligned interval product over one or more parameters — the everyday set-based model. When every parameter declares a direction of badness the worst case is a single corner, read off in constant time. Signature: Box(**params) where each value is (lo, hi) or (lo, hi, direction).

from codesign import Box
b = Box(se=(1.7e6, 2.3e6, "more_is_better"),
        eff=(0.83, 0.97, "more_is_better"))
print(b.param_names())                    # ['se', 'eff']
print(b.worst_case_values(None, {}, {}))
# {'se': 1700000.0, 'eff': 0.83}

Parameters:

• **params — one keyword per uncertain parameter. The value is (lo, hi) (direction undeclared) or (lo, hi, direction) with direction one of "more_is_better", "more_is_worse", "less_is_better", "less_is_worse" (in parameter units).

With \(n\) undeclared parameters all \(2^n\) endpoint corners are probed (cheap for modest \(n\)). Errors: raises ValueError if a value is neither a 2- nor 3-tuple, if an unknown direction token is passed, or if lo > hi (empty interval).

13.4.1.3 Ellipsoid

An \(n\)-D ellipsoid \((p-c)^\top \Sigma^{-1}(p-c) \leq 1\) in parameter space; the covariance \(\Sigma\) captures scale and correlation, so it is a more honest set than a box when parameters vary together.

from codesign import Ellipsoid
e = Ellipsoid(
    center={"se": 2.0e6, "eff": 0.9},
    cov=[[1.0e10, 0.0], [0.0, 2.5e-3]],
    params=["se", "eff"],
    directions={"se": "more_is_better",
                "eff": "more_is_better"},
)
print(e.worst_case_values(None, {}, {}))
# {'se': 1929289.32..., 'eff': 0.86464...}

Parameters:

• center — dict of centre coordinates by parameter name.

• cov — symmetric positive-definite shape matrix \(\Sigma\), rows/columns ordered as params.

• params — parameter names in cov’s order.

• directions=None — optional per-parameter badness tokens; when all are given the worst case is the analytic boundary point \(c + L\,u^\star\), else the boundary is sampled.

• boundary_samples=8 — samples per dimension used when directions are not fully declared.

Errors: raises ValueError if cov’s shape does not match params or if cov is not positive definite; ImportError if numpy is absent.

13.4.1.4 Disk

Two-dimensional convenience: a filled circle, i.e.the isotropic ellipsoid \(\Sigma = r^2 I\). A factory function returning an Ellipsoid, not a class. Signature: Disk(center, radius, params=None, directions=None).

from codesign import Disk, Ellipsoid
d = Disk(center={"x": 0.0, "y": 0.0}, radius=1.0,
         directions={"x": "more_is_worse",
                     "y": "more_is_worse"})
print(isinstance(d, Ellipsoid))               # True
print(d.worst_case_values(None, {}, {}))
# {'x': 0.70710..., 'y': 0.70710...}

Parameters: center (2-key dict), radius (in parameter units), params (inferred from center if None), directions (as for Ellipsoid). Errors: raises ValueError unless exactly two parameters are supplied.

13.4.1.5 Circle

Two-dimensional convenience for the circle boundary only. For monotone modules with declared directions the worst case lies on the boundary, so Circle and Disk return identical worst-case answers; the implementation delegates to Disk. Signature: Circle(center, radius, params=None, directions=None).

from codesign import Circle, Ellipsoid
c = Circle(center={"x": 0.0, "y": 0.0}, radius=2.0,
           directions={"x": "more_is_worse",
                       "y": "more_is_worse"})
print(isinstance(c, Ellipsoid))               # True
print(c.worst_case_values(None, {}, {}))
# {'x': 1.41421..., 'y': 1.41421...}

Parameters: as Disk. Errors: as Disk (raises ValueError unless exactly two parameters are supplied).

13.4.1.6 Independence

The independence copula: each component drawn independently in \(U(0,1)\). The default dependence structure of Stochastic. Signature: Independence(); method sample_uniform(n, d, rng) -> ndarray of shape (n, d).

from codesign.uncertainty import Independence
import numpy as np
rng = np.random.default_rng(0)
u = Independence().sample_uniform(3, 2, rng)
print(u.shape)                                # (3, 2)

Parameters: none. sample_uniform takes the sample count n, the dimension d, and a numpy Generator rng. Errors: none raised directly.

13.4.1.7 GaussianCopula

Gaussian copula with a given correlation matrix: draw \(z \sim N(0, R)\) and apply \(\Phi\) to each coordinate to obtain correlated uniforms. Signature: GaussianCopula(correlation).

from codesign import GaussianCopula
import numpy as np
gc = GaussianCopula(correlation=[[1.0, 0.4],
                                 [0.4, 1.0]])
rng = np.random.default_rng(0)
u = gc.sample_uniform(4, 2, rng)
print(u.shape, repr(gc))
# (4, 2) GaussianCopula(d=2)

Parameters: correlation — a \((d,d)\) symmetric positive-definite matrix with unit diagonal. Errors: the constructor raises ValueError if correlation is not square or not positive definite; sample_uniform raises ValueError if d differs from the copula’s dimension. Sampling additionally requires scipy.

13.4.1.8 Stochastic

Joint distribution over uncertain parameters, built from named scipy.stats frozen marginals plus an optional copula. This is the uncertain_dist declaration consumed by the Monte Carlo path of solve; it has a full treatment in Section 8. Signature: Stochastic(marginals=None, copula=None, **marginals_kw).

from codesign import Stochastic
from scipy import stats
import numpy as np
st = Stochastic(marginals={
    "m": stats.norm(loc=1.0, scale=0.1)})
print(st.param_names())                       # ['m']
print(st.sample(2, np.random.default_rng(0)))
# [{'m': 1.0350...}, {'m': 0.9386...}]

Parameters:

• marginals — dict of parameter name to a frozen scipy distribution (e.g.stats.norm(loc, scale)). Marginals may also be given as keyword arguments (**marginals_kw).

• copula=None — a Copula; defaults to Independence.

Methods: param_names() -> list[str]; sample(n, rng) -> list[dict] of \(n\) parameter dicts (copula uniforms mapped through each marginal’s ppf). Errors: the constructor raises ValueError if no marginal is supplied.

13.4.1.9 UncertaintyResult

Dataclass returned by an uncertainty-aware solve. Which fields are populated depends on the summaries requested via uncertainty=[...]; the rest stay None.

from codesign.uncertainty import UncertaintyResult
r = UncertaintyResult(mean={"mass": 0.55},
                      feasibility_rate=1.0,
                      n_samples_used=1000)
print(r.mean, r.feasibility_rate)
# {'mass': 0.55} 1.0

Fields:

Statistics are computed over feasible samples only. Errors: none raised directly (a plain dataclass).

13.4.1.10 solve_with_uncertainty

Entry point for an uncertainty-aware solve, returning an UncertaintyResult. Ordinarily reached through solve(dp, f, uncertainty=[...]), which dispatches here when the uncertainty argument is present. Signature:

solve_with_uncertainty(dp, functionality, uncertainty,
    n_samples=1000, rng_seed=None, max_iter=200,
    verbose=0)

from codesign import Module, Reals, solve, Box
class Battery(Module):
    F = {"capacity": Reals(unit="J")}
    R = {"mass":     Reals(unit="kg")}
    def __init__(self, se=2.0e6, eff=0.9):
        self.se, self.eff = se, eff
        super().__init__()
    def h(self, f):
        return {"mass": f["capacity"] / (self.se * self.eff)}
bat = Battery()
bat.uncertain_set = Box(se=(1.7e6, 2.3e6, "more_is_better"),
                        eff=(0.83, 0.97, "more_is_better"))
res = solve(bat, {"capacity": 1.0e6}, uncertainty=["worst_case"])
print(round(list(res.worst_case.antichain.points)[0]["mass"], 5))
# 0.70872

Parameters:

• dp, functionality — the design problem and its outer F vector, as for solve.

• uncertainty — list drawn from "worst_case", "mean", "p95", "cvar95", "samples".

• n_samples=1000 — Monte Carlo draws for the stochastic summaries; rng_seed=None seeds the numpy generator.

• max_iter=200, verbose=0 — forwarded to each inner solve.

Errors: raises ValueError for an unknown summary label, when no module on dp carries an uncertain_set/uncertain_dist, or when stochastic summaries are requested but no uncertain_dist is present; ImportError if numpy is absent.

13.4.2 online↩︎

Reference cards for the compositional online learner of module codesign.online — a port of [18] (arXiv:2604.22624). The evaluators maintain history-dependent bounds on each candidate’s inner-solve output so that provably suboptimal candidates are eliminated without ever being solved; solve_online is the budgeted driver. See Section 9 for the modelling discussion and Section 9.3 for the warm-start and picker options. Every listing requires numpy; LinearParametricEvaluator additionally requires scipy.

13.4.2.1 OptimisticEvaluator

Abstract base holding the observation history and the bound contract. The public bound(candidate) returns a pair of dicts (lower, upper) mapping each numeric R component to its current lower and upper bound at the queried feature point; the base implementation is the trivial fallback \((0, +\infty)\), which subclasses tighten. Signature: OptimisticEvaluator(features, r_components).

from codesign.online import OptimisticEvaluator
ev = OptimisticEvaluator(features=["x"],
                         r_components=["cost"])
print(ev.bound({"x": 2.0}))
# ({'cost': 0.0}, {'cost': inf})

Parameters: features — candidate keys used as the feature vector; r_components — resource-port names to bound.

Errors: observe/bound raise KeyError if a candidate is missing a declared feature.

13.4.2.2 MonotonicityEvaluator

Bounds from monotonicity alone: assuming the output is component-wise monotone in the features, an observation at \(f_0\) with antichain-min \(R_0\) lower-bounds every candidate with features \(\geq f_0\) by \(R_0\) and upper-bounds every candidate with features \(\leq f_0\) by \(R_0\). The aggressive choice when the assumption genuinely holds. Signature: MonotonicityEvaluator(features, r_components).

from codesign import (MonotonicityEvaluator, Antichain,
                      Ports, Reals)
R = Ports({"cost": Reals()})
ev = MonotonicityEvaluator(["x"], ["cost"])
ev.observe(0, {"x": 1.0}, Antichain.singleton(R, {"cost": 5.0}))
print(ev.bound({"x": 2.0}))   # ({'cost': 5.0}, {'cost': inf})
print(ev.bound({"x": 0.5}))   # ({'cost': 0.0}, {'cost': 5.0})

Parameters: as OptimisticEvaluator. Behaviour is described, not its implementation: the bound is exact to float comparisons, computed by a vectorised numpy rescan of the whole history when numpy is present, with an equivalent pure-Python fallback for numpy-absent installs. Errors: raises KeyError for a missing feature.

13.4.2.3 LipschitzEvaluator

Bounds from a Lipschitz assumption: with \(|h(c_1)-h(c_2)| \leq L\,\|features(c_1)-features(c_2)\|\), each observation tightens the bound everywhere by a cone of slope \(L\). Signature: LipschitzEvaluator(features, r_components, L).

from codesign import (LipschitzEvaluator, Antichain,
                      Ports, Reals)
R = Ports({"cost": Reals()})
ev = LipschitzEvaluator(["x"], ["cost"], L=2.0)
ev.observe(0, {"x": 1.0}, Antichain.singleton(R, {"cost": 5.0}))
print(ev.bound({"x": 2.0}))   # ({'cost': 3.0}, {'cost': 7.0})

Parameters: L — a single positive float (same constant for every R component) or a dict mapping R component name to its own constant. Behaviour, not implementation: distances use a vectorised numpy path when numpy is present and an equivalent pure-Python fallback otherwise. Errors: raises KeyError for a missing feature, or if L is a dict lacking an R component.

13.4.2.4 LinearParametricEvaluator

The certified optimistic bound of [18] (Section V-C3, eqs.–28 and Lemma V.5). Each resource coordinate is assumed an exact affine function of the features, \(\text{req}_k(c)=\phi(c)^\top\theta_k^*\) with \(\phi(c)=[1,features(c)]\); the evaluator maintains the confidence polytope \(\Theta(H)\) of parameters consistent with every observation and returns, per coordinate, \(\min_{\theta\in\Theta(H)}\phi(i)^\top\theta\) — one LP solved with scipy.optimize.linprog. Because \(\theta_k^*\) is feasible, the optimum is \(\leq \text{req}_k(i)\): a guaranteed lower bound (the optimism guarantee), so a Pareto-optimal candidate is never wrongly eliminated. No upper bound is certified, so the returned upper bound stays at \(+\infty\). Signature:

LinearParametricEvaluator(features, r_components,
    confidence=None, min_obs=3, *, prior_box=None,
    noise_bound=0.0, solver="highs")

from codesign import (LinearParametricEvaluator, Antichain,
                      Ports, Reals)
R = Ports({"cost": Reals()})
ev = LinearParametricEvaluator(["x"], ["cost"],
                               prior_box=(-100.0, 100.0))
for x in (1.0, 2.0, 3.0):            # cost = 2x + 1, affine
    ev.observe(0, {"x": x},
               Antichain.singleton(R, {"cost": 2*x + 1}))
lo, hi = ev.bound({"x": 4.0})        # certified lower bound
print(round(lo["cost"], 3), hi["cost"])   # 9.0 inf

Parameters:

• prior_box=None — the prior box \(\Theta_0\) keeping the LP bounded while the fit is under-determined. None leaves every parameter unbounded (always safe, but trivial until \(\theta\) is pinned down); a (lo, hi) float pair applies to every parameter; a dict gives each R component its own (lo, hi). A finite box must contain the true \(\theta^*\) or the bound may become invalid.

• noise_bound=0.0 — half-width of the observation band \(|\phi(j)^\top\theta - r| \leq\) noise_bound; 0.0 is the paper’s exact/noiseless model, positive values a documented bounded-noise extension.

• min_obs=3 — below this many observations the trivial \((0,+\infty)\) bound is returned.

• solver="highs" — method forwarded to linprog.

• confidence=None — deprecated and ignored, retained only for backward compatibility with the former OLS \(\pm\) confidence-band version; passing it emits a DeprecationWarning. The certified bound has no confidence-band parameter.

Errors: raises ValueError if noise_bound is negative; TypeError if a prior_box entry is not a (lo, hi) pair; ImportError if numpy or scipy is absent. An unbounded, infeasible, or failed LP degrades safely to the trivial lower bound.

13.4.2.5 GaussianProcessEvaluator

Bounds from a zero-mean Gaussian-process surrogate with an RBF kernel, fitted in pure numpy; the bound at a query is \(\text{mean} \pm \texttt{confidence}\cdot\sigma\). The right tool when the response surface has feature interactions or local nonlinearity, where the certified linear evaluator would degrade to the no-information bound. The bounds are informative but not certified. Signature:

GaussianProcessEvaluator(features, r_components,
    length_scale=0.3, sigma_f=1.0, noise=1e-3,
    confidence=2.0, min_obs=3)

from codesign import (GaussianProcessEvaluator, Antichain,
                      Ports, Reals)
R = Ports({"cost": Reals()})
ev = GaussianProcessEvaluator(["x"], ["cost"],
                              length_scale=0.5)
for x in (0.0, 0.5, 1.0):
    ev.observe(0, {"x": x},
               Antichain.singleton(R, {"cost": x}))
lo, hi = ev.bound({"x": 0.25})
print(round(lo["cost"], 3), round(hi["cost"], 3))
# 0.066 0.289

Parameters:

• length_scale=0.3 — RBF length scale (0.3 suits features in roughly \([0,1]\)).

• sigma_f=1.0 — kernel amplitude, rescaled per output by the empirical standard deviation.

• noise=1e-3 — observation-noise jitter stabilising the Cholesky factor.

• confidence=2.0 — multiplier on the predictive standard deviation (\(\approx 95\%\) band under Gaussian residuals).

• min_obs=3 — below this the fallback \((0,+\infty)\) bound is returned.

observe additionally invalidates the fit cache (a lazy refit runs on the next bound). Errors: raises KeyError for a missing feature; ImportError if numpy is absent.

13.4.2.6 OnlineResult

Dataclass returned by solve_online, carrying the surviving antichain plus a full audit of the elimination cascade.

from codesign import OnlineResult, Antichain, Ports, Reals
R = Ports({"cost": Reals()})
r = OnlineResult(
    antichain=Antichain.singleton(R, {"cost": 5.0}),
    n_evaluated=3, n_eliminated=7, n_candidates=10)
print(r.n_evaluated, r.n_eliminated, r.n_candidates)
# 3 7 10

Fields:

Errors: none raised directly (a plain dataclass).

13.4.2.7 solve_online

The budgeted online driver: solve’s elimination-based cousin. For each candidate, candidate_fn builds a fresh DP; the evaluator’s bounds prune provably suboptimal candidates before their inner solve is run. Signature:

solve_online(candidate_fn, functionality, *,
    candidates, evaluator, budget=None, max_iter=200,
    verbose=0, warm_start=None, picker="lcb")

from codesign import (AlgebraicDP, Ports, Reals,
                      LipschitzEvaluator, solve_online)
F, R = Ports({"load": Reals()}), Ports({"cost": Reals()})
def candidate_fn(c):
    return AlgebraicDP(F, R, {"cost":
        lambda f, u=c["price"]: f["load"] * u})
cands = [{"price": p} for p in (2.0, 3.0, 5.0, 8.0)]
ev = LipschitzEvaluator(["price"], ["cost"], L=10.0)
res = solve_online(candidate_fn, {"load": 1.0},
                   candidates=cands, evaluator=ev,
                   warm_start=2, picker="ucb")
print(res.n_evaluated, res.n_eliminated, res.antichain)
# 4 0 Antichain[(cost=2)]

Parameters:

• candidate_fn — candidate_fn(candidate) -> DP.

• functionality — the outer F vector, passed to every inner solve unchanged.

• candidates — list of feature dicts; each must contain every name in evaluator.features.

• evaluator — the OptimisticEvaluator instance (reset and repopulated during the solve).

• budget=None — maximum inner solves; None is unlimited (elimination still prunes).

• max_iter=200, verbose=0 — forwarded to each inner solve; verbose 0/1/2 for silent/summary/trace.

• warm_start=None — seed evaluations before the picker takes over: an int \(n\) triggers greedy farthest-point sampling of \(n\) candidates; a list of indices seeds exactly those.

• picker="lcb" — selection strategy: "lcb" (exploit), "ucb" (add exploration bonus, tune via ("ucb", {"kappa": 1.0})), "random", or a callable (lo, hi, r_components, **kwargs) -> score.

Errors: raises ValueError for an unknown picker name and TypeError for a picker of the wrong type; KeyError if a candidate is missing a declared feature.

13.5 Temporal layers↩︎

The open-loop temporal planners: architecture switching (temporal), scalar dynamic programming (dynamic), and antichain-valued sequential co-design (sequential). See Section 15.

13.5.1 temporal↩︎

The temporal module implements Case 1 of the temporal co-design layers (section 15): architecture switching across a changing environment. A sequence of epochs each poses its own functionality and admits a set of candidate architectures; solve_schedule chooses one architecture per epoch to minimise the total of the per-epoch co-design costs plus the switching costs, by an exact Viterbi (shortest-path) pass over the epoch-by-architecture lattice. The per-epoch cost is an ordinary solve(), so this layer sits on top of static co-design without modifying it. The Architecture defined here is the shared decision object reused by every temporal layer.

13.5.1.1 Architecture

A named candidate configuration that can be active in an epoch or stage; the common decision object of the whole temporal family. It is a dataclass wrapping a design problem.

Architecture(name: str, dp: Any, tags: Mapping = {})

from codesign import AlgebraicDP, Architecture, Ports, Reals, System
s = System("eco")
d = s.provides("demand")
s.requires("cost")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"cost": Reals()}),
      equations={"cost": lambda f: 10.0})).demand >= d
s.constrain("cost", lambda x: x["m.cost"])
eco = Architecture("eco", s.build(), tags={"mode": "eco"})
print(eco.name, eco.tags)          # -> eco {'mode': 'eco'}

• name — str — identifier used in results, plots and switching bookkeeping.

• dp — Any — a design problem accepted by solve() (typically a built System); may be shared across epochs or epoch-specific.

• tags — Mapping — free-form metadata carried through to results; not interpreted by the solver.

Errors: none raised directly.

13.5.1.2 Epoch

One environment regime over which a single architecture is active. The changing environment enters through the per-epoch functionality.

Epoch(name: str, functionality: Mapping,
      candidates: Optional[Sequence[Architecture]] = None,
      duration: float = 1.0)

from codesign import Architecture, Epoch
arch = Architecture("eco", None, tags={"mode": "eco"})
ep = Epoch("glucose", {"mu": 0.8}, candidates=[arch], duration=2.0)
print(ep.name, ep.functionality, ep.duration)
# -> glucose {'mu': 0.8} 2.0

• name — str — epoch identifier.

• functionality — Mapping — the outer functionality demanded in this epoch, passed straight to solve().

• candidates — Sequence[Architecture] or None — architectures admissible here; when None the schedule-level default set is used, so some architectures can be made unavailable in some regimes.

• duration — float — weight on the epoch’s running cost, expressing unequal epoch lengths without rescaling the cost function.

Errors: none raised directly (a missing default set is reported by solve_schedule; see below).

13.5.1.3 EpochResult

The per-epoch outcome recorded inside a solved schedule. A dataclass; its total_cost property sums the running and switching contributions.

from codesign import EpochResult
er = EpochResult(epoch="glucose", architecture="eco",
                 running_cost=8.0, switch_cost=0.5, feasible=True,
                 point={"cost": 8.0})
print(er.total_cost)               # -> 8.5

• epoch — str — the epoch’s name.

• architecture — str — name of the chosen architecture.

• running_cost — float — duration-weighted scalar cost of the chosen design point (inf if infeasible).

• switch_cost — float — switching cost charged at this epoch’s boundary (zero when the architecture is unchanged).

• feasible — bool — whether this epoch was satisfiable by its chosen architecture.

• point — Mapping or None — the chosen resource point.

• tags — Mapping — the chosen architecture’s tags.

• total_cost (property) — running_cost \(+\) switch_cost.

Errors: none raised directly.

13.5.1.4 ScheduleResult

The outcome of solve_schedule: the chosen architecture per epoch with decomposed costs. A dataclass with a schedule property and a compact __repr__.

from codesign import EpochResult, ScheduleResult
e0 = EpochResult("glucose", "eco", 8.0, 0.0, True, {"cost": 8.0})
e1 = EpochResult("acetate", "bal", 6.0, 0.5, True, {"cost": 6.0})
sr = ScheduleResult([e0, e1], total_cost=14.5, feasible=True,
                    n_switches=1)
print(sr.schedule)
print(sr)
# -> ['eco', 'bal']
# -> ScheduleResult(feasible, cost=14.5, switches=1, [eco -> bal])

• epochs — List[EpochResult] — one entry per epoch, in order.

• total_cost — float — running plus switching cost summed over epochs.

• feasible — bool — True iff every epoch was satisfiable by its chosen architecture.

• n_switches — int — number of epoch boundaries at which the active architecture changed.

• schedule (property) — List[str] — the chosen architecture name per epoch, in order.

Errors: none raised directly.

13.5.1.5 solve_schedule

Chooses an architecture per epoch minimising running plus switching cost, by an exact Viterbi pass over the epoch/architecture lattice. With switch_cost == 0 and hysteresis == 0 the result reduces to the epoch-local greedy choice.

solve_schedule(epochs, architectures=None, *, cost_fn,
               switch_cost=0.0, hysteresis=0.0, max_iter=200,
               solve_kwargs=None) -> ScheduleResult

from codesign import (AlgebraicDP, Architecture, Epoch, Ports,
                      Reals, System, solve_schedule)
def mk(name, c):
    s = System(name); d = s.provides("demand"); s.requires("cost")
    s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
          R=Ports({"cost": Reals()}),
          equations={"cost": lambda f, c=c: c})).demand >= d
    s.constrain("cost", lambda x: x["m.cost"]); return s.build()
lo = Architecture("lo", mk("lo", 5.0))
hi = Architecture("hi", mk("hi", 3.0))
eps = [Epoch("e1", {"demand": 1.0}, candidates=[lo]),
       Epoch("e2", {"demand": 1.0}, candidates=[hi])]
r = solve_schedule(eps, cost_fn=lambda p: p["cost"], switch_cost=1.0)
print(r.schedule, r.total_cost, r.n_switches)   # -> ['lo', 'hi'] 9.0 1

• epochs — Sequence[Epoch] — the ordered environment regimes.

• architectures — Sequence[Architecture] or None — default candidate set for any epoch that does not supply its own; required if any epoch leaves candidates as None.

• cost_fn — callable — maps a resource point to a scalar; lower is better.

• switch_cost — float or callable — cost charged when the active architecture changes; a float applies uniformly, a callable f(prev, next) -> float makes it transition-dependent.

• hysteresis — float — extra margin a challenger must beat the incumbent by before a switch is preferred; suppresses chattering between near-equal options.

• max_iter — int — forwarded to solve() for each epoch solve.

• solve_kwargs — Mapping or None — extra keyword arguments forwarded to solve() (for example uncertainty).

Errors: raises ValueError when an epoch supplies no candidates and no default architectures set was passed; the message reads

epoch <name> has no candidates and no default architecture set
was supplied. Either set candidates=[...] on the epoch, or pass
architectures=[...] to solve_schedule().

An infeasible epoch is not an error: the schedule is returned with feasible=False and the offending running cost recorded as inf for inspection.

13.5.2 dynamic↩︎

The dynamic module implements Case 2 with a scalar carried state (section 15): a finite-horizon dynamic program whose per-stage decision is which architecture to instantiate, whose per-stage cost is itself a co-design solve(), and whose carried scalar state (fuel, charge, cumulative wear) is advanced by a user transition. The value function is scalar — at each (stage, state) the single best cost-to-go — found by a backward Bellman recursion. State is discretised onto an explicit StateGrid; a transition that would leave the grid envelope is rejected before snapping, so an over-spent resource is never silently rescued by the nearest in-range node. It reuses the Architecture of the temporal module (13.5.1.1).

13.5.2.1 Stage

One step of the finite-horizon dynamic program: everything needed to score and advance every candidate architecture at this point in time given the carried state. A dataclass.

Stage(name: str,
      functionality: Callable[[float], Mapping],
      transition: Callable[[float, Mapping], float],
      candidates: Optional[Sequence[Architecture]] = None,
      admissible: Optional[Callable[[float], bool]] = None)

from codesign import Architecture, Stage
mode = Architecture("burn", None)
st = Stage("leg1", functionality=lambda s: {"demand": 1.0},
           transition=lambda s, p: s - p["fuel"],
           candidates=[mode], admissible=lambda s: s >= 0.0)
print(st.name, st.transition(10.0, {"fuel": 3.0}))   # -> leg1 7.0

• name — str — stage identifier.

• functionality — callable — functionality(state) -> Mapping, the outer functionality demanded when the carried resource has value state. A state-independent stage ignores its argument.

• transition — callable — transition(state, point) -> new_state, mapping the incoming state and the chosen resource point to the outgoing state; the returned value is snapped to the nearest grid node.

• candidates — Sequence[Architecture] or None — architectures admissible here; falls back to the schedule-level default when None.

• admissible — callable or None — admissible(state) -> bool; states returning False are forbidden (infinite value).

Errors: none raised directly.

13.5.2.2 StageResult

A per-stage record along a rolled-out optimal policy. A dataclass with no public methods; constructed by rollout.

from codesign import StageResult
sr = StageResult(stage="leg1", architecture="burn", state_in=10.0,
                 state_out=7.0, stage_cost=4.0, feasible=True,
                 point={"fuel": 3.0})
print(sr.stage, sr.state_in, sr.state_out, sr.stage_cost)
# -> leg1 10.0 7.0 4.0

• stage — str — the stage’s name.

• architecture — str — chosen architecture name.

• state_in — float — snapped carried state entering the stage.

• state_out — float — snapped carried state leaving the stage.

• stage_cost — float — co-design cost incurred at this stage.

• feasible — bool — whether a finite-cost action existed here.

• point — Mapping or None — the chosen resource point.

• tags — Mapping — the chosen architecture’s tags.

Errors: none raised directly.

13.5.2.3 StateGrid

A one-dimensional discretisation of the carried scalar resource. Tabular dynamic programming needs a finite state set, so a continuous resource is bucketed onto explicit, sorted nodes; transitions landing between nodes are snapped to the nearest. Keeping the grid an explicit object makes the discretisation visible and lets the caller trade accuracy against cost.

StateGrid(nodes: Sequence[float])

from codesign import StateGrid
g = StateGrid.linspace(0.0, 10.0, 6)   # nodes 0,2,4,6,8,10
print(list(g), len(g), g.snap(5.4))
# -> [0.0, 2.0, 4.0, 6.0, 8.0, 10.0] 6 6.0

• nodes — Sequence[float] — the grid node values (stored sorted). Also iterable, with len() giving the node count.

Errors: the constructor raises ValueError on an empty node sequence (“StateGrid needs at least one node, got an empty sequence.…”); linspace raises ValueError when n < 1 (“StateGrid.linspace needs at least one node, got n=<n>.Pass n >= 1 …”). With n == 1 a single node at lo is returned.

13.5.2.4 DynamicPolicy

A solved, state-indexed policy: the best architecture per (stage, state), produced by solve_dynamic. Beyond the single rolled-out path the full table supports closed-loop control, re-queried at the actual realised state when it diverges from the nominal roll-out.

from codesign import (AlgebraicDP, Architecture, Ports, Reals, Stage,
                      StateGrid, System, solve_dynamic)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
burn = Architecture("burn", s.build())
st = Stage("leg", lambda s: {"demand": 1.0},
           lambda s, p: s - p["fuel"], candidates=[burn])
policy = solve_dynamic([st], StateGrid.linspace(0.0, 10.0, 11),
                       cost_fn=lambda p: p["fuel"])
print(policy.cost_to_go(0, 8.0), policy.action_at(0, 8.0))
# -> 2.0 ('burn', {'fuel': 2.0}, 2.0, 6.0)

Errors: none raised directly (an out-of-table state returns inf / None rather than raising).

13.5.2.5 DynamicResult

The outcome of a policy rolled out from an initial state, returned by rollout and solve_and_rollout. A dataclass with a schedule property and a compact __repr__; it also carries the full policy for closed-loop re-querying.

from codesign import (AlgebraicDP, Architecture, Ports, Reals, Stage,
                      StateGrid, System, solve_and_rollout)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
burn = Architecture("burn", s.build())
st = Stage("leg", lambda s: {"demand": 1.0},
           lambda s, p: s - p["fuel"], candidates=[burn])
grid = StateGrid.linspace(0.0, 10.0, 11)
res = solve_and_rollout([st, st], grid, 10.0, cost_fn=lambda p: p["fuel"])
print(res.schedule, res.total_cost, res.feasible)
# -> ['burn', 'burn'] 4.0 True

• stages — List[StageResult] — per-stage record of the optimal policy from the initial state.

• total_cost — float — optimal cost-to-go from the initial state (running plus terminal).

• feasible — bool — True iff a finite-cost policy exists from the initial state.

• policy — DynamicPolicy — the full state-indexed policy table.

• schedule (property) — List[str] — architecture chosen at each stage along the rolled-out path.

Errors: none raised directly.

13.5.2.6 solve_dynamic

Solves the finite-horizon architecture DP by a backward Bellman pass over the (stage, state) lattice, returning a full state-indexed DynamicPolicy (13.5.2.4). At each (stage, state) it solves the co-design problem for each admissible architecture at the stage’s state-dependent functionality, reads the cheapest feasible point, advances the state, and adds the discretised cost-to-go of the successor.

solve_dynamic(stages, grid, *, cost_fn, architectures=None,
              terminal_cost=None, max_iter=200, solve_kwargs=None,
              cache=True) -> DynamicPolicy

from codesign import (AlgebraicDP, Architecture, Ports, Reals, Stage,
                      StateGrid, System, solve_dynamic)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
burn = Architecture("burn", s.build())
st = Stage("leg", lambda s: {"demand": 1.0},
           lambda s, p: s - p["fuel"], candidates=[burn])
grid = StateGrid.linspace(0.0, 10.0, 11)
policy = solve_dynamic([st, st], grid, cost_fn=lambda p: p["fuel"])
print(policy.stage_names, policy.cost_to_go(0, 10.0))
# -> ['leg', 'leg'] 4.0

• stages — Sequence[Stage] — the horizon in forward order (solved backward internally).

• grid — StateGrid — discretisation of the carried scalar resource.

• cost_fn — callable — scalar cost on a resource point; lower is better.

• architectures — Sequence[Architecture] or None — default candidates for stages without their own.

• terminal_cost — callable or None — terminal_cost(state) -> float scoring the state after the final stage; defaults to zero everywhere.

• max_iter — int — forwarded to solve().

• solve_kwargs — Mapping or None — extra keyword arguments forwarded to solve().

• cache — bool — cache co-design solves keyed by (architecture name, functionality); safe when functionality(state) is deterministic.

Errors: raises ValueError when a stage supplies no candidates and no default architectures set was passed; the message reads

stage <name> has no candidates and no default architecture set
was supplied. Either set candidates=[...] on the stage, or pass
architectures=[...] to solve_dynamic().

13.5.2.7 rollout

Rolls a solved policy forward from a concrete initial state, following the policy’s chosen architecture at each stage, advancing the snapped carried state, and assembling a DynamicResult (13.5.2.5).

rollout(policy, stages, initial_state, *, arch_lookup=None)
    -> DynamicResult

from codesign import (AlgebraicDP, Architecture, Ports, Reals, Stage,
                      StateGrid, System, rollout, solve_dynamic)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
burn = Architecture("burn", s.build())
st = Stage("leg", lambda s: {"demand": 1.0},
           lambda s, p: s - p["fuel"], candidates=[burn])
grid = StateGrid.linspace(0.0, 10.0, 11)
policy = solve_dynamic([st, st], grid, cost_fn=lambda p: p["fuel"])
res = rollout(policy, [st, st], 10.0)
print(res.schedule, res.total_cost)   # -> ['burn', 'burn'] 4.0

• policy — DynamicPolicy — a policy from solve_dynamic.

• stages — Sequence[Stage] — the same stage list passed to solve_dynamic.

• initial_state — float — the starting carried value.

• arch_lookup — Mapping or None — maps architecture name to Architecture so tags attach to each record; built automatically from the stage candidate lists when omitted.

Errors: none raised directly. When the policy has no feasible action from the reached state the roll-out stops early, the offending record is marked feasible=False, and DynamicResult.feasible becomes False.

13.5.2.8 solve_and_rollout

A convenience wrapper equivalent to solve_dynamic followed by rollout from initial_state. Prefer the two-step form when the policy will be queried at several initial states or in closed loop.

solve_and_rollout(stages, grid, initial_state, *, cost_fn,
                  architectures=None, terminal_cost=None, max_iter=200,
                  solve_kwargs=None, cache=True) -> DynamicResult

from codesign import (AlgebraicDP, Architecture, Ports, Reals, Stage,
                      StateGrid, System, solve_and_rollout)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
burn = Architecture("burn", s.build())
st = Stage("leg", lambda s: {"demand": 1.0},
           lambda s, p: s - p["fuel"], candidates=[burn])
grid = StateGrid.linspace(0.0, 10.0, 11)
res = solve_and_rollout([st, st], grid, 10.0, cost_fn=lambda p: p["fuel"])
print(res.schedule, res.total_cost, res.feasible)
# -> ['burn', 'burn'] 4.0 True

Parameters: as solve_dynamic (13.5.2.6), plus the positional initial_state (float) rolled out from. Errors: same ValueError as solve_dynamic when candidates are missing.

13.5.3 sequential↩︎

The sequential module is the antichain-valued generalisation of the scalar dynamic layer (section 15.3): a finite-horizon decision process whose stages are co-design problems coupled by a carried state, solved by a backward Bellman recursion whose value at each (stage, state) is a Pareto antichain rather than a scalar. The accumulated resource (the named cost axes on the antichain) is kept deliberately distinct from the carried state: the transition reads any carried quantity off the full solved point, while only the cost axes accumulate on the front. The scalar dynamic layer is the width-one special case. This module also makes three theory results operational — monotone value (Q1, check_monotonicity), front \(=\) reachable frontier (Q2), and exact factorisation at a reset (Q3, detect_resets / factorise_at_resets) — and exposes the precompute-then-DP structure (section 15.5). It reuses the Architecture (13.5.1.1) and StateGrid (13.5.2.3).

13.5.3.1 SeqStage

One stage of a sequential co-design problem. Structurally the antichain- valued analogue of Stage (13.5.2.1): returning the grid’s bottom node from transition makes the stage a reset. A dataclass.

SeqStage(name: str,
         functionality: Callable[[float], Mapping],
         transition: Callable[[float, Mapping], float],
         candidates: Optional[Sequence[Architecture]] = None,
         admissible: Optional[Callable[[float], bool]] = None)

from codesign import Architecture, SeqStage
mode = Architecture("m", None)
st = SeqStage("leg", functionality=lambda x: {"demand": 1.0},
              transition=lambda x, p: x - p["fuel"],
              candidates=[mode], admissible=lambda x: x >= 0.0)
print(st.name, st.transition(6.0, {"fuel": 2.0}))   # -> leg 4.0

• name — str — stage identifier.

• functionality — callable — functionality(state) -> Mapping, the outer functionality demanded when the carried state is state.

• transition — callable — transition(state, point) -> new_state, from the incoming state and a chosen resource point (over the resource ports) to the outgoing state.

• candidates — Sequence[Architecture] or None — architectures admissible here; each one’s co-design solve supplies part of the stage antichain \(h_k(x)\) (the union of their fronts). Falls back to the problem-level default when None.

• admissible — callable or None — admissible(state) -> bool; forbidden states carry the empty (infeasible) value.

Errors: none raised directly.

13.5.3.2 sum_combine

The consumable/accumulating monoid combination \(\oplus\): component-wise sum over the cost axes. Resources accumulate across stages (fuel burned, money spent, wear incurred); the reachable frontier can grow polynomially with the horizon. The default combine for solve_sequential.

from codesign import sum_combine
a = {"cost": 5.0, "co2": 1.0}
b = {"cost": 2.0, "co2": 8.0}
print(sum_combine(a, b))           # -> {'cost': 7.0, 'co2': 9.0}

• a, b — Mapping — two resource points over the same axes; returns a dict of their component-wise sum.

Errors: none raised directly (a KeyError results if b lacks an axis present in a).

13.5.3.3 join_combine

The renewable monoid combination \(\oplus\): component-wise maximum (join). The requirement of two stages is the peak, not the sum (a worker, an oven, a bus reused across stages); the reachable frontier stays bounded uniformly in the horizon, and if every stage resets the problem collapses to a static co-design.

from codesign import join_combine
a = {"cost": 5.0, "co2": 1.0}
b = {"cost": 2.0, "co2": 8.0}
print(join_combine(a, b))          # -> {'cost': 5.0, 'co2': 8.0}

• a, b — Mapping — two resource points over the same axes; returns a dict of their component-wise maximum.

Errors: none raised directly (a KeyError results if b lacks an axis present in a).

13.5.3.4 solve_sequential

Solves the antichain-valued sequential co-design problem by the backward Bellman recursion in the upper-set value space, carrying a full Pareto antichain of cumulative resources at each (stage, state). Returns a SeqResult (13.5.3.5) whose value is \(V_0(\texttt{initial\_state})\). The snippet below uses a single architecture (a width-one front); with several incomparable candidates the value has width greater than one — see dp_over_catalog (13.5.3.12) for a multi-point front.

solve_sequential(stages, grid, *, cost_axes, initial_state,
                 combine=sum_combine, architectures=None,
                 max_iter=200, solve_kwargs=None) -> SeqResult

from codesign import (AlgebraicDP, Architecture, Ports, Reals, SeqStage,
                      StateGrid, System, solve_sequential)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
a = Architecture("a", s.build())
st = SeqStage("leg", lambda x: {"demand": 1.0},
              lambda x, p: x - p["fuel"], candidates=[a],
              admissible=lambda x: x >= -1e-9)
res = solve_sequential([st, st], StateGrid.linspace(0., 6., 13),
                       cost_axes=["fuel"], initial_state=6.)
print(res.width, sorted(p["fuel"] for p in res.value), res.feasible)
# -> 1 [4.0] True

• stages — Sequence[SeqStage] — the horizon in forward order (solved backward internally).

• grid — StateGrid — discretisation of the carried scalar state.

• cost_axes — Sequence[str] — names of the resource ports to accumulate on the antichain (for example ["cost", "co2"]); these define the product poset \(R\).

• initial_state — float — starting carried state; the returned value is \(V_0\) of it.

• combine — callable — the monoid combination \(\oplus\): sum_combine (default) for a consumable resource, join_combine for a renewable one.

• architectures — Sequence[Architecture] or None — default candidates for stages without their own.

• max_iter — int — forwarded to solve().

• solve_kwargs — Mapping or None — extra keyword arguments forwarded to solve().

Errors: raises ValueError when a stage supplies no candidates and no default architectures set was passed; the message reads

stage <name> has no candidates and no default architecture set
was supplied. Either set candidates=[...] on the stage, or pass
architectures=[...] to the sequential solver (solve_sequential /
check_monotonicity).

13.5.3.5 SeqResult

The outcome of solve_sequential (and dp_over_catalog): the whole-horizon Pareto front from the initial state. A dataclass with a compact __repr__. The snippet obtains one from dp_over_catalog to show a genuine multi-point front.

from codesign import StateGrid, dp_over_catalog
cat = [("clean", {"cost": 10.0, "co2": 1.0}),
       ("cheap", {"cost": 2.0, "co2": 8.0})]
res = dp_over_catalog([cat, cat], StateGrid.linspace(0., 4., 5),
                      cost_axes=["cost", "co2"], initial_state=4.0,
                      transition=lambda s, p: s)
print(res.width, res.feasible,
      sorted((p["cost"], p["co2"]) for p in res.value))
# -> 3 True [(4.0, 16.0), (12.0, 9.0), (20.0, 2.0)]

• value — Antichain — the value antichain \(\mathrm{Min}\,V_0(\texttt{initial\_state})\); each point is a Mapping over the cost axes.

• width — int — len(value), the number of incomparable Pareto-optimal totals (the reachable-frontier width).

• feasible — bool — True iff the value antichain is non-empty and free of tops.

• policy — SeqPolicy — the full state-indexed antichain-valued policy table.

Errors: none raised directly.

13.5.3.6 SeqPolicy

The state-indexed antichain-valued policy from the backward pass, held on SeqResult.policy. For each (stage, state node) it stores the value antichain \(V_k(x)\) and, per Pareto point, the realising architecture and chosen stage point, so an optimal choice sequence can be traced forward.

from codesign import StateGrid, dp_over_catalog
cat = [("clean", {"cost": 10.0, "co2": 1.0}),
       ("cheap", {"cost": 2.0, "co2": 8.0})]
res = dp_over_catalog([cat, cat], StateGrid.linspace(0., 4., 5),
                      cost_axes=["cost", "co2"], initial_state=4.0,
                      transition=lambda s, p: s)
pol = res.policy
print(pol.width_at(0, 4.0), pol.stage_names)
# -> 3 ['stage_0', 'stage_1']

Errors: none raised directly.

13.5.3.7 MonotonicityReport

The result of check_monotonicity: whether the three hypotheses of the monotone-value theorem — (H1), the joint (H2), and (H3) — held on the grid, with sampled witnesses when they did not. A dataclass whose monotone_value_guaranteed property is True exactly when all hold (the Q1 theorem then guarantees a monotone value).

from codesign import MonotonicityReport
rep = MonotonicityReport(h1_ok=True, h2_ok=True, h2_joint_ok=True, h3_ok=True)
print(rep.monotone_value_guaranteed, repr(rep))
# -> True MonotonicityReport(H1=ok, H2=ok, H2joint=ok, H3=ok,
#         value_monotone=guaranteed)

• h1_ok — bool — (H1) held: the stage map h_k is consistently oriented at every tested grid pair.

• h2_ok — bool — the state slice of (H2) held: phi_k(.,r) is monotone in the state at every pair.

• h2_joint_ok — bool — the resource slice of (H2) held, so the transition is jointly monotone on X x R (the paper’s (H2-joint)): a larger resource sends the successor no lower in the certified orientation.

• h3_ok — bool — (H3) held: the admissible region is a down-set of X in the certified orientation.

• h1_violations, h2_violations, h2_joint_violations, h3_violations — List[Tuple[str, float, float]] — sampled witnesses (the h2_joint witness is the out-of-order successor pair (stage, phi_x_r, phi_x_rp); the others are (stage, x, x’)), empty when the condition held.

• monotone_value_guaranteed (property) — h1_ok and h2_ok and h2_joint_ok and h3_ok.

Errors: none raised directly.

13.5.3.8 check_monotonicity

Numerically verifies the three hypotheses of the monotone-value theorem on the state grid — (H1) consistent orientation of h_k, the joint (H2) monotonicity of the transition on X x R (both its state slice h2_ok and its resource slice h2_joint_ok), and the (H3) down-set property of admissibility — and returns a MonotonicityReport (13.5.3.7). It is orientation-aware: it accepts the benign “consumable but monotone” orientation (a budget carried as remaining slack) and flags only genuinely non-monotone stages (perishable / fatigue-as-state, with an interior optimum), a transition that is not jointly monotone, or a non-down-set admissible region, all of which fall outside the Q1 guarantee. Spurious violations most often signal a too-coarse grid, since snapping is not order-preserving at bucket boundaries.

check_monotonicity(stages, grid, *, cost_axes, architectures=None,
                   max_iter=200, solve_kwargs=None,
                   max_violations=8) -> MonotonicityReport

from codesign import (AlgebraicDP, Architecture, Ports, Reals, SeqStage,
                      StateGrid, System, check_monotonicity)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
a = Architecture("a", s.build())
st = SeqStage("leg", lambda x: {"demand": 1.0},
              lambda x, p: x - p["fuel"], candidates=[a],
              admissible=lambda x: x >= -1e-9)
rep = check_monotonicity([st, st], StateGrid.linspace(0., 6., 13),
                         cost_axes=["fuel"])
print(rep.h1_ok, rep.h2_ok, rep.h2_joint_ok, rep.h3_ok,
      rep.monotone_value_guaranteed)
# -> True True True True True

• stages, grid, cost_axes, architectures, max_iter, solve_kwargs — as solve_sequential (13.5.3.4).

• max_violations — int — cap on the number of sampled witnesses recorded per hypothesis.

Errors: raises the same ValueError as solve_sequential when a stage has no candidates and no default set (the message names both “solve_sequential / check_monotonicity”).

13.5.3.9 detect_resets

Returns the indices of stages that are resets: a stage is a reset iff its transition lands on the grid’s bottom node for every incoming state and every chosen resource point. At a reset the horizon factorises exactly (Q3), in every regime, since the proof uses only distributivity.

detect_resets(stages, grid, *, cost_axes, architectures=None,
              max_iter=200, solve_kwargs=None) -> List[int]

from codesign import (AlgebraicDP, Architecture, Ports, Reals, SeqStage,
                      StateGrid, System, detect_resets)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
a = Architecture("a", s.build())
run = SeqStage("run", lambda x: {"demand": 1.0},
               lambda x, p: x - p["fuel"], candidates=[a])
reset = SeqStage("reset", lambda x: {"demand": 1.0},
                 lambda x, p: 0.0, candidates=[a])
print(detect_resets([run, reset, run],
                    StateGrid.linspace(0., 6., 13), cost_axes=["fuel"]))
# -> [1]

• stages, grid, cost_axes, architectures, max_iter, solve_kwargs — as solve_sequential (13.5.3.4).

Errors: raises ValueError (naming the stage) when a stage has no candidates and no default architectures was supplied — the same guard as solve_sequential and check_monotonicity, so an under-specified stage fails loudly instead of being silently (and vacuously) reported as a reset.

13.5.3.10 factorise_at_resets

Partitions the horizon into maximal quiescence-free runs given the reset indices, returning inclusive (start, end) ranges. By Q3 the value over the full horizon is the \(\oplus\)-product of the values of these independent sub-problems — the order-theoretic, deterministic analogue of regeneration-point decomposition. Pure bookkeeping: it solves nothing.

factorise_at_resets(stages, resets) -> List[Tuple[int, int]]

from codesign import SeqStage, factorise_at_resets
f = lambda x: {"demand": 1.0}
t = lambda x, p: x
stages = [SeqStage(f"s{i}", f, t) for i in range(4)]
print(factorise_at_resets(stages, resets=[1]))
# -> [(0, 1), (2, 3)]

• stages — Sequence[SeqStage] — the horizon (only its length is used).

• resets — Sequence[int] — reset stage indices, typically from detect_resets (13.5.3.9). A reset closes the run containing it.

Errors: none raised directly.

13.5.3.11 precompute_catalog

Solves every architecture once at a fixed functionality and returns a flat, tagged Pareto catalog — the co-design precomputation step of the Formula 1 seasonal framework (section 15.5). Points that look dominated at this single-stage level are retained if non-dominated in the combined antichain, since they may become optimal once aggregated in an outer DP.

precompute_catalog(architectures, functionality, cost_axes, *,
                   max_iter=200, solve_kwargs=None)
    -> List[Tuple[str, Mapping]]

from codesign import (AlgebraicDP, Architecture, Ports, Reals,
                      System, precompute_catalog)
s = System("m"); d = s.provides("demand"); s.requires("fuel")
s.add("m", AlgebraicDP(F=Ports({"demand": Reals()}),
      R=Ports({"fuel": Reals()}),
      equations={"fuel": lambda f: 2.0})).demand >= d
s.constrain("fuel", lambda x: x["m.fuel"])
a = Architecture("a", s.build())
cat = precompute_catalog([a], {"demand": 1.0}, cost_axes=["fuel"])
print(cat)                         # -> [('a', {'fuel': 2.0})]

• architectures — Sequence[Architecture] — the candidates to solve once.

• functionality — Mapping — the fixed outer functionality at which each architecture is solved.

• cost_axes — Sequence[str] — resource ports accumulated on the catalog’s antichain.

• max_iter, solve_kwargs — forwarded to solve().

Returns a list of (arch_name, point) pairs (each point the full solved resource point, so a carried-state axis remains readable); feed it to dp_over_catalog (13.5.3.12). Errors: none raised directly.

13.5.3.12 dp_over_catalog

Runs the antichain-valued DP that selects from precomputed catalogs: each stage’s action set is a frozen catalog rather than a live co-design solve, so no solve() happens inside the Bellman sweep. Use it when the per-stage co-design is independent of the carried state (the Formula 1 regime); otherwise use solve_sequential. Returns a SeqResult (13.5.3.5).

dp_over_catalog(catalogs, grid, *, cost_axes, initial_state,
                transition, combine=sum_combine,
                admissible=None) -> SeqResult

from codesign import StateGrid, dp_over_catalog, sum_combine
cat = [("clean", {"cost": 10.0, "co2": 1.0}),
       ("cheap", {"cost": 2.0, "co2": 8.0})]
res = dp_over_catalog([cat, cat], StateGrid.linspace(0., 4., 5),
                      cost_axes=["cost", "co2"], initial_state=4.0,
                      transition=lambda s, p: s, combine=sum_combine)
print(res.width, sorted((p["cost"], p["co2"]) for p in res.value))
# -> 3 [(4.0, 16.0), (12.0, 9.0), (20.0, 2.0)]

• catalogs — Sequence[Sequence[Tuple[str, Mapping]]] — one catalog per stage, each a sequence of (arch_name, point) pairs from precompute_catalog.

• grid — StateGrid — discretisation of the carried scalar state.

• cost_axes — Sequence[str] — resource ports accumulated on the antichain.

• initial_state — float — starting carried state.

• transition — callable — transition(state, point) -> new_state.

• combine — callable — monoid combination \(\oplus\) (default sum_combine).

• admissible — callable or None — state admissibility predicate.

Errors: none raised directly.

13.6 State and closed loop↩︎

The vector-state generalisation (state, vector_dp) and the closed-loop online_codesign controller. See Sections 15.4 and 15.6.

13.6.1 state↩︎

The state module supplies the general carried state of the temporal layer (section 15.4): a hashable state-vector value type and a product grid of named axes carrying the component-wise order the monotonicity results are stated over. A one-axis grid reproduces the scalar StateGrid exactly, so the vector layer strictly generalises the scalar carried state rather than parallelling it.

13.6.1.1 StateVec

The hashable value type of a carried state vector: a tuple of (axis_name, value) pairs sorted by name, so two states with the same contents compare and hash equal. It is a type alias, not a class — StateVec = Tuple[Tuple[str, Any], ...] — and is used directly as the key of the dynamic-program table. Build one with make_state and read it with state_get / state_as_dict rather than by indexing the raw tuple.

from codesign import make_state
sv = make_state(fuel=12.0, flag=0)  # build a StateVec
print(sv)          # (('flag', 0), ('fuel', 12.0)): sorted, hashable
table = {sv: 0.0}  # usable directly as a DP-table key
print(table[make_state(flag=0, fuel=12.0)])  # 0.0: order-independent

• Contents — one (str, Any) pair per axis; canonicalised (sorted by name) so keyword order at construction is irrelevant.

Errors: none raised directly (it is a plain tuple).

13.6.1.2 make_state

Build a canonical StateVec from keyword axis values.

make_state(**axis_values) -> StateVec

from codesign import make_state
x = make_state(fuel=12.0, flag=0)
print(x)                     # (('flag', 0), ('fuel', 12.0))
print(isinstance(x, tuple))  # True: hashable, usable as a DP-table key

• **axis_values — each keyword names an axis and gives its value; the result is normalised so equal contents hash and compare equal.

Errors: none raised directly.

13.6.1.3 state_get

Read one axis value from a state vector. Signature: state_get(state: StateVec, axis: str) -> Any.

from codesign import make_state, state_get
x = make_state(fuel=12.0, flag=0)
print(state_get(x, "fuel"))  # 12.0

• state — the state vector to read.

• axis — the axis name whose value is wanted.

Errors: raises KeyError when the axis is absent, with the message "state vector has no axis ’charge’; it carries axes [’flag’, ’fuel’]" — the message names the missing axis and lists the axes the vector actually carries.

13.6.1.4 state_as_dict

Return a plain dict view of a state vector, convenient inside a transition or functionality callback. Signature: state_as_dict(state: StateVec) -> Dict[str, Any].

from codesign import make_state, state_as_dict
x = make_state(fuel=12.0, flag=0)
print(state_as_dict(x))  # {'flag': 0, 'fuel': 12.0}

• state — the state vector to expand.

Errors: none raised directly.

13.6.1.5 Axis

Abstract base class defining the single-axis contract that every concrete axis satisfies and that VectorStateGrid composes into the product order. An axis carries a name and the four operations below; the two concrete axes (ContinuousAxis, DiscreteAxis) are documented as deltas from this contract. The snippet exercises the contract through a concrete axis, since Axis itself is not instantiable.

from codesign import ContinuousAxis, Axis
ax = ContinuousAxis("charge", 0.0, 1.0, 5)  # a concrete Axis
print(isinstance(ax, Axis))   # True
print(list(ax.nodes()))       # [0.0, 0.25, 0.5, 0.75, 1.0]
print(ax.snap(0.3))           # 0.25 (nearest node)
print(ax.in_bounds(1.2))      # False
print(ax.leq(0.25, 0.5))      # True

Errors: the four base methods raise NotImplementedError; use a concrete axis.

13.6.1.6 ContinuousAxis

A bucketed real interval, the vector analogue of the scalar grid and a concrete Axis (13.6.1.5). It places n evenly spaced nodes on [lo, hi]; snap returns the nearest node, in_bounds tests [lo, hi] within a tolerance, and leq applies the interval order, reversed when increasing_is_larger is False. Signature: ContinuousAxis(name, lo, hi, n, increasing_is_larger=True).

from codesign import ContinuousAxis
ax = ContinuousAxis("slack", 0.0, 10.0, 3, increasing_is_larger=False)
print(list(ax.nodes()))  # [0.0, 5.0, 10.0]
print(ax.leq(10.0, 0.0))  # True: reversed, so less slack ranks higher

Errors: raises ValueError when n < 1 (a message stating the axis needs at least one node and that n=1 yields a single node at lo).

13.6.1.7 DiscreteAxis

A finite labelled axis with a caller-supplied order, a concrete Axis (13.6.1.5) for quantities like a regulatory flag or an operating regime. Its nodes are the admissible labels; snap is the identity on admissible labels (labels are used verbatim); and leq follows order (least first), treating labels omitted from order, or all labels when order is None, as mutually incomparable — which the monotonicity guard reads as “no order to exploit on this axis”. Signature: DiscreteAxis(name, values, order=None).

from codesign import DiscreteAxis
flag = DiscreteAxis("penalty", values=[0, 1], order=[0, 1])
print(flag.nodes())      # [0, 1]
print(flag.leq(0, 1))    # True: 0 ranks below 1
print(flag.snap(1))      # 1 (snap is identity on admissible labels)
free = DiscreteAxis("regime", values=["a", "b"])  # no order given
print(free.leq("a", "b"))  # False: distinct labels are incomparable

Errors: raises ValueError when values is empty (a message asking for the admissible labels).

13.6.1.8 VectorStateGrid

A product grid over several named axes, carrying a full state vector. It enumerates the Cartesian product of its axes’ nodes as the DP state set, snaps and bounds-checks a proposed successor axis by axis, and exposes the component-wise product order used by the monotonicity results: a <= b iff a is at or below b on every axis. A single-axis grid reproduces the scalar StateGrid. Signature: VectorStateGrid(axes: Sequence[Axis]).

from codesign import VectorStateGrid, ContinuousAxis, DiscreteAxis
grid = VectorStateGrid([
    ContinuousAxis("charge", 0.0, 1.0, 3),
    DiscreteAxis("flag", values=[0, 1], order=[0, 1]),
])
print(len(grid))         # 6 state nodes (3 x 2)
print(grid.bottom())     # (('charge', 0.0), ('flag', 0))
s = grid.snap({"charge": 0.4, "flag": 1})
print(s)                 # (('charge', 0.5), ('flag', 1))
print(grid.in_bounds({"charge": 2.0, "flag": 0}))  # False
print(grid.leq(grid.bottom(), s))  # True

• axes — one or more Axis objects with distinct names; their order fixes the product.

Errors: raises ValueError on an empty axis sequence, and on duplicate axis names with the message "VectorStateGrid axis names must be unique, but [’wear’] appear more than once in [’wear’, ’wear’]. Give each axis a distinct name."

13.6.2 vector_dp↩︎

The vector_dp module is the general carried-state sequential solver (section 15.4): the antichain-valued Bellman recursion of solve_sequential (13.5.3.4) carrying a full state vector on a VectorStateGrid rather than a single scalar. Its card structure mirrors sequential throughout, with the product order of the vector grid substituted for the scalar order.

13.6.2.1 VecStage

One stage of a vector-state sequential co-design problem: a dataclass bundling the stage’s state-dependent functionality, its carried-state transition, and (optionally) its candidate architectures and admissibility predicate. The vector-state analogue of sequential’s SeqStage.

from codesign import VecStage, state_get
stage = VecStage(
    name="leg",
    functionality=lambda s: {"demand": state_get(s, "wear")},
    transition=lambda s, p: {"wear": dict(s)["wear"] + p["d_wear"]},
)
print(stage.name, stage.candidates, stage.admissible)  # leg None None

Errors: none raised directly (a stage with no candidates and no solver-level default is reported by the solver; see below).

13.6.2.2 VecResult

The outcome of solve_vector_sequential: a dataclass holding the value antichain at the initial state, its width, feasibility, and the full policy for roll-out.

from codesign import (AlgebraicDP, Architecture, Ports, Reals, VecStage,
    VectorStateGrid, ContinuousAxis, solve_vector_sequential)
dp = AlgebraicDP(F=Ports({"demand": Reals()}),
    R=Ports({"cost": Reals(), "d_wear": Reals()}),
    equations={"cost": 1.0, "d_wear": 1.0})
stage = VecStage("leg", lambda s: {"demand": 1.0},
    lambda s, p: {"wear": dict(s)["wear"] + p["d_wear"]},
    candidates=[Architecture("run", dp)])
grid = VectorStateGrid([ContinuousAxis("wear", 0.0, 4.0, 5)])
res = solve_vector_sequential([stage], grid,
    cost_axes=["cost"], initial_state={"wear": 0.0})
print(res.feasible, res.width)  # True 1
print(res.value)                # Antichain[(cost=1)]

Errors: none raised directly.

13.6.2.3 VecPolicy

A state-vector-indexed, antichain-valued policy. It is produced by solve_vector_sequential (read from VecResult.policy) rather than constructed directly, and it answers value and best-action queries at any (stage, state), snapping the queried state onto the grid first.

from codesign import (AlgebraicDP, Architecture, Ports, Reals, VecStage,
    VectorStateGrid, ContinuousAxis, solve_vector_sequential, make_state)
dp = AlgebraicDP(F=Ports({"demand": Reals()}),
    R=Ports({"cost": Reals(), "d_wear": Reals()}),
    equations={"cost": 1.0, "d_wear": 1.0})
stage = VecStage("leg", lambda s: {"demand": 1.0},
    lambda s, p: {"wear": dict(s)["wear"] + p["d_wear"]},
    candidates=[Architecture("run", dp)])
grid = VectorStateGrid([ContinuousAxis("wear", 0.0, 4.0, 5)])
pol = solve_vector_sequential([stage], grid,
    cost_axes=["cost"], initial_state={"wear": 0.0}).policy
s0 = make_state(wear=0.0)
print(pol.width_at(0, s0))                       # 1
print([p["cost"] for p in pol.value_at(0, s0)])  # [1.0]
print(pol.best_action_at(0, s0))  # ('run', {...}, (('wear', 1.0),))

13.6.2.3.1 value_at(k, state)

the value antichain at stage k and (snapped) state; an empty antichain if the node is unreached.

13.6.2.3.2 width_at(k, state)

the width (front size) of value_at(k, state).

13.6.2.3.3 best_action_at(k, state)

one realising (arch, full_point, succ_state) at (k, state), picking the choice behind the cost axis’s minimal-first point (a myopic roll-out); None when no finite choice exists.

Errors: none raised directly.

13.6.2.4 solve_vector_sequential

Solve the vector-state antichain-valued sequential co-design problem. Identical in structure to solve_sequential (13.5.3.4) but carrying a full state vector on a VectorStateGrid; the transition returns a {axis: value} mapping snapped onto the product grid, and out-of-bounds successors are rejected before snapping on every axis.

solve_vector_sequential(stages, grid, *, cost_axes, initial_state,
    combine=sum_combine, architectures=None, max_iter=200,
    solve_kwargs=None)

from codesign import (AlgebraicDP, Architecture, Ports, Reals, VecStage,
    VectorStateGrid, ContinuousAxis, solve_vector_sequential)
dp = AlgebraicDP(F=Ports({"demand": Reals()}),
    R=Ports({"cost": Reals(), "d_wear": Reals()}),
    equations={"cost": 1.0, "d_wear": 1.0})
stage = VecStage("leg", lambda s: {"demand": 1.0},
    lambda s, p: {"wear": dict(s)["wear"] + p["d_wear"]},
    candidates=[Architecture("run", dp)])
grid = VectorStateGrid([ContinuousAxis("wear", 0.0, 4.0, 5)])
res = solve_vector_sequential([stage, stage], grid,
    cost_axes=["cost"], initial_state={"wear": 0.0})
print(res)                              # VecResult(feasible, width=1, ...)
print([p["cost"] for p in res.value])   # [2.0]: two legs, unit cost each

Errors: raises ValueError when a stage has neither candidates nor a solver-level architectures default (the message names the offending stage and both remedies).

13.6.2.5 VectorMonotonicityReport

The result of check_vector_monotonicity: a dataclass recording whether (H1), the joint (H2), and (H3) hold over the product order, with sample violations.

from codesign import VectorMonotonicityReport
rep = VectorMonotonicityReport(h1_ok=True, h2_ok=True,
                               h2_joint_ok=True, h3_ok=True)
print(rep.monotone_value_guaranteed)  # True: all hypotheses hold
print(rep)  # VectorMonotonicityReport(H1=ok, H2=ok, H2joint=ok, H3=ok, ...)

The property monotone_value_guaranteed is h1_ok and h2_ok and h2_joint_ok and h3_ok. Errors: none raised directly.

13.6.2.6 check_vector_monotonicity

Verify (H1), the joint (H2), and (H3) over the vector grid’s product order, the vector-state analogue of the scalar guard. For every ordered comparable pair x <= x’ it checks the stage antichain is consistently oriented (a benign consumable-but-monotone axis passes; a genuinely non-monotone perishable one is flagged), that the transition is jointly monotone on X x R (both its state slice and its resource slice), and that admissibility is a down-set. The product grid can be large, so at most max_pairs comparable pairs per stage are sampled; the report is exact when the grid is small enough that the cap is not reached.

check_vector_monotonicity(stages, grid, *, cost_axes,
    architectures=None, max_iter=200, solve_kwargs=None,
    max_violations=8, max_pairs=4000)

from codesign import (AlgebraicDP, Architecture, Ports, Reals, VecStage,
    VectorStateGrid, ContinuousAxis, check_vector_monotonicity)
dp = AlgebraicDP(F=Ports({"demand": Reals()}),
    R=Ports({"cost": Reals(), "d_wear": Reals()}),
    equations={"cost": 1.0, "d_wear": 1.0})
stage = VecStage("leg", lambda s: {"demand": 1.0},
    lambda s, p: {"wear": dict(s)["wear"] + p["d_wear"]},
    candidates=[Architecture("run", dp)])
grid = VectorStateGrid([ContinuousAxis("wear", 0.0, 4.0, 5)])
rep = check_vector_monotonicity([stage], grid, cost_axes=["cost"])
print(rep)  # H1=ok, H2=ok, H2joint=ok, H3=ok, value_monotone=...
print(rep.monotone_value_guaranteed)  # True

• stages, grid, cost_axes, architectures, max_iter, solve_kwargs — as for the solver above.

• max_violations — cap on recorded witnesses per hypothesis.

• max_pairs — cap on comparable pairs sampled per stage; below it the check is exact.

Errors: raises the same ValueError as the solver when a stage has no candidates and no default architecture set.

13.6.3 online_codesign↩︎

The online_codesign module is the closed-loop, measurement-driven counterpart of the open-loop temporal planners (section 15.6, worked in example 15.11). At each control step it senses the measured state, reads the live requirement and environment, re-solves the co-design over the admissible architectures, applies the cheapest feasible choice by stepping the true plant, and logs the outcome, so divergence between plan and execution is absorbed by the next sensing step.

13.6.3.1 ControlStep

The audit record of one closed-loop control step: a dataclass capturing the step’s measured state, conditions, chosen architecture and point, cost, and feasibility. run_online_codesign emits one per step.

from codesign import ControlStep
step = ControlStep(step=0, measured_state=10.0, requirement={"rate": 2.0},
    environment={}, architecture="node", point={"energy": 2.0, "ops": 3.0},
    cost=3.0, feasible=True)
print(step)  # ControlStep(t=0, arch='node', cost=3, ok)

Errors: none raised directly.

13.6.3.2 OnlineCoDesignResult

The trajectory produced by run_online_codesign: a dataclass holding the per-step log, the total cost of feasible steps, and whether every step was feasible. The schedule property extracts the sequence of chosen architectures.

from codesign import OnlineCoDesignResult, ControlStep
s = ControlStep(0, 10.0, {"rate": 2.0}, {}, "node",
    {"energy": 2.0, "ops": 3.0}, 3.0, True)
res = OnlineCoDesignResult(steps=[s], total_cost=3.0, feasible=True)
print(res.schedule)  # ['node']
print(res)           # OnlineCoDesignResult(feasible, steps=1, ...)

The property schedule returns [s.architecture for s in steps]. Errors: none raised directly.

13.6.3.3 resolve_at

Solve every architecture at a fixed functionality and pick the cheapest feasible point by cost_fn. This is the single myopic co-design solve performed once per control step; it is exposed separately so a caller can drive the re-solve directly. Signature: resolve_at(architectures, functionality, cost_fn, *, max_iter=200, solve_kwargs=None).

from codesign import AlgebraicDP, Architecture, Ports, Reals, resolve_at
dp = AlgebraicDP(F=Ports({"rate": Reals()}),
    R=Ports({"energy": Reals(), "ops": Reals()}),
    equations={"energy": lambda f: f["rate"], "ops": 3.0})
name, point, cost = resolve_at([Architecture("node", dp)],
    {"rate": 2.0}, cost_fn=lambda p: p["ops"])
print(name, cost)       # node 3.0
print(point["energy"])  # 2.0

• architectures — the admissible configurations to solve.

• functionality — the live demand to solve each at.

• cost_fn — scalar cost on a resource point; the cheapest feasible point wins.

• max_iter, solve_kwargs — forwarded to solve.

Returns (arch_name, point, cost), or ("", None, inf) when no architecture is feasible. Errors: none raised directly.

13.6.3.4 run_online_codesign

Run the closed-loop, measurement-driven co-design loop for n_steps steps. Each step senses, reads the live requirement and environment, re-solves via resolve_at, applies the cheapest feasible choice by stepping the plant (the true process, not a nominal model), and logs a ControlStep; an infeasible step holds the state and is recorded without aborting the run.

run_online_codesign(architectures, *, n_steps, sensor, requirement,
    plant, cost_fn, environment=None, initial_state=None,
    max_iter=200, solve_kwargs=None)

from codesign import (AlgebraicDP, Architecture, Ports, Reals,
    run_online_codesign)
dp = AlgebraicDP(F=Ports({"rate": Reals()}),
    R=Ports({"energy": Reals(), "ops": Reals()}),
    equations={"energy": lambda f: f["rate"], "ops": 3.0})
res = run_online_codesign([Architecture("node", dp)], n_steps=3,
    sensor=lambda t, s: s,
    requirement=lambda t, s: {"rate": 2.0},
    plant=lambda t, s, a, p: s - p["energy"],
    cost_fn=lambda p: p["ops"], initial_state=10.0)
print(res)             # OnlineCoDesignResult(feasible, steps=3, ...)
print(res.schedule)    # ['node', 'node', 'node']
print(res.total_cost)  # 9.0

Errors: none raised directly. Fold environment into requirement when it should influence the solve; it is otherwise only recorded.

14 Worked Examples↩︎

The package ships 25 runnable examples under examples/ and the same 25 as executed notebooks under notebooks/. Each is referenced below by its filename. This section walks through examples 1–17; the temporal, vector-state, and online examples 18–25 are covered in Section 15.

Examples 1 to 5 are reproductions of models published by [1], and the section headings name the figure or section of that paper each one comes from. Example 25 replicates the synthetic benchmarks of [18]. The remaining examples are the author’s, though several are modelled on published case studies: the mobility and autonomous-system studies of [7], [9], [11], the perception and compute co-design of [12], and the Formula 1 season study of [21].

14.1 Example 1: the drone (Fig. 48 of [1])↩︎

01_drone.py. The canonical MCDP from [1]: the “hello world” of monotone co-design, and the smallest complete use of the operator API (FunctionDP closed by loop; see 13.2.1.3 and 13.2.3.3).

Model↩︎

A battery powers an actuator that must lift the battery itself plus a fixed payload, against a fixed avionics power draw. The functionalities are the mission endurance (s), an extra_payload (kg) and an extra_power (W); the resource of interest is the battery mass (kg). The physics is a lumped model, \[m \;\geq\; \frac{\bigl(c\,(g(m+m_{\text{pay}}))^2 + P_{\text{ex}}\bigr)\,T}{\alpha},\] with \(\alpha = 1.8\times10^6\) J/kg the cell specific energy. The battery mass \(m\) appears on both sides: it is the loop axis. The inner FunctionDP emits the required mass under two names – battery_mass (fed back by loop, axis battery_mass) and report_mass (mirrored onto the outer \(R\) so the converged value survives on the interface rather than collapsing to a unit poset). All ports are on \(\mathbb{R}_{\geq 0}\) (see 13.1.1.2).

Code↩︎

def build_drone():
    F = Ports({
        "endurance": Reals(unit="s"),
        "extra_payload": Reals(unit="kg"),
        "extra_power": Reals(unit="W"),
        "battery_mass": Reals(unit="kg"),
    })
    R = Ports({
        "battery_mass": Reals(unit="kg"),   # loop axis
        "report_mass": Reals(unit="kg"),    # mirror, stays visible
    })
    inner = FunctionDP(F, R, drone_inner_h, name="drone-inner")
    return loop(inner, axis="battery_mass", name="drone")

Output↩︎

[language={},numbers=none,keywordstyle={},  extendedchars=true,literate={⊤}{{$\top$}}1]
Short, light: endurance=60s, payload=0.10kg, extra_P=1.0W
   iters=9, feasible=True, Antichain[(report_mass=0.0003564 kg)]

Medium, modest: endurance=300s, payload=0.50kg, extra_P=5.0W
   iters=22, feasible=True, Antichain[(report_mass=0.04921 kg)]

Longer mission: endurance=600s, payload=0.50kg, extra_P=5.0W
   iters=41, feasible=True, Antichain[(report_mass=0.1283 kg)]

Marginal: endurance=600s, payload=1.00kg, extra_P=10.0W
   iters=16, feasible=False, Antichain[(report_mass=⊤)]

Infeasible: endurance=1800s, payload=1.00kg, extra_P=10.0W
   iters=8, feasible=False, Antichain[(report_mass=⊤)]

Reading↩︎

The example sweeps five mission profiles, three feasible and two infeasible. The converged battery mass grows steeply with endurance – doubling the mission from 300 s to 600 s roughly triples the mass (0.049 kg to 0.128 kg), because heavier batteries demand more lift power, which demands still more battery, the feedback the loop resolves. Beyond the feasibility boundary the fixed point is \(\top\) (printed as the top element): no finite mass closes the loop. Iteration counts (8–41) track distance from that boundary: cases that diverge to \(\top\) or settle near the origin resolve fastest, while the near-feasible 600 s/0.5 kg case takes longest to climb. Examples 14.6 and 14.7 rebuild this exact model through the MCDP and System builders; the solver itself is solve (13.2.4.1).

14.2 Example 2: integer optimisation (Sec. VI-D)↩︎

02_integer_optimization.py. The classic Sec.VI-D problem from the paper, and the cleanest illustration of multi-point antichains ascending under the Kleene iteration. It exercises Naturals (13.1.1.3), a multi-valued FunctionDP (13.2.1.3) and the antichain algebra (13.1.2.1).

Model↩︎

Find the minimal pair \((x,y)\in\mathbb{N}\times\mathbb{N}\) with \[x + y \;\geq\; \lceil\sqrt{x}\rceil + \lceil\sqrt{y}\rceil + c .\] Both loop variables are bundled into a single composite axis xy (a nested Ports of two Naturals), so one Loop closes the pair. For each candidate the inner relation enumerates every split of the deficit \(c\) across the two coordinates – \(x_{\text{out}}\) ranging over \([\lceil\sqrt{x}\rceil,\,\text{target}-\lceil\sqrt{y}\rceil]\) – so the antichain gains several incomparable points at once. Infinities are mapped to the top pair, giving the iteration a well-defined \(\top\) to rest at when no finite pair works.

Code↩︎

target = sx + sy + c          # sx,sy = ceil(sqrt(x_in)),ceil(sqrt(y_in))
pts = []
for x_out in range(sx, target - sy + 1):
    y_out = target - x_out
    if y_out < sy:
        break
    pts.append({"xy": {"x": x_out, "y": y_out},
                "xy_report": {"x": x_out, "y": y_out}})
return Antichain.from_set(R, pts)
# ... inner wrapped and closed:
inner = FunctionDP(F=F, R=R, h_fn=h, name=f"sqrt_sum(c={c_value})")
return Loop(inner, axis="xy")

Output↩︎

c = 1: iters = 6, feasible = True
   S_0: { (0, 0) }
   S_1: { (0, 1), (1, 0) }
   S_2: { (0, 2), (1, 1), (2, 0) }
   S_3: { (0, 3), (1, 2), (2, 1), (3, 0) }
   S_4: { (0, 3), (2, 2), (3, 0) }
   S_5: { (0, 3), (3, 0) }
   S_6: { (0, 3), (3, 0) }
   M(c=1) = { (0, 3), (3, 0) }

c = 4: iters = 6, feasible = True
   ...
   M(c=4) = { (0, 7), (3, 6), (4, 4), (6, 3), (7, 0) }

Reading↩︎

The trace makes the ascent concrete: the seed \(S_0=\{(0,0)\}\) grows, peaks (four points at \(S_3\) for \(c=1\)), then contracts as the Min operation drops points that become dominated, settling at the antichain \(M(c)\) of minimal solutions. For \(c=4\) it converges in six steps to a five-point front. The example deliberately flags a typo in the paper: the paper’s \(M(1)=\{(1,0),(0,1)\}\) is wrong; the solver’s \(\{(0,3),(3,0)\}\) is correct, since \((1,0)\) gives \(1 \not\geq \lceil\sqrt1\rceil + 0 + 1 = 2\). That the number of points is non-monotone across iterations is the whole point of working in antichains rather than scalars. Example 14.5 plots these very iterates; the same make_looped builder is imported there.

14.3 Example 3: AUV seabed surveying (Sec. VIII)↩︎

03_auv_seabed.py. Example 10 of the paper: a genuinely cyclic co-design where velocity trades against sensor width, both feeding back through mission time. It shows a Loop (13.2.3.3) over a FunctionDP that emits a small Pareto spread, then collapses it with minimize_cost (13.2.4.3).

Model↩︎

An autonomous underwater vehicle sweeps an area \(A\) (m\(^2\)) at velocity \(v\) (m/s) with sensor field of view \(r\) (m). The couplings are intrinsically cyclic: a larger \(v\) shortens the time \(T = kA/(vr)\) but raises actuation power \(P_{\text{act}}=\psi\,v^3\); a wider \(r\) also shortens \(T\) but raises sensing power \(P_{\text{sens}}=\chi\,r\) and cost. Energy is \(E=(P_{\text{act}}+ P_{\text{sens}})\,T\). The design pair \((v,r)\) is bundled into the loop axis design; the outer resources \(T\), \(E\), cost stay visible. At each step the inner relation enumerates a few \((v,r)\) candidates above the current loop value (up to the physical caps \(v\le3\), \(r\le5\)), producing the front; hitting a cap lifts the axis to \(\top\) to signal infeasibility.

Code↩︎

pts = []
for v_try in (v, min(v*1.3, V_MAX), min(v*1.7, V_MAX)):
    for r_try in (r, min(r*1.3, R_MAX), min(r*1.7, R_MAX)):
        if v_try < v_in or r_try < r_in:
            continue
        T_try = K_GEOM * A / (v_try * r_try)
        E_try = (PSI_A*v_try**3 + CHI_A*r_try) * T_try
        pts.append({"design": {"v": v_try, "r": r_try},
                    "T": T_try, "E": E_try,
                    "cost": SENSOR_COST_A * r_try})
return Antichain.from_set(R, pts)

Output↩︎

Area = 100 m^2
   iters=2, feasible=True
   T=1176s, E=29.6kJ, $=100
   T=905s, E=29.5kJ, $=130
   T=692s, E=29.5kJ, $=170
   best by composite cost: T=692s, E=29.5kJ, $=170

Area = 1000 m^2
   iters=2, feasible=True
   T=11765s, E=295.9kJ, $=100
   ...
   best by composite cost: T=6920s, E=295.1kJ, $=170

Reading↩︎

Each scenario yields a three-point Pareto front trading survey time against sensor cost: the fastest sweep (\(T=692\) s) also costs the most ($170, the widest sensor), while the cheapest ($100) is slowest. Energy barely moves across the front – the drag and sensing terms nearly cancel over these \((v,r)\) choices – so it is time-versus-cost that drives the decision. Time scales linearly with \(A\) (100 m\(^2\to\)​10 000 m\(^2\) scales \(T\) by exactly \(100\times\)), confirming the coverage relation. The composite scalar cost (\(1/\text{s}+0.05/\text{kJ}+1/\$\)) is dominated by mission time, so minimize_cost always selects the fastest design. Contrast with the multi-objective handling in Example 14.8, where the front is genuinely two-dimensional.

14.4 Example 4: UncertainDP and ODE_DP↩︎

04_uncertain_and_ode.py. Two short demos of the advanced primitives: UncertainDP (13.2.1.8) for set-bracketed uncertainty, and ODE_DP (13.2.1.7) for deriving a monotone relation from a differential equation. Both are built on AlgebraicDP (13.2.1.2) internally.

Model↩︎

Uncertainty. A battery’s specific energy is known only to lie in \([1.6,2.0]\) MJ/kg. Two AlgebraicDPs bracket it: \(m=\text{cap}/1.6\text{e}6\) (pessimistic, more mass) and \(m=\text{cap}/2.0\text{e}6\) (optimistic). UncertainDP pairs them as lower/upper; with_mode selects which bound to solve. ODE. A heater holds a temperature rise \(\Delta T\) above ambient; Newton cooling gives \(\dot{T}=(P_{\text{in}}-h\Delta T)/C\), so at steady state \(P_{\text{in}}=h\,\Delta T\). ODE_DP integrates the right-hand side \(\dot{x} = h\,\Delta T - x\) to its steady state and extracts the root as the required power – a monotone \(R\) relation recovered from dynamics rather than written in closed form.

Code↩︎

uncertain = UncertainDP(F=F, R=R, lower=optimistic,
                        upper=pessimistic, mode="upper",
                        name="battery_uncertain")
for mode in ("lower", "upper"):
    result = solve(uncertain.with_mode(mode), {"capacity": 3.6e6})

heater = ODE_DP(
    F=F, R=R,
    rhs=lambda x, t, f: H_LOSS * f["delta_T"] - x,
    extract=lambda x: {"power": float(x)},
    mode="steady_state", x0_fn=lambda f: 0.0, name="heater_ode")

Output↩︎

UncertainDP demo: battery sizing under specific-energy uncertainty

   optimistic   (lower): mass = 1.800 kg
   pessimistic  (upper): mass = 2.250 kg

ODE_DP demo: power required to hold a steady temperature rise

   delta_T =    5 K  -->  P_in =   4.0 W   (=  h_loss * delta_T = 4.0 W)
   delta_T =   20 K  -->  P_in =  16.0 W   (=  h_loss * delta_T = 16.0 W)
   delta_T =   50 K  -->  P_in =  40.0 W   (=  h_loss * delta_T = 40.0 W)

Reading↩︎

For 1 kWh (\(3.6\) MJ) the optimistic cell needs 1.80 kg, the pessimistic 2.25 kg; the true mass lies between, and any design that survives the 2.25 kg (upper) case is robust to the chemistry actually delivered. This is the set-based philosophy of Sec. VII made operational: solve both brackets, commit to the pessimistic one. The heater demo is a sanity check – the numerically integrated steady state reproduces the closed form \(P_{\text{in}}=h\,\Delta T\) exactly (\(h=0.8\) W/K), confirming that ODE_DP recovers the algebraic relation a modeller would otherwise hand-derive. See 8 for the theory these two primitives implement; Examples 14.11 and 14.12 push the uncertainty machinery further.

14.5 Example 5: visualising the Kleene ascent↩︎

05_visualize_kleene.py. A visual companion to Example 14.2: it renders each iterate \(S_0, S_1, \ldots\) of the Sec.VI-D problem as a scatter panel, reproducing the structure of Fig. 36 of the paper. It shows how the trace recorded by solve (13.2.4.1, via record_trace=True) feeds a plotting routine.

Model↩︎

The MCDP is unchanged from Example 14.2 – the script imports make_looped from 02_integer_optimization rather than redefining it (the module name starts with a digit, so importlib is used). The novelty is entirely in consuming result.trace: a list of TraceEntry records (13.2.4.1), one per Kleene iterate, each carrying the antichain at that step. Each panel plots the finite points of one \(S_k\) in \(\mathbb{N}\times\mathbb{N}\); \(\top\) points are skipped.

Code↩︎

result = solve(looped, {"c": c_value}, max_iter=50, record_trace=True)
trace = result.trace
for k, entry in enumerate(trace):
    A = entry.antichain
    xs = [p["xy"]["x"] for p in A.points if p["xy"]["x"] != inf]
    ys = [p["xy"]["y"] for p in A.points if p["xy"]["y"] != inf]
    axes[k].scatter(xs, ys, s=60, c="C3", zorder=3)
    axes[k].set_title(f"$S_{{{k}}}$  ({len(A.points)} pts)")
fig.savefig(out_path, dpi=140, bbox_inches="tight")

Output↩︎

The script writes one PNG per value of \(c\in\{1,4,8\}\) (matplotlib in the Agg backend, no display required). Console output is only the progress log:

Plotting Kleene iteration traces into .../outputs ...
   wrote .../outputs/kleene_trace_c1.png
   wrote .../outputs/kleene_trace_c4.png
   wrote .../outputs/kleene_trace_c8.png

Committed copies of the figures live under docs/images/ (Figure 1).

Reading↩︎

The panels make the antichain algebra tangible: points appear along a frontier \(x+y=\text{const}\), the frontier advances outward, and then dominated points fall away under Min until only the incomparable minimal pairs remain. The figure is the exact geometric content of the Sec.VI-D trace printed textually in Example 14.2; reading the two together is the fastest way to build intuition for what a Kleene fixed point over antichains actually looks like. This is the only example in the first eight that touches matplotlib and therefore has a non-trivial import cost (\(\approx 1.2\) s wall-clock, essentially all of it library import).

14.6 Example 6: the drone in MCDPL syntax↩︎

06_drone_mcdpl_syntax.py. The Example 14.1 drone rebuilt through the declarative MCDP builder (13.3.1.1), whose notation mirrors the paper’s mcdp{...} blocks. It demonstrates the highest-level authoring layer: provides/requires/ constraint/loop_on inside a context manager.

Model↩︎

The MCDP is identical to Example 14.1: same physics, same feedback on battery_mass. Only the spelling changes. Inside a with MCDP("drone") as m: block, functionalities are declared with m.provides, resources with m.requires, each resource is given its defining expression by m.constraint, and the recursion is closed by m.loop_on("battery_mass"). As before, report_mass mirrors the loop value onto the visible interface. m.build() compiles the block down to the same FunctionDP-in-a-Loop the operator API produces by hand.

Code↩︎

with MCDP("drone") as m:
    m.provides("endurance", unit="s")
    m.provides("extra_payload", unit="kg")
    m.provides("extra_power", unit="W")
    m.provides("battery_mass", unit="kg")     # loop variable
    m.requires("battery_mass", unit="kg")     # loop axis
    m.requires("report_mass", unit="kg")      # mirror
    m.constraint("battery_mass", battery_mass_eq)
    m.constraint("report_mass", battery_mass_eq)
    m.loop_on("battery_mass")
return m.build()

Output↩︎

[language={},numbers=none,keywordstyle={},  extendedchars=true,literate={⊤}{{$\top$}}1 {×}{{$\times$}}1]
DP(loop(drone, axis=battery_mass): endurance:R+[s]×extra_payload:R+[kg]
   ×extra_power:R+[W] -> A[report_mass:R+[kg]])

Short, light: endurance=60.0, extra_payload=0.1, extra_power=1.0
   iters=9, feasible=True, Antichain[(report_mass=0.0003564 kg)]

Medium, modest: endurance=300.0, extra_payload=0.5, extra_power=5.0
   iters=22, feasible=True, Antichain[(report_mass=0.04921 kg)]

Longer mission: endurance=600.0, extra_payload=0.5, extra_power=5.0
   iters=41, feasible=True, Antichain[(report_mass=0.1283 kg)]

Infeasible: endurance=1800.0, extra_payload=1.0, extra_power=10.0
   iters=8, feasible=False, Antichain[(report_mass=⊤)]

Reading↩︎

The numbers are identical to Example 14.1, to the last digit (0.0003564 kg, 0.04921 kg, 0.1283 kg) and identical iteration counts (9, 22, 41, 8). That is the point: the MCDP builder is a front-end that lowers to exactly the same internal DP as the operator API, so the choice between them is purely one of authoring ergonomics. The printed signature confirms the compiled object is loop(drone, axis=battery_mass) – the same Loop as before. Prefer this layer when a model is naturally a flat list of constraint equations; drop to FunctionDP/loop (13.2.1.3, 13.2.3.3) when a relation needs arbitrary Python. Example 14.7 shows the third route, the System/Module builder.

14.7 Example 7: the drone, modular↩︎

07_drone_modular.py. The Example 14.1 drone a third way: as a System (13.3.2.1) of Module subclasses wired by operator-overloaded >= constraints. This is the layer for genuinely modular models – each subsystem carries its own \(F\)/\(R\) and relation.

Model↩︎

Two Module subclasses declare class-level \(F\)/\(R\) and an h: Battery maps capacity\(\to\)mass via specific energy, Actuator maps lift_force\(\to\)power as \(c\,\text{lift}^2\). The System exposes the outer ports, adds the two modules, and states three inequalities that read like the textbook constraints: battery capacity \(\ge\) power\(\times\)endurance, actuator lift \(\ge g(m_{\text{batt}}+m_{\text{pay}})\), and total mass \(\ge m_{\text{batt}}+ m_{\text{pay}}\). Each module.port is a handle; >= records a constraint rather than evaluating a boolean.

Code↩︎

sys = System("drone")
endurance     = sys.provides("endurance", unit="s")
extra_payload = sys.provides("extra_payload", unit="kg")
extra_power   = sys.provides("extra_power", unit="W")
total_mass    = sys.requires("total_mass", unit="kg")
battery  = sys.add("battery",  Battery())
actuator = sys.add("actuator", Actuator())
battery.capacity    >= (actuator.power + extra_power) * endurance
actuator.lift_force >= G * (battery.mass + extra_payload)
total_mass          >= battery.mass + extra_payload
return sys        # sys.build() compiles to a looped DP

Output↩︎

[language={},numbers=none,keywordstyle={},  extendedchars=true,literate={⊤}{{$\top$}}1]
System('drone'):
  provides: endurance, extra_payload, extra_power
  requires: total_mass
  subsystems:
    battery: (capacity) -> (mass)
    actuator: (lift_force) -> (power)
  constraints:
    battery.capacity >= ((actuator.power + extra_power) * endurance)
    actuator.lift_force >= (9.81 * (battery.mass + extra_payload))
    total_mass >= (battery.mass + extra_payload)

Short, light: endurance=60.0, ...
   iters=17, feasible=True, Antichain[(total_mass=0.1004 kg)]
Medium, modest: ...
   iters=43, feasible=True, Antichain[(total_mass=0.5492 kg)]
Longer mission: ...
   iters=81, feasible=True, Antichain[(total_mass=0.6283 kg)]
Infeasible: ...
   iters=16, feasible=False, Antichain[(total_mass=⊤)]

Reading↩︎

The System reports total_mass (battery + payload), so the numbers differ from Examples 14.1 and 14.6 – which reported only the battery – but the underlying battery masses agree exactly: e.g.the medium case here is \(0.5492 = 0.04921 + 0.5\) kg. Iteration counts are roughly \(2\times\) those of the single-axis model (43 vs 22, 81 vs 41) because the compiled loop now closes two coupled internal axes (battery mass and actuator power) rather than one, so the fixed point takes more sweeps to propagate. The declarative >= form and the lambda form sys.constrain(target, fn) compile to the same constraint list. This is the recommended layer for multi-subsystem designs; Example 14.8 scales it to three modules including a catalog.

14.8 Example 8: a modular vehicle with Pareto tradeoffs↩︎

08_vehicle_modular.py. A System (13.3.2.1) of three subsystems where one is a discrete CatalogDP (13.2.1.5). Because the catalog holds Pareto-incomparable motors, the system result is a genuine multi-point front over total mass and total cost, collapsed by minimize_cost (13.2.4.3).

Model↩︎

Three modules are wired: a CatalogDP motor with seven entries (each a CatalogEntry giving torque at a (mass, cost); some share torque but trade mass against cost, so they are incomparable), a linear Chassis (mass and cost \(\propto\) load) and a Battery (sized by energy). Five >= constraints close the loop: chassis load \(\ge\) payload \(+\) motor \(+\) battery mass, motor torque \(\ge\) a load-dependent demand, battery energy \(\ge\) mission energy, and two aggregations for total_mass and total_cost. The mutual dependence of chassis mass and motor choice makes this cyclic.

Code↩︎

motor   = sys.add("motor",   make_motor_catalog())  # 7-entry CatalogDP
chassis = sys.add("chassis", Chassis())
battery = sys.add("battery", Battery())
chassis.load   >= payload + motor.mass + battery.mass
motor.torque   >= TORQUE_PER_KG * G * (payload + chassis.mass
                                       + battery.mass)
battery.energy >= mission_energy
total_mass >= payload + motor.mass + chassis.mass + battery.mass
total_cost >= motor.cost + chassis.cost + battery.cost

Output↩︎

Small parcel: payload=2.0, mission_energy=200000.0
   iters=4, feasible=True
   Pareto front (2 points):
      total_mass=  4.82 kg,  total_cost=$ 413.00
      total_mass=  5.78 kg,  total_cost=$ 165.00
   cheapest: total_mass=5.78 kg, total_cost=$165.00

Medium load: payload=10.0, mission_energy=1000000.0
   iters=3, feasible=True
   Pareto front (2 points):
      total_mass= 22.49 kg,  total_cost=$ 475.00
      total_mass= 20.41 kg,  total_cost=$ 969.00
   cheapest: total_mass=22.49 kg, total_cost=$475.00

Heavy + long: payload=20.0, mission_energy=5000000.0
   iters=2, feasible=False

Reading↩︎

Two of the three missions produce genuine two-point fronts: a lighter but pricier motor versus a heavier but cheaper one. For the small parcel the choice is stark – the light design is only 0.96 kg lighter yet $248 dearer – so minimize_cost on cost picks the 5.78 kg/$165 option. The heavy-plus-long mission is infeasible and the loop detects it in just two iterations: no motor in the catalog can supply the torque once the chassis has grown to carry a 5 MJ battery, so the fixed point is \(\top\). This is the first-eight showcase of a discrete catalog driving a system-level antichain; the microgrid flagship (Example 14.13) and the full-vehicle study (Example 14.17) build on the same CatalogDP pattern at much larger scale.

14.9 Example 9: a robotic arm with non-trivial topology↩︎

09_robotic_arm.py. Five subsystems wired with eight constraints, showing where the operator-overloaded builder (12) earns its keep: the connection graph is neither series nor parallel nor a single clean loop, and the cyclic structure emerges from the constraint lines rather than from any explicit loop machinery.

Model↩︎

The arm carries two revolute Joint modules (shoulder, elbow), a Sensor, a Controller, and a Battery. Each joint maps functionality \((\text{torque}, \text{speed})\) to resources \((\text{mass}, \text{electric\_power})\) with \(\text{mass} = \rho\,\tau\) and \(\text{electric\_power} = \tau\,\omega / \eta\). The topology is genuinely coupled in three places. The controller’s power draw grows with both the sensor sample rate and the joint command rate; the shoulder must lift the payload and the elbow joint’s own mass at the end of its longer arm, so shoulder.torque reads elbow.mass and closes a mechanical cycle through the elbow; and the battery is sized by the mission-integrated power of every electrical subsystem. The single system-level resource is total_mass over \(\mathbb{R}_{\geq 0}\) (kg).

Code↩︎

The constraint block is the whole point: each line is one physical relationship, and elbow.mass on the right of the shoulder torque constraint is what makes the Kleene iteration iterate.

elbow.torque >= G * payload_mass * elbow_arm_length
# shoulder lifts payload + the elbow joint's own mass: a cycle
shoulder.torque >= G * (payload_mass + elbow.mass) * shoulder_arm_length
sensor.sample_rate >= 2.0 * control_rate           # Nyquist
controller.input_rate   >= 2.0 * control_rate
controller.command_rate >= control_rate
battery.energy >= operating_time * (
    elbow.electric_power + shoulder.electric_power
    + controller.power + sensor.power)
total_mass >= (payload_mass + elbow.mass + shoulder.mass
    + sensor.mass + controller.mass + battery.mass)

Output↩︎

Printing the assembled System echoes the nine constraints back in normalised form, then three mission profiles solve to a single minimal mass each. Every profile converges in four Kleene steps.

Pick-and-place light:
   iters=4, feasible=True
   total_mass = 1.97 kg
Heavier payload:
   iters=4, feasible=True
   total_mass = 7.30 kg
Long-reach precise:
   iters=4, feasible=True
   total_mass = 7.24 kg

No figures are produced.

Reading↩︎

The four-step convergence for every profile shows the elbow-mass cycle is shallow: one pass to size the elbow, one to feed its mass into the shoulder, and the fixed point settles. The heavier-payload and long-reach cases land at almost the same total mass by different routes — one is payload-dominated, the other actuator- and battery-dominated through the higher control rate and longer operating time. This is the cleanest small demonstration of the constraint-graph style; contrast it with the modular vehicle of example 14.8 and the flagship microgrid cycle of example 14.13, both of which use the same System builder (12) at larger scale.

14.10 Example 10: watching the solver work↩︎

10_solver_trace.py. The two-module drone of example 14.7, instrumented to exhibit every observability feature of solve in one script (see the solver observability section, 7.4, and the reference card 13.2.4.1): silent mode, the end-of-solve summary, the per-iteration feed, the structured trace, a custom callback, and the status field.

Model↩︎

A battery (\(\text{mass} = \text{capacity}/1.8\,\text{MJ/kg}\)) and an actuator (\(\text{power} = 10\,F^2\)) are coupled by the usual drone cycle: capacity depends on actuator power and endurance, lift force depends on battery mass. The model is deliberately trivial so that the focus stays on the six calls the __main__ block makes, not on the physics.

Code↩︎

The instrumentation is entirely in the keyword arguments to solve: verbose controls printing, trace collects a list of TraceEntry records (13.2.4.5), and on_iteration receives each entry as it is produced.

r = solve(drone, f, trace=True, max_iter=100)
print(f"collected {len(r.trace)} trace entries")
print(f"first 5 deltas: {[e.delta for e in r.trace[:5]]}")

def my_logger(entry):
    if entry.iteration         print(f"iter {entry.iteration}, |A|={entry.n_points}, "
              f"delta={entry.delta}")

solve(drone, f, on_iteration=my_logger, max_iter=100)

Output↩︎

The nominal solve converges in 43 iterations to the canonical 0.5492 kg. Capping max_iter short and posing an under-resourced mission drives status through its three distinct values — these are the abbreviated forms of the final beat:

# 2. Final summary (verbose=1)
[solve] converged: 43 iters, |A|=1, total=0.3ms, feasible=True

# 6. Status field (max_iter vs converged vs feasible)
   max_iter=3:   status='max_iter', feasible=True,  iters=3
   max_iter=200: status='converged', feasible=True, iters=43
   infeasible:   status='diverged', feasible=False, iters=16

No figures are produced by the script; the companion notebook (10_solver_trace.ipynb) adds a delta-versus-iteration plot showing the characteristic oscillating decay of the Kleene residual to machine precision.

Reading↩︎

The decisive line is the last block. status separates two conditions that the boolean feasible flag conflates: hitting the iteration cap on a solvable problem (’max_iter’, where raising the cap would help) versus a genuinely infeasible mission (’diverged’, where it would not). The trace’s deltas confirm the Kleene iteration is monotone in the fixed-point sense but not monotone step-to-step — the residual zig-zags downward, which is why a single small max_iter can report a non-converged but already-feasible incumbent.

14.11 Example 11: set-based deterministic uncertainty↩︎

11_uncertain_drone.py. The drone of example 14.7 with its battery parametrised by two internal quantities — specific energy and efficiency — that are known only up to an uncertainty set. The question posed is: under the worst-case point of that set, how heavy is the drone? This exercises the set-based half of the uncertainty layer (8), the Box (13.4.1.2) and Ellipsoid (13.4.1.3) brackets.

Model↩︎

The battery mass is \(\text{capacity} / (\text{specific\_energy} \cdot \text{efficiency})\). The nominal parameters \((2.0\,\text{MJ/kg}, 0.90)\) are chosen so their product, 1.8 MJ/kg of delivered energy density, reduces to the canonical drone of examples 14.1, 14.6, and 14.7; the nominal solve converges to the same 0.5492 kg, and both uncertainty sets perturb around that shared reference. Each parameter is declared “more is better,” so worst-case badness is the low-parameter direction. The Box admits every combination in the two independent ranges; the Ellipsoid is a tilted, correlated set that excludes the joint-extreme corner.

Code↩︎

The uncertainty is attached to the module, not the solve; the same system is rebuilt with a different uncertain_set per query, and solve is asked for the worst_case summary.

bat = Battery()
bat.uncertain_set = Box(
    specific_energy=(1.7e6, 2.3e6, "more_is_better"),
    efficiency=(0.83, 0.97, "more_is_better"),
)
drone = make_drone(bat)
r_box = solve(drone, f, uncertainty=["worst_case"])
wc = list(r_box.worst_case.antichain.points)[0]["total_mass"]

Output↩︎

Nominal parameters:  total_mass = 0.5492 kg
Box uncertainty:
   worst-case total_mass = 0.5668 kg
   uncertainty penalty   = +0.0176 kg
Ellipsoid uncertainty (smaller, correlated set):
   worst-case total_mass = 0.5527 kg
   uncertainty penalty   = +0.0035 kg

No figures are produced.

Reading↩︎

The two penalties, \(+0.0176\) kg for the box against \(+0.0035\) kg for the ellipsoid, are the whole lesson. The box worst case sits at the single corner where both parameters are simultaneously at their lowest; the ellipsoid, being correlated and smaller, carves that implausible corner out and returns a far milder worst case. Choosing the uncertainty model is therefore a modelling decision about which joint outcomes are physically credible, not a mere tolerance bound. Example 14.12 takes the same drone stochastically and shows where the worst case sits relative to the mean, p95, and CVaR of a Monte Carlo distribution.

14.12 Example 12: stochastic uncertainty with a Gaussian copula↩︎

12_stochastic_drone.py. The same drone battery as example 14.11, but now the two internal parameters carry marginal distributions tied by a positive correlation — more energy-dense cells tend also to be more efficient. This is the stochastic half of the uncertainty layer (8): a Stochastic distribution (13.4.1.8) whose marginals are glued by a GaussianCopula (13.4.1.7).

Model↩︎

Both parameters take a uniform marginal over the same ranges as the box in example 14.11, coupled by a Gaussian copula at correlation \(0.4\). The nominal product is again 1.8 MJ/kg delivered, so the nominal solve is the canonical 0.5492 kg. A single solve call requests all five summaries at once — worst_case (from the box still attached), mean, p95, cvar95, and the raw samples — over one thousand Monte Carlo draws at a fixed seed.

Code↩︎

bat.uncertain_dist = Stochastic(
    marginals={
        "specific_energy": stats.uniform(loc=1.7e6, scale=0.6e6),
        "efficiency":      stats.uniform(loc=0.83, scale=0.14),
    },
    copula=GaussianCopula(correlation=[[1.0, 0.4],
                                       [0.4, 1.0]]),
)
res = solve(drone, f,
    uncertainty=["worst_case", "mean", "p95", "cvar95", "samples"],
    n_samples=1000, rng_seed=42)

Output↩︎

The five summaries come out in the canonical ordering \(\text{nominal} < \text{mean} < \text{p95} < \text{cvar95} < \text{worst\_case}\), with a feasibility rate of 1.000.

Nominal mass: 0.5492 kg
   mean total_mass   = 0.5506 kg
   p95 total_mass    = 0.5632 kg
   cvar95 total_mass = 0.5647 kg
   worst-case mass   = 0.5668 kg   (Box, deterministic)
Raw samples: 1000 antichains (1000 feasible)

The script also prints an ASCII histogram of the sampled masses; the companion notebook renders it as a matplotlib histogram with vertical markers at each summary. No image files are written by the script.

Reading↩︎

The ordering is the reason to compute all five together. The mean sits barely above nominal because the marginals are centred near the nominal point; p95 and CVaR95 climb into the tail, and the set-based worst case of example 14.11 bounds them all from above. The design engineer reads this as a spectrum of conservatism: size to the mean and half the fleet misses the target, size to CVaR95 and only the worst 5% of parameter draws exceed the mass budget, size to the box worst case and none do — at a heavier vehicle. The copula correlation matters because it concentrates draws along the “both good / both bad” diagonal, shifting the tail summaries relative to an independent model.

14.13 Example 13: the microgrid flagship↩︎

13_microgrid.py. An off-grid cabin must supply a daily energy demand and a peak load without grid power. The example is the flagship case study, exercising catalogue choice, a genuine fixed-point cycle, warm-started sweeps, stochastic uncertainty, and the visualisation suite in one script — the whole System builder (12) and solver (13.2.4.1) stack end to end.

Model↩︎

Four subsystems contribute over \(\mathbb{R}_{\geq 0}\): a solar PV array (cheap but sun-limited), a lithium battery whose chemistry is a discrete parameter (LFP, NMC, LCO, NaIon, each a \((\text{Wh/kg}, \text{USD/kWh}, \text{cycle life})\) triple), a diesel generator (reliable, carbon-heavy), and a mounting frame. The frame is the source of the cycle: its cost and mass scale with the total supported mass, frame included, so the constraint

\[\text{frame.supported\_mass} \;\ge\; \text{solar.mass} + \text{battery.mass} + \text{diesel.mass} + \text{frame.mass}\]

reads the frame’s own output on its right-hand side and forces the Kleene iteration to iterate to a fixed point rather than resolve in one pass.

Code↩︎

The sweep contrasts a cold restart against warm-starting each solve from the previous load’s fixed point via start_from.

warm_iters = 0
prev = None
for L in loads:                      # 50 points, 5..30 kWh
    f = {"daily_load_kwh": float(L),
         "peak_load_kw": 3.0, "backup_hours": 12.0}
    r = solve(dp, f, max_iter=400, start_from=prev)
    warm_iters += r.iterations
    prev = r                         # reuse the fixed point next step

Output↩︎

At the test mission (15 kWh/day, 3 kW peak) all four chemistries solve in 24 Kleene iterations; the warm-started sweep saves roughly 10 % of total steps.

LFP     cost=$   11689  mass= 388.6kg  CO2= 416.1kg/yr  (iters=24)
   NMC     cost=$   14628  mass= 356.2kg  CO2= 416.1kg/yr  (iters=24)
   NaIon   cost=$   11328  mass= 432.8kg  CO2= 416.1kg/yr  (iters=24)

   cold-start total iters: 1205
   warm-start total iters: 1083
   speedup: 1.11x

   nominal cost:  $11689   MC mean cost: $11930
   p95 cost:      $15162   CVaR95 cost:  $16127

The script writes three files to outputs/: microgrid_convergence.png (the Kleene residual trace), microgrid_uncertainty.png (the Monte Carlo cost distribution, x-axis labelled “total cost (USD)”), and microgrid.dot (Graphviz source for the system diagram).

Reading↩︎

The chemistry table is a Pareto snapshot: NaIon is the cheapest but heaviest, NMC the lightest but dearer, LFP the balanced default, and CO\(_2\) is fixed because it is set by the diesel backup, not the chemistry. The 24-iteration count confirms the frame cycle is real — a acyclic model would resolve in one or two passes. The 1.11\(\times\) warm-start speedup is modest but free: consecutive loads have nearby fixed points, so seeding from the previous solution shaves off the early Kleene steps. The stochastic pass shows the design is robust — feasibility stays at 1.000 as sun hours vary — while the tail cost climbs by roughly 40 % from nominal to CVaR95, the premium for sizing against bad-sun years. Compare the discrete-catalogue mechanics here with the online elimination of example 14.14, which prunes a catalogue instead of enumerating it.

14.14 Example 14: online elimination over a robot catalogue↩︎

14_online_fleet.py. A logistics service must hit a target throughput and range by buying robots from a catalogue of 200 candidate types, each described by four features (speed, payload, unit cost, energy per km). Naively this is 200 inner solves; the online solver (15.6, 13.4.2.7) prunes provably suboptimal candidates using an optimistic evaluator and runs inner solves only for the survivors. The example reproduces the spirit of the multi-robot fleet study of [18].

Model↩︎

Each robot’s inner DP is a smooth (fractional-fleet) AlgebraicDP mapping the mission \((\text{throughput}, \text{range})\) to \((\text{total\_cost}, \text{total\_energy})\), with \(\text{total\_cost} = (\text{throughput}/(sp))\,c\) and \(\text{total\_energy} = \text{range}\cdot e \cdot 24\). Three evaluators are compared: a LipschitzEvaluator (13.4.2.2 and neighbours) with per-output constants; a MonotonicityEvaluator on the derived feature \(\text{cost\_per\_capacity} = c/(sp)\), under which total_cost is exactly monotone; and the certified confidence-polytope LinearParametricEvaluator.

Code↩︎

The three flavours differ only in the evaluator object handed to solve_online; the inner-DP factory and mission are shared.

res = solve_online(make_dp, mission,
                   candidates=candidates,
                   evaluator=MonotonicityEvaluator(
                       features=["cost_per_capacity", "energy_per_km"],
                       r_components=["total_cost", "total_energy"]),
                   verbose=0)
print(res.n_evaluated, res.n_eliminated, len(res.antichain))

Output↩︎

The exhaustive baseline finds a 5-point Pareto front. All three evaluators recover every Pareto point; they differ only in how much of the catalogue they had to evaluate.

True Pareto front: 5 non-dominated candidates

Lipschitz (L=300/30):        evaluated 192/200  Pareto recovery: OK
Monotonicity (cost_per_cap): evaluated  13/200  Pareto recovery: OK
LinearParametric (certified):evaluated 122/200  Pareto recovery: OK

The script writes fleet_online.png (three feature-plane panels shading each candidate by status) and fleet_online_convergence.png (incumbent antichain size against inner solves performed) to outputs/.

Reading↩︎

The evaluated counts are the story. Monotonicity with the right derived feature prunes hardest — 13 of 200 — because a single monotone coordinate orders the whole catalogue. Lipschitz is the most general and most conservative, evaluating 192 of 200: with a safe \(L\) it never drops a Pareto point but rarely eliminates one either. The certified LinearParametric evaluator sits between, at 122 of 200, and recovers all five: it maintains the polytope of linear parameter vectors consistent with every observation and lower-bounds each query by one LP over that polytope, so the bound is guaranteed and never wrongly eliminates an optimal candidate. It prunes less aggressively here because total_cost is not affine in the raw features; on a truly linear map it becomes exact and aggressive. This certified bound replaces an earlier OLS \(\pm\) confidence-band heuristic that pruned to about 35 of 200 but wrongly dropped one Pareto point on this seed — the trade is aggressiveness for a guarantee. Example 14.16 pushes the same three evaluators onto a deliberately nonlinear surface and shows the certified evaluator declining to guess.

14.15 Example 15: monoclonal antibody fed-batch co-design↩︎

15_bioprocess.py. A biopharmaceutical company must deliver a monoclonal antibody at a target titer (g/L) and annual demand (kg/year). Four subsystems are coupled cyclically, and two of them are catalogue lookups (13.2.1.5) resolved by the solver. This is the first biotech-upstream application in the library, calibrated to the 2024–2026 bioprocessing literature [22]–[28], and it uses the same System builder (12) as the engineering examples.

Model↩︎

A CellLine module maps demanded titer to the required peak viable-cell density and oxygen uptake via an effective integrated specific productivity \(q_P\). A FeedStrategy module turns the glucose set-point into a metabolic-burden factor with a U-shaped batch-failure penalty around 8 mM. Two CatalogDP lookups pick the smallest sufficient bioreactor (single-use 200 L through stainless 25,000 L) and the cheapest sufficient media (four commercial CD formulations). The cycle the Kleene iteration resolves runs through the peak cell density: higher titer needs higher density, which inflates the metabolic burden through the feed, which inflates the density needed for the same titer. The outer objective is the vector \((\text{COGS}, \text{footprint}, \text{CO}_2\text{ per gram})\).

Code↩︎

Cell line and feed are fixed at construction so the outer loop sweeps them; the bioreactor and media are chosen by the solver. The three outer resources are assembled from ports plus closed-over parameters via sys.constrain.

cell  = sys.add("cell",  CellLine(*cell_line))
feed  = sys.add("feed",  FeedStrategy(glucose_setpoint_mm))
bior  = sys.add("bior",  make_bioreactor_dp())   # CatalogDP
media = sys.add("media", make_media_dp())        # CatalogDP
# metabolically-inflated density must be supported by both catalogues
bior.peak_vcd  >= cell.peak_vcd * feed.metabolic_factor
media.peak_vcd >= cell.peak_vcd * feed.metabolic_factor
sys.constrain("cogs_per_g", cogs_eq)   # ports + annual_demand, yield

Output↩︎

At each of the three mission scales the global Pareto front is a genuine two-point tradeoff between the high-producer workhorse CHO-K1 and the next-generation short-batch CHO-MK, both at the 8 mM feed optimum.

=== titer 5.0 g/L, annual demand 100 kg/year ===
  design            COGS      footprint      CO2
  CHO-K1/glu=8mM   $ 44.04/g   39.9 m^2    214.29 g/g
  CHO-MK/glu=8mM   $ 63.65/g   30.5 m^2    214.29 g/g
  Cheapest design overall: CHO-K1/glu=8mM

No figures are produced.

Reading↩︎

The two Pareto points are not comparable: CHO-K1 wins on cost of goods (low licence fee, long batch) while CHO-MK wins on footprint (high licence fee, short batch, so fewer parallel lines at the same annual output). CHO-K1 is the cheapest overall at every scale, and per-gram COGS, footprint intensity, and CO\(_2\) all fall as demand rises, the signature of amortising catalogue capex over more grams. The 8 mM feed wins at both cell lines because the U-shaped penalty punishes both the starvation and the waste-accumulation extremes. All parameters sit within published ranges (specific productivity from Reinhart et al., metabolic limits from Khattak et al. and Lao & Toth 1997); the values are illustrative, not production-precision. Example 14.16 fixes CHO-K1 and turns to the operating-condition optimisation this model leaves open.

14.16 Example 16: online DOE for the mAb fed-batch process↩︎

16_online_doe.py. The bioprocess of example 14.15 is fixed (CHO-K1, 100 kg/year, 5 g/L) and the question shifts to where to operate within a \(5\times5\times5\times3 = 375\)-point grid of operating conditions — temperature, pH, glucose target, feed-start day. Each candidate stands in for a 10–14-day bioreactor run costing $20,000–$100,000, so the online elimination of example 14.14 (15.6, 13.4.2.7) is now a Design-of-Experiments budget question.

Model↩︎

A closed-form effect model maps each condition vector to an effective productivity and thence to the outer triple \((\text{COGS}, \text{footprint}, \text{CO}_2)\), calibrated from the process-parameter literature (temperature downshift from Yoon et al. and Sou et al., pH from Trummer et al., glucose and feed from Khattak et al. and Gagnon et al.). The surface is deliberately nonlinear — a temperature/pH/glucose U-shape and a power-law batch-failure rate. Two non-online baselines (a 75-run factorial DOE at the pH = 7.1 slice, and a 40-run uniform random sample) are compared against the Lipschitz, Monotonicity, and certified LinearParametric evaluators at a budget of 40 inner solves, scored by how many of the four distinct Pareto classes each recovers.

Code↩︎

Recovery is counted by output value, not candidate key, because many grid points share the same \((\text{COGS}, \text{footprint})\) outcome.

true_classes = {(round(p["cogs_per_g"], 2),
                 round(p["footprint_m2"], 1)) for p in true_pareto}
for name, ev in evaluators:
    res = solve_online(make_dp, {"target_titer": 5.0},
                       candidates=candidates, evaluator=ev, budget=40)
    recovered = {(round(o["cogs_per_g"], 2),
                  round(o["footprint_m2"], 1))
                 for i in res.incumbent_ids for o in [simulate_run(...)]}
    print(name, len(true_classes & recovered), "/", len(true_classes))

Output↩︎

The exhaustive baseline (375 solves) yields four Pareto classes. At budget 40 the tuned Lipschitz evaluator matches the 75-run factorial DOE at recovering 3 of 4; the certified LinearParametric evaluator recovers zero, and correctly so.

strategy                         evals   elim   Pareto recovered
  Factorial DOE (pH=7.1 slice)        75     -    3 / 4 classes (75  Random sample (seed=42)             40     -    3 / 4 classes (75  Lipschitz (normalised features)     40     0    3 / 4 classes (75  Monotonicity (cogs only)            40     0    0 / 4 classes ( 0  LinearParametric                    40     0    0 / 4 classes ( 0
  Monotonicity + 4 corner warm-start runs: recovery 1/4 classes (25

The companion notebook visualises the pick trajectory in the (temperature, glucose) projection.

Reading↩︎

The zero from LinearParametric is the safe-degradation property, not a failure. On the linear catalogue of example 14.14 the same evaluator recovered every Pareto point; here the map is markedly nonlinear, no single affine parameter set fits the observations, the confidence polytope cannot support a non-trivial bound, and the evaluator falls back to the no-information value rather than risk eliminating a candidate it cannot certify as suboptimal. Monotonicity alone is likewise silent because its lower bounds tighten only for candidates above every observation, and the Pareto front lives at the low-feature corner; seeding four hand-picked corner runs lifts it from 0 to 1 class. The generalisable lesson is that the structural prior does the work and must match the problem: a certified global prior is exact on an affine map and deliberately silent off it, a Gaussian process buys reach on smooth nonlinear surfaces at the price of an uncertified bound, and pure local bounds are always safe but need denser observations or a warm-start to concentrate.

14.17 Example 17: full-vehicle co-design across ICE, hybrid, and EV↩︎

17_car_codesign.py. The largest example in the package (3351 lines) and the first to model three architecturally distinct variants in one design study. A passenger car is decomposed into per-architecture MCDP modules — 24 for ICE, 23 for hybrid, 18 for EV — each with its own F and R ports, assembled with the System builder (12) and the catalogue DP (13.2.1.5) over engines, transmissions, motors, and battery packs. Because the file is so large, only the EV mass/energy loop and the sweep skeleton are quoted here; the three build_*_car bodies follow the same template. Component masses, efficiencies, and costs are calibrated from the automotive engineering literature [29]–[39].

Model↩︎

ICE, hybrid, and EV share a common chassis core (body, front/rear suspension, front/rear brakes, steering, tires, wheels) and auxiliary core (HVAC, interior, safety, lighting/infotainment); they differ in the powertrain modules. Two coupled cycles close inside every architecture and are resolved by the Kleene iteration. The first is the mass spiral: every load-bearing module reads the design curb weight as an F input, and the macro aggregation sums every module weight back into curb_weight, so heavier vehicle drives heavier components drives heavier vehicle. The second is the energy-storage loop: fuel-tank size (ICE/HEV) or battery capacity (EV/HEV) is sized by consumption \(\times\) range, and the storage’s own mass re-enters the curb weight. For EVs the pack is 20–35 % of curb weight, so this loop dominates and convergence takes 30–50 iterations against 15–25 for ICE. The macro objective is the vector \((\text{cost}, \text{curb weight}, \text{CO}_2)\), with fuel or energy consumption, maintenance, and durability also reported.

Code↩︎

Inside build_ev_car, one weight expression is fed to every load-bearing module (the mass spiral) and the battery is sized from the mission range and a consumption lambda (the energy loop).

total_weight_expr = (body_m.weight + fr_m.weight + rr_m.weight
    + bf_m.weight + br_m.weight + st_m.weight + ti_m.weight
    + wh_m.weight + mot_m.weight + bat_m.weight + pe_m.weight
    + obc_m.weight + red_m.weight + btm_m.weight + DRIVER_MASS_KG)
for m in (fr_m, rr_m, bf_m, br_m, st_m, ti_m, wh_m, li_m):
    m.design_mass >= total_weight_expr          # mass spiral
sys.constrain("hv_battery.target_range_km",
              lambda x, r=mission["target_range_km"]: r)
sys.constrain("hv_battery.energy_kWh_per_100km",
              ev_energy_kWh_per_100km)           # closes energy loop

Each architecture has a sweep function enumerating its catalogue combinations; main runs all three across four missions and keeps the cheapest feasible design per pair.

def sweep_ev(mission, *, verbose=False):
    results = []
    for body in _eligible_bodies(mission):
        for motor in EV_MOTORS:
            for battery in EV_BATTERIES:
                dp = build_ev_car(mission=mission, body=body,
                                  motor=motor, battery=battery, ...)
                res = solve(dp, dict(mission), max_iter=200)
                if res.feasible and res.antichain.points:
                    results.append((label, dict(list(
                        res.antichain.points)[0])))
    return results

Output↩︎

The script prints a per-mission table and a 10-year-TCO summary. The trimmed summary below reproduces the headline claims. For the Family Daily mission (5 passengers, 700 km range, 9 s 0–100): the cheapest ICE is $35,523 at 196 g/km, the cheapest HEV $40,272 at 91 g/km with the best 10-year TCO at $56,427, and the cheapest EV $64,666 at 58 g/km but at 2434 kg curb weight. The 7-passenger 800 km EV (Suburban Utility) has zero feasible designs.

Mission             Arch    Cost   Weight     CO2   10y TCO
Family Daily        ICE   $35,523   1710kg  196g/km $64,998
Family Daily        HEV   $40,272   1719kg   91g/km $56,427
Family Daily        EV    $64,666   2434kg   58g/km $76,052
Suburban Utility    HEV   $42,510   1910kg  113g/km $61,066
Suburban Utility    EV        -        -       -       -   (0 feasible)

The full module graphs are shipped as docs/diagrams/car_{ice,hev,ev}.svg; the draw_system renderer shades the mass-spiral strongly-connected component in amber.

Reading↩︎

The three architectures trade cost against CO\(_2\) monotonically — ICE cheapest to build and dirtiest, EV dearest and cleanest, HEV in between — but the 10-year TCO reorders them: the HEV’s low fuel bill gives it the best total cost of ownership on the Family Daily mission, undercutting both the cheaper-to-buy ICE and the far dearer EV. The 2434 kg EV curb weight is the mass spiral and energy loop compounding: a 100 kWh 800V pack is needed for 700 km, and the pack’s own mass inflates every load-bearing component. The Suburban Utility EV infeasibility — 7 passengers over 800 km — is the model correctly reproducing the real-world absence of such a product against 2024 pack-level energy density (160–180 Wh/kg): the required pack mass never reaches a fixed point below the chassis limit. Calibration cites Genta & Morello, the Bosch handbook, Heywood, Hofmann, and the IEA EV Outlook among others; the numbers are “within published ranges for 2020–2024 production vehicles,” and the contribution is the compositional analysis, not the point values. This example is the capstone of the static-solve examples; the temporal, sequential, and online-co-design regimes begin with the following section (15).

15 Temporal, Vector-State, and Online Co-Design↩︎

The examples up to 17 solve a co-design problem once: a fixed functionality is posed and the Kleene iteration returns the minimal resource antichain. Many engineering systems are not static. The best architecture changes as the environment changes; a resource is spent or worn across a sequence of stages; the design must be re-decided online as measurements arrive. This section documents the six layers the package adds on top of the ordinary solve() to handle these regimes, and the theory that makes them exact. All six reuse a common decision object, the Architecture (a named wrapper around a System), so they compose with the whole rest of the framework.

Unlike the rest of the manual, this section does not document published theory. The static co-design machinery underneath it is Censi’s [1]; the temporal, sequential, vector-state, and receding-horizon layers built on top of it are the author’s own and are being prepared for separate publication. They have not been peer reviewed. The closest published work is the hierarchical seasonal co-design of [21], which was developed independently.

15.1 Case 1: architecture switching (temporal)↩︎

The simplest temporal problem is switching: a system whose best architecture changes because the environment changes, with a cost to switch. The environment is a sequence of epochs, each posing its own functionality and admitting a subset of candidate architectures. solve_schedule chooses one architecture per epoch to minimise the total of the per-epoch co-design costs plus the switching costs, by an exact Viterbi (shortest-path) pass over the epoch-by-architecture lattice. The switching cost may be a uniform scalar or a per-transition callable, and a hysteresis parameter lets a brief unfavourable epoch be ridden out on the incumbent rather than switched into and back out of. The per-epoch cost is a full solve(), so switching sits on top of static co-design rather than replacing it.

15.2 Case 2, scalar value: dynamic programming (dynamic)↩︎

When a resource is carried between stages, the problem becomes a dynamic program. solve_dynamic runs a finite-horizon backward Bellman recursion whose per-stage decision is which architecture to instantiate, whose per-stage cost is itself a co-design solve(), and whose carried scalar state (fuel, charge, cumulative wear) is advanced by a user transition. The value function is scalar: at each (stage, state) it is the single best cost-to-go. The state is discretised onto an explicit StateGrid; a transition that would leave the grid envelope is rejected before snapping, so an over-spent resource (a negative fuel level) is never silently rescued by the nearest in-range node. The solver returns a state-indexed DynamicPolicy, queryable in closed loop at off-nominal states, and rollout threads it forward from a concrete initial state.

15.3 Case 2, antichain value: sequential co-design (sequential)↩︎

The multi-objective generalisation replaces the scalar value with a Pareto antichain. Following the sequential co-design theory, fix an ordered commutative resource monoid \((R, \le, \oplus, 0)\) and work in the value space of upper sets of \(R\) ordered by reverse inclusion (a larger feasible set is a lower resource threshold, hence “easier”). A sequential co-design problem over stages \(0,\dots,N\) and a state poset \((X, \le_X, \bot_X)\) assigns to each stage \(k\) a state-parametrised minimal-resource map \(h_k : X \to \mathcal{U}R\) and a transition \(\phi_k : X \times R \to X\). The value is the backward recursion \(V_{N+1} \equiv R\) and \[V_k(x) \;=\; \mathrm{Min} \bigcup_{r \in h_k(x)} \Bigl( \uparrow\! r \;\oplus\; V_{k+1}\bigl(\phi_k(x, r)\bigr) \Bigr),\] the antichain-valued Bellman operator, implemented by Antichain.union_min. In the package the accumulated resource (the named cost axes on the antichain) is kept deliberately distinct from the carried state \(x\): the transition reads any carried quantity off the full solved point, while only the cost axes accumulate on the front. The scalar dynamic layer is the width-one special case (a single cost axis gives a singleton antichain).

Three results are made operational. (Q1) Monotone value. If, for every stage, (H1) \(h_k\) is monotone in the carried state (more carried state makes the stage no easier), (H2) \(\phi_k\) is jointly monotone on \(X\times R\) (monotone in the state and in the chosen resource point, so the successor built from a smaller state and a smaller resource stays comparable), and (H3) admissibility is a down-set of \(X\) (a state below an admissible one is admissible, so the out-of-bounds guard never deletes the successor the proof constructs), then every \(V_k\) is monotone. (H1) is load-bearing and necessary; the joint slice of (H2) and (H3) are the two obligations the informal statement previously omitted, and both hold for the canonical consumable-sum and renewable-join transitions. check_monotonicity verifies all three on the state grid; it is orientation-aware, accepting the benign “consumable but monotone” orientation (a budget carried as remaining slack) and flagging only genuinely non-monotone (perishable or fatigue-as-state) stages, a transition that is not jointly monotone, or a non-down-set admissible region, all of which fall outside the guarantee. (Q2) Front \(=\) reachable frontier. \(\mathrm{Min}\,V_k(x)\) equals the antichain of minimal reachable cumulative resources over feasible choice sequences; there is no tail-pruning gap, and the front size is the width of the reachable set, which grows polynomially in the horizon for a summed objective on a fixed number of axes (sum_combine) and stays bounded for a renewable one (join_combine). (Q3) Exact factorisation. At a reset stage (a quiescent transition landing on \(\bot_X\) regardless of input), detect_resets and factorise_at_resets split the horizon into a \(\oplus\)-product of independent sub-problems, the order-theoretic deterministic analogue of regeneration-point decomposition; the proof uses only distributivity, so it holds in every regime.

15.4 The general carried state: vectors (state, vector_dp)↩︎

Realistic problems carry structured state, not a lone scalar. The Formula 1 seasonal co-design of Neumann, Zardini, and colleagues [21] carries a state vector of two battery-wear levels plus a discrete regulatory-penalty flag; a self-reconfiguring robot carries the accumulated wear of each drive module plus a shared energy budget. The state module provides a VectorStateGrid: a product grid over named axes, each either a ContinuousAxis (a bucketed real interval, snapped and bounds-checked like the scalar grid) or a DiscreteAxis (a finite labelled set with a caller-supplied order, or mutually incomparable when no order is given). solve_vector_sequential runs the same antichain-valued Bellman recursion carrying the full state vector, with the transition returning an \(\{\text{axis} : \text{value}\}\) mapping snapped onto the product grid and out-of-bounds successors rejected on every axis. The monotonicity results carry over verbatim with the component-wise product order substituted for the scalar order, and check_vector_monotonicity verifies (H1), the joint (H2), and (H3) over that order. A single-axis grid reproduces the scalar sequential result exactly, so the vector layer strictly generalises rather than parallels it.

A subtle correctness point governs the backward pass. The stage antichain must not be Pareto-reduced on the cost axes before the transition is applied: a point dominated on cost may be the only feasible choice from a constrained carried state (a higher-cost morphology that spares a worn module), and its consequence for the carried state is invisible in the cost projection. The implementation therefore enumerates all solved points and lets union_min prune only after the transition. For the same reason the pre-transition deduplication keys candidate points on their full resource point (all ports), not on the cost projection: two points that share a cost but differ on a carried axis (equal cost, different fuel) are incomparable in the full resource poset and must both survive, so the \(\mathrm{Min}\) lives in \(R\) rather than its projection. This is what makes the reconfigurable-robot example (21) feasible.

15.5 Precompute-then-DP: the Formula 1 structure↩︎

There are two distinct ways to combine co-design with dynamic programming, and they should not be conflated. The sequential solver above re-solves the co-design at every (stage, state), which is necessary when the per-stage functionality depends on the carried state. The Formula 1 seasonal framework [21] uses the cheaper opposite order: the co-design layer is run once to produce track-dependent Pareto-optimal performance-versus-wear mappings, which are then frozen into a catalog of implementations that an outer finite-horizon dynamic program selects among. No co-design solve happens inside the Bellman sweep; the DP only indexes the precomputed catalog. The package exposes this explicitly: precompute_catalog solves the architectures once at a fixed functionality and returns the tagged Pareto catalog (retaining points that look dominated at the single-stage level, since they may become optimal once aggregated in the DP, exactly the F1 paper’s observation), and dp_over_catalog runs the antichain-valued DP selecting catalog indices. Precompute-then-DP agrees exactly with solve_sequential when the co-design is state-independent, and is preferred there for cost; when the co-design genuinely depends on the carried state, only the re-solving sequential form is correct.

15.6 Online feedback co-design (online_codesign)↩︎

The layers above are open-loop planners: the whole horizon is solved in advance. The closed-loop counterpart re-decides the design online as measurements arrive. run_online_codesign implements a receding-horizon, measurement-in-the-loop controller: at each step it (i) senses the measured state of the running process, (ii) reads the live requirement and environmental conditions, (iii) re-solves the co-design at those conditions over the admissible architectures (resolve_at picks the cheapest feasible point), (iv) applies the chosen configuration by stepping the true plant, and (v) logs a ControlStep. Because the plant is stepped with the true process rather than a nominal transition, any divergence between plan and execution is absorbed by the next sensing step, which is what distinguishes feedback from open-loop replay. This is the co-design instance of control co-design (CCD, [40]) in its nested form, here the myopic variant (re-solve a single static co-design each step); the model is assumed known, so measurements update the carried state and the conditions but not the co-design model itself. Two natural extensions are a receding-horizon lookahead that plans several steps ahead with the vector DP and commits only the first, and online learning of the model from measurements (coupling in the online-learning layer of section 9); both are left for future increments.

15.7 Example 18: metabolic architecture switching↩︎

18_metabolic_switching.py. The first temporal example (Case 1 of section 15): the best architecture changes over time because the environment changes, and the object solved is a schedule of which architecture to run in each epoch, charged for switching. It exercises solve_schedule (13.5.1.5) over the Epoch (13.5.1.2) and Architecture (13.5.1.1) objects of section 15.

Model↩︎

An organism alternates between a glycolytic architecture (cheap per unit flux, capped growth rate, no overhead; models fast growth on glucose) and a gluconeogenic one (a larger per-flux burden and a fixed glyoxylate-shunt overhead, but a higher ceiling; models slower, costlier growth on acetate). Each architecture is a small co-design System: the outer functionality is the demanded growth rate \(\mu\in\mathbb{R}_{\geq 0}\) (1/h), and the outer resource is a dimensionless burden summed from an enzyme module (burden linear in \(\mu\), or \(\infty\) above \(\mu_{\max}\)) and a maintenance module (the fixed overhead). The environment is five epochs. Only the substrate-matching architecture is admissible per epoch, except a contested mixed epoch, flanked by acetate on both sides, where both are candidates, so the scheduler has a genuine choice. Serving mixed on the cheaper glycolytic pathway forces two extra switches (acetate \(\to\) glucose-type and back); riding it out on the incumbent costs none. The trade is resolved by the per-switch re-acclimation cost, the diauxic-shift lag.

Code↩︎

epochs = [
    Epoch("glucose_1", {"mu": 0.8}, candidates=[GLYCOLYTIC]),
    Epoch("acetate_1", {"mu": 0.5}, candidates=[GLUCONEOGENIC]),
    Epoch("mixed",     {"mu": 0.6},
          candidates=[GLYCOLYTIC, GLUCONEOGENIC]),
    Epoch("acetate_2", {"mu": 0.5}, candidates=[GLUCONEOGENIC]),
    Epoch("glucose_2", {"mu": 0.8}, candidates=[GLYCOLYTIC]),
]
# Same environment, two switching economies:
sched_lo = solve_schedule(epochs, cost_fn=burden_cost, switch_cost=0.05)
sched_hi = solve_schedule(epochs, cost_fn=burden_cost, switch_cost=0.8)

Output↩︎

Low re-acclimation cost (0.05 per switch):
  mixed      -> glycolytic     burden= 0.600  (+switch 0.05)
  total burden = 5.800, switches = 4, feasible = True

High re-acclimation cost (0.8 per switch):
  mixed      -> gluconeogenic  burden= 1.920
  total burden = 8.520, switches = 2, feasible = True

Mixed-epoch pathway: glycolytic (cheap switching)
                     -> gluconeogenic (costly switching)

Reading↩︎

Identical environment, different switching economics, qualitatively different metabolic program. Under a low switch cost the scheduler takes four switches and serves mixed on the locally cheaper glycolytic pathway (total burden \(5.800\)); under a high switch cost the two round-trip switches no longer pay, so the organism holds the incumbent gluconeogenic pathway through mixed (two switches, total burden \(8.520\)). The mixed epoch’s pathway is the visible pivot, and it is exactly the temporal coupling that a static co-design cannot see.

15.8 Example 19: rover module activation↩︎

19_rover_modules.py. Case 2 of section 15, with a genuinely carried scalar resource. A planetary rover activates and deactivates operating modes across mission phases on a battery that depletes and recharges from solar. It exercises solve_dynamic (13.5.2.6) over a StateGrid (13.5.2.3), returning a DynamicPolicy (13.5.2.4) rolled out from two starting charges.

Model↩︎

Each mode (drive, science, comms, survival) is a small co-design System that, given a phase capability demand, sizes a shared power bus and returns two outer resources: energy_Wh (subtracted from the battery) and a cost proxy that adds an opportunity penalty \(\mathrm{BEST}-\mathrm{value}\) for the mission value not delivered, kept non-negative as the resource poset requires. The outer object is a finite-horizon dynamic program over six phases: the per-phase decision is the mode; the per-phase cost is read from that mode’s co-design solve; and the battery state of charge \(s\in\mathbb{R}_{\geq 0}\) (Wh) is carried by the transition \(s'=\min(\text{cap},\,s-\text{spent}+\text{solar})\), admissible while \(s\ge 0\). A heavy-drawing mode buys capability now at the cost of charge later, so the optimal policy threads mode activations through the evolving state of charge, standard spacecraft power-mode practice.

Code↩︎

policy = solve_dynamic(
    stages, grid,
    cost_fn=mission_cost,          # point -> energy + opportunity penalty
    terminal_cost=terminal_reward, # small credit for leftover charge
)
# One policy, valid for any initial charge on the grid; roll it out twice.
full = rollout(policy, stages, BATTERY_CAPACITY_WH)   # 300 Wh
low  = rollout(policy, stages, 90.0)                  # depleted start

Output↩︎

Start at full charge (300 Wh):
  phase_1  -> science   soc_in= 300.0 Wh  ...  soc_out= 245.0 Wh
  ...
  phase_5  -> science   soc_in=  80.0 Wh  ...  soc_out=  25.0 Wh
  phase_6  -> comms     soc_in=  25.0 Wh  ...  soc_out=  10.0 Wh
  total cost = 727.50, feasible = True

Start at low charge (90 Wh):
  phase_1  -> comms     soc_in=  90.0 Wh  ...  soc_out=  75.0 Wh
  ...  (comms held every phase)  ...
  total cost = 918.00, feasible = True

Reading↩︎

The same policy, read at different starting charges, produces different module-activation schedules. From a full battery the rover runs high-value science until reserve falls, then steps down to comms for the last phase; starting depleted it can never afford science’s draw and holds comms throughout, so it sheds load exactly as a power-constrained mission would. The higher total cost from the depleted start (\(918.00\) against \(727.50\)) is the accumulated opportunity penalty of the deferred objectives, the price the DP pays to keep the battery non-negative.

15.9 Example 20: antichain-valued sequential co-design↩︎

20_sequential_codesign.py. The sequential layer of section 15.3: the value at each stage and state is a whole Pareto front of cumulative resource totals, not one number. Its antichain-valued Bellman recursion is solve_sequential (13.5.3.4), and the \((H1)/(H2\text{-joint})/(H3)\) guard is check_monotonicity (13.5.3.8).

Model↩︎

A multi-leg survey mission of four legs. At each leg the operator picks one of three incomparable modes, spanning a cost/CO\(_2\) trade (eco: low CO\(_2\), high cost; balanced; rapid: low cost, high CO\(_2\)), each drawing down a shared energy budget carried between legs. Because cost and CO\(_2\) are incommensurable, the value of a whole-mission plan is not a single number but a Pareto front of cumulative \((\text{cost},\text{CO}_2)\) totals, accumulated on two fixed axes by sum_combine, while the carried energy is a separate state axis read off the chosen point (admissible while non-negative). The framework guarantees the value front equals the reachable frontier: every returned point is achievable by some feasible mode sequence, and no non-dominated achievable total is missed.

Code↩︎

res = solve_sequential(
    stages, grid,
    cost_axes=["cost", "co2"],       # the two accumulated resource axes
    initial_state=ENERGY_CAPACITY,   # carried scalar budget (30 units)
    combine=sum_combine,             # totals add across legs
)
rep = check_monotonicity(stages, grid, cost_axes=["cost", "co2"])
# rep.monotone_value_guaranteed is True: (H1), joint (H2) and (H3) hold.

Output↩︎

Whole-mission Pareto front of (cost, CO2) totals (9 points):
  cost=   8.0   co2=  36.0        <- all rapid: cheapest, dirtiest
  cost=  12.0   co2=  31.0
  ...
  cost=  20.0   co2=  21.0        <- interior mixes of modes
  ...
  cost=  40.0   co2=   4.0        <- all eco: cleanest, costliest

Monotonicity guard: MonotonicityReport(H1=ok, H2=ok, H2joint=ok,
                    H3=ok, value_monotone=guaranteed)

Reading↩︎

The solver returns the exact reachable frontier of whole-mission plans: nine incomparable points running from all-rapid (\(8.0\), \(36.0\)) through all-eco (\(40.0\), \(4.0\)), with interior mixes between. Each point is a different per-leg mode schedule. The front’s size is the width of the reachable set; with a summed objective on two fixed axes it grows polynomially, here linearly, in the horizon, not exponentially, which is why it stays small and legible. The passing \((H1)/(H2\text{-joint})/(H3)\) guard certifies the value is monotone in the carried budget: more energy never shrinks the achievable front.

15.10 Example 21: reconfigurable robot with a state vector↩︎

21_reconfigurable_robot.py. The vector-state layer of section 15.4: the carried state is not a scalar but a vector. It exercises solve_vector_sequential (13.6.2.4) over a VectorStateGrid (13.6.1.8), guarded by check_vector_monotonicity over the product order.

Model↩︎

A field robot reconfigures between three morphologies across five legs: tracked (low energy, wears the track module fast), wheeled (wears the wheel module fast, spares the track), and hybrid (splits the load, wears both a little). Each morphology is a co-design System emitting \((\text{energy},\text{ops})\) plus the per-module wear increments it applies this leg. The carried state is a three-axis vector: per-module wear for the two independent drive modules plus a shared energy budget, advanced by \(\text{wear}'=\text{wear}+\Delta\text{wear}\) and \(\text{energy}'=\min(\text{cap},\text{energy}-\text{draw}+\text{solar})\), admissible while each wear stays under its limit and energy stays non-negative. Cost (energy and an operations penalty) accumulates on two axes; the value is therefore an antichain over the cost axes parametrised by the state vector, which a scalar DP could not represent.

Code↩︎

grid = VectorStateGrid([
    ContinuousAxis("track_wear", 0.0, WEAR_LIMIT, 11),
    ContinuousAxis("wheel_wear", 0.0, WEAR_LIMIT, 11),
    ContinuousAxis("energy", 0.0, ENERGY_CAPACITY, 13),
])
res = solve_vector_sequential(
    stages, grid, cost_axes=["energy", "ops"],
    initial_state={"track_wear": 0.0, "wheel_wear": 0.0,
                   "energy": ENERGY_CAPACITY},
    combine=sum_combine,
)

Output↩︎

grid size: 1573 state nodes

Mission Pareto front of (energy, ops) totals (5 points):
  energy=  13.0   ops=  17.0
  energy=  14.0   ops=  16.0
  energy=  15.0   ops=  15.0
  energy=  16.0   ops=  14.0
  energy=  17.0   ops=  13.0

Vector monotonicity guard: VectorMonotonicityReport(
    H1=ok, H2=ok, H2joint=ok, H3=ok, value_monotone=guaranteed)

Reading↩︎

The product grid of \(11\times 11\times 13 = 1573\) state nodes carries the full three-axis state. The mission front is a five-point \((\text{energy},\text{ops})\) staircase, each point a different morphology schedule: the low-energy end leans on the wheeled morphology, and the low-ops end mixes morphologies to spread wear so neither drive module hits its limit and forces an expensive fallback. This is the structured multi-component state of the Formula 1 seasonal problem (15.12) in a robotics setting. The passing vector guard certifies the value is monotone over the product order: more spare wear budget and more energy never shrink the achievable front.

15.11 Example 22: online feedback co-design↩︎

22_online_feedback_codesign.py. The closed-loop layer of section 15.6: no offline plan is followed; each step is solved against what was actually measured. It exercises run_online_codesign (13.6.3.4), whose audit trail is a list of ControlStep (13.6.3.1) records.

Model↩︎

A solar-powered environmental sensor node runs in the field over twelve control steps. At each step it senses its true battery charge, reads the live data-rate demand and environment, re-solves its co-design at those measured conditions, applies the cheapest feasible configuration, and steps the true plant. The node has three configurations, each a co-design System sizing radio and compute against a demanded rate and returning \((\text{energy},\text{ops})\): low_power (always feasible), nominal, and high_rate (feasible only at a high enough rate and when charge allows). The measured charge enters the solve as functionality, so an energy-hungry configuration whose draw would exceed the available charge is simply infeasible, which is how the feedback path gates the co-design. This is control co-design (CCD, [40]) in its nested, myopic form: re-solve one static co-design each step. A storm at steps 4–7 raises the demand and the day/night cycle changes solar recharge; the loop adapts to both.

Code↩︎

result = run_online_codesign(
    CONFIGS, n_steps=N_STEPS,
    sensor=sensor,             # reads true charge back (closes the loop)
    requirement=requirement,   # {"rate": demand, "charge": measured_soc}
    environment=environment,   # solar input this step
    plant=plant,               # steps the true battery, clips to [0, cap]
    cost_fn=cost_fn,           # minimise ops cost
    initial_state=BATTERY_CAPACITY,
)

Output↩︎

Closed-loop audit trail:
   t  soc_in  rate  solar  config     energy   ops
   3   100.0     3     14  nominal      11.0   3.0
   4   100.0     8     12  high_rate    34.0   1.0
   5    78.0     9     10  high_rate    37.0   1.0
   6    51.0     9      8  high_rate    37.0   1.0
   7    22.0     7      8  nominal      19.0   3.0
   8    11.0     3     10  nominal      11.0   3.0
   9    10.0     3     12  low_power     5.0   5.0
  10    17.0     3     14  nominal      11.0   3.0
total ops cost = 32.0, feasible = True

Reading↩︎

The node escalates to high_rate through the storm (steps 4–7) when the demanded rate justifies it. The closed loop then reads the depleting night-time charge and the re-solve gates the hungry configuration out: at step 9 charge has fallen to \(10.0\) and only low_power is feasible, before solar recovers and the node climbs back to nominal. No offline plan is consulted, and because the plant is stepped with the true process rather than a nominal transition, any divergence between plan and execution is absorbed by the next sensing step, which is what distinguishes feedback from open-loop replay. Two natural extensions, receding-horizon lookahead with the vector DP and online learning of the model (section 9), are left for later increments.

15.12 Example 23: hierarchical co-design for a Formula 1 season↩︎

23_formula1_season.py. A faithful reproduction, in the framework’s vocabulary, of the seasonal co-design of Neumann, Habermacher, Fieni, Cerofolini, Zardini, and Onder [21], and the canonical instance of the precompute-then-DP structure of section 15.5. Layer one freezes a CatalogDP (13.2.1.5) with precompute_catalog (13.5.3.11); layer two is a scalar-maximisation season DP built on solve_dynamic (13.5.2.6). A heavy runner.

Model↩︎

The two layers are explicit. The race-level co-design (layer one) is a genuine MCDP solve: for each track, battery size, and discretised incoming battery age, a CatalogDP over energy-deployment strategies emits a Pareto front of \((\text{race\_time},\text{wear})\), since deploying more electrical energy lowers race time but ages the battery faster, and an aged pack cannot reach the low race times a fresh one can. precompute_catalog freezes those fronts once. The season-level dynamic program (layer two) is a scalar-maximisation MDP carrying the state vector \((w_{b,1}, w_{b,2}, x_{\mathrm{ex}})\), the fractional wear of the two regulation-permitted battery units plus a flag recording whether a replacement penalty has been incurred. At each race the controls are which unit to run, which frozen deployment strategy to select, and whether to install a fresh unit (a grid penalty of ten places for the first replacement, five for each later one, per the modelled FIA rules). Race time maps to a finishing position through an empirical time-gap model and a probabilistic grid-start correction scaled by the track’s overtaking difficulty, and the finishing-position distribution is integrated against the FIA points table to give expected championship points; the DP maximises the season total by backward induction. No co-design solve happens inside the season sweep, which is what distinguishes this from the re-solving sequential DP.

Code↩︎

# Layer 1: freeze one (race_time, wear) front per track x battery x age.
catalog = precompute_catalog(
    tracks, battery_sizes, age_buckets, deploy_strategies)
# 112 fronts = 8 tracks x 2 batteries x 7 age buckets.

# Layer 2: scalar-max season DP over the frozen catalog, carrying
# (w_b1, w_b2, x_ex); no co-design solve inside the season sweep.
season = solve_season(catalog, races, points_table)

Output↩︎

112 race fronts precomputed (8 tracks x 2 batteries x 7 ages).
Optimal season expected points: 134.3
Optimal per-race policy:
  race         unit  batt      deploy     time  repl grid E[pts]
  Monza           1   4MJ deploy_1.00   5024.0   yes   13   15.7
  ...
  Zandvoort       1   4MJ deploy_1.00   5191.2    no    3   15.9
Finding 1 (local penalty, global gain): 2 replacement(s).
Finding 2: forward total = 134.3, reversed total = 134.8
  forward replaces at ['Monza','Spa']; reversed at ['Spa','Monza'].
Rendered paper-analogue figures (outputs/f1_paper_fig1..3.png)

Reading↩︎

Two of the paper’s findings are reproduced. First, the optimal policy accepts a local penalty for a global gain: it takes a battery replacement at one race, sacrificing points that day (a worse grid slot), to keep a fresh, low-wear pack available for aggressive deployment across the remaining races, raising the season total to \(134.3\). Second, race order shifts the optimal policy: because the grid-penalty cost is track-dependent, reversing the calendar moves the optimal replacement onto a cheaper track, so the attainable total is near-invariant (\(134.8\) against \(134.3\)) while the replacement set reorders, the temporal coupling the paper emphasises. The season DP is validated against exhaustive brute-force enumeration in tests/test_formula1.py; the DP optimum matches exactly. The example and notebook also regenerate the paper’s three key figures in the same format for visual comparison (the race-time-versus-wear fronts, the grid-start position penalty, and the finishing-position density); these now write to the outputs/ directory as f1_paper_fig1, _fig2, and _fig3. They reproduce the paper’s figure structure and qualitative behaviour from the framework’s own co-design solves, not the paper’s specific numbers, since the paper’s fronts come from an optimal-control lap simulation and a fitted position model that are not reproduced here. The sense in which the framework “does the same thing” is that it produces the same structured artefacts feeding the same season-level dynamic program.

15.13 Example 24: catalog-driven car co-design from one table↩︎

24_car_catalog_codesign.py. The co-design counterpart of the hand-wired car study (14.17): here a single architecture table plus one build_architecture() function replaces the three separate build_ice/build_hybrid/build_ev builders, and every discrete powertrain choice becomes a CatalogDP (13.2.1.5) entry rather than a bespoke module. It is the whole-vehicle instance of catalog-driven assembly (13.5.3.12) over a System (13.3.2.1). It is a heavy runner, though it completes in about one second.

Model↩︎

A passenger vehicle is decomposed into 22 subsystems (powertrain, chassis, interior and auxiliary). A single 12-row ARCHITECTURE_CATALOG drives the whole study: each row pre-selects the discrete powertrain choices (engine, transmission, e-motor topology, boost, battery strategy, charger, drivetrain) for one point on the modern technology spectrum, from pure ICE through mild, full, and plug-in hybrid, range-extender EV, to battery-electric, while the parametric modules size themselves from mission demand. Body, suspension, and tyres stay full CatalogDPs, swept one row at a time by solve_architecture so that each assembled System is single-valued and its constraint network carries only monotone quantities. The outer mission is a compact C-segment target (5 seats, 370 L, 170 km/h, 500 km range, 0–100 in \(11.5\) s); the outer resources are production cost, curb weight, unified primary energy per 100 km, fuel per 100 km, tailpipe CO\(_2\), and durability. Two coupled feedback cycles are closed by Kleene iteration on the monotone network: the mass spiral (every load-bearing module reads the total curb mass, which their own masses in turn determine) and the energy/battery spiral (consumption depends on mass, and for a BEV the battery is sized to range \(\times\) consumption, feeding back into mass). The unified energy metric is primary energy input, which is why an ICE lands near 48 kWh/100 km while a BEV lands near 15: the combustion path discards most of the fuel’s chemical energy as heat.

The data-provenance policy is explicit. A SOURCES block at the top of the file footnotes every value that was added or spot-checked against public data (2024–2025): battery pack price (BloombergNEF 2024), pack gravimetric density (a 2025 World Electric Vehicle Journal review), fuel lower-heating values and CO\(_2\) factors (UBC and Waterloo notes), engine dry masses, and the Corolla-bracketed mission spec. Values not individually footnoted are kept from the scaffold and lie within the cited published ranges. They are illustrative, not OEM-specific: the MCDP framework is what is being validated; the numbers serve the framework.

Code↩︎

# One 12-row table drives the whole study (name, indices, params):
ARCHITECTURE_CATALOG = [
  ("ICE_economy",    1, 0, 0, 1.4, 1.00,  0.0,  0.0, F, F, F),
  ("HEV_full",       2, 3, 2, 1.0, 0.60,  2.0,  0.0, T, F, F),
  ("BEV_long_range", 0, 4, 4, 1.0, 0.00, 80.0,150.0, T, F, F),
  # ... 9 more rows spanning MHEV / PHEV / REEV / diesel ...
]

def build_architecture(arch, mission, *, body_row, susp_row, tire_row):
    sys = System(arch["name"])
    sys.add("engine", make_engine_dp([eng]))   # 1-entry CatalogDP slice
    # ... 21 more modules ...
    # Mass spiral: load-bearing modules read the converged curb mass.
    sys.constrain("susp_front.axle_load_kg", lambda x: 0.55*total_mass(x))
    sys.constrain("brakes_front.vehicle_mass_kg", total_mass)
    # Energy/battery spiral: BEV pack sized to range x consumption.
    sys.constrain("battery.target_battery_kWh", battery_target_kWh)
    return sys.build()   # Kleene iteration closes both cycles

Output↩︎

--- Per-architecture feasibility (cheapest feasible design) ---
  Architecture      cls  feas   cost$  mass kg  kWh/100  L/100  CO2
  ICE_economy       ICE  yes   17,447      930     47.8    5.3   123
  Diesel_long_range ICE  yes   24,400     1081     40.6    4.1   109
  HEV_full          FHEV yes   26,945     1315     42.2    4.1    94
  PHEV_small        PHEV yes   36,162     1537     34.9    3.1    71
  REEV              REEV yes   37,182     1613     21.6    1.1    26
  BEV_long_range    BEV  yes   44,572     1788     14.8    0.0     0
  BEV_AWD_perf      BEV  yes   68,076     2517     17.0    0.0     0

--- Pareto front over (cost, energy/100km, mass) ---
  ICE_economy    $ 17,447   47.8 kWh/100    930 kg
  Diesel_long..  $ 24,400   40.6 kWh/100   1081 kg
  PHEV_small     $ 36,162   34.9 kWh/100   1537 kg
  REEV           $ 37,182   21.6 kWh/100   1613 kg
  BEV_long_range $ 44,572   14.8 kWh/100   1788 kg

Reading↩︎

All twelve architectures are feasible for the compact mission, and the cheapest-feasible sweep recovers the expected technology ordering: the ICE economy car is lightest (\(930\) kg) and cheapest (\(\$17{,}447\)) but thirstiest (\(47.8\) kWh/100 km, \(123\) g/km CO\(_2\)), while the battery-electric long-range car is the most efficient (\(14.8\) kWh/100 km, zero tailpipe CO\(_2\)) at the highest mass (\(1788\) kg) and near the highest cost. Between them the hybrids and the range-extender trade cost against energy and CO\(_2\) monotonically, and the mass climb from ICE to BEV is the visible signature of the energy/battery spiral: a BEV’s pack is sized to range \(\times\) consumption, and the extra pack mass raises consumption, which the Kleene loop resolves to a fixed point. The five-point Pareto front over (cost, energy, mass) keeps only the non-dominated members, the same structured artefact the temporal examples consume, produced here from a single declarative table rather than three hand-written builders.

15.14 Example 25: online co-design paper benchmarks↩︎

25_online_paper_benchmarks.py. A statistical replication of the synthetic benchmarks of Alharbi, Dahleh, and Zardini, “Compositional Online Learning for Multi-Objective System Co-Design” [18], the paper behind the online-learning layer of section 9. It drives the library’s optimistic evaluators, MonotonicityEvaluator (13.4.2.2) and LipschitzEvaluator (13.4.2.3), whose bound() reproduces the paper’s eqs. (23)–(25); the outer sampling loop is written out here rather than delegated to solve_online (13.4.2.7) because the library scans a candidate list with an upper-confidence picker whereas the paper samples a base measure with rejection. A heavy runner (about 6 s with the EA panel).

Model↩︎

The paper studies online multi-objective decision-making in monotone co-design and recovers the target-feasible antichain of non-dominated resources using as few expensive evaluations as possible. Its engine is Algorithm 1 (Rejection Sampler with Optimistic Evaluators): draw a candidate from a low-discrepancy Halton base measure, compute history-dependent optimistic bounds, and skip it without evaluation when either (13) its optimistic resource is already dominated by the incumbent antichain, or (14) its optimistic functionality can never meet the target. Only survivors consume the budget. Two synthetic families of Section VII-B are replicated: E1/Monotone (\(g:[0,1]^3\!\to\![0,1]^2\), a weighted sum of thresholded indicators; the monotone join bound) and E2/Lipschitz (\(g:[0,1]^4\!\to\![0,1]^2\), a triangle wave through an \(L{=}2\) linear map; the Lipschitz cone bound). The score is the cumulative hypervolume difference (HVD), \(\sum_t [\mathrm{HV}(A_{\mathrm{ref}}) - \mathrm{HV}(\hat{A}_t)]\), of the incumbent antichain against a reference point \((1,1)\); lower is better and HVD\({}\to 0\) on convergence. Baselines are the same Halton proposals with no elimination, uniform random draws over the shared discretised pool, and an optional pymoo EA panel (NSGA-III / MOEA/D / RVEA).

Two deviations from the paper are documented in the source and preserved here. First, taken literally the monotone family is degenerate (an unconstrained monotone \(g\) is minimised at the box corner, so there is nothing to learn); the example therefore makes the target binding, recovering \(\mathrm{FixFunMinRes}(f)\) over the feasible upper set, which exercises both elimination conditions and leaves the Lipschitz family untouched. Second, families E3 (intermodal mobility) and E4 (heterogeneous multi-robot) are not replicated: their expensive blocks are external models not present in this repository. Because the paper’s seeds, grid, and atom count are unpublished, the replication is statistical, not bit-exact; only the relative claims are validated.

Code↩︎

# Paper Algorithm 1: rejection sampler with optimistic evaluators.
while accepted < budget and ptr < len(pool):
    i = ptr; ptr += 1
    if front and rng.random() >= delta:            # else force-accept
        lo, hi = evaluator.bound(cand[i])          # optimistic bounds
        if constrained and (hi[RC[0]] < f[0] or hi[RC[1]] < f[1]):
            continue        # (14): optimistic functionality misses target
        if _dominated_by_incumbent(lo, front):
            continue        # (13): optimistic resource already dominated
    r = (float(res[i, 0]), float(res[i, 1]))        # accept: query block
    evaluator.observe(i, cand[i], _singleton_antichain(r))
    accepted += 1
    if feas[i]:
        front = pareto_insert(front, r)
    cum += hv_star - hv2d(front, ref)               # cumulative HVD

Output↩︎

TABLE I  -- Monotone problems (E1), FixFunMinRes(f)   [M1..M4 shown]
  Ours        9.51e+00  7.94e+00  4.85e+00  1.50e+01
  Halton      1.55e+01  1.46e+01  7.29e+00  2.63e+01
  MOEA/D      1.62e+01  1.50e+01  1.78e+01  2.30e+01

TABLE II -- Lipschitz problems (E2), pure minimization  [L1..L4 shown]
  Ours        1.07e+01  2.56e+01  1.26e+01  9.25e+00
  Halton      1.15e+01  2.75e+01  1.29e+01  1.05e+01
  MOEA/D      9.00e+00  2.98e+01  8.71e+00  1.30e+01

VALIDATION
(a) Ours < Halton on ALL monotone instances
    ours : Ours beats Halton on 8/8  -> PASS
(b) Ours < Halton on MOST Lipschitz instances
    ours : Ours beats Halton on 8/8 (need >= 6) -> PASS
(c) [EA panel] paper: Ours best-or-2nd; MOEA/D blowups on Lipschitz.
    ours : best-or-2nd monotone 8/8, lipschitz 3/8
    ours : Lipschitz MOEA/D max/med=3.3x   Ours max/med=2.3x
    -> NOT REPRODUCED at higher statistical strength (WP-P3).

VERDICT: (a) PASS, (b) PASS  [elimination-vs-Halton reproduce]
         (c) NOT REPRODUCED -- EA panel.

Reading↩︎

The replication is deliberately honest about what carries and what does not. Claims (a) and (b) reproduce robustly: the optimistic elimination beats the identical Halton proposals with no elimination on all eight instances of both families (here \(8/8\) and \(8/8\); even stronger than the paper’s \(7/8\) under independent re-runs at \(\mathrm{RUNS}{=}20\), \(\mathrm{ITERS}{=}500\), two master seeds and the full EA panel), and Ours is best-of-all on every monotone instance. These are the fair apples-to-apples comparison, since all three methods draw from the same shared discrete pool and differ only by the elimination step. Claim (c) did not reproduce: on the smooth Lipschitz family the continuous EAs (MOEA/D in particular) match or beat Ours, and MOEA/D shows no blowups here (max/median \(\approx 3.3\times\), not the paper’s \(\approx 20\times\) on L8), so it is neither poor nor high-variance in this harness. The leading hypotheses are that the EA panel optimises the continuous box while Ours, Halton and Random are confined to the fixed discrete pool (an asymmetry favouring EAs on smooth maps, with EA HVD clipped at zero against a discrete reference they can beat), that the budget is far below paper scale, and that the EAs are untuned (fixed \(\text{pop}=20\)). The example prints NOT REPRODUCED rather than a fragile PASS, and mirrors the reasoning in its own Validation status docstring; only the elimination-vs-Halton claims are counted as reproduced. Two log-scale HVD-versus-iteration figures are written to outputs/ (25_online_monotone.png and 25_online_lipschitz.png).

16 Realistic Problem Domains for Temporal Co-Design↩︎

The temporal, vector-state, and online layers of section 15 were built with concrete application domains in mind. This section surveys realistic problems, drawn from the recent literature, that map onto these layers, as candidates for future worked examples and as evidence that the abstractions are not invented in a vacuum. Each paragraph names the co-design structure (which layer applies) alongside the domain.

16.1 Self-reconfiguring and modular robots↩︎

Self-reconfiguring modular robots deliberately change their morphology, by rearranging the connectivity of their modules, to adapt to terrain, task, or damage [41]: a worm shape for a pipe, legs for rough ground, a wheel for flat terrain. This is directly the temporal switching and vector-state setting: the morphology is the architecture, the terrain or task sequence is the environment, and per-module wear plus shared energy is the carried state vector (example 21 is a first cut). Recent work makes the coupling sharper. Terrain-aware morphology search in dynamic environments [42] chooses a morphology per scene under an energy budget; multi-task space applications map task requirements to modular-unit morphology through an inverse-mapping design step [43]; and reconfigurable cobots such as the CONCERT platform and freeform strut-node systems (FreeSN, SMORES) expose a catalog of admissible configurations with real reconfiguration costs, the switching-cost term of solve_schedule. Floor-cleaning tiling robots (hTrihex, hTetro) reconfigure to cover obstacle-strewn areas at lower energy than fixed-shape coverage planning [44], an energy-versus- coverage Pareto per configuration.

16.2 Programmable self-assembly↩︎

Programmable self-assembly designs building blocks with tunable binding interactions so they assemble into a target structure. The design problem is multi-objective, trading assembly yield, assembly speed, and economy (number of distinct particle species), which is exactly a Pareto-front co-design. Hübl and Goodrich show that simultaneously optimising assembly time and yield through the binding-energy and concentration design space can speed assembly by orders of magnitude without lowering yield, and that economical, semi-addressable designs (reusing building blocks) can beat fully-addressable ones on all three axes [45], [46]. The temporal angle is the assembly protocol: a staged schedule of temperature or concentration set-points is a sequence of co-design decisions coupled by the carried state of intermediate-species concentrations, a natural fit for the vector-state sequential DP, with Markov-state-model protocol optimisation [47] as the dynamic backbone.

16.3 Multi-objective protein design↩︎

Protein design inherently trades competing objectives, stability, function, solubility, immunogenicity, and manufacturability, and Pareto-based multi-objective sampling produces candidate sets spanning these tradeoffs for informed downstream selection [48]. This is a co-design antichain over design objectives. Recent methods make the front explicit: ST-PARM aligns inference-time generation to a controllable stability-versus-function tradeoff; MosPro and NSGA-2 variants produce Pareto-optimal sequence sets; MAProt negotiates a structure-versus-function consensus across agents. The temporal and online angle is closed-loop directed evolution: frameworks such as EVOLVEpro couple a generative or predictive model with iterative wet-lab feedback, proposing a batch, measuring it, and updating, which is precisely the online feedback co-design loop of section 15.6 with the measured functionality replacing the modelled one (the model- learning extension of that loop).

16.4 Control co-design (CCD)↩︎

Control co-design, the simultaneous optimisation of a physical plant and its controller, is an established engineering discipline, described by Garcia-Sanz as a game changer [40], with two standard paradigms: nested (an outer plant loop wrapping an inner control solve) and simultaneous [49]. The nested paradigm is exactly the structure of the online co-design loop. Concrete instances that map onto the temporal and online layers include robust control co-design with a receding-horizon MPC inner loop [50]; plant-and-controller optimisation for energy systems with MPC [51]; control co-design of spacecraft trajectory and rocket-engine design in NASA’s OpenMDAO/GMAT tooling; digital-twin control co-design of active suspensions via reinforcement learning; and adaptive MPC with online recursive-least-squares identification for load-frequency control under renewable disturbances, the archetype of the model-learning extension. These supply both the vocabulary (plant, controller, nested vs simultaneous) and the benchmark problems (active suspension, hybrid-electric vehicle, airborne wind energy) for positioning the online layer.

16.5 Summary of the mapping↩︎

Across these domains the same three structural questions recur, and the package answers each: which architecture to run as conditions change (switching, temporal); how a resource or wear vector carried across stages shapes the optimal schedule (vector-state DP, vector_dp); and how to re-decide the design online from measurements (feedback, online_codesign). The reconfigurable robot (21) and adaptive sensor node (22) are first worked examples in the robotics and control-co-design cells of this table; self-assembly and protein design are identified as high-value future examples whose multi-objective, staged-protocol, and closed-loop-evolution structure the existing layers already support.

17 Modelling Guidelines and Pitfalls↩︎

A few patterns that come up repeatedly when authoring MCDPs.

17.0.0.1 Expose loop variables you want to inspect.

When using the operator API directly, Loop projects its axis out of the outer \(R\). To inspect the converged loop value, include it in the inner \(R\) under a different name: for example, both battery_mass as the loop axis and report_mass mirrored for outer-\(R\) visibility. The System builder handles this transparently.

17.0.0.2 Cap physical maxima.

When a design variable has a physical ceiling (\(v_{\max}\), \(r_{\max}\)), make \(h\) return a \(\top\)-valued antichain whenever the iterate exceeds it. The Kleene ascent will then converge to infeasibility rather than oscillating or diverging.

17.0.0.3 Use FunctionDP or CatalogDP for breadth.

AlgebraicDP always returns a singleton antichain. For a genuine Pareto front, use FunctionDP (enumerate the tradeoffs explicitly) or CatalogDP (provide multiple incomparable entries).

17.0.0.4 Scalarisation is downstream.

The MCDP solver returns the Pareto front. Choosing which point to ship is a separate concern: apply minimize_cost with a scalar objective, or inspect the front and pick by hand.

17.0.0.5 Catalog Cartesian product can blow up.

When several subsystems return multi-valued antichains, the inner \(h\) of the System takes their Cartesian product. For typical models (one or two catalog subsystems among many algebraic ones), this is not a bottleneck; but if many subsystems are catalog-driven, the inner antichain may grow large between \(\mathsf{Min}\) prunings.

17.0.0.6 Numerical precision in feedback loops.

Antichain.eq compares points exactly. For contractive floating-point feedback, the iteration may oscillate by an ulp near the fixed point and consume many iterations before converging. The divergence cap and max_iter guard against runaway, but choosing max_iter generously (200 to 500) is wise for nontrivial models.

18 Related Work and Further Reading↩︎

This section maps the literature this library implements. It is a reading guide, not a survey. The order is roughly chronological within each theme.

18.0.0.1 Founding theory.

The monotone co-design framework is due to Andrea Censi. The monograph [1] is the primary reference: posets, antichains, design problems as monotone relations, the three composition operators, the Kleene solver, and the approximation theory. It has been revised several times on arXiv and, as far as we can establish, has never appeared in a journal, so cite the preprint. [3] is the short conference statement of the same ideas. [19] is the standard background on lattices and order.

18.0.0.2 Cyclic constraints.

[4] isolates the class of co-design problems whose constraint graph contains feedback loops and gives the solution method. This is the theory behind Loop (Section 6) and behind every cyclic example in the manual.

18.0.0.3 Uncertainty.

[5] treats uncertainty in monotone co-design and is the source of the upper/lower bracketing that UncertainDP implements. More recent work makes uncertainty composable and parametric [13], gives it a categorical footing in symmetric monoidal categories [14], and handles distributional uncertainty with adaptive decision-making [15]. The library implements the bracketing and a set-based and stochastic layer on top of it; it does not implement the categorical constructions.

18.0.0.4 Beyond monotonicity.

[20] studies co-design problems that violate the monotonicity assumption. Nothing in this library addresses that case.

18.0.0.5 Mobility systems.

The application of co-design to future mobility begins with [6] and continues through [7] on autonomous-vehicle-enabled mobility and [8] on tooling for impact assessment. These are the model for the multi-modal and fleet-scale examples in Section 14.

18.0.0.6 Autonomous systems.

[9] runs co-design from hardware selection through to control synthesis. [10] gives a structured approach to co-designing embodied intelligence. [11] does task-driven modular co-design of vehicle control systems, and is the closest published relative of the System builder (Section 12).

18.0.0.7 Perception and compute.

[12] co-designs perception and decision-making subject to a compute budget, on autonomous vehicles. This is the reference for any model in which sensing quality, compute, power, and task performance are traded against each other.

18.0.0.8 Recent extensions.

[16] generalises design problems to quantale-enriched categories, aiming at quantitative heterogeneous design. [17] attacks scalability through linear design problems. Neither is implemented here.

18.0.0.9 Online learning.

[18] gives the compositional online learning method for multi-objective co-design. The codesign.online module (Section 9) is a port of it, including the certified optimistic bound.

18.0.0.10 Sequential and hierarchical applications.

[21] combines monotone co-design with sequential optimisation over a Formula 1 season. It is the published work closest in spirit to the temporal layers of Section 15, and it was written independently of them.

18.0.0.11 Categorical background.

[52] is a book in preparation on applied category theory for engineering, which covers the categorical setting that design problems live in.

18.0.0.12 What is not from this literature.

The temporal, sequential, vector-state, and receding-horizon layers of Section 15, and the problem domains surveyed in Section 16, are the author’s own work. Architecture switching by Viterbi over epochs, the vector-state backward Bellman recursion, the factorisation at reset points, and the receding-horizon online loop are not taken from the co-design papers above. They are being prepared for separate publication and have not been published or peer reviewed. Readers should treat them as software documentation rather than as established results, and should not cite them as though a paper existed. Where those layers make contact with the wider control-co-design literature, the connection is noted in place: the nested and simultaneous paradigms [49] and control co-design as a discipline [40].

19 Limitations and Future Work↩︎

Version 0.2.1 covers the algorithmic core: posets, antichains, six primitive DP types, three composition operators, the Kleene solver, and the two builders. The following are intentional omissions, all tractable extensions:

• MCDPL surface syntax. The paper’s mcdp { ... } text format of [2] is not parsed. The MCDP builder gives the same shape in Python.

• Non-finitely-representable antichains. The library works on finite, exactly-represented antichains. For relations whose Pareto front is continuous, the only support is the bracket pattern of UncertainDP. Approximation kernels (Sec. VII of [1]) are not yet implemented.

• Visualisation of the design graph. Example 5 visualises the Kleene ascent over \(\mathbb{N}\times \mathbb{N}\), but there is no general renderer for the operator tree or for the subsystem network of a System.

• Symbolic/automatic differentiation of \(h\). All relations are evaluated numerically.

• Distributed or parallel solving. The Kleene iteration is sequential.

Contributions, particularly on these items, are welcome.

20 API Reference Summary↩︎

The public API exported by import codesign is summarised below.

20.0.0.1 Posets.

Reals, Naturals, Ports, Discrete, abstract base Poset.

20.0.0.2 Antichains.

Antichain.

20.0.0.3 Design problems.

DesignProblem (abstract base), AlgebraicDP, FunctionDP, CatalogDP, CatalogEntry, ConstraintDP, ODE_DP, UncertainDP.

20.0.0.4 Composition.

Series, Parallel, Loop (classes); series, par, loop (function aliases).

20.0.0.5 Primitives.

adder, multiplier, scale, constant, identity.

20.0.0.6 Solver.

solve, kleene_loop, minimize_cost, SolveResult, TraceEntry. Warm-start via the start_from= argument on solve.

20.0.0.7 Builders.

MCDP, System.

20.0.0.8 Class-based modules.

Module (declarative DesignProblem base class).

20.0.0.9 Constraint DSL.

Port, Expr, ModuleHandle (returned by System.add), and the helper functions sqrt, exp, log for use inside constraint expressions.

20.0.0.10 Uncertainty.

Set-based: UncertaintySet (abstract), Box, Ellipsoid, Disk, Circle. Stochastic: Stochastic, plus copulas Independence and GaussianCopula. Result: UncertaintyResult. Entry point: solve_with_uncertainty (the same machinery is also reachable through solve(dp, f, uncertainty=[...])).

20.0.0.11 Online learning.

OptimisticEvaluator (abstract); MonotonicityEvaluator, LipschitzEvaluator, LinearParametricEvaluator. Result: OnlineResult. Entry point: solve_online.

20.0.0.12 Visualisation.

Imported as from codesign import viz. Functions: plot_antichain, plot_convergence, plot_uncertainty, to_dot.

20.0.0.13 Module-level.

__version__ (the string "0.2.1").

References↩︎

The bibliography is generated from references.bib in the same directory. Entries in the co-design group were checked against a primary record before being written: the arXiv API entry for preprints, Crossref for anything carrying a DOI. Author lists are printed in full.

The bioprocess sources are the parameter-calibration basis for examples 15 and 16, and the automotive sources for example 17. All calibration values in this manual are illustrative, drawn from the published ranges in those sources, and are not OEM- or product-specific.

How to cite↩︎

If you use codesign-mcdp in academic work, please cite it. A machine-readable CITATION.cff file at the repository root lets GitHub’s “Cite this repository” button generate a citation automatically. In BibTeX:

@software{briat_codesign_mcdp,
  author  = {Briat, Corentin},
  title   = {codesign-mcdp: A Python Library for
             Monotone Co-Design Problems},
  year    = {2026},
  version = {0.2.1},
  url     = {https://github.com/cbriat/codesign-mcdp}
}

Citing the software is not a substitute for citing the theory. Please cite [1] for the monotone co-design framework itself, and the paper behind whichever layer you use: [5] for uncertainty, [18] for the online-learning layer, and the relevant application paper from the co-design group below.

[52]

References↩︎