diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs
index c4a290e..60d192a 100644
--- a/crates/engine/src/cmd/mod.rs
+++ b/crates/engine/src/cmd/mod.rs
@@ -14,6 +14,7 @@ pub mod pivot;
 pub mod print;
 pub mod review;
 pub mod sheet;
+pub mod solver;
 pub mod spelling;
 pub mod view;
 pub mod whatif;
@@ -118,6 +119,7 @@ pub fn command_specs() -> &'static [CommandSpec] {
         v.extend(extra::specs());
         v.extend(draw::specs());
         v.extend(spelling::specs());
+        v.extend(solver::specs());
         v.extend(pivot::specs());
         v.extend(whatif::specs());
         v.extend(print::specs());
diff --git a/crates/engine/src/cmd/solver.rs b/crates/engine/src/cmd/solver.rs
new file mode 100644
index 0000000..5759f54
--- /dev/null
+++ b/crates/engine/src/cmd/solver.rs
@@ -0,0 +1,1487 @@
+//! Solver (Data › Analysis › Solver): finds values for variable cells that maximise, minimise
+//! or hit a value in an objective cell, subject to constraints.
+//!
+//! The workbook is the model: every evaluation sets the variable cells on a copy of the
+//! workbook and recalculates their dependents ([`Trial`]), so any formula works. Two methods:
+//!
+//! - **Simplex LP**: the model is linearised from `n + 1` evaluations and checked for linearity
+//!   at one more point; then a dense two-phase simplex (Dantzig's rule, Bland's rule when it
+//!   stalls) solves it exactly. Integer and binary constraints use depth-first branch and bound.
+//! - **GRG Nonlinear** (the name Excel's dialog uses): GridCraft's method is *not* a generalised
+//!   reduced gradient. It is an augmented-Lagrangian method (Powell–Hestenes–Rockafellar) whose
+//!   inner problems are solved by projected gradient descent with Barzilai–Borwein steps and an
+//!   Armijo line search, on central-difference gradients. Variable bounds are kept by
+//!   projection. It finds a local optimum of a smooth model; integer constraints use a small
+//!   branch and bound on top.
+//!
+//! Work is bounded by an evaluation budget (no clock: `Instant` traps on wasm). The model is
+//! saved as hidden sheet-scoped names (`solver_opt`, `solver_adj`, `solver_lhs1`…), so it
+//! travels with the workbook through XLSX.
+
+use gridcraft_core::{CellRef, RangeRef, Value};
+use gridcraft_model::{Cell, DefinedName, Workbook};
+use serde_json::{Value as Json, json};
+
+use super::whatif::Trial;
+use super::*;
+
+const CMD: &str = "data.solver";
+/// Excel's limit on decision variables.
+const MAX_VARS: usize = 200;
+/// Most scalar constraints (after expanding ranges).
+const MAX_CONSTRAINTS: usize = 1000;
+const DEFAULT_EVALS: u64 = 20_000;
+const MAX_EVALS: u64 = 200_000;
+/// Most simplex pivots × tableau cells per solve.
+const LP_WORK: u64 = 2_000_000_000;
+const MAX_NODES_LP: usize = 10_000;
+const MAX_NODES_NLP: usize = 200;
+/// Constraint precision (Excel's default 0.000001).
+const PRECISION: f64 = 1e-6;
+const EPS: f64 = 1e-9;
+
+pub fn specs() -> Vec<CommandSpec> {
+    vec![cmd!(
+        "data.solver",
+        "Solver…",
+        ["Data", "Analysis"],
+        None,
+        "{objective: \"B10\" (formula cell), goal?: max|min|value (default max), value?: number (for goal value), variables: \"B2:B4,D1\" (≤200 constant cells, same sheet), constraints?: [\"B5 <= 10\", \"B2:B4 >= C2:C4\", \"B2 int\", \"B3 bin\", {lhs: \"B6\", op: \"=\", rhs: 3 | \"C6\"}], nonNegative?: bool (default true: variables without a lower bound get 0), method?: grg|simplex (default grg), keep?: bool (default true; false restores the original values), saveModel?: bool (default true: hidden solver_* names), maxEvaluations?: 20000, integerOptimality?: percent (default 0), saved?: bool (solve the sheet's saved model; other params override)} → {status: optimal|converged|iterationLimit|infeasible|unbounded|notLinear, message, objective, values: {cell: value}, evaluations, changed}. No `objective` → {model: the saved model} (the menu opens the dialog).",
+        has_doc,
+        run
+    )]
+}
+
+// ---------------------------------------------------------------- model
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Rel {
+    Le,
+    Eq,
+    Ge,
+    Int,
+    Bin,
+}
+
+impl Rel {
+    fn parse(s: &str) -> Option<Rel> {
+        match s.trim().to_ascii_lowercase().as_str() {
+            "<=" | "≤" | "=<" | "le" => Some(Rel::Le),
+            ">=" | "≥" | "=>" | "ge" => Some(Rel::Ge),
+            "=" | "==" | "eq" => Some(Rel::Eq),
+            "int" | "integer" => Some(Rel::Int),
+            "bin" | "binary" => Some(Rel::Bin),
+            _ => None,
+        }
+    }
+    fn text(self) -> &'static str {
+        match self {
+            Rel::Le => "<=",
+            Rel::Eq => "=",
+            Rel::Ge => ">=",
+            Rel::Int => "int",
+            Rel::Bin => "bin",
+        }
+    }
+    /// The relation code of the saved model (`solver_relN`), as in Solver's VBA `SolverAdd`.
+    fn code(self) -> u8 {
+        match self {
+            Rel::Le => 1,
+            Rel::Eq => 2,
+            Rel::Ge => 3,
+            Rel::Int => 4,
+            Rel::Bin => 5,
+        }
+    }
+    fn from_code(c: &str) -> Option<Rel> {
+        match c.trim() {
+            "1" => Some(Rel::Le),
+            "2" => Some(Rel::Eq),
+            "3" => Some(Rel::Ge),
+            "4" => Some(Rel::Int),
+            "5" => Some(Rel::Bin),
+            _ => None,
+        }
+    }
+}
+
+#[derive(Clone, Debug)]
+enum Rhs {
+    None,
+    Const(f64),
+    Cells(Vec<CellRef>),
+}
+
+#[derive(Clone, Debug)]
+struct Constraint {
+    lhs: Vec<CellRef>,
+    lhs_text: String,
+    rel: Rel,
+    rhs: Rhs,
+    rhs_text: String,
+}
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Goal {
+    Max,
+    Min,
+    Value(f64),
+}
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Method {
+    Simplex,
+    Grg,
+}
+
+#[derive(Clone, Debug)]
+struct Model {
+    sheet: usize,
+    objective: CellRef,
+    goal: Goal,
+    vars: Vec<CellRef>,
+    vars_text: String,
+    cons: Vec<Constraint>,
+    non_negative: bool,
+    method: Method,
+}
+
+/// Cells of `"B2:B4,D1"` on one sheet (row-major per area, in order), at most `cap`.
+fn cells_of(wb: &Workbook, t: &str, default_sheet: usize, cap: usize) -> Result<(usize, Vec<CellRef>)> {
+    let mut sheet = None;
+    let mut out = Vec::new();
+    for part in t.split(',').map(str::trim).filter(|x| !x.is_empty()) {
+        let (sh, body) = split_sheet(part);
+        let si = match sh {
+            Some(n) => wb.sheet_index(&n).ok_or_else(|| bad(CMD, format!("no sheet named `{n}`")))?,
+            None => default_sheet,
+        };
+        if *sheet.get_or_insert(si) != si {
+            return Err(bad(CMD, "all cells of a reference must be on one sheet"));
+        }
+        let r = RangeRef::parse(body.trim()).ok_or_else(|| bad(CMD, format!("not a reference: `{part}`")))?;
+        if r.count() > cap as u64 || out.len() + r.count() as usize > cap {
+            return Err(bad(CMD, format!("too many cells (at most {cap})")));
+        }
+        out.extend(r.iter());
+    }
+    if out.is_empty() {
+        return Err(bad(CMD, format!("not a reference: `{t}`")));
+    }
+    Ok((sheet.unwrap_or(default_sheet), out))
+}
+
+fn num_text(t: &str) -> Option<f64> {
+    gridcraft_core::parse::parse_number_text(t.trim()).filter(|v| v.is_finite())
+}
+
+/// Splits `"B5 <= 10"` / `"B2 int"` into (lhs, relation, rhs).
+fn split_constraint(t: &str) -> Option<(String, Rel, String)> {
+    for op in ["<=", ">=", "=<", "=>", "≤", "≥", "=="] {
+        if let Some((a, b)) = t.split_once(op) {
+            return Some((a.trim().into(), Rel::parse(op)?, b.trim().into()));
+        }
+    }
+    if let Some((a, b)) = t.split_once('=') {
+        return Some((a.trim().into(), Rel::Eq, b.trim().into()));
+    }
+    let lower = t.trim().to_ascii_lowercase();
+    for (suffix, rel) in [("integer", Rel::Int), ("binary", Rel::Bin), ("int", Rel::Int), ("bin", Rel::Bin)] {
+        if let Some(a) = lower.strip_suffix(suffix) {
+            let a = a.trim();
+            return Some((t.trim().get(..a.len())?.trim().to_string(), rel, String::new()));
+        }
+    }
+    None
+}
+
+fn parse_constraint(wb: &Workbook, sheet: usize, c: &Json) -> Result<Constraint> {
+    let (lhs_text, rel, rhs_json) = match c {
+        Json::String(t) => {
+            let (l, r, rhs) = split_constraint(t).ok_or_else(|| bad(CMD, format!("not a constraint: `{t}`")))?;
+            (l, r, Json::String(rhs))
+        }
+        Json::Object(_) => {
+            let lhs = str_param(c, "lhs").or_else(|| str_param(c, "cell")).or_else(|| str_param(c, "range")).unwrap_or("").to_string();
+            let op = str_param(c, "op").or_else(|| str_param(c, "relation")).unwrap_or("");
+            let rel = Rel::parse(op).ok_or_else(|| bad(CMD, format!("unknown constraint relation `{op}`")))?;
+            (lhs, rel, c.get("rhs").or_else(|| c.get("value")).cloned().unwrap_or(Json::Null))
+        }
+        _ => return Err(bad(CMD, "a constraint is a string like \"B5 <= 10\" or {lhs, op, rhs}")),
+    };
+    let (lsi, lhs) = cells_of(wb, &lhs_text, sheet, MAX_CONSTRAINTS)?;
+    if lsi != sheet {
+        return Err(EngineError::Other("Constraint cells must be on the objective cell's sheet.".into()));
+    }
+    let (rhs, rhs_text) = match rel {
+        Rel::Int | Rel::Bin => (Rhs::None, String::new()),
+        _ => match &rhs_json {
+            Json::Number(n) => {
+                let v = n.as_f64().filter(|v| v.is_finite()).ok_or_else(|| bad(CMD, "the constraint's right side must be a number"))?;
+                (Rhs::Const(v), v.to_string())
+            }
+            Json::String(t) => {
+                let t = t.trim().trim_start_matches('=').trim();
+                if let Some(v) = num_text(t) {
+                    (Rhs::Const(v), t.to_string())
+                } else {
+                    let (rsi, cells) = cells_of(wb, t, sheet, MAX_CONSTRAINTS)?;
+                    if rsi != sheet {
+                        return Err(EngineError::Other("Constraint cells must be on the objective cell's sheet.".into()));
+                    }
+                    if cells.len() != 1 && cells.len() != lhs.len() {
+                        return Err(EngineError::Other("The constraint's two sides must be the same size (or the right side one cell).".into()));
+                    }
+                    (Rhs::Cells(cells), t.to_string())
+                }
+            }
+            _ => return Err(bad(CMD, "missing the constraint's right side")),
+        },
+    };
+    Ok(Constraint { lhs, lhs_text, rel, rhs, rhs_text })
+}
+
+fn parse_model(s: &Session, p: &Json) -> Result<Model> {
+    let d = s.doc()?;
+    let wb = &d.wb;
+    let obj_text = str_param(p, "objective").or_else(|| str_param(p, "setObjective")).or_else(|| str_param(p, "set")).unwrap_or("");
+    let (sheet, obj) = cells_of(wb, obj_text, wb.active_sheet, 1).map_err(|_| EngineError::Other("The objective must be a single cell.".into()))?;
+    let objective = *obj.first().ok_or_else(|| EngineError::Other("The objective must be a single cell.".into()))?;
+    let sh = wb.sheet(sheet).ok_or(EngineError::NoDocument)?;
+    if sh.cell(objective).is_none_or(|c| c.formula.is_none()) {
+        return Err(EngineError::Other("The objective cell must contain a formula.".into()));
+    }
+    let goal = match str_param(p, "goal").or_else(|| str_param(p, "to")).unwrap_or("max").to_ascii_lowercase().as_str() {
+        "max" | "maximize" | "maximum" => Goal::Max,
+        "min" | "minimize" | "minimum" => Goal::Min,
+        "value" | "valueof" | "value of" => {
+            let v = f64_param(p, "value")
+                .or_else(|| str_param(p, "value").and_then(num_text))
+                .filter(|v| v.is_finite())
+                .ok_or_else(|| bad(CMD, "goal `value` needs a number in `value`"))?;
+            Goal::Value(v)
+        }
+        g => return Err(bad(CMD, format!("unknown goal `{g}` (max, min or value)"))),
+    };
+    let vars_text = str_param(p, "variables").or_else(|| str_param(p, "changing")).or_else(|| str_param(p, "byChanging")).unwrap_or("").to_string();
+    let (vsi, mut vars) = cells_of(wb, &vars_text, sheet, MAX_VARS).map_err(|e| match e {
+        EngineError::BadParams { msg, .. } if msg.starts_with("too many") => {
+            EngineError::Other(format!("Too many variable cells (at most {MAX_VARS})."))
+        }
+        _ => EngineError::Other("The variable cells must be a reference like B2:B4.".into()),
+    })?;
+    if vsi != sheet {
+        return Err(EngineError::Other("The variable cells must be on the objective cell's sheet.".into()));
+    }
+    let mut seen = std::collections::HashSet::new();
+    vars.retain(|c| seen.insert(*c));
+    if vars.iter().any(|c| sh.cell(*c).is_some_and(|c| c.formula.is_some())) {
+        return Err(EngineError::Other("Variable cells must hold constants, not formulas.".into()));
+    }
+    let mut cons = Vec::new();
+    let mut scalar = 0usize;
+    if let Some(list) = p.get("constraints") {
+        let list = list.as_array().ok_or_else(|| bad(CMD, "`constraints` must be an array"))?;
+        if list.len() > MAX_CONSTRAINTS {
+            return Err(EngineError::Other(format!("Too many constraints (at most {MAX_CONSTRAINTS}).")));
+        }
+        for c in list {
+            let c = parse_constraint(wb, sheet, c)?;
+            scalar = scalar.saturating_add(c.lhs.len());
+            if scalar > MAX_CONSTRAINTS {
+                return Err(EngineError::Other(format!("Too many constraints (at most {MAX_CONSTRAINTS} cells).")));
+            }
+            if matches!(c.rel, Rel::Int | Rel::Bin) && c.lhs.iter().any(|x| !vars.contains(x)) {
+                return Err(EngineError::Other("Integer and binary constraints must be on variable cells.".into()));
+            }
+            cons.push(c);
+        }
+    }
+    let non_negative = bool_param(p, "nonNegative").unwrap_or(true);
+    let method = match str_param(p, "method").unwrap_or("grg").to_ascii_lowercase().replace([' ', '_', '-'], "").as_str() {
+        "grg" | "grgnonlinear" | "nonlinear" => Method::Grg,
+        "simplex" | "simplexlp" | "lp" | "linear" => Method::Simplex,
+        "evolutionary" => return Err(EngineError::Other("The Evolutionary method isn't available; use GRG Nonlinear or Simplex LP.".into())),
+        m => return Err(bad(CMD, format!("unknown method `{m}` (grg or simplex)"))),
+    };
+    Ok(Model { sheet, objective, goal, vars, vars_text, cons, non_negative, method })
+}
+
+// ---------------------------------------------------------------- saved model
+
+fn ref_text(sheet_name: &str, cells: &str) -> String {
+    let q = crate::pivot::quote_sheet(sheet_name);
+    cells
+        .split(',')
+        .map(str::trim)
+        .filter(|x| !x.is_empty())
+        .filter_map(|part| {
+            let r = RangeRef::parse(split_sheet(part).1.trim())?;
+            let abs = |c: CellRef| format!("${}${}", gridcraft_core::col_to_letters(c.col), c.row as u64 + 1);
+            Some(if r.is_single() { format!("{q}!{}", abs(r.start)) } else { format!("{q}!{}:{}", abs(r.start), abs(r.end)) })
+        })
+        .collect::<Vec<_>>()
+        .join(",")
+}
+
+/// The model as hidden names scoped to its sheet (replacing any saved before).
+fn save_model(wb: &mut Workbook, m: &Model) {
+    let Some(sheet_name) = wb.sheet(m.sheet).map(|s| s.name.clone()) else { return };
+    wb.names.retain(|n| !(n.scope == Some(m.sheet) && n.name.to_ascii_lowercase().starts_with("solver_")));
+    let mut put = |name: String, formula: String| {
+        wb.names.push(DefinedName { name, scope: Some(m.sheet), formula, comment: String::new(), hidden: true });
+    };
+    put("solver_opt".into(), ref_text(&sheet_name, &m.objective.a1()));
+    let (typ, val) = match m.goal {
+        Goal::Max => (1, 0.0),
+        Goal::Min => (2, 0.0),
+        Goal::Value(v) => (3, v),
+    };
+    put("solver_typ".into(), typ.to_string());
+    put("solver_val".into(), val.to_string());
+    put("solver_adj".into(), ref_text(&sheet_name, &m.vars_text));
+    for (i, c) in m.cons.iter().enumerate() {
+        let k = i + 1;
+        put(format!("solver_lhs{k}"), ref_text(&sheet_name, &c.lhs_text));
+        put(format!("solver_rel{k}"), c.rel.code().to_string());
+        let rhs = match &c.rhs {
+            Rhs::None => {
+                if c.rel == Rel::Bin {
+                    "\"binary\"".to_string()
+                } else {
+                    "\"integer\"".to_string()
+                }
+            }
+            Rhs::Const(v) => v.to_string(),
+            Rhs::Cells(_) => ref_text(&sheet_name, &c.rhs_text),
+        };
+        put(format!("solver_rhs{k}"), rhs);
+    }
+    put("solver_num".into(), m.cons.len().to_string());
+    put("solver_neg".into(), if m.non_negative { "1" } else { "2" }.into());
+    put("solver_eng".into(), if m.method == Method::Simplex { "2" } else { "1" }.into());
+    put("solver_ver".into(), "3".into());
+}
+
+/// The saved model of a sheet as `data.solver` params, if any.
+pub fn saved_model(wb: &Workbook, sheet: usize) -> Option<Json> {
+    let get = |n: &str| {
+        wb.names
+            .iter()
+            .find(|d| d.name.eq_ignore_ascii_case(n) && (d.scope == Some(sheet) || d.scope.is_none()))
+            .map(|d| d.formula.trim().trim_start_matches('=').to_string())
+    };
+    let local = |t: String| -> String { t.split(',').map(|p| split_sheet(p.trim()).1.replace('$', "")).collect::<Vec<_>>().join(",") };
+    let objective = local(get("solver_opt")?);
+    if objective.is_empty() {
+        return None;
+    }
+    let goal = match get("solver_typ").as_deref().map(str::trim) {
+        Some("2") => "min",
+        Some("3") => "value",
+        _ => "max",
+    };
+    let value = get("solver_val").and_then(|v| num_text(&v)).unwrap_or(0.0);
+    let count = get("solver_num").and_then(|v| num_text(&v)).unwrap_or(0.0).clamp(0.0, MAX_CONSTRAINTS as f64) as usize;
+    let mut constraints = Vec::new();
+    for k in 1..=count {
+        let (Some(lhs), Some(rel)) = (get(&format!("solver_lhs{k}")), get(&format!("solver_rel{k}")).and_then(|r| Rel::from_code(&r))) else {
+            continue;
+        };
+        let lhs = local(lhs);
+        constraints.push(match rel {
+            Rel::Int | Rel::Bin => json!(format!("{lhs} {}", rel.text())),
+            _ => {
+                let rhs = get(&format!("solver_rhs{k}")).unwrap_or_default();
+                let rhs = if num_text(&rhs).is_some() { rhs } else { local(rhs) };
+                json!(format!("{lhs} {} {rhs}", rel.text()))
+            }
+        });
+    }
+    Some(json!({
+        "objective": objective,
+        "goal": goal,
+        "value": value,
+        "variables": get("solver_adj").map(local).unwrap_or_default(),
+        "constraints": constraints,
+        "nonNegative": get("solver_neg").as_deref().map(str::trim) != Some("2"),
+        "method": if get("solver_eng").as_deref().map(str::trim) == Some("2") { "simplex" } else { "grg" },
+    }))
+}
+
+// ---------------------------------------------------------------- evaluation
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Sense {
+    Le,
+    Ge,
+    Eq,
+}
+
+/// A scalar constraint `lhs (sense) rhs`.
+#[derive(Clone, Debug)]
+struct Row {
+    lhs: CellRef,
+    sense: Sense,
+    rhs: Option<CellRef>,
+    rhs_const: f64,
+}
+
+/// The values the solvers see at a point: the objective and `lhs − rhs` per row.
+#[derive(Clone, Debug)]
+struct Point {
+    f: f64,
+    g: Vec<f64>,
+}
+
+struct Evaluator {
+    trial: Trial,
+    sheet: usize,
+    vars: Vec<CellRef>,
+    objective: CellRef,
+    rows: Vec<Row>,
+    evals: u64,
+    budget: u64,
+}
+
+fn as_num(v: &Value) -> f64 {
+    match v {
+        Value::Number(n) => *n,
+        Value::Empty => 0.0,
+        Value::Bool(b) => f64::from(u8::from(*b)),
+        _ => f64::NAN,
+    }
+}
+
+impl Evaluator {
+    fn exhausted(&self) -> bool {
+        self.evals >= self.budget
+    }
+
+    /// Evaluates the model at `x`; `None` once the budget is spent.
+    fn eval(&mut self, x: &[f64]) -> Option<Point> {
+        if self.exhausted() {
+            return None;
+        }
+        self.evals += 1;
+        let sets: Vec<(usize, CellRef, Value)> = self.vars.iter().zip(x).map(|(c, v)| (self.sheet, *c, Value::number(*v))).collect();
+        let mut reads = Vec::with_capacity(1 + 2 * self.rows.len());
+        reads.push((self.sheet, self.objective));
+        for r in &self.rows {
+            reads.push((self.sheet, r.lhs));
+            if let Some(c) = r.rhs {
+                reads.push((self.sheet, c));
+            }
+        }
+        let vals = self.trial.run(&sets, &reads);
+        let mut it = vals.iter().map(as_num);
+        let f = it.next().unwrap_or(f64::NAN);
+        let mut g = Vec::with_capacity(self.rows.len());
+        for r in &self.rows {
+            let l = it.next().unwrap_or(f64::NAN);
+            let rv = if r.rhs.is_some() { it.next().unwrap_or(f64::NAN) } else { r.rhs_const };
+            g.push(l - rv);
+        }
+        Some(Point { f, g })
+    }
+}
+
+// ---------------------------------------------------------------- problem
+
+struct Problem {
+    lo: Vec<f64>,
+    hi: Vec<f64>,
+    integer: Vec<bool>,
+    rows: Vec<Row>,
+}
+
+/// Splits the constraints into variable bounds and general rows.
+fn build_problem(m: &Model) -> Problem {
+    let n = m.vars.len();
+    let mut lo = vec![f64::NEG_INFINITY; n];
+    let mut hi = vec![f64::INFINITY; n];
+    let mut has_lower = vec![false; n];
+    let mut integer = vec![false; n];
+    let mut rows = Vec::new();
+    let var_index = |c: &CellRef| m.vars.iter().position(|v| v == c);
+    for c in &m.cons {
+        for (k, cell) in c.lhs.iter().enumerate() {
+            let j = var_index(cell);
+            match c.rel {
+                Rel::Int | Rel::Bin => {
+                    if let Some(j) = j {
+                        if let Some(b) = integer.get_mut(j) {
+                            *b = true;
+                        }
+                        if c.rel == Rel::Bin
+                            && let (Some(l), Some(h), Some(hl)) = (lo.get_mut(j), hi.get_mut(j), has_lower.get_mut(j))
+                        {
+                            *l = l.max(0.0);
+                            *h = h.min(1.0);
+                            *hl = true;
+                        }
+                    }
+                }
+                rel => {
+                    let sense = match rel {
+                        Rel::Le => Sense::Le,
+                        Rel::Ge => Sense::Ge,
+                        _ => Sense::Eq,
+                    };
+                    let (rhs, rhs_const) = match &c.rhs {
+                        Rhs::Const(v) => (None, *v),
+                        Rhs::Cells(cells) => (cells.get(k).or_else(|| cells.first()).copied(), 0.0),
+                        Rhs::None => (None, 0.0),
+                    };
+                    match (j, rhs) {
+                        // A constant bound on a variable cell is a box bound.
+                        (Some(j), None) => {
+                            if let (Some(l), Some(h), Some(hl)) = (lo.get_mut(j), hi.get_mut(j), has_lower.get_mut(j)) {
+                                if sense != Sense::Ge {
+                                    *h = h.min(rhs_const);
+                                }
+                                if sense != Sense::Le {
+                                    *l = l.max(rhs_const);
+                                    *hl = true;
+                                }
+                            }
+                        }
+                        _ => rows.push(Row { lhs: *cell, sense, rhs, rhs_const }),
+                    }
+                }
+            }
+        }
+    }
+    if m.non_negative {
+        for (l, hl) in lo.iter_mut().zip(&has_lower) {
+            if !hl {
+                *l = l.max(0.0);
+            }
+        }
+    }
+    Problem { lo, hi, integer, rows }
+}
+
+/// What a relaxation (or the whole solve) ended with.
+#[derive(Clone, Debug)]
+enum Outcome {
+    /// x, objective to minimise, and whether the solver converged (vs. ran out of budget).
+    Solved(Vec<f64>, f64, bool),
+    Infeasible,
+    Unbounded,
+    Limit,
+}
+
+/// Depth-first branch and bound over a relaxation solver; returns the best integer point.
+fn branch_and_bound(
+    lo: &[f64],
+    hi: &[f64],
+    integer: &[bool],
+    max_nodes: usize,
+    gap: f64,
+    relax: &mut dyn FnMut(&[f64], &[f64]) -> Outcome,
+) -> (Outcome, usize) {
+    let mut stack: Vec<(Vec<f64>, Vec<f64>)> = vec![(lo.to_vec(), hi.to_vec())];
+    let mut best: Option<(Vec<f64>, f64)> = None;
+    let mut nodes = 0usize;
+    let mut limited = false;
+    let mut root = true;
+    while let Some((l, h)) = stack.pop() {
+        if nodes >= max_nodes {
+            limited = true;
+            break;
+        }
+        nodes += 1;
+        let out = relax(&l, &h);
+        let was_root = std::mem::replace(&mut root, false);
+        let (x, obj, converged) = match out {
+            Outcome::Solved(x, obj, c) => (x, obj, c),
+            Outcome::Unbounded if was_root => return (Outcome::Unbounded, nodes),
+            Outcome::Limit => {
+                limited = true;
+                break;
+            }
+            _ => continue,
+        };
+        if !converged {
+            limited = true;
+        }
+        if let Some((_, b)) = &best
+            && obj >= b - gap * b.abs().max(1e-12) - 1e-9 * (1.0 + b.abs())
+        {
+            continue;
+        }
+        // The most fractional integer variable.
+        let frac = x
+            .iter()
+            .zip(integer)
+            .enumerate()
+            .filter(|(_, (_, int))| **int)
+            .map(|(j, (v, _))| (j, (v - v.round()).abs()))
+            .filter(|(_, f)| *f > 1e-6)
+            .max_by(|a, b| a.1.total_cmp(&b.1));
+        match frac {
+            None => {
+                let snapped: Vec<f64> = x.iter().zip(integer).map(|(v, int)| if *int { v.round() } else { *v }).collect();
+                best = Some((snapped, obj));
+            }
+            Some((j, _)) => {
+                let v = x.get(j).copied().unwrap_or(0.0);
+                let (mut down_h, mut up_l) = (h.clone(), l.clone());
+                if let Some(e) = down_h.get_mut(j) {
+                    *e = e.min(v.floor());
+                }
+                if let Some(e) = up_l.get_mut(j) {
+                    *e = e.max(v.ceil());
+                }
+                // Explore the nearer side first (pushed last).
+                if v - v.floor() < 0.5 {
+                    stack.push((up_l, h));
+                    stack.push((l, down_h));
+                } else {
+                    stack.push((l, down_h));
+                    stack.push((up_l, h));
+                }
+            }
+        }
+    }
+    match best {
+        Some((x, obj)) => (Outcome::Solved(x, obj, !limited), nodes),
+        None if limited => (Outcome::Limit, nodes),
+        None => (Outcome::Infeasible, nodes),
+    }
+}
+
+// ---------------------------------------------------------------- simplex
+
+/// Pivots the tableau on (r, j), updating the objective row.
+fn pivot(t: &mut [Vec<f64>], obj: &mut [f64], r: usize, j: usize) {
+    let Some(prow) = t.get_mut(r) else { return };
+    let pv = prow.get(j).copied().unwrap_or(0.0);
+    if pv == 0.0 {
+        return;
+    }
+    for v in prow.iter_mut() {
+        *v /= pv;
+    }
+    let prow = prow.clone();
+    for (i, row) in t.iter_mut().enumerate() {
+        if i == r {
+            continue;
+        }
+        let f = row.get(j).copied().unwrap_or(0.0);
+        if f != 0.0 {
+            for (v, p) in row.iter_mut().zip(&prow) {
+                *v -= f * p;
+            }
+        }
+    }
+    let f = obj.get(j).copied().unwrap_or(0.0);
+    if f != 0.0 {
+        for (v, p) in obj.iter_mut().zip(&prow) {
+            *v -= f * p;
+        }
+    }
+}
+
+enum Step {
+    Optimal,
+    Unbounded,
+    Limit,
+}
+
+/// Minimises with the reduced costs in `obj` (last entry: −objective). Columns not `allowed`
+/// never enter.
+fn simplex(t: &mut [Vec<f64>], basis: &mut [usize], obj: &mut [f64], allowed: &[bool], budget: &mut u64) -> Step {
+    let n = obj.len().saturating_sub(1);
+    let cost = (t.len() as u64 + 1).saturating_mul(n as u64 + 1).max(1);
+    let mut stalled = 0u32;
+    loop {
+        if *budget < cost {
+            return Step::Limit;
+        }
+        *budget -= cost;
+        let bland = stalled > 50;
+        let mut enter = None;
+        let mut best = -EPS;
+        for (j, d) in obj.iter().take(n).enumerate() {
+            if !allowed.get(j).copied().unwrap_or(false) || *d >= -EPS {
+                continue;
+            }
+            if bland {
+                enter = Some(j);
+                break;
+            }
+            if *d < best {
+                best = *d;
+                enter = Some(j);
+            }
+        }
+        let Some(j) = enter else { return Step::Optimal };
+        let mut leave: Option<(usize, f64)> = None;
+        for (i, row) in t.iter().enumerate() {
+            let a = row.get(j).copied().unwrap_or(0.0);
+            if a <= EPS {
+                continue;
+            }
+            let ratio = row.get(n).copied().unwrap_or(0.0) / a;
+            let better = match leave {
+                None => true,
+                Some((li, lr)) => ratio < lr - 1e-12 || (ratio <= lr + 1e-12 && basis.get(i) < basis.get(li)),
+            };
+            if better {
+                leave = Some((i, ratio));
+            }
+        }
+        let Some((r, ratio)) = leave else { return Step::Unbounded };
+        stalled = if ratio <= EPS { stalled + 1 } else { 0 };
+        pivot(t, obj, r, j);
+        if let Some(b) = basis.get_mut(r) {
+            *b = j;
+        }
+    }
+}
+
+/// Minimises `c·x` subject to `rows` (coefficients, sense, rhs) and `lo ≤ x ≤ hi`.
+fn solve_lp(c: &[f64], rows: &[(Vec<f64>, Sense, f64)], lo: &[f64], hi: &[f64], budget: &mut u64) -> Outcome {
+    let n = c.len();
+    if lo.iter().zip(hi).any(|(l, h)| l > h) {
+        return Outcome::Infeasible;
+    }
+    // x_j = offset + Σ sign·y_k with y ≥ 0.
+    struct Map {
+        var: usize,
+        sign: f64,
+    }
+    let mut cols: Vec<Map> = Vec::new();
+    let mut offset = vec![0.0; n];
+    let mut extra: Vec<(Vec<(usize, f64)>, Sense, f64)> = Vec::new();
+    for j in 0..n {
+        let (l, h) = (lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), hi.get(j).copied().unwrap_or(f64::INFINITY));
+        if l.is_finite() {
+            if let Some(o) = offset.get_mut(j) {
+                *o = l;
+            }
+            cols.push(Map { var: j, sign: 1.0 });
+            if h.is_finite() {
+                extra.push((vec![(cols.len() - 1, 1.0)], Sense::Le, h - l));
+            }
+        } else if h.is_finite() {
+            if let Some(o) = offset.get_mut(j) {
+                *o = h;
+            }
+            cols.push(Map { var: j, sign: -1.0 });
+        } else {
+            cols.push(Map { var: j, sign: 1.0 });
+            cols.push(Map { var: j, sign: -1.0 });
+        }
+    }
+    let ny = cols.len();
+    // Rows over y: (coefficients, sense, rhs).
+    let mut yrows: Vec<(Vec<f64>, Sense, f64)> = Vec::with_capacity(rows.len() + extra.len());
+    for (a, sense, b) in rows {
+        let mut coef = vec![0.0; ny];
+        for (k, m) in cols.iter().enumerate() {
+            if let Some(v) = coef.get_mut(k) {
+                *v = a.get(m.var).copied().unwrap_or(0.0) * m.sign;
+            }
+        }
+        let shift: f64 = a.iter().zip(&offset).map(|(x, o)| x * o).sum();
+        yrows.push((coef, *sense, b - shift));
+    }
+    for (entries, sense, b) in extra {
+        let mut coef = vec![0.0; ny];
+        for (k, v) in entries {
+            if let Some(e) = coef.get_mut(k) {
+                *e = v;
+            }
+        }
+        yrows.push((coef, sense, b));
+    }
+    let m = yrows.len();
+    let n_slack = yrows.iter().filter(|r| r.1 != Sense::Eq).count();
+    // Rows with b < 0 are negated so every right side is ≥ 0.
+    let normalized: Vec<(Vec<f64>, Sense, f64)> = yrows
+        .into_iter()
+        .map(|(a, s, b)| {
+            if b < 0.0 {
+                let s = match s {
+                    Sense::Le => Sense::Ge,
+                    Sense::Ge => Sense::Le,
+                    Sense::Eq => Sense::Eq,
+                };
+                (a.iter().map(|v| -v).collect(), s, -b)
+            } else {
+                (a, s, b)
+            }
+        })
+        .collect();
+    let n_art = normalized.iter().filter(|r| r.1 != Sense::Le).count();
+    let width = ny + n_slack + n_art;
+    let mut t: Vec<Vec<f64>> = Vec::with_capacity(m);
+    let mut basis = vec![0usize; m];
+    let (mut si, mut ai) = (ny, ny + n_slack);
+    for (i, (a, s, b)) in normalized.iter().enumerate() {
+        let mut row = vec![0.0; width + 1];
+        for (v, x) in row.iter_mut().zip(a) {
+            *v = *x;
+        }
+        match s {
+            Sense::Le => {
+                if let Some(v) = row.get_mut(si) {
+                    *v = 1.0;
+                }
+                if let Some(bi) = basis.get_mut(i) {
+                    *bi = si;
+                }
+                si += 1;
+            }
+            Sense::Ge => {
+                if let Some(v) = row.get_mut(si) {
+                    *v = -1.0;
+                }
+                si += 1;
+                if let Some(v) = row.get_mut(ai) {
+                    *v = 1.0;
+                }
+                if let Some(bi) = basis.get_mut(i) {
+                    *bi = ai;
+                }
+                ai += 1;
+            }
+            Sense::Eq => {
+                if let Some(v) = row.get_mut(ai) {
+                    *v = 1.0;
+                }
+                if let Some(bi) = basis.get_mut(i) {
+                    *bi = ai;
+                }
+                ai += 1;
+            }
+        }
+        if let Some(v) = row.get_mut(width) {
+            *v = *b;
+        }
+        t.push(row);
+    }
+    let is_art = |j: usize| j >= ny + n_slack && j < width;
+    // Phase 1: minimise the sum of the artificials.
+    let mut obj = vec![0.0; width + 1];
+    for j in (ny + n_slack)..width {
+        if let Some(v) = obj.get_mut(j) {
+            *v = 1.0;
+        }
+    }
+    for (i, row) in t.iter().enumerate() {
+        if basis.get(i).is_some_and(|b| is_art(*b)) {
+            for (v, x) in obj.iter_mut().zip(row) {
+                *v -= x;
+            }
+        }
+    }
+    let all: Vec<bool> = vec![true; width];
+    if n_art > 0 {
+        match simplex(&mut t, &mut basis, &mut obj, &all, budget) {
+            Step::Limit => return Outcome::Limit,
+            Step::Unbounded | Step::Optimal => {}
+        }
+        let scale = normalized.iter().map(|r| r.2.abs()).fold(1.0f64, f64::max);
+        if -obj.get(width).copied().unwrap_or(0.0) > 1e-7 * scale {
+            return Outcome::Infeasible;
+        }
+        // Drive zero-valued artificials out of the basis where possible.
+        for i in 0..m {
+            if !basis.get(i).is_some_and(|b| is_art(*b)) {
+                continue;
+            }
+            let j = t.get(i).and_then(|row| (0..ny + n_slack).find(|&j| row.get(j).is_some_and(|v| v.abs() > EPS)));
+            if let Some(j) = j {
+                pivot(&mut t, &mut obj, i, j);
+                if let Some(b) = basis.get_mut(i) {
+                    *b = j;
+                }
+            }
+        }
+    }
+    // Phase 2: the real costs; artificials may not re-enter.
+    let mut cost = vec![0.0; width + 1];
+    for (k, mp) in cols.iter().enumerate() {
+        if let Some(v) = cost.get_mut(k) {
+            *v = c.get(mp.var).copied().unwrap_or(0.0) * mp.sign;
+        }
+    }
+    let mut obj = cost.clone();
+    for (i, row) in t.iter().enumerate() {
+        let cb = basis.get(i).and_then(|b| cost.get(*b)).copied().unwrap_or(0.0);
+        if cb != 0.0 {
+            for (v, x) in obj.iter_mut().zip(row) {
+                *v -= cb * x;
+            }
+        }
+    }
+    let allowed: Vec<bool> = (0..width).map(|j| !is_art(j)).collect();
+    match simplex(&mut t, &mut basis, &mut obj, &allowed, budget) {
+        Step::Limit => return Outcome::Limit,
+        Step::Unbounded => return Outcome::Unbounded,
+        Step::Optimal => {}
+    }
+    let mut y = vec![0.0; width];
+    for (i, row) in t.iter().enumerate() {
+        if let (Some(b), Some(v)) = (basis.get(i), row.get(width))
+            && let Some(e) = y.get_mut(*b)
+        {
+            *e = *v;
+        }
+    }
+    let mut x = offset;
+    for (k, mp) in cols.iter().enumerate() {
+        if let Some(e) = x.get_mut(mp.var) {
+            *e += mp.sign * y.get(k).copied().unwrap_or(0.0);
+        }
+    }
+    let val: f64 = c.iter().zip(&x).map(|(a, b)| a * b).sum();
+    Outcome::Solved(x, val, true)
+}
+
+/// Linearises the model at `x0` (n + 1 evaluations) and checks it at one more point. Returns the
+/// objective gradient and constant, and per row the gradient and constant of `lhs − rhs`.
+#[allow(clippy::type_complexity)]
+fn linearize(ev: &mut Evaluator, x0: &[f64]) -> std::result::Result<((Vec<f64>, f64), Vec<(Vec<f64>, f64)>), &'static str> {
+    let n = x0.len();
+    let p0 = ev.eval(x0).ok_or("budget")?;
+    if !p0.f.is_finite() || p0.g.iter().any(|v| !v.is_finite()) {
+        return Err("error");
+    }
+    let mut cf = vec![0.0; n];
+    let mut ga = vec![vec![0.0; n]; p0.g.len()];
+    for j in 0..n {
+        let mut x = x0.to_vec();
+        if let Some(v) = x.get_mut(j) {
+            *v += 1.0;
+        }
+        let p = ev.eval(&x).ok_or("budget")?;
+        if let Some(v) = cf.get_mut(j) {
+            *v = p.f - p0.f;
+        }
+        for (row, (gi, g0)) in ga.iter_mut().zip(p.g.iter().zip(&p0.g)) {
+            if let Some(v) = row.get_mut(j) {
+                *v = gi - g0;
+            }
+        }
+    }
+    let dot = |a: &[f64], x: &[f64]| a.iter().zip(x).map(|(p, q)| p * q).sum::<f64>();
+    let f_const = p0.f - dot(&cf, x0);
+    let rows: Vec<(Vec<f64>, f64)> = ga.into_iter().zip(&p0.g).map(|(a, g0)| (a.clone(), g0 - dot(&a, x0))).collect();
+    // Linearity check at a spread-out point.
+    let xt: Vec<f64> =
+        x0.iter().enumerate().map(|(j, v)| v + if j % 2 == 0 { 1.7 + 0.3 * (j % 5) as f64 } else { -0.6 - 0.2 * (j % 3) as f64 }).collect();
+    let pt = ev.eval(&xt).ok_or("budget")?;
+    let close = |actual: f64, pred: f64, terms: f64| actual.is_finite() && (actual - pred).abs() <= 1e-6 * (1.0 + actual.abs() + terms);
+    let terms = |a: &[f64]| a.iter().zip(&xt).map(|(p, q)| (p * q).abs()).sum::<f64>();
+    if !close(pt.f, f_const + dot(&cf, &xt), terms(&cf)) {
+        return Err("nonlinear");
+    }
+    for ((a, k), g) in rows.iter().zip(&pt.g) {
+        if !close(*g, k + dot(a, &xt), terms(a)) {
+            return Err("nonlinear");
+        }
+    }
+    Ok(((cf, f_const), rows))
+}
+
+// ---------------------------------------------------------------- nonlinear
+
+/// Augmented-Lagrangian value of a point: objective `phi` plus the constraint terms.
+fn aug(phi: f64, cons: &[f64], senses: &[Sense], lam: &[f64], mu: f64) -> f64 {
+    let mut v = phi;
+    for ((c, s), l) in cons.iter().zip(senses).zip(lam) {
+        v += match s {
+            Sense::Eq => l * c + 0.5 * mu * c * c,
+            _ => {
+                let t = (l + mu * c).max(0.0);
+                (t * t - l * l) / (2.0 * mu)
+            }
+        };
+    }
+    v
+}
+
+/// Constraint values in the `≤ 0` / `= 0` form, and the objective to minimise.
+fn shape(p: &Point, rows: &[Row], goal: Goal) -> (f64, Vec<f64>) {
+    let phi = match goal {
+        Goal::Max => -p.f,
+        Goal::Min => p.f,
+        Goal::Value(v) => (p.f - v) * (p.f - v),
+    };
+    let cons = p.g.iter().zip(rows).map(|(g, r)| if r.sense == Sense::Ge { -g } else { *g }).collect();
+    (phi, cons)
+}
+
+fn violation(cons: &[f64], senses: &[Sense]) -> f64 {
+    cons.iter().zip(senses).map(|(c, s)| if *s == Sense::Eq { c.abs() } else { c.max(0.0) }).fold(0.0, f64::max)
+}
+
+fn project(x: &mut [f64], lo: &[f64], hi: &[f64]) {
+    for ((v, l), h) in x.iter_mut().zip(lo).zip(hi) {
+        if *v < *l {
+            *v = *l;
+        }
+        if *v > *h {
+            *v = *h;
+        }
+    }
+}
+
+/// Local minimisation of the model within bounds `lo..hi` from `x0`.
+fn nlp(ev: &mut Evaluator, x0: &[f64], lo: &[f64], hi: &[f64], goal: Goal) -> Outcome {
+    let n = x0.len();
+    let rows = ev.rows.clone();
+    let senses: Vec<Sense> = rows.iter().map(|r| r.sense).collect();
+    let mut x = x0.to_vec();
+    project(&mut x, lo, hi);
+    let mut lam = vec![0.0; rows.len()];
+    let mut mu = 10.0;
+    let Some(p) = ev.eval(&x) else { return Outcome::Limit };
+    let (mut phi, mut cons) = shape(&p, &rows, goal);
+    if !phi.is_finite() {
+        return Outcome::Infeasible;
+    }
+    let phi0 = phi;
+    let mut prev_viol = violation(&cons, &senses);
+    let mut converged = false;
+    let mut out_of_budget = false;
+    let mut carried: Option<f64> = None;
+    for _outer in 0..40 {
+        // Inner: projected gradient on the augmented Lagrangian.
+        let lval = |phi: f64, cons: &[f64], lam: &[f64], mu: f64| aug(phi, cons, &senses, lam, mu);
+        let grad = |ev: &mut Evaluator, x: &[f64], lam: &[f64], mu: f64| -> Option<Vec<f64>> {
+            let mut gr = vec![0.0; n];
+            for j in 0..n {
+                let xj = x.get(j).copied().unwrap_or(0.0);
+                let h = 1e-6 * xj.abs().max(1.0);
+                let (l, u) = (lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), hi.get(j).copied().unwrap_or(f64::INFINITY));
+                let (a, b) = (if xj - h >= l { xj - h } else { xj }, if xj + h <= u { xj + h } else { xj });
+                if b <= a {
+                    continue;
+                }
+                let mut xa = x.to_vec();
+                let mut xb = x.to_vec();
+                if let Some(v) = xa.get_mut(j) {
+                    *v = a;
+                }
+                if let Some(v) = xb.get_mut(j) {
+                    *v = b;
+                }
+                let (pa, ca) = shape(&ev.eval(&xa)?, &rows, goal);
+                let (pb, cb) = shape(&ev.eval(&xb)?, &rows, goal);
+                let d = (lval(pb, &cb, lam, mu) - lval(pa, &ca, lam, mu)) / (b - a);
+                if let Some(v) = gr.get_mut(j) {
+                    *v = if d.is_finite() { d } else { 0.0 };
+                }
+            }
+            Some(gr)
+        };
+        let mut l_cur = lval(phi, &cons, &lam, mu);
+        let Some(mut g) = grad(ev, &x, &lam, mu) else {
+            out_of_budget = true;
+            break;
+        };
+        // The step length carries over between outer iterations once known.
+        let gnorm = g.iter().map(|v| v.abs()).fold(0.0, f64::max);
+        let mut alpha = carried.unwrap_or(if gnorm > 0.0 { 1.0 / gnorm } else { 1.0 });
+        let mut inner_done = false;
+        let mut quiet = 0;
+        for _inner in 0..500 {
+            let mut accepted = None;
+            let mut a = alpha;
+            for _ in 0..40 {
+                let mut xn: Vec<f64> = x.iter().zip(&g).map(|(v, d)| v - a * d).collect();
+                project(&mut xn, lo, hi);
+                let Some(pn) = ev.eval(&xn) else { break };
+                let (phn, cn) = shape(&pn, &rows, goal);
+                let ln = lval(phn, &cn, &lam, mu);
+                let decrease: f64 = g.iter().zip(x.iter().zip(&xn)).map(|(d, (u, v))| d * (u - v)).sum();
+                if ln.is_finite() && ln <= l_cur - 1e-4 * decrease {
+                    accepted = Some((xn, phn, cn, ln));
+                    break;
+                }
+                a *= 0.25;
+            }
+            if ev.exhausted() {
+                out_of_budget = true;
+            }
+            let Some((xn, phn, cn, ln)) = accepted else {
+                inner_done = true;
+                break;
+            };
+            // Variables or objective running off towards infinity: unbounded.
+            if phn.abs() > 1e15 * (1.0 + phi0.abs()) || xn.iter().zip(x0).any(|(v, s)| v.abs() > 1e15 * (1.0 + s.abs())) {
+                return Outcome::Unbounded;
+            }
+            let s: Vec<f64> = xn.iter().zip(&x).map(|(a, b)| a - b).collect();
+            let step = s.iter().map(|v| v.abs()).fold(0.0, f64::max);
+            let xmax = xn.iter().map(|v| v.abs()).fold(1.0, f64::max);
+            let small = step <= 1e-10 * xmax || (l_cur - ln).abs() <= 1e-13 * (1.0 + l_cur.abs());
+            quiet = if small { quiet + 1 } else { 0 };
+            x = xn;
+            phi = phn;
+            cons = cn;
+            l_cur = ln;
+            if quiet >= 3 {
+                inner_done = true;
+                break;
+            }
+            let Some(gn) = grad(ev, &x, &lam, mu) else {
+                out_of_budget = true;
+                break;
+            };
+            // Barzilai–Borwein step.
+            let yv: Vec<f64> = gn.iter().zip(&g).map(|(a, b)| a - b).collect();
+            let sy: f64 = s.iter().zip(&yv).map(|(a, b)| a * b).sum();
+            let ss: f64 = s.iter().map(|v| v * v).sum();
+            alpha = if sy > 0.0 { (ss / sy).clamp(1e-12, 1e12) } else { (a * 4.0).min(1e12) };
+            g = gn;
+            if g.iter().all(|v| *v == 0.0) {
+                inner_done = true;
+                break;
+            }
+        }
+        carried = Some(alpha);
+        if out_of_budget {
+            break;
+        }
+        let viol = violation(&cons, &senses);
+        let scale = rows.iter().map(|r| r.rhs_const.abs()).fold(1.0, f64::max);
+        if viol <= PRECISION * scale && inner_done {
+            converged = true;
+            if rows.is_empty() {
+                break;
+            }
+        }
+        // Multiplier and penalty updates.
+        for ((l, c), s) in lam.iter_mut().zip(&cons).zip(&senses) {
+            *l = match s {
+                Sense::Eq => *l + mu * c,
+                _ => (*l + mu * c).max(0.0),
+            };
+        }
+        if viol > 0.25 * prev_viol {
+            mu = (mu * 10.0).min(1e10);
+        }
+        if converged && viol <= PRECISION * scale * 1e-2 {
+            break;
+        }
+        prev_viol = viol;
+        converged = false;
+    }
+    let viol = violation(&cons, &senses);
+    let scale = rows.iter().map(|r| r.rhs_const.abs()).fold(1.0, f64::max);
+    if viol > PRECISION * scale * 10.0 {
+        return if out_of_budget { Outcome::Limit } else { Outcome::Infeasible };
+    }
+    Outcome::Solved(x, phi, (converged || viol <= PRECISION * scale) && !out_of_budget)
+}
+
+// ---------------------------------------------------------------- command
+
+fn run(s: &mut Session, p: &Json) -> Result<Json> {
+    let d = s.doc()?;
+    let active = d.wb.active_sheet;
+    let mut params = p.clone();
+    if bool_param(p, "saved").unwrap_or(false) {
+        let saved = saved_model(&d.wb, active).ok_or_else(|| EngineError::Other("This sheet has no saved Solver model.".into()))?;
+        if let (Some(base), Some(over)) = (saved.as_object(), p.as_object()) {
+            let mut merged = base.clone();
+            for (k, v) in over {
+                merged.insert(k.clone(), v.clone());
+            }
+            params = Json::Object(merged);
+        }
+    }
+    let has_objective = ["objective", "setObjective", "set"].iter().any(|k| str_param(&params, k).is_some_and(|t| !t.trim().is_empty()));
+    if !has_objective {
+        let model = saved_model(&d.wb, active);
+        s.ui_requests.push(crate::UiRequest::Dialog("solver".into(), json!({})));
+        return Ok(json!({"model": model}));
+    }
+    let model = parse_model(s, &params)?;
+    let budget = params
+        .get("maxEvaluations")
+        .and_then(Json::as_f64)
+        .filter(|v| v.is_finite() && *v >= 1.0)
+        .map_or(DEFAULT_EVALS, |v| (v as u64).min(MAX_EVALS));
+    let gap = params.get("integerOptimality").and_then(Json::as_f64).filter(|v| v.is_finite()).map_or(0.0, |v| v.clamp(0.0, 100.0) / 100.0);
+    let keep = bool_param(&params, "keep").unwrap_or(true);
+    let save = bool_param(&params, "saveModel").unwrap_or(true);
+    let prob = build_problem(&model);
+    let d = s.doc()?;
+    let sh = d.wb.sheet(model.sheet).ok_or(EngineError::NoDocument)?;
+    let x0: Vec<f64> = model.vars.iter().map(|c| as_num(&sh.value(*c))).map(|v| if v.is_finite() { v } else { 0.0 }).collect();
+    let mut ev = Evaluator {
+        trial: Trial::new(d),
+        sheet: model.sheet,
+        vars: model.vars.clone(),
+        objective: model.objective,
+        rows: prob.rows.clone(),
+        evals: 0,
+        budget,
+    };
+    let any_int = prob.integer.iter().any(|b| *b);
+    let mut not_linear = false;
+    let outcome = if prob.lo.iter().zip(&prob.hi).any(|(l, h)| l > h) {
+        Outcome::Infeasible
+    } else {
+        match model.method {
+            Method::Simplex => match linearize(&mut ev, &x0) {
+                Err("nonlinear") | Err("error") => {
+                    not_linear = true;
+                    Outcome::Infeasible
+                }
+                Err(_) => Outcome::Limit,
+                Ok(((cf, f_const), lin)) => {
+                    let mut rows: Vec<(Vec<f64>, Sense, f64)> = lin.iter().zip(&prob.rows).map(|((a, k), r)| (a.clone(), r.sense, -k)).collect();
+                    let c: Vec<f64> = match model.goal {
+                        Goal::Max => cf.iter().map(|v| -v).collect(),
+                        Goal::Min => cf.clone(),
+                        Goal::Value(v) => {
+                            rows.push((cf.clone(), Sense::Eq, v - f_const));
+                            vec![0.0; cf.len()]
+                        }
+                    };
+                    let size = (rows.len() as u64 + 2 * model.vars.len() as u64 + 1) * (4 * model.vars.len() as u64 + 2 * rows.len() as u64 + 1);
+                    let mut work = LP_WORK.max(size);
+                    let mut relax = |l: &[f64], h: &[f64]| solve_lp(&c, &rows, l, h, &mut work);
+                    if any_int {
+                        branch_and_bound(&prob.lo, &prob.hi, &prob.integer, MAX_NODES_LP, gap, &mut relax).0
+                    } else {
+                        relax(&prob.lo, &prob.hi)
+                    }
+                }
+            },
+            Method::Grg => {
+                let goal = model.goal;
+                if any_int {
+                    let mut relax = |l: &[f64], h: &[f64]| {
+                        let mut start = x0.clone();
+                        project(&mut start, l, h);
+                        nlp(&mut ev, &start, l, h, goal)
+                    };
+                    branch_and_bound(&prob.lo, &prob.hi, &prob.integer, MAX_NODES_NLP, gap, &mut relax).0
+                } else {
+                    nlp(&mut ev, &x0, &prob.lo, &prob.hi, goal)
+                }
+            }
+        }
+    };
+    let (status, message, solution) = match outcome {
+        _ if not_linear => ("notLinear", "The model isn't linear in the variable cells, so Simplex LP can't solve it. Try GRG Nonlinear.", None),
+        Outcome::Solved(x, _, true) => (
+            if model.method == Method::Simplex { "optimal" } else { "converged" },
+            if model.method == Method::Simplex {
+                "Solver found an optimal solution. All constraints are satisfied."
+            } else {
+                "Solver converged to a locally optimal solution. All constraints are satisfied."
+            },
+            Some(x),
+        ),
+        Outcome::Solved(x, _, false) => {
+            ("iterationLimit", "Solver stopped at its evaluation limit; the best solution found satisfies the constraints.", Some(x))
+        }
+        Outcome::Limit => ("iterationLimit", "Solver stopped at its evaluation limit without finding a feasible solution.", None),
+        Outcome::Infeasible => ("infeasible", "Solver could not find a solution that satisfies all constraints.", None),
+        Outcome::Unbounded => ("unbounded", "The objective can grow without limit: the problem is unbounded.", None),
+    };
+    let evaluations = ev.evals;
+    // The objective's value at the solution.
+    let objective = solution.as_ref().and_then(|x| {
+        ev.budget = ev.budget.saturating_add(1);
+        ev.eval(x).map(|p| p.f)
+    });
+    let mut values = serde_json::Map::new();
+    if let Some(x) = &solution {
+        for (c, v) in model.vars.iter().zip(x) {
+            values.insert(c.a1(), json!(v));
+        }
+    }
+    let changed = keep && solution.is_some();
+    if changed || save {
+        let sheet = model.sheet;
+        let vars = model.vars.clone();
+        let m = model.clone();
+        edit(s, |cx| {
+            if changed && let Some(x) = &solution {
+                let sh = cx.sheet_mut(sheet)?;
+                for (c, v) in vars.iter().zip(x) {
+                    let style = sh.style_id(*c);
+                    sh.set_cell(*c, Cell { value: Value::number(*v), formula: None, style });
+                }
+                for c in &vars {
+                    cx.touch(sheet, *c);
+                }
+            }
+            if save {
+                save_model(&mut cx.wb, &m);
+            }
+            Ok(())
+        })?;
+    }
+    Ok(json!({
+        "status": status,
+        "message": message,
+        "objective": objective,
+        "values": values,
+        "evaluations": evaluations,
+        "changed": changed,
+    }))
+}
+
+#[cfg(test)]
+mod tests {
+    use gridcraft_core::CellRef;
+    use serde_json::json;
+
+    use crate::Session;
+
+    fn s() -> Session {
+        let mut s = Session::new();
+        s.new_workbook();
+        s
+    }
+    fn n(s: &Session, a: &str) -> f64 {
+        let d = s.doc().unwrap();
+        d.wb.active().unwrap().value(CellRef::parse(a).unwrap()).as_f64().unwrap_or(f64::NAN)
+    }
+    fn set(s: &mut Session, cells: &[(&str, &str)]) {
+        for (c, v) in cells {
+            s.execute("cell.set", json!({"cell": c, "input": v})).unwrap();
+        }
+    }
+
+    #[test]
+    fn simplex_solves_the_wyndor_glass_lp() {
+        // Hillier & Lieberman's Wyndor Glass problem: max 3x + 5y, x ≤ 4, 2y ≤ 12, 3x + 2y ≤ 18
+        // → x = 2, y = 6, Z = 36.
+        let mut s = s();
+        set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=3*A1+5*A2"), ("B2", "=2*A2"), ("B3", "=3*A1+2*A2")]);
+        let r = s
+            .execute(
+                "data.solver",
+                json!({"objective": "B1", "goal": "max", "variables": "A1:A2", "constraints": ["A1 <= 4", "B2 <= 12", {"lhs": "B3", "op": "<=", "rhs": 18}], "method": "simplex"}),
+            )
+            .unwrap();
+        assert_eq!(r["status"], "optimal", "{r}");
+        assert!((n(&s, "A1") - 2.0).abs() < 1e-9 && (n(&s, "A2") - 6.0).abs() < 1e-9);
+        assert!((n(&s, "B1") - 36.0).abs() < 1e-9);
+        // The model is saved as hidden names, survives XLSX and solves again from them.
+        let names: Vec<String> = s.doc().unwrap().wb.names.iter().map(|n| n.name.clone()).collect();
+        assert!(names.iter().any(|n| n == "solver_opt") && names.iter().any(|n| n == "solver_lhs3"));
+        let b64 = s.execute("file.saveBytes", json!({"format": "xlsx"})).unwrap()["base64"].as_str().unwrap().to_string();
+        s.execute("file.open", json!({"name": "m.xlsx", "base64": b64})).unwrap();
+        set(&mut s, &[("A1", "0"), ("A2", "0")]);
+        let m = s.execute("data.solver", json!({})).unwrap();
+        assert_eq!(m["model"]["objective"], "B1");
+        assert_eq!(m["model"]["method"], "simplex");
+        let r = s.execute("data.solver", json!({"saved": true})).unwrap();
+        assert_eq!(r["status"], "optimal");
+        assert!((n(&s, "B1") - 36.0).abs() < 1e-9);
+        // Minimise with keep: false leaves the cells alone.
+        set(&mut s, &[("A1", "1"), ("A2", "1")]);
+        let r = s.execute("data.solver", json!({"saved": true, "goal": "min", "keep": false})).unwrap();
+        assert_eq!(r["objective"], 0.0);
+        assert_eq!((n(&s, "A1"), n(&s, "A2")), (1.0, 1.0));
+    }
+
+    #[test]
+    fn integer_and_binary_constraints_use_branch_and_bound() {
+        // max 5x + 4y, 6x + 4y ≤ 24, x + 2y ≤ 6: the LP optimum (3, 1.5) = 21; integers give
+        // (4, 0) = 20.
+        let mut s = s();
+        set(&mut s, &[("B1", "=5*A1+4*A2"), ("B2", "=6*A1+4*A2"), ("B3", "=A1+2*A2")]);
+        let base = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["B2 <= 24", "B3 <= 6"], "method": "simplex"});
+        let r = s.execute("data.solver", base.clone()).unwrap();
+        assert!((r["objective"].as_f64().unwrap() - 21.0).abs() < 1e-9);
+        let mut p = base.clone();
+        p["constraints"].as_array_mut().unwrap().push(json!("A1:A2 int"));
+        let r = s.execute("data.solver", p).unwrap();
+        assert_eq!(r["status"], "optimal");
+        assert_eq!((n(&s, "A1"), n(&s, "A2"), n(&s, "B1")), (4.0, 0.0, 20.0));
+        // 0/1 knapsack: values 10, 13, 7, weights 5, 6, 4, capacity 10 → items 2 and 3 (20).
+        set(&mut s, &[("D1", "=10*C1+13*C2+7*C3"), ("D2", "=5*C1+6*C2+4*C3")]);
+        let r = s
+            .execute("data.solver", json!({"objective": "D1", "variables": "C1:C3", "constraints": ["D2 <= 10", "C1:C3 bin"], "method": "simplex"}))
+            .unwrap();
+        assert_eq!(r["status"], "optimal");
+        assert_eq!((n(&s, "C1"), n(&s, "C2"), n(&s, "C3"), n(&s, "D1")), (0.0, 1.0, 1.0, 20.0));
+    }
+
+    #[test]
+    fn grg_finds_smooth_optima() {
+        // Largest rectangle with perimeter at most 20: a 5 × 5 square, area 25.
+        let mut s = s();
+        set(&mut s, &[("A1", "1"), ("A2", "1"), ("B1", "=A1*A2"), ("B2", "=2*A1+2*A2")]);
+        let r = s.execute("data.solver", json!({"objective": "B1", "variables": "A1:A2", "constraints": ["B2 <= 20"]})).unwrap();
+        assert_eq!(r["status"], "converged", "{r}");
+        assert!((n(&s, "A1") - 5.0).abs() < 1e-3 && (n(&s, "A2") - 5.0).abs() < 1e-3, "{r}");
+        // Unconstrained minimum of (x − 3)² + (y + 1)² at (3, −1) once negatives are allowed.
+        set(&mut s, &[("C1", "0"), ("C2", "0"), ("D1", "=(C1-3)^2+(C2+1)^2")]);
+        s.execute("data.solver", json!({"objective": "D1", "goal": "min", "variables": "C1:C2", "nonNegative": false})).unwrap();
+        assert!((n(&s, "C1") - 3.0).abs() < 1e-4 && (n(&s, "C2") + 1.0).abs() < 1e-4);
+        // Value Of: x³ = 8 → x = 2.
+        set(&mut s, &[("E1", "1"), ("F1", "=E1^3")]);
+        s.execute("data.solver", json!({"objective": "F1", "goal": "value", "value": 8, "variables": "E1"})).unwrap();
+        assert!((n(&s, "E1") - 2.0).abs() < 1e-4);
+    }
+
+    #[test]
+    fn infeasible_unbounded_and_nonlinear_models_are_reported() {
+        let mut s = s();
+        set(&mut s, &[("A1", "1"), ("B1", "=2*A1")]);
+        for method in ["simplex", "grg"] {
+            let r = s
+                .execute("data.solver", json!({"objective": "B1", "variables": "A1", "constraints": ["A1 >= 5", "B1 <= 6"], "method": method}))
+                .unwrap();
+            assert_eq!(r["status"], "infeasible", "{method}: {r}");
+            assert_eq!(r["changed"], false);
+            assert_eq!(n(&s, "A1"), 1.0, "an infeasible model leaves the cells alone");
+        }
+        let r = s.execute("data.solver", json!({"objective": "B1", "variables": "A1", "method": "simplex"})).unwrap();
+        assert_eq!(r["status"], "unbounded");
+        let r = s.execute("data.solver", json!({"objective": "B1", "variables": "A1", "method": "grg", "maxEvaluations": 5000})).unwrap();
+        assert!(r["status"] == "unbounded" || r["status"] == "iterationLimit", "{r}");
+        set(&mut s, &[("A1", "1"), ("C1", "=A1^2")]);
+        let r = s.execute("data.solver", json!({"objective": "C1", "goal": "min", "variables": "A1", "method": "simplex"})).unwrap();
+        assert_eq!(r["status"], "notLinear");
+    }
+
+    #[test]
+    fn hostile_models_are_refused_without_panicking() {
+        let mut s = s();
+        let r = s.execute("data.solver", json!({})).unwrap();
+        assert!(r["model"].is_null());
+        assert!(s.take_ui_requests().iter().any(|q| matches!(q, crate::UiRequest::Dialog(n, _) if n == "solver")));
+        set(&mut s, &[("B1", "=SUM(A1:A1000)"), ("C1", "5")]);
+        for p in [
+            json!({"objective": "B1", "variables": "A1:A1000"}),
+            json!({"objective": "B1", "variables": "A1:XFD1048576"}),
+            json!({"objective": "C1", "variables": "A1"}),
+            json!({"objective": "B1:B2", "variables": "A1"}),
+            json!({"objective": "B1", "variables": "B1"}),
+            json!({"objective": "B1", "variables": "A1", "goal": "sideways"}),
+            json!({"objective": "B1", "variables": "A1", "goal": "value"}),
+            json!({"objective": "B1", "variables": "A1", "method": "evolutionary"}),
+            json!({"objective": "B1", "variables": "A1", "constraints": ["A1 ?? 3"]}),
+            json!({"objective": "B1", "variables": "A1", "constraints": "A1 <= 3"}),
+            json!({"objective": "B1", "variables": "A1", "constraints": ["A1:A3 <= C1:C2"]}),
+            json!({"objective": "B1", "variables": "A1", "constraints": ["C1 int"]}),
+            json!({"objective": "B1", "variables": "A1", "constraints": ["A1:A1048576 <= 1"]}),
+            json!({"objective": "B1", "variables": "A1", "constraints": vec![json!("A1 <= 1"); 1001]}),
+            json!({"objective": "Nope!B1", "variables": "A1"}),
+        ] {
+            assert!(s.execute("data.solver", p.clone()).is_err(), "{p}");
+        }
+        // Huge evaluation requests are capped, and a model erroring everywhere doesn't panic.
+        set(&mut s, &[("D1", "=1/0+A1")]);
+        let r = s.execute("data.solver", json!({"objective": "D1", "variables": "A1", "maxEvaluations": 1e300})).unwrap();
+        assert_ne!(r["status"], "converged");
+    }
+}
diff --git a/crates/ui-egui/src/dialogs.rs b/crates/ui-egui/src/dialogs.rs
index 7dd3acb..d67bc1a 100644
--- a/crates/ui-egui/src/dialogs.rs
+++ b/crates/ui-egui/src/dialogs.rs
@@ -317,6 +317,23 @@ impl Dialog {
                 json!({"header": true}),
             ),
             "goalSeek" => Dialog::custom("goalSeek", "Goal Seek", json!({"set": active, "to": "0", "changing": ""})),
+            "solver" => {
+                // Starts from the sheet's saved model, else the active cell as objective.
+                let saved = app.session.active().and_then(|doc| gridcraft_engine::cmd::solver::saved_model(&doc.wb, doc.wb.active_sheet));
+                let mut values = json!({"objective": active, "goal": "max", "value": "0", "variables": "", "constraints": [], "nonNegative": true, "method": "grg"});
+                if let (Some(m), Some(v)) = (saved.as_ref().and_then(Json::as_object), values.as_object_mut()) {
+                    for (k, x) in m {
+                        v.insert(k.clone(), if k == "value" { json!(x.to_string()) } else { x.clone() });
+                    }
+                }
+                if let Some(v) = values.as_object_mut() {
+                    v.insert("newOp".into(), json!("<="));
+                }
+                let mut d = Dialog::custom("solver", "Solver Parameters", values);
+                d.command = Some("data.solver");
+                d
+            }
+            "solverResults" => Dialog::custom("solverResults", "Solver Results", params.clone()),
             "subtotal" => Dialog::form(
                 "subtotal",
                 "Subtotal",
@@ -659,6 +676,7 @@ pub fn show(app: &mut SheetApp, ctx: &egui::Context) {
         "formatCells" => 560.0,
         "about" => 660.0,
         "insertFunction" | "commandSearch" | "nameManager" | "manageRules" | "journal" | "agents" => 520.0,
+        "solver" => 460.0,
         _ => 380.0,
     };
     let window = egui::Window::new(d.title.clone())
@@ -681,6 +699,8 @@ pub fn show(app: &mut SheetApp, ctx: &egui::Context) {
                 "manageRules" => manage_rules(app, ui),
                 "commandSearch" => command_search(app, ui, &mut d, &mut confirm),
                 "goalSeek" => goal_seek(app, ui, &mut d, &mut confirm),
+                "solver" => solver(ui, &mut d, &mut confirm, &mut open),
+                "solverResults" => solver_results(app, ui, &mut d, &mut open),
                 "spelling" => spelling(app, ui, &mut d, &mut open),
                 "agents" => {
                     ui.label("GridCraft is fully drivable by agents. Every menu item, button and gesture is a command:");
@@ -877,7 +897,12 @@ pub fn show(app: &mut SheetApp, ctx: &egui::Context) {
     if confirm && let Some(cmd) = d.command {
         let params = build_params(&d);
         match app.run(cmd, params) {
-            Ok(_) => open = false,
+            Ok(r) => {
+                open = false;
+                if d.name == "solver" {
+                    app.open_dialog("solverResults", r);
+                }
+            }
             Err(e) => app.message = Some((d.title.clone(), crate::clean_error(&e))),
         }
     } else if confirm {
@@ -1607,6 +1632,144 @@ fn goal_seek(app: &mut SheetApp, ui: &mut egui::Ui, d: &mut Dialog, confirm: &mu
     }
 }
 
+/// Data › Solver: the Solver Parameters form. Solve runs `data.solver` through the dialog's
+/// command; the result opens Solver Results.
+fn solver(ui: &mut egui::Ui, d: &mut Dialog, confirm: &mut bool, open: &mut bool) {
+    let text = |d: &Dialog, k: &str| d.values.get(k).map(|v| v.as_str().map(str::to_string).unwrap_or_else(|| v.to_string())).unwrap_or_default();
+    let edit = |ui: &mut egui::Ui, d: &mut Dialog, k: &str, width: f32| {
+        let mut t = text(d, k);
+        if ui.add(egui::TextEdit::singleline(&mut t).desired_width(width)).changed() {
+            d.values.insert(k.into(), json!(t));
+        }
+    };
+    ui.horizontal(|ui| {
+        ui.label("Set Objective:");
+        edit(ui, d, "objective", 200.0);
+    });
+    let goal = text(d, "goal");
+    ui.horizontal(|ui| {
+        ui.label("To:");
+        for (g, label) in [("max", "Max"), ("min", "Min"), ("value", "Value Of:")] {
+            if ui.radio(goal == g, label).clicked() {
+                d.values.insert("goal".into(), json!(g));
+            }
+        }
+        edit(ui, d, "value", 80.0);
+    });
+    ui.label("By Changing Variable Cells:");
+    edit(ui, d, "variables", 420.0);
+    ui.label("Subject to the Constraints:");
+    let list: Vec<String> = d
+        .values
+        .get("constraints")
+        .and_then(Json::as_array)
+        .map(|a| a.iter().map(|c| c.as_str().map(str::to_string).unwrap_or_else(|| c.to_string())).collect())
+        .unwrap_or_default();
+    egui::Frame::group(ui.style()).show(ui, |ui| {
+        ui.set_min_size(vec2(420.0, 90.0));
+        egui::ScrollArea::vertical().max_height(110.0).show(ui, |ui| {
+            for (i, c) in list.iter().enumerate() {
+                if ui.selectable_label(d.list_index == i, c).clicked() {
+                    d.list_index = i;
+                }
+            }
+        });
+    });
+    ui.horizontal(|ui| {
+        edit(ui, d, "newLhs", 110.0);
+        let op = text(d, "newOp");
+        egui::ComboBox::from_id_salt("solverOp").selected_text(op.clone()).width(60.0).show_ui(ui, |ui| {
+            for o in ["<=", "=", ">=", "int", "bin"] {
+                if ui.selectable_label(op == o, o).clicked() {
+                    d.values.insert("newOp".into(), json!(o));
+                }
+            }
+        });
+        if !matches!(op.as_str(), "int" | "bin") {
+            edit(ui, d, "newRhs", 110.0);
+        }
+        if ui.button("Add").clicked() {
+            let lhs = text(d, "newLhs");
+            if !lhs.trim().is_empty() {
+                let rhs = if matches!(op.as_str(), "int" | "bin") { String::new() } else { format!(" {}", text(d, "newRhs").trim()) };
+                let mut list = list.clone();
+                list.push(format!("{} {op}{rhs}", lhs.trim()));
+                d.values.insert("constraints".into(), json!(list));
+                d.values.insert("newLhs".into(), json!(""));
+                d.values.insert("newRhs".into(), json!(""));
+            }
+        }
+        if ui.add_enabled(d.list_index < list.len(), egui::Button::new("Delete")).clicked() {
+            let mut list = list.clone();
+            if d.list_index < list.len() {
+                list.remove(d.list_index);
+            }
+            d.values.insert("constraints".into(), json!(list));
+        }
+    });
+    let mut nonneg = d.values.get("nonNegative").and_then(Json::as_bool).unwrap_or(true);
+    if ui.checkbox(&mut nonneg, "Make Unconstrained Variables Non-Negative").changed() {
+        d.values.insert("nonNegative".into(), json!(nonneg));
+    }
+    let method = text(d, "method");
+    ui.horizontal(|ui| {
+        ui.label("Select a Solving Method:");
+        let label = |m: &str| if m == "simplex" { "Simplex LP" } else { "GRG Nonlinear" };
+        egui::ComboBox::from_id_salt("solverMethod").selected_text(label(&method)).show_ui(ui, |ui| {
+            for m in ["grg", "simplex"] {
+                if ui.selectable_label(method == m, label(m)).clicked() {
+                    d.values.insert("method".into(), json!(m));
+                }
+            }
+        });
+    });
+    ui.label(
+        egui::RichText::new(if method == "simplex" {
+            "Simplex LP: for linear models; integer and binary constraints by branch and bound."
+        } else {
+            "GRG Nonlinear: a local optimum of a smooth model, from the current values."
+        })
+        .small(),
+    );
+    ui.add_space(8.0);
+    ui.with_layout(egui::Layout::right_to_left(egui::Align::Center), |ui| {
+        let solve =
+            ui.add(egui::Button::new(egui::RichText::new("  Solve  ").color(Color32::WHITE)).fill(crate::theme::Tokens::get(ui.ctx()).accent));
+        if solve.clicked() && !text(d, "objective").trim().is_empty() {
+            *confirm = true;
+        }
+        if ui.button(" Close ").clicked() {
+            *open = false;
+        }
+    });
+}
+
+/// Solver Results: the outcome, and keeping or restoring the variable cells.
+fn solver_results(app: &mut SheetApp, ui: &mut egui::Ui, d: &mut Dialog, open: &mut bool) {
+    let message = d.values.get("message").and_then(Json::as_str).unwrap_or("").to_string();
+    ui.label(message);
+    let changed = d.values.get("changed").and_then(Json::as_bool).unwrap_or(false);
+    let restore = d.values.get("restore").and_then(Json::as_bool).unwrap_or(false);
+    if changed {
+        ui.add_space(4.0);
+        if ui.radio(!restore, "Keep Solver Solution").clicked() {
+            d.values.insert("restore".into(), json!(false));
+        }
+        if ui.radio(restore, "Restore Original Values").clicked() {
+            d.values.insert("restore".into(), json!(true));
+        }
+    }
+    let mut confirm = false;
+    ok_cancel(ui, &mut confirm, open);
+    if confirm {
+        // The solution (and the saved model) is one undo step.
+        if changed && restore {
+            app.run_or_alert("edit.undo", json!({}));
+        }
+        *open = false;
+    }
+}
+
 fn spelling(app: &mut SheetApp, ui: &mut egui::Ui, d: &mut Dialog, open: &mut bool) {
     let r = d.result.clone().unwrap_or(Json::Null);
     if let Some(e) = r.get("error").and_then(Json::as_str) {
diff --git a/crates/ui-egui/src/i18n.rs b/crates/ui-egui/src/i18n.rs
index 1e7b212..36a20e2 100644
--- a/crates/ui-egui/src/i18n.rs
+++ b/crates/ui-egui/src/i18n.rs
@@ -313,6 +313,7 @@ const JAPANESE: &[(&str, &str)] = &[
     ("Group", "グループ化"),
     ("Ungroup", "グループ解除"),
     ("Subtotal", "小計"),
+    ("Solver", "ソルバー"),
     ("What-If\nAnalysis", "What-If 分析"),
     ("Data Table…", "データ テーブル…"),
     ("Scenario Manager…", "シナリオ マネージャー…"),
diff --git a/crates/ui-egui/src/i18n/ko.rs b/crates/ui-egui/src/i18n/ko.rs
index c37ae3c..13fdc2f 100644
--- a/crates/ui-egui/src/i18n/ko.rs
+++ b/crates/ui-egui/src/i18n/ko.rs
@@ -164,6 +164,7 @@ pub(super) fn translate(text: &str) -> Option<&'static str> {
         "Group" => "그룹",
         "Ungroup" => "그룹 해제",
         "Subtotal" => "부분합",
+        "Solver" => "해 찾기",
         "What-If\nAnalysis" => "가상 분석",
         "Data Table…" => "데이터 표…",
         "Scenario Manager…" => "시나리오 관리자…",
diff --git a/crates/ui-egui/src/i18n/pt.rs b/crates/ui-egui/src/i18n/pt.rs
index 1bdc2c3..b0b30e0 100644
--- a/crates/ui-egui/src/i18n/pt.rs
+++ b/crates/ui-egui/src/i18n/pt.rs
@@ -165,6 +165,7 @@ pub(super) fn translate(text: &str) -> Option<&'static str> {
         "Group" => "Agrupar",
         "Ungroup" => "Desagrupar",
         "Subtotal" => "Subtotal",
+        "Solver" => "Solver",
         "What-If\nAnalysis" => "Análise de\nHipóteses",
         "Data Table…" => "Tabela de Dados…",
         "Scenario Manager…" => "Gerenciador de Cenários…",
diff --git a/crates/ui-egui/src/i18n/ru.rs b/crates/ui-egui/src/i18n/ru.rs
index 8d54856..aaaadbb 100644
--- a/crates/ui-egui/src/i18n/ru.rs
+++ b/crates/ui-egui/src/i18n/ru.rs
@@ -164,6 +164,7 @@ pub(super) fn translate(text: &str) -> Option<&'static str> {
         "Group" => "Группировать",
         "Ungroup" => "Разгруппировать",
         "Subtotal" => "Промежуточный итог",
+        "Solver" => "Поиск решения",
         "What-If\nAnalysis" => "Анализ «что\nесли»",
         "Data Table…" => "Таблица данных…",
         "Scenario Manager…" => "Диспетчер сценариев…",
diff --git a/crates/ui-egui/src/i18n/zh.rs b/crates/ui-egui/src/i18n/zh.rs
index 9f04ada..462462c 100644
--- a/crates/ui-egui/src/i18n/zh.rs
+++ b/crates/ui-egui/src/i18n/zh.rs
@@ -164,6 +164,7 @@ pub(super) fn translate(text: &str) -> Option<&'static str> {
         "Group" => "组合",
         "Ungroup" => "取消组合",
         "Subtotal" => "分类汇总",
+        "Solver" => "规划求解",
         "What-If\nAnalysis" => "假设分析",
         "Data Table…" => "数据表…",
         "Scenario Manager…" => "方案管理器…",
diff --git a/crates/ui-egui/src/ribbon.rs b/crates/ui-egui/src/ribbon.rs
index 264479a..70b24a3 100644
--- a/crates/ui-egui/src/ribbon.rs
+++ b/crates/ui-egui/src/ribbon.rs
@@ -1580,6 +1580,10 @@ fn data(app: &mut SheetApp, ui: &mut Ui) {
             app.open_dialog("subtotal", json!({}));
         }
     });
+    sep(ui);
+    if big_button(ui, Icon::Function, "Solver", "Solver", false).clicked() {
+        app.open_dialog("solver", json!({}));
+    }
 }
 
 fn review(app: &mut SheetApp, ui: &mut Ui) {
diff --git a/crates/ui-egui/tests/interaction.rs b/crates/ui-egui/tests/interaction.rs
index 84f20a6..21f4f04 100644
--- a/crates/ui-egui/tests/interaction.rs
+++ b/crates/ui-egui/tests/interaction.rs
@@ -228,6 +228,7 @@ fn dialogs_open_and_close() {
         "spelling",
         "goTo",
         "about",
+        "solver",
     ] {
         h.state_mut().open_dialog(name, json!({}));
         h.run_steps(2);
@@ -237,6 +238,37 @@ fn dialogs_open_and_close() {
     }
 }
 
+#[test]
+fn solver_dialog_loads_the_saved_model_and_results_can_restore() {
+    let mut s = blank();
+    for (c, v) in [("A1", "0"), ("A2", "0"), ("B1", "=3*A1+5*A2"), ("B2", "=2*A2"), ("B3", "=3*A1+2*A2")] {
+        s.execute("cell.set", json!({"cell": c, "input": v})).unwrap();
+    }
+    let r = s
+        .execute(
+            "data.solver",
+            json!({"objective": "B1", "variables": "A1:A2", "constraints": ["A1 <= 4", "B2 <= 12", "B3 <= 18"], "method": "simplex"}),
+        )
+        .unwrap();
+    assert_eq!(s.active().unwrap().wb.active().unwrap().value(CellRef::parse("B1").unwrap()), Value::Number(36.0));
+    let mut h = harness(s);
+    h.state_mut().open_dialog("solver", json!({}));
+    h.run_steps(2);
+    let d = h.state().dialog.as_ref().unwrap();
+    assert_eq!(d.values.get("objective"), Some(&json!("B1")));
+    assert_eq!(d.values.get("method"), Some(&json!("simplex")));
+    assert_eq!(d.values.get("constraints").and_then(|c| c.as_array()).map(Vec::len), Some(3));
+    h.state_mut().dialog = None;
+    // Solver Results › Restore Original Values undoes the solution.
+    h.state_mut().open_dialog("solverResults", r);
+    h.state_mut().dialog.as_mut().unwrap().values.insert("restore".into(), json!(true));
+    h.run_steps(2);
+    key(&mut h, Key::Enter, Modifiers::NONE);
+    assert!(h.state().dialog.is_none());
+    assert_eq!(value(&h, "A1"), Value::Number(0.0));
+    assert_eq!(value(&h, "B1"), Value::Number(0.0));
+}
+
 #[test]
 fn about_renders_every_tab() {
     let mut h = harness(blank());
diff --git a/docs/gaps.md b/docs/gaps.md
index aa752b8..ef1e08d 100644
--- a/docs/gaps.md
+++ b/docs/gaps.md
@@ -98,8 +98,10 @@ No alpha blockers: every core workflow passes the [alpha gate](roadmap.md#alpha-
   data model (`data.fromWeb`, `data.relationships` missing). **Estimate:** 60–100 h.
 
 ### 13. Analysis tools: Solver, Analysis ToolPak, Forecast Sheet (F)
-- `data.solver`, `data.analysisToolPak`, `data.forecastSheet` missing; FORECAST.ETS family missing.
-  **Estimate:** 22–35 h.
+- `data.analysisToolPak`, `data.forecastSheet` missing; FORECAST.ETS family missing. Solver
+  (`data.solver`) has Simplex LP with integer/binary branch and bound and a local nonlinear method
+  (augmented Lagrangian, not GRG); Evolutionary, AllDifferent constraints, Solver options and
+  reports remain. **Estimate:** 16–28 h.
 
 ### 14. Page Layout view and print polish (U, F)
 - No paper view with rulers and in-place header/footer editing; web print exports only (#141, PR
diff --git a/docs/parity-checklist.md b/docs/parity-checklist.md
index e98835d..623b2b6 100644
--- a/docs/parity-checklist.md
+++ b/docs/parity-checklist.md
@@ -4,7 +4,7 @@
 
 Each entry is an Excel ribbon or menu command (feature names only). It counts as implemented when GridCraft registers a command with the same id. This measures breadth only; the overall assessment (depth, fidelity, effort) is in [target-app-parity.md](target-app-parity.md).
 
-**Overall: 258 / 290 commands implemented (89.0%).**
+**Overall: 259 / 290 commands implemented (89.3%).**
 
 | Tab | Implemented | Total | Parity |
 |---|---|---|---|
@@ -14,7 +14,7 @@ Each entry is an Excel ribbon or menu command (feature names only). It counts as
 | Draw | 4 | 4 | 100% |
 | Page Layout | 23 | 23 | 100% |
 | Formulas | 23 | 23 | 100% |
-| Data | 20 | 28 | 71% |
+| Data | 21 | 28 | 75% |
 | Review | 9 | 18 | 50% |
 | View | 18 | 20 | 90% |
 | Automate | 3 | 3 | 100% |
@@ -240,7 +240,7 @@ Implemented:
 - `formulas.calculateNow` Calculate Now (Calculation)
 - `formulas.calculateSheet` Calculate Sheet (Calculation)
 
-## Data (20 / 28, 71%)
+## Data (21 / 28, 75%)
 
 Implemented:
 
@@ -264,6 +264,7 @@ Implemented:
 - `data.subtotal` Subtotal (Outline)
 - `data.showDetail` Show Detail (Outline)
 - `data.hideDetail` Hide Detail (Outline)
+- `data.solver` Solver (Analysis)
 
 Missing:
 
@@ -273,7 +274,6 @@ Missing:
 - `data.relationships` Relationships (Data Tools)
 - `data.forecastSheet` Forecast Sheet (Forecast)
 - `data.analyzeData` Analyze Data (Analysis)
-- `data.solver` Solver (Analysis)
 - `data.analysisToolPak` Data Analysis (Analysis)
 
 ## Review (9 / 18, 50%)
diff --git a/docs/target-app-parity.md b/docs/target-app-parity.md
index 2e5ed01..b3ec3a6 100644
--- a/docs/target-app-parity.md
+++ b/docs/target-app-parity.md
@@ -87,7 +87,7 @@ performance 8, stability 8, localization 4, hardware 2, platforms 2, ecosystem 2
 | PivotTables | 7 | ~65% | ~45% | 25–40 | Insert/Recommended, field list, three layouts, date grouping, show-values-as, refresh, XLSX round-trip. Missing: slicers, timelines, calculated fields/items, PivotCharts, value/label filters, number grouping, GETPIVOTDATA, data-model pivots. Round-trip never checked in Excel |
 | Charts and sparklines | 8 | ~50% | ~30% | 40–60 | 25 chart kinds, all 2-D. The chart model holds title, legend position, a data-labels flag, gridlines flag and axis titles; no axis scaling, trendlines, error bars, per-point formatting, chart sheets, 3-D, surface, map, Pareto. Several kinds don't survive XLSX round trip (open PRs #158, #161). [chart-parity.md](chart-parity.md) |
 | Page layout, print, PDF | 5 | ~80% | ~55% | 12–20 | Page setup, print titles, breaks, Page Break Preview, headers/footers, PDF export and printing. No true Page Layout view; the PDF writer uses only the 14 standard fonts without embedding (no CJK or other non-Latin text in PDF); web print only exports (#141, #143) |
-| What-if and analysis | 4 | ~50% | ~40% | 22–35 | Goal Seek, Scenario Manager, data tables. Missing: Solver, Analysis ToolPak, Forecast Sheet, Analyze Data |
+| What-if and analysis | 4 | ~50% | ~40% | 22–35 | Goal Seek, Scenario Manager, data tables, Solver (Simplex LP and a local nonlinear method; no Evolutionary or reports). Missing: Analysis ToolPak, Forecast Sheet, Analyze Data |
 | Get & Transform, external data, data model | 5 | ~15% | ~10% | 60–100 | Import of CSV/TSV, JSON, HTML. No Power Query editor, no From Web, database or folder sources, no connections or refresh, no relationships or data model, no external workbook links (dropped with a warning) |
 | Insert objects | 3 | ~55% | ~45% | 25–40 | Pictures, 8 shape kinds (Excel has ~170), text boxes, own icon set, ink with Ink to Shape, hyperlinks, checkboxes. Missing: SmartArt, WordArt, equations, signature line, objects, 3-D models, screenshot, pictures in cells (PR #105) |
 | Review, comments, protection, collaboration | 4 | ~55% | ~40% | 35–55 | Threaded comments, notes, spelling (system dictionary), thesaurus, sheet and workbook protection. Missing: comment navigation, Show Changes, Allow Edit Ranges, notes-to-comments conversion, accessibility checker, co-authoring, opening password-encrypted files (refused with a message since #17) |
diff --git a/docs/ui-parity.md b/docs/ui-parity.md
index e85bca1..71b6602 100644
--- a/docs/ui-parity.md
+++ b/docs/ui-parity.md
@@ -65,7 +65,7 @@ shortcut lists for Windows and Mac, then fill the gaps: 6–10 h.
 
 Format Cells, Insert Function/Formula Builder, Find and Replace, Go To and Go To Special, Sort,
 Name Manager, Data Validation, Paste Special, Conditional Formatting manager, Page Setup, Goal
-Seek, Scenario Manager and others exist. Missing or partial: Options, Solver, Data Analysis,
+Seek, Scenario Manager, Solver and others exist. Missing or partial: Options, Data Analysis,
 Forecast Sheet, Power Query editor, Insert Object, Equation editor, Protect/Allow Edit Ranges,
 Workbook Links, Accessibility checker.