# Nexus — complete documentation

> Generated from the same markdown that is published as WIKI.md in each
> repository and rendered at https://vishalram24.github.io/nexus/docs/. Every code example here has been
> executed against a clean install.

Install: `pip install nexus-energy` (pulls in nexus-opt automatically).

---

# nexus-energy — User Guide

**Energy-system optimisation in Python — that also answers backwards.**

Build and optimise an energy-system digital twin the way you do today, at exact
parity with PyPSA and GenX. Then ask the question they can't: *which of my
inputs is wrong?*

`nexus-energy` is for engineers and researchers who build energy-system models
in Python — capacity expansion, economic dispatch, unit commitment, optimal
power flow — and who need the model to agree with a system that actually
exists.

> **Every code block in this guide has been executed against a clean install.**
> Where a printed result is shown, that is the actual output.

---

## Contents

**Getting started**

1. [What this library does](#1-what-this-library-does)
2. [Install and verify](#2-install-and-verify)
3. [Your first model](#3-your-first-model)
4. [Adding time](#4-adding-time)
5. [The five building blocks](#5-the-five-building-blocks)
6. [Reading results](#6-reading-results)

**Doing real work**

7. [Storage](#7-storage)
8. [Unit commitment](#8-unit-commitment)
9. [Capacity expansion](#9-capacity-expansion)
10. [Networks and transmission](#10-networks-and-transmission)
11. [Policy constraints](#11-policy-constraints)
12. [Coming from PyPSA](#12-coming-from-pypsa)
13. [Multi-stage planning](#13-multi-stage-planning)

**The part that is new**

14. [Differentiable dispatch and inverse calibration](#14-differentiable-dispatch-and-inverse-calibration)
15. [MPC and auto-calibration](#15-mpc-and-auto-calibration)

**Making it fast**

16. [Solver controls](#16-solver-controls)
17. [Temporal aggregation](#17-temporal-aggregation)
18. [Decomposition at scale](#18-decomposition-at-scale)

**Reference**

19. [Uncertainty — stochastic and robust](#19-uncertainty--stochastic-and-robust)
20. [AC power flow](#20-ac-power-flow)
21. [ML-guided solving](#21-ml-guided-solving)
22. [The component library](#22-the-component-library)
23. [Benchmarks](#23-benchmarks)
24. [Honest scope and known limits](#24-honest-scope-and-known-limits)
25. [Troubleshooting](#25-troubleshooting)
26. [API index](#26-api-index)

---

## 1. What this library does

Every production energy tool solves the problem **forward**: given these costs
and efficiencies, here is the optimal plan. That is well served already, and
`nexus-energy` does it too — at exact parity with PyPSA and GenX, several times
faster.

The expensive daily question is the **reverse**. Your twin's output does not
match the real system — which input is wrong, and by how much? Today that is
answered by brute force: perturb a parameter, re-solve, repeat. Morris, Sobol
and MGA sweeps cost hours to days, scale with every parameter, and return a
confident number even when the data cannot identify it.

`nexus-energy` differentiates through the optimisation itself, so the answer
comes from the model's own structure — in a handful of solves — and it says so
when the data is silent. That is [section 14](#14-differentiable-dispatch-and-inverse-calibration).

You model by placing components on a network graph:

```
                  ┌──────────────────────┐
                  │     EnergySystem     │
                  └──────────┬───────────┘
         ┌───────────────────┼───────────────────┐
┌────────▼────────┐ ┌────────▼────────┐ ┌────────▼────────┐
│       Bus       │ │    Generator    │ │     Storage     │
│  (electricity,  │ │  (solar, gas,   │ │ (battery, pump, │
│   H2, heat, CO2)│ │   hydro, …)     │ │  thermal, V2G)  │
└────────┬────────┘ └─────────────────┘ └─────────────────┘
         │
         ├───────────────────┐
┌────────▼────────┐ ┌────────▼────────┐
│      Load       │ │      Link       │
│  (demand side)  │ │ (line, HP, HVDC,│
└─────────────────┘ │  electrolyser)  │
                    └─────────────────┘
```

---

## 2. Install and verify

```bash
pip install nexus-energy
```

This pulls in [`nexus-opt`](https://github.com/VishalRam24/nexus-opt), the Rust
solver core, automatically. No separate step, no Rust toolchain. Python 3.9+.

Verify:

```python
import nexus_energy as ne

sys = ne.EnergySystem("check")
elec = sys.add_bus("elec")
sys.add_generator("g", bus=elec, capacity=100, marginal_cost=10)
sys.add_load("d", bus=elec, amount=50)
print(sys.optimise().total_cost)      # 500.0
```

### Development install

```bash
git clone https://github.com/VishalRam24/nexus-energy
cd nexus-energy
uv sync
uv run pytest
```

Benchmark scripts that compare against pandapower, PowerModels.jl or GenX
expect those reference installations outside this repo and skip automatically
when absent. Point them elsewhere with `NEXUS_PANDAPOWER_DIR` /
`NEXUS_POWERMODELS_DIR`.

---

## 3. Your first model

One bus, two generators, one load. Solar is free, gas costs 50/MWh, demand is
300 MW.

```python
import nexus_energy as ne

sys = ne.EnergySystem("my_system")
elec = sys.add_bus("elec", carrier="electricity")

sys.add_generator("solar", bus=elec, capacity=500, marginal_cost=0)
sys.add_generator("gas",   bus=elec, capacity=200, marginal_cost=50)
sys.add_load("demand", bus=elec, amount=300)

result = sys.optimise()
print(result.status, result.total_cost)
print(result.generator_dispatch)
```

```
optimal 0.0
{'solar': array([300.]), 'gas': array([0.])}
```

Solar covers the whole load, so the cost is zero and gas stays off. That is the
merit order, and it is the entire idea of economic dispatch.

Note that `generator_dispatch` values are **arrays**, not scalars — this model
has one timestep, so each array has length 1.

---

## 4. Adding time

Real models have profiles. `set_timesteps(n, dt)` sets the number of snapshots
and the hours each one spans; loads and availabilities then take arrays.

```python
import numpy as np
import nexus_energy as ne

sys = ne.EnergySystem("with_time")
elec = sys.add_bus("elec")
sys.set_timesteps(6, dt=1.0)

sys.add_generator("solar", bus=elec, capacity=400, marginal_cost=0,
                  carrier_factor=np.array([0, .2, .8, 1., .5, 0.]))
sys.add_generator("gas", bus=elec, capacity=500, marginal_cost=60)
sys.add_load("demand", bus=elec,
             amount=np.array([200, 250, 300, 350, 300, 220.]))

r = sys.optimise()
print(r.status, r.total_cost)
print("solar:", r.generator_dispatch["solar"])
print("gas  :", r.generator_dispatch["gas"])
print("price:", r.bus_shadow_prices["elec"])
```

```
optimal 41400.0
solar: [  0.  80. 300. 350. 200.   0.]
gas  : [200. 170.   0.   0. 100. 220.]
price: [60. 60.  0.  0. 60. 60.]
```

Three things to read out of that:

- **`carrier_factor`** is the per-timestep availability in `[0, 1]`. Solar's
  usable output is `capacity × carrier_factor[t]`, so at `t=3` its ceiling is
  400 MW and it serves the whole 350 MW load.
- **Solar is curtailed at `t=3`**: it could produce 400 but only 350 is needed.
- **`bus_shadow_prices`** is the marginal price of energy at that bus — the
  cost of one more MWh of demand. It is 60 when gas is on the margin and 0
  when free solar is, which is exactly the market-clearing price.

`dt` matters for anything integrated over time: storage state of charge, energy
costs, emissions. With `dt=0.25` each snapshot is 15 minutes.

---

## 5. The five building blocks

### Bus

A connection point. Everything attaches to a bus, and each bus enforces energy
balance at every timestep.

```python
elec = sys.add_bus("elec", carrier="electricity")
h2   = sys.add_bus("h2",   carrier="hydrogen")
```

The carrier sets default units (MWh for electricity, tCO₂ for carbon) and lets
the library check that you are not wiring hydrogen into an electricity balance.

### Generator

`add_generator(name, bus, capacity, marginal_cost=0.0, **kwargs)`

| Parameter | Unit | Meaning |
|---|---|---|
| `capacity` | MW | nameplate output |
| `marginal_cost` | $/MWh | variable cost of production |
| `capital_cost` | **$/MW/year** | annualised capex — only used when `extendable=True` |
| `carrier_factor` | array in `[0,1]` | per-timestep availability |
| `p_min` | MW | minimum stable output |
| `emission_factor` | tCO₂/MWh | for carbon pricing and caps |
| `extendable` | bool | make capacity a decision variable |
| `max_capacity`, `min_capacity` | MW | bounds when extendable |
| `committable` | bool | binary on/off unit commitment |
| `ramp_up`, `ramp_down` | MW/timestep | ramp limits |
| `min_up_time`, `min_down_time` | timesteps | dwell constraints |
| `startup_cost`, `shutdown_cost` | $ | per transition |
| `must_run` | bool | always on |
| `tech` | str | technology tag, used by policy carve-outs |

### Load

`add_load(name, bus, amount)` — `amount` is a scalar or a per-timestep array.
Load is inelastic by default: it must be served.

### Storage

`add_storage(name, bus, power_capacity, energy_capacity, **kwargs)` — see
[section 7](#7-storage).

### Link

`add_link(name, bus_from, bus_to, capacity, efficiency=1.0, **kwargs)`

A link moves or converts energy between two buses. A transmission line is a
link with `efficiency` slightly below 1; an electrolyser is a link from an
electricity bus to a hydrogen bus with `efficiency` around 0.7; a heat pump is
a link with `efficiency` above 1 (its coefficient of performance).

---

## 6. Reading results

`optimise()` returns an `OptimisationResult`:

| Field | Type | Contents |
|---|---|---|
| `status` | str | `"optimal"`, `"infeasible"`, … — **check this first** |
| `total_cost` | float | objective value |
| `solve_time` | float | seconds |
| `generator_dispatch` | `{name: array}` | MW per timestep |
| `storage_charge` / `storage_discharge` | `{name: array}` | MW per timestep |
| `storage_soc` | `{name: array}` | MWh per timestep |
| `link_flow` | `{name: array}` | MW per timestep |
| `bus_shadow_prices` | `{name: array}` | $/MWh — marginal price at each bus |
| `unit_status` | `{name: array}` | commitment `u[t]`; only when `committable=True` |
| `capacity_additions` | `{name: float}` | MW built; only for extendable components |
| `storage_soc_duals` | `{name: array}` | $/MWh marginal value of stored energy (LP solves, full-mode storages) |
| `cap_dual` | `{name: float}` | marginal operational value of capacity (Benders) |

Inspect the model itself before solving:

```python
print(sys.summary())
print(sys.n_buses, sys.n_components, sys.n_timesteps)   # properties, not methods
```

```
EnergySystem: sum
  Buses: 1
  Generators: 1
  Storages: 0
  Loads: 1
  Links: 0
  Timesteps: 1 (dt=1.0h)
```

---

## 7. Storage

```python
sys.add_storage(
    "battery",
    bus=elec,
    power_capacity=100.0,        # MW  — max charge/discharge rate
    energy_capacity=400.0,       # MWh — 4-hour battery
    efficiency_charge=0.92,
    efficiency_discharge=0.92,
    self_discharge=0.001,        # fraction lost per timestep
    soc_initial=0.5,             # fraction of energy_capacity
    cyclic=True,                 # soc(0) == soc(T)
)
```

Worked example — a battery shifting solar into the evening:

```python
import numpy as np, nexus_energy as ne

sys = ne.EnergySystem("storage_demo")
elec = sys.add_bus("elec")
sys.set_timesteps(8)

sys.add_generator("solar", bus=elec, capacity=400, marginal_cost=0,
                  carrier_factor=np.array([0, 0, .6, 1., 1., .4, 0, 0.]))
sys.add_generator("gas", bus=elec, capacity=400, marginal_cost=80)
sys.add_storage("battery", bus=elec, power_capacity=100.0, energy_capacity=400.0,
                efficiency_charge=0.92, efficiency_discharge=0.92,
                self_discharge=0.001, soc_initial=0.5, cyclic=True)
sys.add_load("demand", bus=elec, amount=np.full(8, 250.0))

r = sys.optimise()
print(r.status, r.total_cost)
print("soc:", np.round(r.storage_soc["battery"], 1))
```

```
optimal 74561.3
soc: [200.   91.1  80.1 172.1 263.9 200.6 200.4 200.2]
```

The battery discharges into the dark early hours, refills through the solar
peak, and returns to its starting level because `cyclic=True`.

### Key parameters

| Parameter | Default | Meaning |
|---|---|---|
| `efficiency_charge` / `efficiency_discharge` | 0.95 | one-way efficiencies |
| `self_discharge` | 0.0 | fraction lost per timestep |
| `soc_min`, `soc_max`, `soc_initial` | 0, 1, 0.5 | fractions of `energy_capacity` |
| `cyclic` | `True` | enforce `soc(0) == soc(T)` |
| `cyclic_level` | `"fixed"` | `"fixed"` also pins the level to `soc_initial`; `"free"` lets the optimiser choose it (PyPSA's convention) |
| `no_simultaneous` | `False` | binary lock preventing simultaneous charge and discharge |
| `marginal_cost` | 0.0 | $/MWh on the discharge leg |
| `extendable` | `False` | make power and energy capacity decision variables |
| `max_hours` / `min_hours` | `None` | bind energy capacity to a duration window |
| `inflow`, `spill_to` | `None` | hydro reservoirs and cascades |
| `pump_capacity` / `turbine_capacity` | `None` | asymmetric rates for pumped hydro |
| `long_duration` | `False` | inter-period SOC carry-over under representative periods |

> **`no_simultaneous` costs one binary per timestep per storage.** You only
> need it when round-trip efficiency is exactly 1 (otherwise losses already
> make simultaneous charge/discharge uneconomic) or when the physical
> constraint is required for its own sake. Leave it off for pure LP dispatch.

### Store vs full formulation

`storage_model` selects the LP variable shape: `"auto"` (default), `"store"` or
`"full"`. The `"store"` model compiles a single energy-state variable `e[t]` per
timestep instead of the three-variable `(charge, discharge, soc)` form — a
smaller polyhedron that helps HiGHS find integer-feasible vertices at the root
in models with committable links.

`"auto"` picks `"store"` only when round-trip efficiency is ≈ 1, charge and
discharge capacities are symmetric, and none of `inflow`, `spill_to`,
`ramp_cost`, `no_simultaneous`, `availability`, `long_duration` or a non-zero
marginal cost is in use. Result extraction derives `charge[t]` / `discharge[t]`
from the signed `Δe/dt`, so downstream code sees the same arrays either way.

---

## 8. Unit commitment

Thermal plants cannot be switched on and off freely. Setting `committable=True`
introduces binary on/off variables and turns the problem into a MILP.

```python
import numpy as np, nexus_energy as ne

sys = ne.EnergySystem("uc")
elec = sys.add_bus("elec")
sys.set_timesteps(8)

sys.add_generator(
    "coal", bus=elec, capacity=100.0,
    p_min=20.0,                # minimum stable output when on
    marginal_cost=30.0,
    startup_cost=5000.0,
    shutdown_cost=1000.0,
    min_up_time=3,             # timesteps
    min_down_time=2,
    ramp_up=40.0,              # MW per timestep
    ramp_down=40.0,
    committable=True,          # <- activates binary commitment
)
sys.add_generator("peaker", bus=elec, capacity=100.0, marginal_cost=200.0)
sys.add_load("demand", bus=elec, amount=np.array([10, 50, 90, 80, 20, 10, 80, 90.]))

r = sys.optimise()
print(r.status, r.total_cost)
print("u[t]:", r.unit_status["coal"])
```

```
optimal 49400.0
u[t]: [ 0.  1.  1.  1. -0.  0.  1.  1.]
```

Coal is off for the first cheap hour, runs for three (satisfying
`min_up_time=3`), shuts down through the trough, and restarts for the evening
peak — paying `startup_cost` twice rather than idling at `p_min` throughout.

> **Parameter names to know.** These differ from some other tools: it is
> `committable` (not `unit_commitment`), `p_min` (not
> `min_stable_generation`), and `ramp_up` / `ramp_down` (not `ramp_limit_up` /
> `ramp_limit_down`).

### Formulation notes

For committable links and non-clustered committable generators, the
`start_up` / `shut_down` indicators are declared **binary** and the state
transition compiles as two inequalities (`v[t] ≥ u[t] − u[t−1]`,
`w[t] ≥ u[t−1] − u[t]`) plus four upper-bound cuts (`v[t] ≤ u[t]`,
`v[t] ≤ 1 − u[t−1]`, `w[t] ≤ 1 − u[t]`, `w[t] ≤ u[t−1]`) — the PyPSA-style
LP-tight 3-bin form. The `v + w ≤ 1` mutex is implied and is not emitted as a
separate row. Clustered generators keep the equality transition form.

**Multi-state (hot / warm / cold) starts.** `start_up_segments` takes a list of
`(min_off_timesteps, startup_cost)` tuples — ascending in `min_off_timesteps`,
non-decreasing in cost, first entry `0`. When set it overrides the flat
`startup_cost` and uses the tight Morales-España (2013) start-type formulation.
Binary UC only; clustered UC falls back to the flat cost.

**Clustered UC.** `clustered=True` with `n_units=N` lumps N identical units
together with continuous `u, v, w ∈ [0, N]` — the GenX `UCommit: 2` analogue.
Far cheaper than N sets of binaries.

**Piecewise heat rates.** `heat_rate_segments` takes
`(p_MW_breakpoint, marginal_cost)` points defining a convex increasing
piecewise-linear fuel cost, overriding the flat `marginal_cost`.

---

## 9. Capacity expansion

Set `extendable=True` and give a `capital_cost`, and capacity becomes a
decision variable.

```python
import numpy as np, nexus_energy as ne

sys = ne.EnergySystem("expand")
elec = sys.add_bus("elec")
sys.set_timesteps(24)
sys.set_snapshot_weights(np.full(24, 365.0))     # one representative day × 365

cf = np.clip(np.sin(np.linspace(0, np.pi, 24)), 0, None)

sys.add_generator("solar", bus=elec, capacity=0.0, marginal_cost=0.0,
                  capital_cost=60_000.0,          # $/MW/YEAR
                  extendable=True, max_capacity=2000.0,
                  carrier_factor=cf)
sys.add_generator("gas", bus=elec, capacity=600.0, marginal_cost=90.0)
sys.add_load("demand", bus=elec, amount=np.full(24, 400.0))

r = sys.optimise()
print(r.status, r.total_cost)
print(r.capacity_additions)
```

```
optimal 110626907.4
{'solar': 769.8}
```

> ⚠ **The single most common beginner mistake: `capital_cost` is `$/MW/year`,
> but your model may only span a few hours.** If you model 6 hours and price
> capex annually, no plant will ever be worth building — the model sees a whole
> year of capital cost recovered by six hours of fuel savings, and builds
> nothing. `set_snapshot_weights()` is the fix: weight each snapshot by the
> number of hours of the year it represents, so the two sides of the trade-off
> are on the same clock. In the example above, 24 snapshots × 365 = 8760 hours.

Related parameters: `min_capacity` (a floor, PyPSA's `p_nom_min`), `fixed_om`
($/MW/year paid on built capacity whether or not it runs), `integer_investment`
with `unit_size` (build in discrete units), and `capex_segments` for piecewise
economies of scale.

Storage is extendable the same way, through `capital_cost_power` ($/MW/year)
and `capital_cost_energy` ($/MWh/year).

---

## 10. Networks and transmission

```python
import numpy as np, nexus_energy as ne

sys = ne.EnergySystem("network")
north = sys.add_bus("north")
south = sys.add_bus("south")
sys.set_timesteps(4)

sys.add_generator("wind", bus=north, capacity=400, marginal_cost=0)
sys.add_generator("gas",  bus=south, capacity=400, marginal_cost=70)
sys.add_load("city", bus=south, amount=np.full(4, 300.0))

sys.add_link("interconnector", bus_from=north, bus_to=south,
             capacity=200.0, efficiency=0.98)

r = sys.optimise()
print(r.status, r.total_cost)
print("flow:", r.link_flow["interconnector"])
```

```
optimal 29400.0
flow: [200. 200. 200. 200.]
```

The line runs at its 200 MW limit in every hour — free northern wind displaces
southern gas up to the point the wire allows, and the remaining 104 MW comes
from gas. A congested interconnector is exactly what you would expect the
shadow prices at the two buses to diverge across.

### Flow formulations

`Link.model_type` selects the physics: `"transport"` (default — a pipe with a
capacity), `"dc_opf"` (linearised DC power flow with voltage angles), `"ptdf"`
(power transfer distribution factors) or `"switched"`.

Internally each link resolves to one of two LP variable shapes:

- **`signed`** — PyPSA-style symmetric bidirectional transport: one
  `f[t] ∈ [−cap, cap]` variable per timestep, no forward/reverse mutex row.
  Halves the per-line variable count.
- **`fwd_rev`** — one non-negative flow variable, plus a reverse variable and a
  `fwd + rev ≤ cap` mutex when bidirectional.

Eligibility for `signed` is deliberately conservative: any feature applying a
coefficient to the *absolute* flow magnitude (losses, marginal cost, ramp cost,
ramp limits, multi-output efficiencies, CO₂ output, UC/switching, linepack)
keeps `fwd_rev` so the LP stays linear. Extendable bidirectional links also
stay `fwd_rev` — signed flow under a capacity *variable* needs a two-sided
coupling, and it measures slower.

The optimum is bit-identical either way. `optimise(link_formulation=…)` pins
the shape for A/B measurement; on a fixed-capacity transport LP, `signed`
measured −19 % variables and constraints and −17 % wall clock.

Other network capabilities: `set_n_minus_1()` for N-1 security,
`set_contingency_reserve()` for single-largest-unit reserve, HVDC via link
parameters, and the AC formulations in [section 20](#20-ac-power-flow).

---

## 11. Policy constraints

### Carbon

```python
sys.set_co2_price(80.0)            # $/tCO2 added to marginal costs
sys.set_emission_limit(600.0)      # tCO2 cap over the horizon
sys.set_co2_rate_cap(...)          # tCO2/MWh intensity cap
sys.set_co2_zone_cap(...)          # per-zone caps
```

Both levers on the same two-generator system:

```
with price 80/tCO2 : total 104400.0   coal 0.0     gas 1200.0
with a 600 t cap   : total  58800.0   coal 240.0   gas  960.0
```

A price of 80 makes coal (0.9 tCO₂/MWh) strictly worse than gas and pushes it
out entirely; a quantity cap lets coal run right up to the limit. Same physical
system, different instrument, different answer — which is the point of modelling
both.

**Pooled multi-zone cap groups** are one constraint over a *set* of buses (the
GenX `Cap_Zone` analogue), which is tighter than independent per-bus caps
because it forbids inter-zone emission averaging:

```python
sys.set_co2_cap_group(
    [bus_a, bus_b, bus_c], limit=0.05,   # tCO2/MWh
    is_rate=True,
    storage_losses_on_rhs=True,          # GenX CO2Cap=2 RHS term
    loss_accounting="net",               # "net" | "dissipation"
)
```

### Clean energy shares

```python
sys.set_rps(0.4, qualifying_techs=["wind"])       # 40% from tagged techs
sys.set_ces(0.6, scores={"wind": 1.0, "gas": 0.4})  # weighted clean standard
```

Both accept `slack_penalty=` to make the target a priced soft constraint rather
than a hard one. `qualifying_techs` matches against each generator's `tech`
tag, so remember to set `tech=` when you add the generator.

### Other instruments

`set_itc()` and `set_ptc()` (investment and production tax credits),
`set_hourly_matching()` (24/7 carbon-free matching), `set_capacity_bucket()`
(technology carve-outs), `set_reserve_margin()`, `set_spinning_reserve()`,
`set_regulation_reserve()`, `set_fuel_supply_limit()`, `set_outage()`,
`set_shared_capacity()`.

---

## 12. Coming from PyPSA

`from_pypsa(network)` converts a PyPSA `Network` into an `EnergySystem`.

```python
import pypsa
from nexus_energy.pypsa_compat import from_pypsa

n = pypsa.Network()
n.set_snapshots(range(4))
n.add("Bus", "elec")
n.add("Generator", "solar", bus="elec", p_nom=400, marginal_cost=0,
      p_max_pu=[0.0, 0.4, 0.9, 0.2])
n.add("Generator", "gas", bus="elec", p_nom=400, marginal_cost=60)
n.add("Load", "demand", bus="elec", p_set=[200, 250, 300, 280])

system = from_pypsa(n)                          # AC lines auto-route to DC-OPF
system = from_pypsa(n, line_model="transport")  # required by the diff layer

r = system.optimise()
print(r.status, r.total_cost)
print(r.generator_dispatch)
```

```
optimal 29400.0
{'solar': array([  0., 160., 300.,  80.]), 'gas': array([200.,  90.,   0., 200.])}
```

This is a **one-way adapter** that reads the network's dataframes. The library
does not depend on PyPSA and never calls it to solve. It exists so a benchmark
can put the identical network in front of both solvers, and so an existing
PyPSA workflow has an on-ramp.

### Two parity-critical conventions

- **Static (scalar) `p_max_pu` is honoured.** When a generator has no
  `p_max_pu` time series, a constant capacity-factor array is synthesised from
  the static column (for example PyPSA-Eur's nuclear at `p_max_pu = 0.781`), so
  `p[t] ≤ capacity · p_max_pu` is enforced exactly as PyPSA does. Dropping
  constant de-ratings used to over-dispatch baseload and under-price the
  system. Exports round-trip the factor back into `generators_t.p_max_pu`.
- **Cyclic storage level is free.** StorageUnits and Stores import with
  `cyclic_level="free"` — only continuity `soc(0) = soc(T)` is enforced and the
  optimiser picks the level, matching PyPSA's `cyclic_state_of_charge` /
  `e_cyclic`. A fixed-level pin over-constrains, and for extendable storage
  (whose energy capacity starts at 0) it forces start and end empty.

Both of these were found *because* the inputs were held identical across the
two solvers.

---

## 13. Multi-stage planning

`MultiStageSystem` chains per-year `EnergySystem` snapshots into one
investment problem with vintage tracking.

```python
import numpy as np, nexus_energy as ne

def stage(year, demand, hours=24):
    s = ne.EnergySystem(f"y{year}")
    b = s.add_bus("elec")
    s.set_timesteps(hours)
    cf = np.clip(np.sin(np.linspace(0, np.pi, hours)), 0, None)
    s.add_generator("solar", bus=b, capacity=0.0, marginal_cost=0.0,
                    capital_cost=500.0,        # see the note below
                    extendable=True, max_capacity=1000.0, carrier_factor=cf)
    s.add_generator("gas", bus=b, capacity=800.0, marginal_cost=95.0)
    s.add_load("demand", bus=b, amount=np.full(hours, demand))
    return s

ms = ne.MultiStageSystem("plan")
ms.add_stage(2030, stage(2030, 300.0))
ms.add_stage(2040, stage(2040, 420.0))

res = ms.optimise()
print(res.status, res.total_cost, res.years)
print("new_builds:     ", res.new_builds)
print("capacity_active:", res.capacity_active)
```

```
optimal 690637.1 [2030, 2040]
new_builds:      {'solar': [665.5, 0.0], 'gas': [0.0, 0.0]}
capacity_active: {'solar': [665.5, 665.5], 'gas': [800.0, 800.0]}
```

Solar is built once in 2030 and is still active in 2040 — that is the vintage
tracking. `MultiStageResult` also carries `stage_dispatch`, `stage_link_flow`,
`stage_storage_soc`, `storage_new_power` / `storage_new_energy` and the link
equivalents.

> **`MultiStageSystem` costs stage dispatch snapshot by snapshot and does not
> apply `set_snapshot_weights`.** Express `capital_cost` on the same basis as
> the modelled horizon, or scale the horizon to a year. The `capital_cost=500`
> above is per MW over the modelled 24 hours, not per year.

Vintage-related generator parameters: `build_year`, `lifetime_years`,
`retire_at_year`, `build_lead_years`, and `retrofit_of` for fuel switching
bounded by retiring host capacity.

Pass `myopic=True` to `optimise()` to solve each stage in sequence with no
foresight, which is the standard comparison against perfect-foresight planning.

---

## 14. Differentiable dispatch and inverse calibration

This is the part that does not exist elsewhere.

The dispatch problem is set up as a **differentiable layer**: analytic
`∂dispatch/∂parameter` Jacobians obtained by implicit differentiation of the
KKT system at the optimum. No sampling, no finite differences.

**Honest scope up front:** this is parameter learning on the inner **LP/QP**.
Every layer uses a mandatory strict-convexity ridge that shifts dispatch versus
the true LP by `O(ridge)` — disclose it in results. It is pure numpy; torch is
optional.

### Recovering a hidden CO₂ price

The headline use: you observe a system's dispatch and want to know what carbon
price explains it.

```python
from nexus_energy.pypsa_compat import from_pypsa
from nexus_energy.diff_bridge import fit_co2_price

system = from_pypsa(network, line_model="transport")
fit = fit_co2_price(system, observed_dispatch, price_bounds=(0.0, 1000.0))

print(fit.price, fit.n_solves, fit.converged)
```

`fit_co2_price` runs in two stages: a coarse **forward-only bracket** — dispatch
is piecewise-linear in price, so the loss has flat pieces that a pure gradient
method is born stuck on — followed by **safeguarded Gauss-Newton** with a
bisection fallback. It returns a `CO2FitResult(price, history, n_solves,
converged)`.

The price enters through `mc_eff = mc + price · emission`, so by the chain rule
`∂dispatch/∂price = d_dispatch_d_mc @ emission`, exposed directly as
`d_dispatch_d_co2_price`.

To inspect the bridge rather than fit through it:

```python
from nexus_energy.diff_bridge import multibus_problem_from_system

bridge = multibus_problem_from_system(system, co2_price=80.0, ridge=1e-2)
# bridge.problem is ready for solve_multibus_dispatch_with_sensitivities;
# bridge.emission and bridge.mc_base satisfy mc = mc_base + price * emission
```

**The bridge fails loudly rather than silently approximating.** It raises on:
extendable components (every capacity must be fixed), DC-OPF links (re-import
with `line_model="transport"`), storage, lossy or costed or multi-output links,
and any of UC / `p_min` / `must_run` / piecewise heat rates.

### Layer inventory

| Layer | Scope | Differentiable inputs |
|---|---|---|
| `solve_dispatch_with_sensitivities` | single-bus, single-period, closed form | mc, capacity, demand |
| `EconomicDispatchLayer` | stateful wrapper of the above | same |
| `MultiBusDispatchProblem` + `solve_multibus_dispatch_with_sensitivities` | multi-bus, multi-period transport QP | mc, capacity, demand, line_limit |
| `StorageDispatchProblem` + `solve_storage_dispatch_with_sensitivities` | single-bus multi-period + one storage | mc, capacity, demand, soc_init, η_c, η_d |
| `MultiBusStorageProblem` + `solve_multibus_storage_dispatch_with_sensitivities` | multi-bus + bus-attached storages, one stacked QP | mc, η_c, η_d, soc_init |
| `SmoothCommitmentLayer` / `fit_commitment_threshold` | sigmoid-smoothed commitment | threshold |
| `CapacityExpansionProblem` / `CapacityExpansionLayer` / `fit_component_params` | differentiable capacity expansion (design gradients) | capex, mc, demand, … |
| `TorchDispatchLayer` | cvxpylayers hook for problems beyond the numpy paths | — |

### Selective Jacobians

Jacobian blocks are expensive. Request only what you will use:

```python
sol = solve_multibus_dispatch_with_sensitivities(
    problem,
    jacobians=("mc",),      # () means forward-only
)
# sol.dispatch (G,T), sol.flows (L,T)
# sol.d_dispatch_d_mc (G·T, G), d_dispatch_d_capacity (G·T, G),
# d_dispatch_d_demand (G·T, B·T), d_dispatch_d_linelimit (G·T, L)
```

Skipped blocks come back as zeros. Pass `()` for the price a sampling baseline
pays per draw; pass `("mc",)` when only cost-side gradients are chained, as in
CO₂-price calibration.

`MultiBusDispatchProblem` also carries `availability` — `(G, T)` in `[0,1]`,
giving per-period derating `p[g,t] ≤ capacity[g] · availability[g,t]` for VRE
capacity factors and outages — and `line_min`, `(L,)` lower flow bounds, where
`None` keeps symmetric `−line_limit` and `0.0` entries make links
unidirectional.

### Storage efficiency Jacobians

`StorageDispatchSolution` exposes the constraint-matrix sensitivities
`d_dispatch_d_charge_eff` and `d_dispatch_d_discharge_eff`, which need the
equality duals recovered from the frozen-active-set KKT. These are what an
auto-calibration loop fits from battery telemetry.

**Dispatch-only telemetry identifies only the η_c · η_d product.** The SOC-trace
rows `d_soc_d_charge_eff` / `d_soc_d_discharge_eff` split charge from discharge
efficiency once observed SOC (BMS logs) is added to the residual. This is a
concrete example of the identifiability question the library is built to answer
honestly.

> `MultiBusStorageProblem` is solved as **one stacked QP over all periods**
> because SOC couples them. Keep the window small — representative days or
> weeks, T ≤ ~48. The dense stacked KKT is `O(((G+L+3S)·T)³)`.

---

## 15. MPC and auto-calibration

### Warm-started rolling horizons

`PersistentDispatchSession` does rolling-horizon LP re-solves **without
rebuilding the model**. `build()` solves once through the normal `optimise()`
path while capturing the assembled model, then loads it into a persistent HiGHS
instance; `advance()` pushes parameter changes straight into HiGHS columns and
rows and re-solves from the retained simplex basis.

```python
import numpy as np
from nexus_energy.mpc import PersistentDispatchSession

sess = PersistentDispatchSession(system)
base = sess.build()

res = sess.advance(
    demand={"elec": np.full(6, 350.0)},   # bus -> (T,) TOTAL demand
    cf={"wind": new_cf},                  # generator -> (T,) availability in [0,1]
    mc={"gas": 75.0},                     # generator -> scalar or (T,) marginal cost
    soc_init={"battery": 120.0},          # storage -> start SOC in MWh
)
print(sess.n_resolves, sess.n_rebuilds, sess.last_iterations)
```

A window re-solve costs a handful of simplex iterations instead of a rebuild
plus a cold solve.

**Honest scope** — `build()` raises otherwise: pure LP dispatch, so no
committable, integer-investment, switchable or piecewise-capex components;
fixed capacities (an extendable generator puts availability into a
`p ≤ cap_var · cf` matrix row this path cannot touch); and storage start-SOC
carry-over requires `soc_initial_free=True`, which makes the start SOC a
pinnable column. Anything out of scope at `advance()` time falls back to a full
rebuild — correct, just slower — and is counted in `n_rebuilds`.

### Keeping the model honest from telemetry

`autocal.fit_params` generalises the scalar CO₂-price recovery to an
m-dimensional parameter vector with a damped (scaled Levenberg-Marquardt)
Gauss-Newton driver. `AutoCalibrator` wraps it into a moving-horizon tracker so
a running MPC keeps correcting its own model.

```python
from nexus_energy.autocal import fit_params, AutoCalibrator

report = fit_params(
    make_solution,          # theta-dict -> SOLVED diff-layer solution
    observed,               # telemetry, flattened like the residual
    params={"eta_c": (0.95, 0.5, 1.0), "mc_gas": (50.0, 10.0, 200.0)},
    jacobian_fn=jac,        # (solution, theta) -> (len(r), m) analytic columns
    gate_threshold=1e-3,    # identifiability gate
    max_rel_step=0.25,      # per-run trust region
)
print(report.changed(), report.frozen)    # frozen[name] is True => data-silent

cal = AutoCalibrator(solution_factory, jac, observed_fn, params,
                     update_every=4, noise_std=0.5, smooth=0.3)
rep = cal.step(window)          # None on skipped cycles
cal.lock_param("eta_c")         # operator override
cal.set_param("mc_gas", 55.0, lock=True)
print(cal.believed)             # current point beliefs
```

Two design rules make this trustworthy in operations:

- **The identifiability gate** freezes any parameter whose scale-aware Jacobian
  column carries no signal in the window. It is flagged `data_silent` and never
  nudged by noise. This is the "it says so when the data is silent" property,
  made operational.
- **The per-cycle slew limit** (`max_rel_step`) bounds how far beliefs can move
  in one cycle, so a single bad window cannot wreck the model.

Noise-aware controls: `noise_std` enables EMA smoothing across windows
(`smooth` gain) and outright **rejection of outlier windows** whose post-fit
loss exceeds `outlier_zscore² · (½nσ²)` — a sensor fault or unmodelled event
leaves beliefs untouched. `residual_fn` switches the fitted observable between
dispatch shares, SOC traces and so on.

Any parameter chain-rulable from the exposed Jacobian blocks works: marginal
costs, fuel price and efficiency through `mc = fuel/η`, CO₂ price through
emissions, storage efficiencies through the efficiency blocks.

---

## 16. Solver controls

`EnergySystem.optimise()` exposes the full solver-control surface. Everything
here is accuracy-preserving.

| Kwarg | Default | Meaning |
|---|---|---|
| `lp_backend` | `None` | `auto` / `simplex` / `ipm` / `ipm_fast` / `pdlp` / `gpu`. `None` resolves: explicit argument > `.nexus_solver.json` sidecar > built-in default `"ipm_fast"`. Pure-LP only; MILP is never forced onto IPM. |
| `solver_method` | `None` | raw HiGHS `solver` option; always wins over `lp_backend`. `"ipm"` auto-forces `run_crossover="on"` so the result stays a bit-exact vertex. |
| `run_crossover` | `None` | `"on"` / `"off"` / `"choose"` — crossover recovers an exact vertex and duals |
| `parallel` | `None` | HiGHS parallel `"on"` / `"off"` / `"choose"`; auto-set `"on"` when `threads > 1` |
| `threads` | `None` | worker threads |
| `time_limit`, `gap` | `None` | stopping criteria |
| `scale_cleanup` | `True` | snap finite bounds/RHS with `0 < abs(v) < 1e-9` to exactly 0 |
| `simplex_scale_strategy` | `None` | HiGHS matrix-equilibration strategy, int 0–5 |
| `eliminate_redundant` | `True` | exact pre-solve pass dropping provably non-binding rows |
| `link_formulation` | `"auto"` | see [section 10](#10-networks-and-transmission) |
| `ramp_cost_formulation` | `"split"` | `"signed"` replaces each ramp-cost up/down aux pair with one aux variable `r ≥ ±Δ`; both price `abs(Δ)` identically at the optimum |
| `mip_strategy` | `"auto"` | see below |
| `uc_fix_schedule` | `None` | `{name: u-array}` pins `u[t]` for committable generators and links |
| `warm_start` | `None` | previous result forwarded to HiGHS `setSolution()` |
| `basis` | `None` | simplex basis hot-start |
| `myopic` | — | `MultiStageSystem.optimise()` only: no inter-stage foresight |

### `mip_strategy`

- **`"auto"`** (default) — for MIPs with ≤ 5000 estimated binaries, solve the LP
  relaxation first; if every binary lands within `1e-4` of {0,1} then the LP
  optimum *is* the MIP optimum and is returned directly (one LP, zero gap).
  Otherwise fall through to HiGHS MIP.
- **`"mip_only"`** — go straight to the MIP solver.
- **`"lp_first"`** — force the LP-first path at any size, on a vertex-producing
  backend (an interior point would spuriously look fractional). Exact by the
  LP-relaxation theorem; the LP cost is paid twice when fractional, hence
  opt-in.
- **`"fix_and_certify"`** — `lp_first` plus: on a fractional vertex, fix every
  `u[t]` that *is* integral there, solve the residual MIP, and return it when
  its cost is within `gap` of the relaxation bound — a valid optimality
  certificate. Falls back to the full MIP otherwise.

> **`ipm_fast` caveat.** It returns an interior (non-vertex) point, so
> degenerate duals and dispatch can differ from simplex, and basis warm-start
> is a no-op. Use `lp_backend="simplex"` when vertex duals or warm-start
> matter. A non-optimal `ipm_fast` result auto-falls back to simplex.

### Auto-tuning the LP method

Static shape heuristics proved unreliable — an early "≥ 50k columns ⇒ IPM" rule
regressed a real benchmark — so the tuner is **empirical**. It races `simplex`,
`ipm` and `ipm_fast` on reduced-horizon proxies of your system, checks the
winner is stable across two probe sizes (a clean simplex→IPM size-crossover
counts as stable), checks objective parity within each size, and writes a
`.nexus_solver.json` sidecar that `optimise()` then reads automatically.

```python
from nexus_energy.solver_tuner import tune_solver, recommend_lp_method

res = tune_solver(system)              # races, prints a report, writes the sidecar
print(res.recommended)                 # e.g. "ipm_fast"

# report-only / finer control:
res = recommend_lp_method(system, time_cap=120.0, threads=8,
                          need_duals=False, probe_hours=[365, 730],
                          verbose=True)
```

Library-agnostic — it works on any `from_pypsa` import. MILP systems
short-circuit to `simplex`, since IPM cannot solve branch-and-bound node LPs.
`need_duals=True` disqualifies `ipm_fast`. The sidecar is advisory: an explicit
`lp_backend=` always wins, and deleting the file reverts. Validated:
PyPSA-Eur → `ipm_fast`, CINDER → `simplex`.

---

## 17. Temporal aggregation

A full 8760-hour planning model is heavy. Cluster the year into representative
days instead.

```python
import nexus_energy as ne

rep = ne.aggregate_to_representative_days(
    timeseries={"demand": demand_8760, "solar_cf": solar_8760},
    n_days=12,
    hours_per_day=24,
    seed=42,
    extreme_periods=[("demand", "max"), ("solar_cf", "min")],
)

result = ne.apply_representative_days(
    system, rep,
    timeseries_map={"demand": "demand", "solar": "solar_cf"},
)
```

`extreme_periods` is important and easy to skip: pure k-medoids clustering
discards the peak-load and dark-doldrums days that determine system adequacy.
Naming them explicitly keeps them in the reduced set.

`aggregate_with_feature_embedding(timeseries, n_days=…, features=…)` embeds
load, capacity factors and ramp statistics in a feature space before running
k-medoids, which keeps the day count low while preserving the correlation
structure between wind and solar dropouts.

**Certified reduction bounds.** `certify_reduction` returns a `CertifiedBound`
— a provable two-sided bound on the cost error a temporal aggregation
introduces, turning "the representative days look close" into a certificate.
See `certified_reduction_demo` for a worked example.

**Variable resolution.** `ResolutionPlan`, `adaptive_resolution_plan`,
`multi_resolution_hierarchy` and `apply_adaptive_resolution` build non-uniform
timestep plans; `set_snapshot_durations()` applies per-snapshot hour lengths
directly, so a coarse merged block correctly moves `power × duration` of energy.

**Long-duration storage.** Set `long_duration=True` on a storage and, when the
system has representative periods plus a chronological mapping,
SOC = `soc_intra[t]` + `soc_inter[original_day]`, tracking carry-over across
the real calendar (Kotzur 2018 inter-period superposition).

---

## 18. Decomposition at scale

### Certified temporal decomposition

`optimise_temporal_certified` solves a long-horizon system as K temporal blocks
**with a global optimality certificate** — the same `(UB − LB)/UB ≤ gap`
contract a monolithic MIP solver gives, at an order of magnitude less wall time.
MILP wall clock scales roughly `T^2.5`, so cutting T is worth a great deal.

```python
from nexus_energy.temporal_certified import optimise_temporal_certified

res = optimise_temporal_certified(
    system_factory,            # callable(t0, t1) -> EnergySystem over [t0, t1)
    total_steps=8760,
    n_blocks=8,
    gap=1e-2,
    lb_blocks=4,               # fewer, longer LB blocks => tighter bound
    lb_rounds=3,               # subgradient ascent rounds on boundary prices
    lb_workers=4,              # process-parallel LB blocks (factory must pickle)
    ub_boundary="both",        # "prices" | "floors" | "both" | "none"
    reachability_envelopes=True,
)
print(res)   # TemporalCertifiedResult(status='certified', objective=…, gap=…)
```

- **Lower bound** — Lagrangian dual decomposition over interior SOC boundaries.
  Each block is a *relaxation* (free interior start SOC, telescoping ±λ boundary
  prices), so the sum of block dual bounds is ≤ the full optimum for arbitrary
  λ. `lb_rounds` refines λ by projected subgradient and the best round is
  reported.
- **Upper bound** — sequential block stitching, handing each block its
  predecessor's terminal SOC. `ub_boundary="prices"` prices terminal energy at
  −λ and subtracts the payment back out, so the UB is the stitched trajectory's
  *true* full-model cost.
- **Reachability envelopes** cap each LB block's free start and terminal SOC by
  what the full problem could physically have stored by that boundary — pure
  bound tightening, never validity-affecting.

Supported scope, guarded with `ValueError`: non-cyclic, non-LDS, non-extendable
storages; committable links with min up/down ≤ 1; no hard ramp limits; no
committable generators; uniform snapshot weights and durations.

Measured on the CINDER MILP: +0.91 % from parity at 58 s (6.1× faster than
monolithic), and beating the incumbent at 87 s.

### Benders

For networks too large for one machine, `BendersDecomposer` splits investment
(master) from dispatch (sub-problems) and solves the sub-problems in parallel.

```python
decomposer = ne.BendersDecomposer(
    system=sys,
    periods=[(0, 2190), (2190, 4380), (4380, 6570), (6570, 8760)],
    stabilisation="level",     # note the British spelling
    max_iter=30,
    tol=1e-3,
    n_jobs=4,
)
result = decomposer.solve()
```

Stabilisation prevents the sub-gradient oscillation that makes textbook Benders
crawl. Other drivers: `solve_with_temporal_benders`,
`solve_with_spatial_benders`, `solve_with_nested_benders`,
`solve_with_dantzig_wolfe`, `solve_with_column_generation`, and
`recommend_decomposition` to choose among them.

`rolling_horizon_solve(system_factory, total_timesteps, window_size, overlap)`
threads the simplex basis between windows for hot-started re-solves.

> Note: a `StageProblem` class exists in both `decomposition` and `stochastic`.
> The top-level aliases are `NestedStageProblem` and `SDDiPStageProblem`.

---

## 19. Uncertainty — stochastic and robust

### Two-stage stochastic programming

```python
scenarios = [
    ne.Scenario(name="low",  probability=0.25, demand_factor=0.8),
    ne.Scenario(name="base", probability=0.50, demand_factor=1.0),
    ne.Scenario(name="high", probability=0.25, demand_factor=1.3),
]

result = ne.solve_stochastic(sys, scenarios,
                             risk_measure="expected",   # or "cvar"
                             cvar_alpha=0.05,
                             method="benders")
```

`Scenario` fields are `name`, `probability`, `demand_factor`,
`carrier_factor_scale`, `fuel_cost_factor` and a free-form `overrides` dict.

### Robust optimisation

```python
uncertainty = ne.BudgetUncertaintySet(
    demand_up=0.15,      # demand can be 15% higher
    cf_down=0.30,        # capacity factors can be 30% lower
    fuel_cost_up=0.20,
    budget=2.0,          # Bertsimas-Sim uncertainty budget
)
result = ne.solve_robust(sys, uncertainty)
```

The `budget` is how many deviations may occur simultaneously — the classic
Bertsimas-Sim Γ. At `budget=0` you recover the deterministic problem; at the
full dimension you get the box-worst case.

### Chance constraints

```python
cc = ne.ChanceConstraint(name="reserve", alpha=0.05, threshold=0.0)
result = ne.solve_saa_chance_constrained(sys, scenarios, alpha=0.05,
                                         reserve_margin=0.15)
```

`alpha` is the violation probability, so `alpha=0.05` is a 95 % reliability
level.

### The extended toolkit

`solve_sddip` (multi-stage SDDiP), `solve_general_chance_constrained`,
`solve_wasserstein_dro` (distributionally robust over Wasserstein balls),
`solve_risk_averse_benders`, `generate_forced_outage_scenarios`,
`generate_demand_scenarios`, `generate_renewable_scenarios`,
`generate_moment_matching_scenarios`, `reduce_scenarios` and
`reduce_scenarios_wasserstein` (optimal-transport scenario reduction).

---

## 20. AC power flow

Linear transport models are insufficient for transmission studies that care
about voltage and reactive power.

### SOCP relaxation

The second-order cone (Jabr) relaxation is convex, so its optimum is a global
bound on the true AC optimum.

```python
res = ne.solve_socp_opf(
    system,
    snapshot=0,
    enable_obbt=True,        # optimisation-based bound tightening
    obbt_iters=3,
    obbt_tol=1e-4,
    cos_envelope_pieces=8,
    enforce_cycle_closure=False,
    enforce_tight_qc=False,
)
```

OBBT iteratively shrinks bounds on voltage magnitudes and angle differences,
tightening the relaxation on networks where plain Jabr is loose.

### Polar AC-OPF

```python
res = ne.solve_ac_opf_polar(system, snapshot=0, slack_bus="bus1")
```

The true non-convex formulation. It reflects real reactive power, at the cost
of being susceptible to local optima.

| Use | Choose |
|---|---|
| transmission **planning** | SOCP + OBBT — convex, gives a defensible envelope |
| operations **analysis** | polar AC-OPF — real non-linear physics |

Conic extensions: `solve_socp_opf_expansion` (capacity expansion on the Jabr
relaxation), `solve_socp_opf_multi` (multi-snapshot), `add_weymouth_pipe` (gas
network Weymouth relaxation) and `add_head_dependent_hydro`.

---

## 21. ML-guided solving

To accelerate rolling-horizon MPC, learned predictors warm-start the UC
binaries so the MILP does not start cold.

```python
sys_feat = ne.extract_system_features(system)
t_feat   = ne.extract_timestep_features(system)

prediction = ne.predict_unit_commitment(...)
schedule   = ne.warm_start_from_prediction(prediction,
                                           confidence_threshold=0.7,
                                           cold_start_fallback=True)
result = system.optimise(uc_fix_schedule=schedule)
```

Predictors: `MeritOrderPredictor` (cheap heuristic baseline),
`HistoricalNeighborPredictor` (k-NN over past solves), `GNNPredictor` (graph
neural network, torch optional), all behind `UCWarmstartPredictor`.

`confidence_threshold` is the safety dial — only predictions above it are
pinned, and `AdaptiveThresholdController` / `solve_with_adaptive_warmstart`
tune it automatically from observed hit rate. Also available: `LearnedVarFixer`
/ `apply_varfix`, `RLVarFixer` / `solve_with_rl_search`, and
`learned_representative_periods` / `feature_embedding_periods`.

---

## 22. The component library

223 components across 15 sectors, each available at several **fidelity** levels
— fidelity being how much physics the component carries. The same LFP cell is
available as:

| Level | Model | What it adds |
|---|---|---|
| `F0a` | round-trip efficiency curve | one lookup versus C-rate; cheapest |
| `F1a` | state of charge | energy in, energy out, SoC over time |
| `F1b` | + thermal | cell temperature moves the efficiency |
| `F1c` | + degradation | capacity fade over cycles |
| `F2a` | equivalent circuit (1RC) | real voltage dynamics |
| `F2d` | single-particle model | electrochemical detail |

**Pick the cheapest level that still answers the question — per component, not
per study.** A capacity screen can run everything at F0 while the one asset
under investigation runs at F2. That is the point of the axis.

Sector split: thermal 41 · power electronics 32 · batteries 26 · solar 18 ·
hydrogen 17 · conventional 13 · hydro and marine 12 · biomass 11 · carbon
capture 11 · gas systems 9 · thermoelectric 8 · desalination 7 · wind 6 ·
mechanical storage 6 · geothermal 6.

F0–F2 are built for all 223. **F3–F6 (distributed physics, AI surrogates,
PINNs) are scaffolded but not yet implemented.**

Access templates through the registry, and compose them with carrier checking:

```python
ne.registry.list_components()
template = ne.registry.get("lfp_cell_f1a")

from nexus_energy.components.composition import Subsystem, FIDELITY_LEVELS
```

`Subsystem` raises `CarrierMismatchError` when you wire a hydrogen output into
an electricity input.

---

## 23. Benchmarks

Full detail in [`COMPARISON_SCORECARD.md`](COMPARISON_SCORECARD.md); stored
results under [`benchmarks/results/`](benchmarks/results), so every number can
be re-run rather than taken on trust.

| Case | Objective vs reference | Wall clock |
|---|---|---|
| PyPSA-Eur capacity expansion (10 bus, 2190 h, real profiles) | **−0.000 %** exact parity | 233 s vs 810 s — **3.47× faster** |
| GenX `1_three_zones_ucommit2` | **−0.02 %** exact | 7.7 s vs 36.1 s — **4.7× faster** |
| pandapower AC-OPF case9 / case14 | +0.0007 % / +0.0792 % | 60.5× / 42.2× |
| PowerModels.jl SOCWR, 3-bus radial | 3.67e-5 | 7.57× |
| CINDER LP | parity | 147 s vs 190 s — 1.3× faster |
| CINDER MILP | MIP gap 0.82 % | 330 s vs 190 s — **1.7× slower** |

Both solvers see the identical network — that is what `from_pypsa` is for. Two
real bugs in this library were found *because* the inputs were held identical
(a static scalar `p_max_pu` being ignored, and cyclic storage being
over-pinned).

### Where the speed comes from

1. **Rust constraint assembly.** Vectorised Python is the modelling speed limit
   in comparable libraries. Sparse CSC constraints are streamed to the solver
   from a PyO3 Rust kernel instead.
2. **Tighter-by-default UC formulations.** LP-tight 3-bin commitment reduces
   branch-and-bound work with no user-visible API change.
3. **Empirical LP-method selection.** The `solver_tuner` races simplex and IPM
   variants on reduced-horizon proxies of *your* problem and locks the winner
   in a sidecar; `ipm_fast` is the lean default for well-conditioned expansion
   LPs, while staircase MILPs keep warm-startable simplex.
4. **Feature-guided clustering.** Embedding load, capacity factors and ramp
   statistics before k-medoids keeps representative-day counts low while
   preserving the extreme periods that determine adequacy.
5. **In-place re-solves.** `PersistentDispatchSession` replaces rebuild plus
   cold solve with a parameter update and a hot-started simplex re-solve.

### An open problem, stated as a bug

Two GenX cases (`rate_co2` at −42.7 %, `mincapreq` at −3.90 %) look like large
wins and are **not**. Feeding the capacity nexus chose back into a real GenX
solve showed GenX's own cost for the nexus solution (5.823e9) essentially
matches GenX's optimum (5.808e9), while nexus *reports* 5.582e9 — a **~4.3 %
OPEX under-count bug here**, not a cheaper optimum. Lead suspect is
transmission-loss modelling (GenX piecewise-linear versus a linear `loss=%`).

The pattern: exact match when the build is determined, opex under-count when
multi-zone renewables are dispatched over transmission. The exact-parity rows
above are unaffected — none of them involve that case.

---

## 24. Honest scope and known limits

- **Not dynamic or EMT simulation.** No transients, no swing equation. Time is
  snapshots and rolling horizons coupled by algebraic constraints — a different
  category from Simulink and Modelica.
- **No integer (UC-MILP) differentiability.** Future work, not a claim.
- **Calibration is LP/QP class.** The small ridge term required to make the
  gradients well defined shifts economics by roughly 1–7 percentage points and
  is disclosed per result.
- **The speed headline is forward-only.** The calibration solve uses a denser
  path.
- **F3–F6 component fidelities are scaffolded, not implemented.**
- **`lp_backend="gpu"` requires cuOpt and a CUDA device**, and falls back to
  the CPU path with a printed warning when either is missing. It is not the
  source of the measured speed numbers above.
- **A first `optimise(threads=N)` call in a fresh process can return
  `status="unknown"`** — this comes from the HiGHS global-scheduler reset in
  the solver core. Check `status` before reading `total_cost`.

Third-party solvers and frameworks appear only as benchmark comparison rows.
They are never wrapped inside the library. `external_solvers` provides an
LP-export bridge towards Gurobi / CPLEX / SCIP / Mosek / Xpress for comparison
runs; `optimise()` itself rejects external solver names with a pointer to it.

---

## 25. Troubleshooting

**`status` is `"infeasible"`.** Demand cannot be met. Check that generation
capacity times availability covers peak load at every timestep, that links have
enough capacity to reach the load, and that a hard policy constraint
(`set_emission_limit`, `set_rps`) is not impossible to satisfy. Policy setters
accept `slack_penalty=` to convert a hard constraint into a priced soft one,
which turns an infeasibility into a diagnosable cost.

**Nothing gets built in a capacity-expansion model.** `capital_cost` is
`$/MW/year`. Weight your snapshots — see the warning in
[section 9](#9-capacity-expansion).

**`total_cost` is 0.0.** Everything was served by zero-marginal-cost
generation. That is usually correct, not a bug.

**Dispatch shows tiny negative numbers** like `-1e-11`. Solver tolerance.
Round before display.

**Shadow prices look wrong or degenerate.** `ipm_fast` (the default LP backend)
returns an interior point. Use `lp_backend="simplex"` when you need vertex
duals.

**A MILP is slow.** Try `mip_strategy="lp_first"`, cut the horizon with
representative periods ([section 17](#17-temporal-aggregation)), use
`clustered=True` for identical thermal units, or decompose
([section 18](#18-decomposition-at-scale)).

**The differentiable bridge raises.** By design — it never silently
approximates. Re-import with `from_pypsa(n, line_model="transport")`, fix all
capacities, and remove storage, UC and lossy links. See
[section 14](#14-differentiable-dispatch-and-inverse-calibration).

**A calibrated parameter never moves.** It was frozen by the identifiability
gate: the data in that window carries no signal about it. Check
`report.frozen`. This is the library working as intended, not a failure.

---

## 26. API index

### EnergySystem

**Building:** `add_bus` · `add_carrier` · `add_generator` · `add_load` ·
`add_storage` · `add_link`

**Time:** `set_timesteps` · `set_snapshot_weights` · `set_snapshot_durations` ·
`set_chronological_mapping`

**Policy:** `set_co2_price` · `set_emission_limit` · `set_co2_rate_cap` ·
`set_co2_zone_cap` · `set_co2_cap_group` · `set_rps` · `set_ces` · `set_itc` ·
`set_ptc` · `set_hourly_matching` · `set_capacity_bucket` ·
`set_fuel_supply_limit`

**Reliability:** `set_n_minus_1` · `set_reserve_margin` ·
`set_contingency_reserve` · `set_spinning_reserve` · `set_regulation_reserve` ·
`set_outage` · `set_shared_capacity`

**Solving and inspection:** `optimise` · `summary` · `n_buses` ·
`n_components` · `n_timesteps` *(the last three are properties)*

### Planning

`MultiStageSystem` (`add_stage`, `optimise`) · `MultiStageResult` · `annuity`

### Temporal

`aggregate_to_representative_days` · `aggregate_with_feature_embedding` ·
`apply_representative_days` · `RepresentativePeriods` ·
`representative_period_error` · `k_medoids` · `rolling_horizon_solve` ·
`ResolutionPlan` · `adaptive_resolution_plan` · `multi_resolution_hierarchy` ·
`apply_adaptive_resolution` · `certify_reduction` · `CertifiedBound`

### Decomposition

`BendersDecomposer` · `solve_with_temporal_benders` ·
`solve_with_spatial_benders` · `solve_with_nested_benders` ·
`solve_with_dantzig_wolfe` · `solve_with_column_generation` ·
`temporal_decomposition` · `recommend_decomposition` ·
`temporal_certified.optimise_temporal_certified`

### Uncertainty

`Scenario` · `solve_stochastic` · `solve_stochastic_ph` · `solve_robust` ·
`BudgetUncertaintySet` · `ChanceConstraint` · `solve_saa_chance_constrained` ·
`solve_general_chance_constrained` · `solve_sddip` · `solve_wasserstein_dro` ·
`solve_risk_averse_benders` · `evaluate_plan` · `reduce_scenarios` ·
`reduce_scenarios_wasserstein` · `generate_demand_scenarios` ·
`generate_renewable_scenarios` · `generate_forced_outage_scenarios` ·
`generate_moment_matching_scenarios` · `cvar_change_of_measure`

### Power flow

`solve_socp_opf` · `solve_socp_opf_multi` · `solve_socp_opf_expansion` ·
`solve_ac_opf_polar` · `obbt_tighten` · `add_weymouth_pipe` ·
`add_head_dependent_hydro`

### Differentiable

`solve_dispatch_with_sensitivities` · `EconomicDispatchLayer` ·
`MultiBusDispatchProblem` · `solve_multibus_dispatch_with_sensitivities` ·
`StorageDispatchProblem` · `solve_storage_dispatch_with_sensitivities` ·
`MultiBusStorageProblem` · `SmoothCommitmentLayer` ·
`fit_commitment_threshold` · `fit_demand_elasticity` ·
`CapacityExpansionProblem` · `CapacityExpansionLayer` ·
`solve_capacity_expansion_with_sensitivities` · `fit_component_params` ·
`TorchDispatchLayer`

### Modules

`nexus_energy.pypsa_compat` — `from_pypsa`
`nexus_energy.diff_bridge` — `multibus_problem_from_system` ·
`d_dispatch_d_co2_price` · `fit_co2_price`
`nexus_energy.autocal` — `fit_params` · `AutoCalibrator`
`nexus_energy.mpc` — `PersistentDispatchSession`
`nexus_energy.solver_tuner` — `tune_solver` · `recommend_lp_method`
`nexus_energy.temporal_certified` — `optimise_temporal_certified`
`nexus_energy.external_solvers` · `nexus_energy.io_tables`

### ML-guided

`extract_system_features` · `extract_timestep_features` ·
`UCWarmstartPredictor` · `MeritOrderPredictor` ·
`HistoricalNeighborPredictor` · `GNNPredictor` · `predict_unit_commitment` ·
`warm_start_from_prediction` · `LearnedVarFixer` · `apply_varfix` ·
`AdaptiveThresholdController` · `solve_with_adaptive_warmstart` · `RLVarFixer` ·
`solve_with_rl_search` · `learned_representative_periods` ·
`feature_embedding_periods`

### Sector coupling and components

`create_power_to_hydrogen` · `create_heat_system` · `create_power_to_gas` ·
`create_temperature_heat_network` · `create_multi_carrier_system` ·
`ComponentTemplate` · `ComponentRegistry` · `registry` · `add_component` ·
`Subsystem` · `CarrierMismatchError`

---

## Licence

MIT. See `LICENSE`.

---

# nexus-opt — User Guide

**One modelling API. Rust core. Every problem class.**

`nexus-opt` is a Python optimisation library with a Rust core. You write the
model once — variables, constraints, an objective — and the same API covers
linear, mixed-integer, quadratic, conic, nonlinear, black-box, multi-objective
and stochastic programming. Changing problem class does not mean changing
library.

If you want the energy-system modelling layer rather than the solver, start at
[nexus-energy](https://github.com/VishalRam24/nexus-energy) — it is built on
this core and installs it for you.

> **Every code block in this guide has been executed against a clean install.**
> Where a printed result is shown, that is the actual output.

---

## Contents

1. [Install and verify](#1-install-and-verify)
2. [Your first model](#2-your-first-model)
3. [Building models](#3-building-models) — variables, expressions, constraints, objectives
4. [Reading results](#4-reading-results)
5. [Problem classes](#5-problem-classes) — LP, MILP, QP, conic, nonlinear
6. [Black-box optimisation](#6-black-box-optimisation) — DE, PSO, SA, CMA-ES
7. [Multi-objective, stochastic and combinatorial](#7-multi-objective-stochastic-and-combinatorial)
8. [Diagnostics](#8-diagnostics) — why a model is infeasible, why a variable took its value
9. [Solver options](#9-solver-options)
10. [Repeated re-solves with PersistentHighs](#10-repeated-re-solves-with-persistenthighs)
11. [Scaling out](#11-scaling-out)
12. [Import, export and validation](#12-import-export-and-validation)
13. [Performance guide](#13-performance-guide)
14. [Architecture](#14-architecture)
15. [Known issues and honest scope](#15-known-issues-and-honest-scope)
16. [API index](#16-api-index)

---

## 1. Install and verify

```bash
pip install nexus-opt
```

Prebuilt wheels cover Linux (glibc and musl; x86_64 / aarch64 / i686), macOS
(arm64 and x86_64) and Windows (x64 and x86). The crate builds as an **abi3**
extension, so one wheel per platform serves every CPython from 3.9 up, and no
Rust toolchain is needed to install.

Verify the install:

```python
import nexus_opt as nx

m = nx.Model()
x = m.variable("x", lower=0)
m.add(x >= 5)
m.minimize(x)
print(m.solve().objective)      # 5.0
```

> **Import `nexus_opt`, not `nexus`.** A short `nexus` alias package exists in
> the source tree for development, but it is not shipped in the wheel. On a
> `pip install`, `import nexus` will fail.

### Building from source

```bash
git clone https://github.com/VishalRam24/nexus-opt
cd nexus-opt
maturin build --release
pip install target/wheels/*.whl
```

Needs a Rust toolchain. Run the tests with `uv run pytest`.

---

## 2. Your first model

A two-variable linear program: meet a demand of 10 units as cheaply as
possible, where `x` costs 2 per unit but is capped at 6, and `y` costs 3.

```python
import nexus_opt as nx

m = nx.Model()
x = m.variable("x", lower=0)
y = m.variable("y", lower=0)

m.add(x + y >= 10)      # meet demand
m.add(x <= 6)           # x is capacity-limited
m.minimize(2 * x + 3 * y)

r = m.solve()
print(r.status, r.objective, r.value(x), r.value(y))
```

```
optimal 24.0 6.0 4.0
```

That is the whole shape of the library: build a `Model`, declare variables,
`add()` constraints written with ordinary Python comparison operators, state an
objective, call `solve()`, read values off the result.

There is no separate "compile" step and no solver-specific syntax. Constraints
compile straight into Rust structures as you write them.

---

## 3. Building models

### Variables

| Call | Kind | Notes |
|---|---|---|
| `m.variable(name, lower=…, upper=…)` | continuous | bounds default to `(-inf, +inf)` |
| `m.integer(name, lower=…, upper=…)` | integer | |
| `m.binary(name)` | binary | equivalent to an integer in `[0, 1]` |
| `m.variables(prefix, count, lower=…, upper=…)` | list of continuous | bulk creation |
| `m.integers(prefix, count, lower=…, upper=…)` | list of integers | bulk creation |

```python
m = nx.Model("kinds")
c = m.variable("c", lower=0.0, upper=10.0)
i = m.integer("i", lower=0, upper=5)
b = m.binary("b")

xs = m.variables("x", 4, lower=0.0, upper=3.0)   # x0 … x3
```

### Expressions

Variables support `+`, `-`, `*`, `/` and `**` with ordinary Python operators.
Multiplying two variables produces a `QuadExpr`; `sin`, `cos` and `sqrt`
produce an `NlExpr`.

```python
linear    = 2 * x + 3 * y - 4
quadratic = w1**2 + 0.5 * (w1 * w2)
nonlinear = nx.sin(x) + nx.sqrt(y)
```

> **Gotcha — parenthesise scalar × variable × variable.** Python evaluates
> `0.05 * w1 * w2` left to right as `(0.05 * w1) * w2`, which is a linear
> expression multiplied by a variable, and raises
> `ValueError: Expr can be multiplied by a scalar or NlExpr`. Write
> `0.05 * (w1 * w2)` instead.

Use `nx.sum()` rather than Python's built-in `sum()` for collections — it
reduces the expression inside Rust in one pass instead of building `N`
intermediate Python objects:

```python
xs = [m.variable(f"x{i}", lower=0, upper=1) for i in range(5)]
m.add(nx.sum(xs) >= 2.5)
m.minimize(nx.sum(xs))
```

### Constraints

`m.add(expr, name=None)` accepts `<=`, `>=` and `==` comparisons. Naming a
constraint is optional but is what makes the diagnostics in
[section 8](#8-diagnostics) readable.

```python
m.add(2 * c + 3 * i - b <= 12.0, name="resource")
m.add(c + i >= 2.0, name="demand")
```

### Objectives

```python
m.minimize(expr)
m.maximize(expr)
```

---

## 4. Reading results

`solve()` returns a `Result`:

| Attribute | Meaning |
|---|---|
| `r.status` | `"optimal"`, `"infeasible"`, `"unbounded"`, `"unknown"` |
| `r.is_optimal` | convenience boolean |
| `r.objective` | objective value at the optimum |
| `r.value(var)` | value of one variable |
| `r.primals` | all primal values, in variable-creation order |
| `r.var_names_list` | variable names, in the same order |

```python
r = m.solve()
if r.is_optimal:
    print(r.objective, r.value(x))
    print(dict(zip(r.var_names_list, r.primals)))
```

**Always check `r.status` before reading `r.objective`.** On a non-optimal
solve, `objective` may be `None`.

Duals and shadow prices come from `r.sensitivity()` — see
[section 8](#8-diagnostics).

---

## 5. Problem classes

The solver is chosen automatically from the model's structure. You can override
it with `solve(solver=…)`, but you rarely need to.

| Model contains | Routed to |
|---|---|
| linear only | HiGHS (simplex / IPM) |
| any integer or binary variable | HiGHS MIP |
| quadratic objective terms | HiGHS QP |
| any cone constraint | Clarabel |
| `sin` / `cos` / `sqrt` / division | Ipopt |

### Mixed-integer linear (MILP)

```python
m = nx.Model("Knapsack")
x = m.variable("x", lower=0.0, upper=10.0)
y = m.integer("y", lower=0, upper=5)
z = m.binary("z")

m.add(2 * x + 3 * y - z <= 12.0, name="resource_limit")
m.add(x + y >= 2.0, name="min_demand")
m.maximize(5 * x + 6 * y + 2 * z)

r = m.solve()
print(r.status, r.objective, r.value(x), r.value(y), r.value(z))
```

```
optimal 34.5 6.5 0.0 1.0
```

### Quadratic (QP)

Quadratic terms appear automatically when variables are multiplied together or
squared.

```python
m = nx.Model("QP_Portfolio")
w1 = m.variable("w1", lower=0, upper=1)
w2 = m.variable("w2", lower=0, upper=1)

m.add(w1 + w2 == 1.0, name="allocation")
m.minimize(0.1 * w1**2 + 0.2 * w2**2 + 0.05 * (w1 * w2))
#                                     ^^^^^^^^^^^^^^^^^ note the parentheses

r = m.solve()
```

### Conic (SOCP / PSD)

A model containing any cone auto-routes to Clarabel — no `solver=` argument.

```python
m = nx.Model("SOCP_Demo")
t  = m.variable("t", lower=0)
x1 = m.variable("x1")
x2 = m.variable("x2")

m.add_soc_cone([t, x1, x2])     # imposes t >= ||(x1, x2)||_2
m.add(x1 >= 3.0)
m.add(x2 >= 4.0)
m.minimize(t)

print(m.solve().objective)      # 5.0  — the 3-4-5 triangle
```

`add_soc_cone` takes `[head, x_1, …, x_k]` and needs at least two variables.

PSD cones take the lower-triangular entries of a `d × d` symmetric matrix —
exactly `d(d+1)/2` variables:

```python
m.add_psd_cone([s11, s21, s22])
```

`quad_form(model, x_vars, P)` builds `xᵀPx` for a PSD matrix `P`, validating
positive-semidefiniteness and raising `ValueError` otherwise:

```python
from nexus_opt import quad_form

expr = quad_form(m, [x1, x2], [[2.0, 0.5], [0.5, 1.0]])
m.minimize(expr)
```

### Nonlinear (NLP)

`sin`, `cos`, `sqrt` and division are first-class expression operators, not
wrappers. Nonlinear constraints live in a parallel container, so the LP/QP
presolve path stays untouched and bit-identical.

```python
m = nx.Model("NLP_Model")
x = m.variable("x", lower=0.1, upper=5.0)
y = m.variable("y", lower=0.1, upper=5.0)

m.add(x / y <= 2.0, name="ratio_bound")
m.minimize(nx.sin(x) + nx.cos(y) + nx.sqrt(x))

r = m.solve()
```

---

## 6. Black-box optimisation

For objectives that are non-convex, non-differentiable, or wrap a simulation,
the metaheuristics take a plain Python callable and a list of bounds. They run
compiled parallel populations in Rust.

```python
# Differential evolution — a good default for rugged landscapes
r = nx.solve_de(
    func=lambda v: (v[0] - 3)**2 + (v[1] - 4)**2,
    bounds=[(-10.0, 10.0), (-10.0, 10.0)],
    pop_size=50,
    max_generations=500,
    f=0.5,          # mutation scale
    cr=0.7,         # crossover rate
    seed=42,
)
print(r.x, r.objective)     # [3.0000, 4.0000] 1.2e-10
```

The others share the same shape:

```python
nx.solve_pso(func, bounds, pop_size=30, max_iterations=500)     # particle swarm
nx.solve_sa(func, bounds, max_iterations=10000)                 # simulated annealing
nx.solve_cmaes(func, bounds, variant="full", restarts=2)        # CMA-ES
```

| Algorithm | Reach for it when |
|---|---|
| `solve_de` | general-purpose global search; the default choice |
| `solve_pso` | flat valleys, cheap objective evaluations |
| `solve_sa` | deep narrow canyons, single-state search |
| `solve_cmaes` | ill-conditioned continuous spaces; `variant="full"` for full covariance with IPOP restarts |

**If you do not know which will win**, race them:

```python
r = nx.solve_portfolio(lambda v: (v[0] - 2.5)**2, [(-10, 10)], max_evals=2000)
```

`solve_portfolio` runs the algorithms concurrently against the same objective
and returns the best result.

---

## 7. Multi-objective, stochastic and combinatorial

### Pareto frontiers

`pareto_frontier` keeps a weighted-sum path for two objectives and routes to
native NSGA-II / NSGA-III when `method` is passed explicitly or when there are
three or more objectives (NSGA-III by default at four or more).

```python
front = nx.pareto_frontier(
    objectives=[f1, f2, f3],        # callables: list[float] -> float
    bounds=[(-5, 5)] * 4,
    n_points=30,
    method="nsga3",                 # "nsga2" | "nsga3" — optional
    generations=300,
    das_dennis_divisions=12,        # NSGA-III reference-point density
    pop_size=150,                   # defaults to n_points * 5
    seed=42,
)
```

### Scenario-based stochastic programming

Minimise expected cost across discrete scenarios:

```python
def cost(x, scenario):
    order, demand = x[0], scenario["demand"]
    return 1.0 * max(0, order - demand) + 5.0 * max(0, demand - order)

r = nx.stochastic_solve(
    cost,
    bounds=[(0, 20)],
    scenarios=[{"demand": d} for d in (8, 10, 12, 14)],
    pop_size=40,
    max_generations=300,
)
print(r.x)      # [14.0] — shortage costs 5x holding, so hedge high
```

### Travelling salesman

`solve_tsp_lazy` adds DFJ subtour-elimination constraints lazily and returns a
`TspResult`:

```python
from nexus_opt import solve_tsp_lazy

res = solve_tsp_lazy(distance_matrix)
print(res.tour)
```

---

## 8. Diagnostics

Production models fail. These three tools are why this library exists in its
current shape.

### Why is my model infeasible?

When `status == "infeasible"`, `infeasibility_report()` isolates the
conflicting subset of constraints and bounds.

```python
m = nx.Model("Infeasible_Setup")
x = m.variable("x", lower=10)
y = m.variable("y", lower=5)
m.add(x + y <= 8, name="demand_cap")     # impossible: x + y >= 15
m.minimize(x + y)

r = m.solve()
if r.status == "infeasible":
    report = r.infeasibility_report()
    print(report["constraint_names"])     # ['demand_cap']
    print(report["explanation"])
    print(report["suggestions"])
```

Name your constraints. An unnamed constraint shows up as `c0`, `c1`, … and the
report becomes much harder to act on.

### Why did this variable take that value?

```python
m = nx.Model("Binding_Setup")
x = m.variable("x", lower=0, upper=100)
m.add(x <= 50, name="grid_limit")
m.minimize(-x)

r = m.solve()
print(r.why(x))
# x = 50.000000 is determined by binding constraint(s): grid_limit

print(r.why_detail(x)["binding_constraints"])
# [(0, 'grid_limit')]
```

`why()` returns a sentence; `why_detail()` returns a dict for automated
reporting.

### Shadow prices

```python
r = m.solve()
sens = r.sensitivity()
```

`sensitivity()` returns four keys:

| Key | Contents |
|---|---|
| `shadow_prices` | `(index, name, dual)` per constraint |
| `most_impactful` | the five constraints with the largest absolute dual |
| `active_constraints` | `(index, name)` for constraints binding at the optimum |
| `variable_bounds` | `(name, status)` — `"between"`, at lower, or at upper |

```python
for idx, name, price in sens["shadow_prices"]:
    print(f"{name}: {price}")
# floor: 2.0
```

---

## 9. Solver options

All of these are lossless — they change how the solve is performed, not what
the optimum is.

| Kwarg | Type / values | Effect |
|---|---|---|
| `solver` | `"highs"`, `"ipopt"`, … | force a specific solver; normally inferred |
| `time_limit` | seconds | stop and report the incumbent |
| `gap` | float | MIP relative gap tolerance |
| `verbose` | bool | solver log to stdout |
| `presolve` | bool, default `True` | run presolve |
| `solver_method` | HiGHS `solver` string | `"simplex"`, `"ipm"`, `"pdlp"`, `"choose"` |
| `run_crossover` | `"on"` / `"off"` / `"choose"` | IPM crossover; `"on"` recovers an exact vertex and duals |
| `parallel` | `"on"` / `"off"` / `"choose"` | HiGHS parallel dual simplex — same algorithm, bit-exact |
| `scale_cleanup` | bool | snap finite bounds/RHS with `0 < abs(v) < 1e-9` to exactly 0 |
| `simplex_scale_strategy` | int 0–5 | HiGHS matrix equilibration |
| `eliminate_redundant` | bool | drop rows whose attainable activity interval sits strictly inside their limits — provably non-binding |
| `warm_start` | a `Result`, or a raw list of floats | forwarded to HiGHS `setSolution()` for MIP warm-start |
| `basis` | basis object | simplex hot-start |
| `threads` | int | worker threads — **see the warning below** |

```python
r = m.solve(time_limit=30.0, presolve=True, solver_method="simplex")
```

> ⚠ **Known bug — the first `solve(threads=N)` in a process can return
> `status="unknown"`.** HiGHS builds its global task scheduler once per
> process; the library tears it down with `resetGlobalScheduler` when the
> thread count changes, and the solve immediately following that reset can come
> back with no result. Subsequent solves are fine. Until this is fixed, either
> leave `threads` unset, or issue one throwaway solve before the one whose
> answer you need, and check `r.status` either way.

### Warm-starting a MIP

```python
r1 = m.solve()
r2 = m.solve(warm_start=r1)     # or warm_start=[…]  (length = variable count)
```

A raw array must have length equal to the model's variable count.

### MIP correctness guard

HiGHS 1.14.0 mis-presolves *degenerate mixed continuous+integer* MIPs — it
returns a wrong "optimal" at gap 0. This reproduces through a one-shot LP-file
read, so it is a solver bug rather than an API issue.

The mitigation is built into `solve()`: for MIPs routed to HiGHS with
`num_vars <= 4000`, the LP relaxation is solved first (LP presolve is correct);
if its optimum is integer-feasible then it *is* the MIP optimum and is returned
directly, bypassing the buggy path. Large MIPs with fractional relaxations fall
through unchanged, so there is no speed cost. `presolve="off"` is always
correct. Guard test: `tests/test_milp_presolve_known_issue.py`.

---

## 10. Repeated re-solves with PersistentHighs

If you are solving a *sequence* of related problems — a rolling horizon, a
calibration loop, a parameter sweep — rebuilding the model each time is waste.
`PersistentHighs` keeps a loaded HiGHS instance alive so each re-solve
hot-starts from the retained simplex basis.

```python
import nexus_opt as nx

m = nx.Model()
x = m.variable("x", lower=0, upper=100)
y = m.variable("y", lower=0, upper=100)
m.add(x + y >= 10, name="demand")
m.minimize(2 * x + 3 * y)

ph = nx.PersistentHighs.from_model(model=m, verbose=False, time_limit=60.0)

r1 = ph.resolve()
print(r1.objective)                  # 20.0

ph.update_col_costs([0], [9.0])      # x now costs 9 instead of 2
r2 = ph.resolve()
print(r2.objective)                  # 30.0
```

In-place update methods (indices are column / row positions):

| Method | Updates |
|---|---|
| `update_col_bounds(idx, lower, upper)` | variable bounds |
| `update_row_bounds(idx, lower, upper)` | constraint RHS |
| `update_col_costs(idx, costs)` | objective coefficients |
| `update_matrix_coeffs(rows, cols, values)` | matrix entries — expert API |

Mirrors of the currently-loaded state, so callers can compute exact deltas
without re-deriving the model:

```python
row_lo, row_hi = ph.row_bounds()
col_lo, col_hi = ph.col_bounds()
costs = ph.col_costs()
print(ph.num_cols(), ph.num_rows())
```

Two things to know:

- **Presolve is disabled** inside the persistent instance. It has to be —
  stable row and column indexing is the whole point.
- **`var_names_list` comes back empty** on a `resolve()` result. Read values
  positionally from `r.primals`, using the column order of the model you built.

`update_matrix_coeffs` is an expert API: you are responsible for keeping any
mirrored higher-level model consistent with it.

---

## 11. Scaling out

### Parallel differential evolution

```python
r = nx.solve_parallel_de(func, bounds, pop_size=40, max_generations=100)
```

### Consensus ADMM

ADMM solves large decoupled problems across sub-agents without any agent
sharing its private objective parameters. Each sub-problem is a callable
`(rho, target) -> list[float]`, and they converge on a shared consensus
variable `z`.

```python
def agent1(rho, target):        # wants z near 3
    return [(2 * 3 + rho * target[0]) / (2 + rho)]

def agent2(rho, target):        # wants z near 7
    return [(2 * 7 + rho * target[0]) / (2 + rho)]

r = nx.solve_admm(
    subproblems=[agent1, agent2],   # note: subproblems, not sub_problems
    dim=1,                          # dimensionality of the consensus variable
    rho=1.0,                        # augmented-Lagrangian penalty
    max_iter=200,
)
print(r["z"])       # ~5.0 — the midpoint
```

---

## 12. Import, export and validation

### Structural summary and validation

Before solving, check the shape of what you built:

```python
print(m.summary())
# Model 'export': 1 vars (1 cont, 0 int, 0 bin), 1 constraints (0 <=, 1 >=, 0 ==), …

for w in m.validate():
    print(f"[{w['severity']}] {w['code']} — {w['message']}")
```

`validate()` flags dangling variables, unconstrained spaces and extreme
coefficients. An empty list means nothing suspicious was found.

### LP export

```python
lp_string = m.to_lp()       # standard CPLEX LP format
```

### Python code generation

`to_python()` emits a self-contained script that rebuilds the identical model —
useful for a reproducible bug report or for handing a colleague the exact
mathematical layout.

```python
with open("reproduced_model.py", "w") as f:
    f.write(m.to_python())
```

---

## 13. Performance guide

1. **Use `nx.sum()`, not Python's `sum()`.** The built-in performs sequential
   pairwise additions, creating `O(N)` intermediate Python objects and crossing
   the PyO3 boundary each time. `nx.sum(vars)` reduces the collection inside
   Rust in one pass.
2. **Reuse the model for repeated solves.** `PersistentHighs` (section 10)
   replaces rebuild-and-cold-solve with an in-place parameter update and a
   hot-started re-solve. On a rolling horizon this is the single largest win
   available.
3. **Leave `presolve` on.** It defaults to `True`, and `eliminate_redundant`
   runs an exact nexus-side row-elimination pass alongside it. Both are
   optimum-preserving.
4. **Let the solver be inferred.** Forcing `solver=` bypasses the routing that
   picks Clarabel for cones and Ipopt for nonlinear expressions.
5. **Name your constraints.** No runtime cost, and it is the difference between
   a usable and a useless infeasibility report.

---

## 14. Architecture

```
┌────────────────────────────────────────────────────────┐
│                      Python API                        │
│  (model building, operator overloading, diagnostics)   │
└──────────────────────────┬─────────────────────────────┘
                           │ PyO3 bindings
┌──────────────────────────▼─────────────────────────────┐
│                       Rust core                        │
│  (flat sparse arrays, expression compiler, parallel    │
│   constraint assembly via Rayon)                       │
└──────────────────────────┬─────────────────────────────┘
                           │ FFI
┌──────────────────────────▼─────────────────────────────┐
│                     Solver backends                    │
│              (HiGHS, Clarabel, Ipopt)                  │
└────────────────────────────────────────────────────────┘
```

Traditional Python optimisation libraries spend most of their time *building*
the model — millions of constraints assembled through nested Python loops.
`nexus-opt` keeps the modelling layer and sparse representation in Rust. Python
expressions built through operator overloading compile directly into Rust
structures (`LinExpr`, `QuadExpr`, `NlExpr`), and at `solve()` time the sparse
constraint matrices are assembled in parallel and streamed into the backend.

Third-party solvers appear in `benchmarks/` as comparison rows only. They are
never wrapped inside the library.

---

## 15. Known issues and honest scope

- **`solve(threads=N)`** can return `status="unknown"` on the first call in a
  process — see the warning in [section 9](#9-solver-options).
- **No GPU backend.** `Model.solve()` accepts no `gpu` argument. Earlier
  documentation described a cuPDLP-C / cuOpt path here; it is not in this
  library. (`nexus-energy` exposes an `lp_backend="gpu"` option that probes for
  cuOpt and falls back to CPU when it is absent.)
- **`sensitivity()` returns no reduced costs.** The four keys it does return
  are listed in [section 8](#8-diagnostics).
- **`PersistentHighs` results carry no variable names** — read `r.primals`
  positionally.
- **`src/presolve.rs` is unused.** `eliminate_redundant` is the presolve pass
  actually wired into the solve path.

### Benchmarks

`benchmarks/` holds the comparison harness and stored results
(`benchmarks/_results/*.jsonl`) against scipy, cvxpy, PuLP, Pyomo, Optuna and
Nevergrad. Numbers are re-runnable rather than quoted.

---

## 16. API index

Everything exported from the top-level `nexus_opt` namespace.

### Core types

| Name | What it is |
|---|---|
| `Model` | the model container |
| `Var`, `Expr`, `QuadExpr`, `NlExpr` | expression types |
| `Constraint` | a compiled constraint |
| `Result` | solve result |
| `BlackBoxResult` | metaheuristic result (`.x`, `.objective`) |
| `TspResult` | TSP result (`.tour`) |
| `PersistentHighs` | in-place re-solve session |

### Model methods

`variable` · `integer` · `binary` · `variables` · `integers` · `add` ·
`add_soc_cone` · `add_psd_cone` · `minimize` · `maximize` · `solve` ·
`summary` · `validate` · `to_lp` · `to_python`

### Result methods

`status` · `is_optimal` · `objective` · `value` · `primals` ·
`var_names_list` · `sensitivity` · `why` · `why_detail` ·
`infeasibility_report`

### Functions

| Name | Purpose |
|---|---|
| `sum`, `minimize`, `maximize` | expression helpers |
| `sin`, `cos`, `sqrt` | nonlinear operators |
| `quad_form(model, xs, P)` | build `xᵀPx`, validating PSD |
| `solve_de`, `solve_pso`, `solve_sa`, `solve_cmaes` | metaheuristics |
| `solve_portfolio` | race the metaheuristics |
| `solve_parallel_de` | parallel differential evolution |
| `solve_admm` | consensus ADMM |
| `pareto_frontier` | multi-objective frontier (NSGA-II / III) |
| `stochastic_solve` | scenario-based expectation minimisation |
| `solve_tsp_lazy` | TSP with lazy subtour elimination |

---

## Licence

MIT. The design and layout of the project site follow Notus React (Creative
Tim, MIT); see `LICENSE` for retained notices.