mirror of
https://github.com/storytold/gridcraft.git
synced 2026-11-10 12:29:58 -06:00
Merge origin/main (3379f12a) into the Python cells branch
Main gained #470 (floating objects) and #501 (Solver Sensitivity report). Conflicts: - cmd/edit.rs, edit/paste_special.rs, sheet/tabs.rs: #470 adds rewrite_cell_links for chart references beside this branch's rewrite_python_refs; every structural edit, sheet rename and sheet delete runs both, the Python one last where it can refuse the edit (the edit then changes nothing). - cmd/sheet.rs: rewrite() keeps this branch's Result with main's doc comment. - model/src/sheet.rs and engine tests.rs: both sides' additions kept. - ROADMAP.md: header combined, both sides' rows kept. Main's new code writes no cells, so no Python handling was needed beyond these.
This commit is contained in:
prxnvt
@@ -1,10 +1,14 @@
|
||||
use super::model::Goal;
|
||||
use super::problem::Evaluator;
|
||||
use super::{CellRef, EPS};
|
||||
use super::{CellRef, EPS, LP_WORK};
|
||||
|
||||
// ---------------------------------------------------------------- evaluation
|
||||
|
||||
type BoundRow = (Vec<(usize, f64)>, Sense, f64);
|
||||
|
||||
/// A row of an LP: coefficients, sense and right side.
|
||||
pub(super) type LpRow = (Vec<f64>, Sense, f64);
|
||||
|
||||
/// Objective (gradient, constant) and per-row (gradient, constant).
|
||||
pub(super) type Linearized = ((Vec<f64>, f64), Vec<(Vec<f64>, f64)>);
|
||||
|
||||
@@ -34,9 +38,10 @@ pub(super) enum Outcome {
|
||||
Limit,
|
||||
}
|
||||
|
||||
struct Map {
|
||||
var: usize,
|
||||
sign: f64,
|
||||
/// A column of the simplex tableau standing for a variable: `x_var = offset + Σ sign·y`.
|
||||
pub(super) struct Map {
|
||||
pub(super) var: usize,
|
||||
pub(super) sign: f64,
|
||||
}
|
||||
|
||||
struct VariableMap {
|
||||
@@ -45,6 +50,23 @@ struct VariableMap {
|
||||
bound_rows: Vec<BoundRow>,
|
||||
}
|
||||
|
||||
/// The final tableau of an optimal LP, for the Sensitivity report.
|
||||
pub(super) struct FinalTableau {
|
||||
/// Rows `B⁻¹A | B⁻¹b`.
|
||||
pub(super) t: Vec<Vec<f64>>,
|
||||
pub(super) basis: Vec<usize>,
|
||||
/// The phase-2 reduced costs.
|
||||
pub(super) obj: Vec<f64>,
|
||||
pub(super) width: usize,
|
||||
/// The first artificial column.
|
||||
pub(super) first_art: usize,
|
||||
pub(super) cols: Vec<Map>,
|
||||
/// Per input row: the column that started as its unit column, and whether the row was
|
||||
/// negated to make its right side non-negative.
|
||||
pub(super) ident: Vec<usize>,
|
||||
pub(super) flipped: Vec<bool>,
|
||||
}
|
||||
|
||||
struct Tableau {
|
||||
rows: Vec<Vec<f64>>,
|
||||
basis: Vec<usize>,
|
||||
@@ -82,7 +104,9 @@ fn map_variables(n: usize, lo: &[f64], hi: &[f64]) -> VariableMap {
|
||||
VariableMap { columns, offset, bound_rows }
|
||||
}
|
||||
|
||||
fn normalize_rows(rows: &[(Vec<f64>, Sense, f64)], bound_rows: Vec<BoundRow>, columns: &[Map], offset: &[f64]) -> Vec<(Vec<f64>, Sense, f64)> {
|
||||
/// The rows over the tableau's columns, each with a non-negative right side, and per row
|
||||
/// whether it was negated to get there.
|
||||
fn normalize_rows(rows: &[(Vec<f64>, Sense, f64)], bound_rows: Vec<BoundRow>, columns: &[Map], offset: &[f64]) -> (Vec<LpRow>, Vec<bool>) {
|
||||
let width = columns.len();
|
||||
let mut transformed = Vec::with_capacity(rows.len() + bound_rows.len());
|
||||
for (coefficients, sense, rhs) in rows {
|
||||
@@ -105,7 +129,8 @@ fn normalize_rows(rows: &[(Vec<f64>, Sense, f64)], bound_rows: Vec<BoundRow>, co
|
||||
transformed.push((row, sense, rhs));
|
||||
}
|
||||
// Rows with a negative right side are negated so every right side is >= 0.
|
||||
transformed
|
||||
let flipped = transformed.iter().map(|row| row.2 < 0.0).collect();
|
||||
let normalized = transformed
|
||||
.into_iter()
|
||||
.map(|(row, sense, rhs)| {
|
||||
if rhs < 0.0 {
|
||||
@@ -119,7 +144,8 @@ fn normalize_rows(rows: &[(Vec<f64>, Sense, f64)], bound_rows: Vec<BoundRow>, co
|
||||
(row, sense, rhs)
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
.collect();
|
||||
(normalized, flipped)
|
||||
}
|
||||
|
||||
fn build_tableau(normalized: &[(Vec<f64>, Sense, f64)], decision_cols: usize, slack_cols: usize, artificial_cols: usize) -> Tableau {
|
||||
@@ -388,20 +414,33 @@ fn simplex(t: &mut [Vec<f64>], basis: &mut [usize], obj: &mut [f64], allowed: &[
|
||||
|
||||
/// Minimises `c·x` subject to `rows` (coefficients, sense, rhs) and `lo ≤ x ≤ hi`.
|
||||
pub(super) fn solve_lp(c: &[f64], rows: &[(Vec<f64>, Sense, f64)], lo: &[f64], hi: &[f64], budget: &mut u64) -> Outcome {
|
||||
solve_lp_tableau(c, rows, lo, hi, budget).0
|
||||
}
|
||||
|
||||
/// [`solve_lp`], with the final tableau when it is optimal.
|
||||
pub(super) fn solve_lp_tableau(
|
||||
c: &[f64],
|
||||
rows: &[(Vec<f64>, Sense, f64)],
|
||||
lo: &[f64],
|
||||
hi: &[f64],
|
||||
budget: &mut u64,
|
||||
) -> (Outcome, Option<FinalTableau>) {
|
||||
if lo.iter().zip(hi).any(|(l, h)| l > h) {
|
||||
return Outcome::Infeasible;
|
||||
return (Outcome::Infeasible, None);
|
||||
}
|
||||
// x_j = offset + Σ sign·y_k with y ≥ 0.
|
||||
let mapping = map_variables(c.len(), lo, hi);
|
||||
let VariableMap { columns: cols, offset, bound_rows: extra } = mapping;
|
||||
let ny = cols.len();
|
||||
// Rows over y: (coefficients, sense, rhs).
|
||||
let normalized = normalize_rows(rows, extra, &cols, &offset);
|
||||
let (normalized, flipped) = normalize_rows(rows, extra, &cols, &offset);
|
||||
let n_slack = normalized.iter().filter(|row| row.1 != Sense::Eq).count();
|
||||
let n_art = normalized.iter().filter(|row| row.1 != Sense::Le).count();
|
||||
let mut tableau = build_tableau(&normalized, ny, n_slack, n_art);
|
||||
// Each row's starting basic column is its unit column.
|
||||
let ident = tableau.basis.clone();
|
||||
if let Err(outcome) = phase_one(&mut tableau, &normalized, budget) {
|
||||
return outcome;
|
||||
return (outcome, None);
|
||||
}
|
||||
let Tableau { rows: mut t, mut basis, width, decision_cols: ny, slack_cols: n_slack, .. } = tableau;
|
||||
let is_art = |column: usize| column >= ny + n_slack && column < width;
|
||||
@@ -423,8 +462,8 @@ pub(super) fn solve_lp(c: &[f64], rows: &[(Vec<f64>, Sense, f64)], lo: &[f64], h
|
||||
}
|
||||
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::Limit => return (Outcome::Limit, None),
|
||||
Step::Unbounded => return (Outcome::Unbounded, None),
|
||||
Step::Optimal => {}
|
||||
}
|
||||
let mut y = vec![0.0; width];
|
||||
@@ -442,7 +481,43 @@ pub(super) fn solve_lp(c: &[f64], rows: &[(Vec<f64>, Sense, f64)], lo: &[f64], h
|
||||
}
|
||||
}
|
||||
let val: f64 = c.iter().zip(&x).map(|(a, b)| a * b).sum();
|
||||
Outcome::Solved(x, val, true)
|
||||
let tableau = FinalTableau { t, basis, obj, width, first_art: ny + n_slack, cols, ident, flipped };
|
||||
(Outcome::Solved(x, val, true), Some(tableau))
|
||||
}
|
||||
|
||||
/// The LP to minimise for `goal`, from the linearised model: the costs, and the rows
|
||||
/// (coefficients, sense, right side) in the order of `model_rows`.
|
||||
pub(super) fn linear_program(goal: Goal, cf: &[f64], f_const: f64, lin: &[(Vec<f64>, f64)], model_rows: &[Row]) -> (Vec<f64>, Vec<LpRow>) {
|
||||
let mut rows: Vec<LpRow> = lin.iter().zip(model_rows).map(|((a, k), r)| (a.clone(), r.sense, -k)).collect();
|
||||
let c = match goal {
|
||||
Goal::Max => cf.iter().map(|v| -v).collect(),
|
||||
Goal::Min => cf.to_vec(),
|
||||
Goal::Value(v) => {
|
||||
rows.push((cf.to_vec(), Sense::Eq, v - f_const));
|
||||
vec![0.0; cf.len()]
|
||||
}
|
||||
};
|
||||
(c, rows)
|
||||
}
|
||||
|
||||
/// The linear program Simplex LP solves: the model linearised at `x0`, its objective
|
||||
/// coefficients, the costs and rows for the goal, and the work limit for its size.
|
||||
pub(super) struct SimplexProgram {
|
||||
pub(super) cf: Vec<f64>,
|
||||
pub(super) c: Vec<f64>,
|
||||
pub(super) rows: Vec<LpRow>,
|
||||
pub(super) work: u64,
|
||||
}
|
||||
|
||||
/// Builds the [`SimplexProgram`] at `x0`. The solve and the Sensitivity report both build it
|
||||
/// here from the start values, so the report's solve repeats the solve's pivots and returns the
|
||||
/// same vertex: linearising at another point changes the coefficients in their last bits, and
|
||||
/// that can break a tie between pivots the other way.
|
||||
pub(super) fn simplex_program(ev: &mut Evaluator, goal: Goal, x0: &[f64]) -> std::result::Result<SimplexProgram, &'static str> {
|
||||
let ((cf, f_const), lin) = linearize(ev, x0)?;
|
||||
let (c, rows) = linear_program(goal, &cf, f_const, &lin, &ev.rows);
|
||||
let size = (rows.len() as u64 + 2 * x0.len() as u64 + 1) * (4 * x0.len() as u64 + 2 * rows.len() as u64 + 1);
|
||||
Ok(SimplexProgram { cf, c, rows, work: LP_WORK.max(size) })
|
||||
}
|
||||
|
||||
/// Linearises the model at `x0` (n + 1 evaluations) and checks it at one more point. Returns the
|
||||
|
||||
@@ -306,6 +306,32 @@ pub(super) fn parse_model(s: &Session, p: &Json) -> Result<Model> {
|
||||
Ok(Model { sheet, objective, goal, vars, vars_text, cons, non_negative, method, evo: evo_options(p), max_evals, gap })
|
||||
}
|
||||
|
||||
/// Whether a `data.solver` result status is a solution the Sensitivity report can describe: an
|
||||
/// optimum (Simplex LP) or a converged local optimum (GRG Nonlinear). A run that stopped at its
|
||||
/// evaluation limit, or an Evolutionary run, has none.
|
||||
pub fn sensitivity_status(status: &str) -> bool {
|
||||
matches!(status, "optimal" | "converged")
|
||||
}
|
||||
|
||||
/// Whether `data.solver` params (such as the `model` of its result) describe a model that can
|
||||
/// have a Sensitivity report: Simplex LP or GRG Nonlinear, a Max or Min objective and no
|
||||
/// integer or binary constraints.
|
||||
pub fn sensitivity_applies(p: &Json) -> bool {
|
||||
let method = str_param(p, "method").unwrap_or("grg").to_ascii_lowercase().replace([' ', '_', '-'], "");
|
||||
let goal = str_param(p, "goal").or_else(|| str_param(p, "to")).unwrap_or("max").to_ascii_lowercase();
|
||||
let integer = p.get("constraints").and_then(Json::as_array).is_some_and(|list| {
|
||||
list.iter().any(|c| {
|
||||
let rel = match c {
|
||||
Json::String(t) => split_constraint(t).map(|(_, rel, _)| rel),
|
||||
Json::Object(_) => str_param(c, "op").or_else(|| str_param(c, "relation")).and_then(Rel::parse),
|
||||
_ => None,
|
||||
};
|
||||
matches!(rel, Some(Rel::Int | Rel::Bin))
|
||||
})
|
||||
});
|
||||
!integer && !matches!(method.as_str(), "evolutionary" | "evo" | "genetic") && !matches!(goal.as_str(), "value" | "valueof" | "value of")
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------- saved model
|
||||
|
||||
/// `"B2:B4,D1"` as absolute references, each qualified by `sheet_name` when given.
|
||||
|
||||
@@ -1,9 +1,10 @@
|
||||
use super::lp::Outcome;
|
||||
use super::model::{Constraint, Goal, Method, Model, Rel, Rhs, parse_model, saved_model};
|
||||
use super::problem::{Evaluator, Problem, as_num, build_problem};
|
||||
use super::sensitivity::{sensitivity_refusal, sensitivity_report};
|
||||
use super::{
|
||||
Cell, CellRef, Ctx, DocState, EngineError, Json, PRECISION, REPORT_CMD, Result, Session, Sheet, Trial, Value, bad, edit, json, optimise,
|
||||
str_param,
|
||||
Cell, CellRef, Ctx, DocState, EngineError, Json, PRECISION, REPORT_CMD, Result, Session, Sheet, Trial, Value, bad, bool_param, edit, json,
|
||||
optimise, str_param,
|
||||
};
|
||||
|
||||
// ---------------------------------------------------------------- reports
|
||||
@@ -11,6 +12,7 @@ use super::{
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum ReportKind {
|
||||
Answer,
|
||||
Sensitivity,
|
||||
Limits,
|
||||
}
|
||||
|
||||
@@ -18,6 +20,7 @@ impl ReportKind {
|
||||
fn title(self) -> &'static str {
|
||||
match self {
|
||||
ReportKind::Answer => "Answer",
|
||||
ReportKind::Sensitivity => "Sensitivity",
|
||||
ReportKind::Limits => "Limits",
|
||||
}
|
||||
}
|
||||
@@ -26,13 +29,14 @@ impl ReportKind {
|
||||
/// The `reports` list (`None` when absent).
|
||||
pub(super) fn report_kinds(cmd: &'static str, p: &Json) -> Result<Option<Vec<ReportKind>>> {
|
||||
let Some(list) = p.get("reports") else { return Ok(None) };
|
||||
let list = list.as_array().ok_or_else(|| bad(cmd, "`reports` must be an array of answer and limits"))?;
|
||||
let list = list.as_array().ok_or_else(|| bad(cmd, "`reports` must be an array of answer, sensitivity and limits"))?;
|
||||
let mut kinds = Vec::new();
|
||||
for k in list {
|
||||
let kind = match k.as_str().map(|t| t.trim().to_ascii_lowercase()).as_deref() {
|
||||
Some("answer") => ReportKind::Answer,
|
||||
Some("sensitivity") => ReportKind::Sensitivity,
|
||||
Some("limits") => ReportKind::Limits,
|
||||
_ => return Err(bad(cmd, format!("unknown report {k} (answer or limits)"))),
|
||||
_ => return Err(bad(cmd, format!("unknown report {k} (answer, sensitivity or limits)"))),
|
||||
};
|
||||
if !kinds.contains(&kind) {
|
||||
kinds.push(kind);
|
||||
@@ -42,14 +46,14 @@ pub(super) fn report_kinds(cmd: &'static str, p: &Json) -> Result<Option<Vec<Rep
|
||||
}
|
||||
|
||||
/// A report's cells: (row, column, value, bold).
|
||||
type Lines = Vec<(u32, u32, Value, bool)>;
|
||||
pub(super) type Lines = Vec<(u32, u32, Value, bool)>;
|
||||
|
||||
fn abs_ref(c: CellRef) -> String {
|
||||
pub(super) fn abs_ref(c: CellRef) -> String {
|
||||
format!("${}${}", gridcraft_core::col_to_letters(c.col), u64::from(c.row) + 1)
|
||||
}
|
||||
|
||||
/// A cell's name in the reports: the nearest text to its left and the nearest text above it.
|
||||
fn cell_name(sh: &Sheet, c: CellRef) -> String {
|
||||
pub(super) fn cell_name(sh: &Sheet, c: CellRef) -> String {
|
||||
let text = |c: CellRef| match sh.value(c) {
|
||||
Value::Text(t) if !t.trim().is_empty() => Some(t.to_string()),
|
||||
_ => None,
|
||||
@@ -68,7 +72,7 @@ fn method_name(m: Method) -> &'static str {
|
||||
}
|
||||
|
||||
/// Header labels in bold from column 1, then the rows below them; returns the next free row.
|
||||
fn table(out: &mut Lines, row: u32, headers: &[&str], rows: Vec<Vec<Value>>) -> u32 {
|
||||
pub(super) fn table(out: &mut Lines, row: u32, headers: &[&str], rows: Vec<Vec<Value>>) -> u32 {
|
||||
for (k, h) in headers.iter().enumerate() {
|
||||
out.push((row, 1 + k as u32, Value::text(*h), true));
|
||||
}
|
||||
@@ -92,7 +96,7 @@ pub(super) struct ReportInput<'a> {
|
||||
pub(super) budget: u64,
|
||||
}
|
||||
|
||||
fn title(out: &mut Lines, kind: ReportKind, sheet: &str) -> u32 {
|
||||
pub(super) fn title(out: &mut Lines, kind: ReportKind, sheet: &str) -> u32 {
|
||||
out.push((0, 0, Value::text(format!("GridCraft {} Report", kind.title())), true));
|
||||
out.push((1, 0, Value::text(format!("Worksheet: {sheet}")), false));
|
||||
3
|
||||
@@ -255,18 +259,67 @@ fn limits_report(d: &DocState, sh: &Sheet, inp: &ReportInput) -> Result<Lines> {
|
||||
}
|
||||
|
||||
/// Builds every report before the solution or saved model can change.
|
||||
pub(super) fn prepare_reports(d: &DocState, kinds: &[ReportKind], inp: &ReportInput) -> Result<Vec<(String, Lines)>> {
|
||||
/// What happens to a Sensitivity report the model allows but the solution doesn't: `Refuse`
|
||||
/// fails the call before anything changes; `Skip` leaves the report out, with its reason, and
|
||||
/// writes the others. `optimum` is false when the solve stopped before an optimum, so the
|
||||
/// report isn't attempted.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum NotOptimal {
|
||||
Refuse,
|
||||
Skip { optimum: bool },
|
||||
}
|
||||
|
||||
/// The reports left out under [`NotOptimal::Skip`], as `{report, reason}`.
|
||||
pub(super) type Skipped = Vec<Json>;
|
||||
|
||||
pub(super) fn prepare_reports(
|
||||
d: &DocState,
|
||||
kinds: &[ReportKind],
|
||||
inp: &ReportInput,
|
||||
not_optimal: NotOptimal,
|
||||
) -> Result<(Vec<(String, Lines)>, Skipped)> {
|
||||
if kinds.is_empty() {
|
||||
return Ok(Vec::new());
|
||||
return Ok((Vec::new(), Vec::new()));
|
||||
}
|
||||
if d.wb.protected_structure {
|
||||
return Err(EngineError::Other("The workbook is protected and cannot be changed.".into()));
|
||||
// A request that must write a sheet is refused before any report is built (a Limits report
|
||||
// runs two optimisations per variable); only a Sensitivity report that may be skipped waits
|
||||
// for the check after the loop.
|
||||
let must_write = kinds.iter().any(|k| *k != ReportKind::Sensitivity || not_optimal == NotOptimal::Refuse);
|
||||
if must_write && d.wb.protected_structure {
|
||||
return Err(protected());
|
||||
}
|
||||
// A model that has no Sensitivity report (integer, Evolutionary, Value Of) is refused in
|
||||
// both modes before any report is built: a Limits report listed first would otherwise run
|
||||
// its two optimisations per variable for an error. Only the solution's own refusal (not
|
||||
// optimal) waits for the loop, where `Skip` may leave the report out.
|
||||
if kinds.contains(&ReportKind::Sensitivity)
|
||||
&& let Some(why) = sensitivity_refusal(inp.model, inp.prob)
|
||||
{
|
||||
return Err(EngineError::Other(why.into()));
|
||||
}
|
||||
let sh = d.wb.sheet(inp.model.sheet).ok_or(EngineError::NoDocument)?;
|
||||
let mut built = Vec::new();
|
||||
let mut skipped = Vec::new();
|
||||
for k in kinds {
|
||||
let lines = match k {
|
||||
ReportKind::Answer => answer_report(d, sh, inp),
|
||||
ReportKind::Sensitivity => match not_optimal {
|
||||
NotOptimal::Skip { optimum } => {
|
||||
let got = if optimum {
|
||||
sensitivity_report(d, sh, inp).map_err(|e| e.to_string())
|
||||
} else {
|
||||
Err("Solver did not reach an optimal solution, so the Sensitivity report was not written.".to_string())
|
||||
};
|
||||
match got {
|
||||
Ok(lines) => lines,
|
||||
Err(reason) => {
|
||||
skipped.push(json!({"report": "sensitivity", "reason": reason}));
|
||||
continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
NotOptimal::Refuse => sensitivity_report(d, sh, inp)?,
|
||||
},
|
||||
ReportKind::Limits => limits_report(d, sh, inp)?,
|
||||
};
|
||||
let name = (1..=10_000u32)
|
||||
@@ -275,7 +328,16 @@ pub(super) fn prepare_reports(d: &DocState, kinds: &[ReportKind], inp: &ReportIn
|
||||
.ok_or_else(|| EngineError::Other("Too many report sheets.".into()))?;
|
||||
built.push((name, lines));
|
||||
}
|
||||
Ok(built)
|
||||
// Only a report sheet changes the workbook's structure, so a request whose reports were all
|
||||
// skipped still succeeds in a protected workbook.
|
||||
if !built.is_empty() && d.wb.protected_structure {
|
||||
return Err(protected());
|
||||
}
|
||||
Ok((built, skipped))
|
||||
}
|
||||
|
||||
fn protected() -> EngineError {
|
||||
EngineError::Other("The workbook is protected and cannot be changed.".into())
|
||||
}
|
||||
|
||||
/// Inserts prepared reports before the model's sheet, which stays active.
|
||||
@@ -292,7 +354,7 @@ pub(super) fn write_reports(cx: &mut Ctx, at: usize, built: Vec<(String, Lines)>
|
||||
let style = if b { bold } else { gridcraft_model::StyleId::DEFAULT };
|
||||
sheet.set_cell(CellRef::new(r, c), Cell { value: v, formula: None, style, py: None });
|
||||
}
|
||||
for (c, w) in [(0u32, 14.0f32), (1, 110.0), (2, 140.0), (3, 120.0), (4, 120.0), (5, 120.0), (6, 120.0), (7, 120.0)] {
|
||||
for (c, w) in [(0u32, 14.0f32), (1, 110.0), (2, 140.0), (3, 120.0), (4, 120.0), (5, 155.0), (6, 140.0), (7, 140.0)] {
|
||||
sheet.cols.insert(c, gridcraft_model::LineInfo { size: Some(w), ..Default::default() });
|
||||
}
|
||||
cx.wb.sheets.insert(at + i, std::sync::Arc::new(sheet));
|
||||
@@ -321,9 +383,9 @@ pub(super) fn insert_report_selections(s: &mut Session, at: usize, count: usize)
|
||||
|
||||
/// `data.solverReport`: reports for a solution of the sheet's saved model.
|
||||
pub(super) fn report(s: &mut Session, p: &Json) -> Result<Json> {
|
||||
let kinds = report_kinds(REPORT_CMD, p)?.ok_or_else(|| bad(REPORT_CMD, "missing `reports` (answer, limits)"))?;
|
||||
let kinds = report_kinds(REPORT_CMD, p)?.ok_or_else(|| bad(REPORT_CMD, "missing `reports` (answer, sensitivity, limits)"))?;
|
||||
if kinds.is_empty() {
|
||||
return Err(bad(REPORT_CMD, "`reports` is empty (answer, limits)"));
|
||||
return Err(bad(REPORT_CMD, "`reports` is empty (answer, sensitivity, limits)"));
|
||||
}
|
||||
let d = s.doc()?;
|
||||
// The given model as is, else the saved model with any params given here over it.
|
||||
@@ -365,8 +427,14 @@ pub(super) fn report(s: &mut Session, p: &Json) -> Result<Json> {
|
||||
let budget = model.max_evals;
|
||||
let message = str_param(p, "message").unwrap_or("").chars().take(500).collect();
|
||||
let inp = ReportInput { model: &model, prob: &prob, original, solution, message, budget };
|
||||
let built = prepare_reports(s.doc()?, &kinds, &inp)?;
|
||||
let sheets = edit(s, |cx| Ok(write_reports(cx, model.sheet, built)))?;
|
||||
let not_optimal = if bool_param(p, "skipNotOptimal").unwrap_or(false) { NotOptimal::Skip { optimum: true } } else { NotOptimal::Refuse };
|
||||
let (built, skipped) = prepare_reports(s.doc()?, &kinds, &inp, not_optimal)?;
|
||||
// Nothing left to write (every report skipped): no edit, so no Undo step.
|
||||
let sheets = if built.is_empty() { Vec::new() } else { edit(s, |cx| Ok(write_reports(cx, model.sheet, built)))? };
|
||||
insert_report_selections(s, model.sheet, sheets.len());
|
||||
Ok(json!({"sheets": sheets}))
|
||||
let mut out = json!({"sheets": sheets});
|
||||
if !skipped.is_empty() {
|
||||
out["skippedReports"] = json!(skipped);
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
@@ -0,0 +1,614 @@
|
||||
//! The Sensitivity report. Simplex LP: each variable's reduced cost and the range of its
|
||||
//! objective coefficient over which the basis stays optimal, and each constraint's shadow price
|
||||
//! and the range of its right side over which the shadow price holds, from the final tableau of
|
||||
//! the LP linearised at the solution. GRG Nonlinear: the reduced gradients and the Lagrange
|
||||
//! multipliers of the binding constraints at the solution, from difference gradients.
|
||||
|
||||
use super::lp::{FinalTableau, LpRow, Outcome, Row, Sense, SimplexProgram, simplex_program, solve_lp_tableau};
|
||||
use super::model::{Goal, Method, Model};
|
||||
use super::problem::{Evaluator, Point, Problem, as_num};
|
||||
use super::reports::{Lines, ReportInput, ReportKind, abs_ref, cell_name, table, title};
|
||||
use super::{CellRef, DocState, EPS, EngineError, Result, Sheet, Trial, Value};
|
||||
|
||||
/// Excel prints an unlimited allowable change as 1E+30.
|
||||
const NO_LIMIT: f64 = 1e30;
|
||||
|
||||
/// Most arithmetic the Lagrange multipliers of a nonlinear model may take (about a second): a
|
||||
/// pass over k binding rows and f free variables costs about 3·k²·f (the normal equations and
|
||||
/// each multiplier's noise), and each pass drops at most one row, so a model at the constraint
|
||||
/// limit could otherwise take minutes.
|
||||
#[cfg(not(test))]
|
||||
const MULTIPLIER_WORK: u64 = 1_000_000_000;
|
||||
/// A low limit in tests, so the refusal is reachable with a small model.
|
||||
#[cfg(test)]
|
||||
const MULTIPLIER_WORK: u64 = 20_000;
|
||||
|
||||
/// Why a model has no Sensitivity report (`None` when it has one).
|
||||
pub(super) fn sensitivity_refusal(m: &Model, prob: &Problem) -> Option<&'static str> {
|
||||
if prob.integer.iter().any(|b| *b) {
|
||||
Some("The Sensitivity report doesn't apply to models with integer or binary constraints.")
|
||||
} else if m.method == Method::Evolutionary {
|
||||
Some("The Sensitivity report is available for Simplex LP and GRG Nonlinear.")
|
||||
} else if matches!(m.goal, Goal::Value(_)) {
|
||||
Some("The Sensitivity report needs an objective to maximise or minimise.")
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}
|
||||
|
||||
/// A sensitivity figure: values within `noise` of zero (rounding) print as 0, unlimited ones
|
||||
/// as 1E+30. The noise level is local to the figure: it comes from the terms the figure is a
|
||||
/// sum of (the rounding error of a sum scales with the sum of its terms' sizes), never from
|
||||
/// the model's largest coefficient or gradient component, so a small price or gradient keeps
|
||||
/// its value beside a large unrelated one, and never from the objective's value, so adding a
|
||||
/// constant to the objective changes nothing the report shows.
|
||||
fn sens_num(v: f64, noise: f64) -> Value {
|
||||
if !v.is_finite() || v >= NO_LIMIT {
|
||||
Value::number(NO_LIMIT)
|
||||
} else if v.abs() <= noise {
|
||||
Value::number(0.0)
|
||||
} else {
|
||||
Value::number(v)
|
||||
}
|
||||
}
|
||||
|
||||
fn num(v: f64) -> Value {
|
||||
if v.is_finite() { Value::number(v) } else { Value::Error(gridcraft_core::CellError::Value) }
|
||||
}
|
||||
|
||||
/// The smallest ratio `num / den` over the `den` above zero (the step to the first bound).
|
||||
fn min_ratio(pairs: impl Iterator<Item = (f64, f64)>) -> f64 {
|
||||
pairs.filter(|(_, den)| *den > EPS).map(|(num, den)| num.max(0.0) / den).fold(f64::INFINITY, f64::min)
|
||||
}
|
||||
|
||||
/// Solves the square system `a·x = b` by Gaussian elimination with partial pivoting; unknowns
|
||||
/// of a singular system that have no pivot are 0.
|
||||
fn solve_dense(mut a: Vec<Vec<f64>>, mut b: Vec<f64>) -> Vec<f64> {
|
||||
let n = b.len();
|
||||
let scale = a.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())).max(f64::MIN_POSITIVE);
|
||||
let mut pivots: Vec<Option<usize>> = vec![None; n];
|
||||
let mut row = 0;
|
||||
for col in 0..n {
|
||||
let Some(p) = (row..n).max_by(|i, k| {
|
||||
let v = |r: &usize| a.get(*r).and_then(|x| x.get(col)).map_or(0.0, |v| v.abs());
|
||||
v(i).total_cmp(&v(k))
|
||||
}) else {
|
||||
break;
|
||||
};
|
||||
if a.get(p).and_then(|x| x.get(col)).is_none_or(|v| v.abs() <= 1e-12 * scale) {
|
||||
continue;
|
||||
}
|
||||
a.swap(row, p);
|
||||
b.swap(row, p);
|
||||
let (Some(prow), Some(pb)) = (a.get(row).cloned(), b.get(row).copied()) else { break };
|
||||
let pv = prow.get(col).copied().unwrap_or(1.0);
|
||||
for (i, (r, bi)) in a.iter_mut().zip(b.iter_mut()).enumerate() {
|
||||
if i == row {
|
||||
continue;
|
||||
}
|
||||
let f = r.get(col).copied().unwrap_or(0.0) / pv;
|
||||
if f != 0.0 {
|
||||
for (v, p) in r.iter_mut().zip(&prow) {
|
||||
*v -= f * p;
|
||||
}
|
||||
*bi -= f * pb;
|
||||
}
|
||||
}
|
||||
if let Some(e) = pivots.get_mut(col) {
|
||||
*e = Some(row);
|
||||
}
|
||||
row += 1;
|
||||
}
|
||||
pivots
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(col, p)| match p {
|
||||
Some(r) => b.get(*r).copied().unwrap_or(0.0) / a.get(*r).and_then(|x| x.get(col)).copied().unwrap_or(1.0),
|
||||
None => 0.0,
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn out_of_evaluations() -> EngineError {
|
||||
EngineError::Other("Solver ran out of evaluations.".into())
|
||||
}
|
||||
|
||||
fn not_optimal() -> EngineError {
|
||||
EngineError::Other("The Sensitivity report needs an optimal solution: these values are not optimal for the model. Solve it first.".into())
|
||||
}
|
||||
|
||||
/// What both report bodies share: the model, the solution, its objective value and the rows'
|
||||
/// final values and right sides.
|
||||
struct At<'a> {
|
||||
inp: &'a ReportInput<'a>,
|
||||
sh: &'a Sheet,
|
||||
point: Point,
|
||||
finals: Vec<(f64, f64)>,
|
||||
}
|
||||
|
||||
impl At<'_> {
|
||||
fn var_head(&self, j: usize, c: CellRef) -> Vec<Value> {
|
||||
vec![Value::text(abs_ref(c)), Value::text(cell_name(self.sh, c)), num(self.inp.solution.get(j).copied().unwrap_or(f64::NAN))]
|
||||
}
|
||||
fn row_head(&self, i: usize, r: &Row) -> Vec<Value> {
|
||||
let fin = self.finals.get(i).map_or(f64::NAN, |f| f.0);
|
||||
vec![Value::text(abs_ref(r.lhs)), Value::text(cell_name(self.sh, r.lhs)), num(fin)]
|
||||
}
|
||||
}
|
||||
|
||||
/// The variable rows, the constraint rows and their headers.
|
||||
type Body = (Vec<Vec<Value>>, Vec<Vec<Value>>, Vec<&'static str>, Vec<&'static str>);
|
||||
|
||||
pub(super) fn sensitivity_report(d: &DocState, sh: &Sheet, inp: &ReportInput) -> Result<Lines> {
|
||||
let m = inp.model;
|
||||
if let Some(why) = sensitivity_refusal(m, inp.prob) {
|
||||
return Err(EngineError::Other(why.into()));
|
||||
}
|
||||
let n = m.vars.len();
|
||||
let x = &inp.solution;
|
||||
let mut ev = Evaluator {
|
||||
trial: Trial::new(d),
|
||||
sheet: m.sheet,
|
||||
vars: m.vars.clone(),
|
||||
objective: m.objective,
|
||||
rows: inp.prob.rows.clone(),
|
||||
evals: 0,
|
||||
// Room for the gradients (GRG) or two linearisations (Simplex LP) at the least.
|
||||
budget: inp.budget.max(3 * n as u64 + 4),
|
||||
target: None,
|
||||
};
|
||||
let point = ev.eval(x).ok_or_else(out_of_evaluations)?;
|
||||
// The rows' final values and right sides at the solution.
|
||||
let mut cells = Vec::new();
|
||||
for r in &inp.prob.rows {
|
||||
cells.push(r.lhs);
|
||||
cells.extend(r.rhs);
|
||||
}
|
||||
let mut trial = Trial::new(d);
|
||||
let sets: Vec<(usize, CellRef, Value)> = m.vars.iter().zip(x).map(|(c, v)| (m.sheet, *c, Value::number(*v))).collect();
|
||||
let reads: Vec<(usize, CellRef)> = cells.iter().map(|c| (m.sheet, *c)).collect();
|
||||
let mut vals = trial.run(&sets, &reads).into_iter().map(|v| as_num(&v));
|
||||
let finals: Vec<(f64, f64)> = inp
|
||||
.prob
|
||||
.rows
|
||||
.iter()
|
||||
.map(|r| {
|
||||
let l = vals.next().unwrap_or(f64::NAN);
|
||||
(l, if r.rhs.is_some() { vals.next().unwrap_or(f64::NAN) } else { r.rhs_const })
|
||||
})
|
||||
.collect();
|
||||
let at = At { inp, sh, point, finals };
|
||||
let (var_rows, con_rows, var_headers, con_headers) =
|
||||
if m.method == Method::Simplex { simplex_body(&mut ev, &at)? } else { nonlinear_body(&mut ev, &at)? };
|
||||
let mut out = Lines::new();
|
||||
let r = title(&mut out, ReportKind::Sensitivity, &sh.name);
|
||||
out.push((r, 0, Value::text("Variable Cells"), true));
|
||||
let r = table(&mut out, r + 1, &var_headers, var_rows);
|
||||
out.push((r, 0, Value::text("Constraints"), true));
|
||||
table(&mut out, r + 1, &con_headers, con_rows);
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Simplex LP: reduced costs, shadow prices and their ranges from the final tableau.
|
||||
fn simplex_body(ev: &mut Evaluator, at: &At) -> Result<Body> {
|
||||
let (inp, m) = (at.inp, at.inp.model);
|
||||
// The program the solve built, from the same start values (data.solver's x0, which Solver
|
||||
// Results passes back as `original`), so the solve below repeats its pivots. Start values
|
||||
// the model can't be evaluated at (given only to the report) fall back to the solution.
|
||||
let start: Vec<f64> = inp.original.iter().map(|v| if v.is_finite() { *v } else { 0.0 }).collect();
|
||||
let built = match simplex_program(ev, m.goal, &start) {
|
||||
Err("nonlinear" | "error") if start != inp.solution => simplex_program(ev, m.goal, &inp.solution),
|
||||
other => other,
|
||||
};
|
||||
let SimplexProgram { cf, c, rows, mut work } = built.map_err(|e| {
|
||||
if e == "budget" {
|
||||
out_of_evaluations()
|
||||
} else {
|
||||
EngineError::Other("The model isn't linear in the variable cells, so Simplex LP can't report its sensitivity.".into())
|
||||
}
|
||||
})?;
|
||||
let (outcome, tab) = solve_lp_tableau(&c, &rows, &inp.prob.lo, &inp.prob.hi, &mut work);
|
||||
let (Outcome::Solved(vertex, best, _), Some(tab)) = (outcome, tab) else {
|
||||
return Err(EngineError::Other("Solver could not find an optimal basis for the Sensitivity report.".into()));
|
||||
};
|
||||
// The report describes the optimum, so the values it shows must be one.
|
||||
if !optimal_point(&c, &rows, (&inp.prob.lo, &inp.prob.hi), &inp.solution, best) {
|
||||
return Err(not_optimal());
|
||||
}
|
||||
// The prices and ranges come from the tableau of `vertex`, so the values shown must be that
|
||||
// vertex. Another optimum with the same objective (when the optimum is not unique) has other
|
||||
// ranges, so it is refused too, although it is optimal.
|
||||
let near = |a: f64, b: f64| (a - b).abs() <= super::PRECISION * (1.0 + a.abs().max(b.abs()));
|
||||
if vertex.len() != inp.solution.len() || !vertex.iter().zip(&inp.solution).all(|(a, b)| near(*a, *b)) {
|
||||
return Err(EngineError::Other(
|
||||
"The Sensitivity report needs the values at the Simplex LP solution: these values are another optimal solution, whose ranges differ. Solve the model with Simplex LP first."
|
||||
.into(),
|
||||
));
|
||||
}
|
||||
// Minimising `sign · objective`.
|
||||
let sign = if m.goal == Goal::Max { -1.0 } else { 1.0 };
|
||||
// A shadow price is the dual yᵢ = Σᵣ c_B(r)·B⁻¹(r, i): the row's unit column of the final
|
||||
// tableau holds B⁻¹'s column i, and c_B(r) is the cost of the column basic in row r (a
|
||||
// slack or artificial costs nothing). The rounding error of a sum scales with the sum of
|
||||
// its terms' sizes, so a price below 1E-9 of Σᵣ |c_B(r)·B⁻¹(r, i)| is noise. The bound is
|
||||
// the row's own: a small price keeps its value beside a large coefficient that plays no
|
||||
// part in it, and a constant in the objective (no coefficient) changes no bound. Allowable
|
||||
// changes to a coefficient or a right side are ratios the tableau's own pivot tolerance
|
||||
// (`EPS`) already bounds.
|
||||
let cost_size = |b: usize| tab.cols.get(b).map_or(0.0, |mp| cf.get(mp.var).copied().unwrap_or(0.0).abs());
|
||||
let price_terms = |col: usize| (0..tab.t.len()).map(|r| tab.basis.get(r).map_or(0.0, |b| cost_size(*b)) * entry(&tab, r, col).abs()).sum::<f64>();
|
||||
let (shadow, terms): (Vec<f64>, Vec<f64>) = (0..rows.len())
|
||||
.map(|i| {
|
||||
let col = tab.ident.get(i).copied().unwrap_or(usize::MAX);
|
||||
(sign * flip(&tab, i) * -tab.obj.get(col).copied().unwrap_or(0.0), price_terms(col))
|
||||
})
|
||||
.unzip();
|
||||
let noise = |terms: f64| 1e-9 * terms.max(f64::MIN_POSITIVE);
|
||||
let con_rows = ev
|
||||
.rows
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, r)| {
|
||||
// The basic values move along the row's unit column as its right side changes.
|
||||
let col = tab.ident.get(i).copied().unwrap_or(usize::MAX);
|
||||
let dir: Vec<(f64, f64)> =
|
||||
(0..tab.t.len()).filter(|r| real_basic(&tab, *r)).map(|r| (entry(&tab, r, tab.width), flip(&tab, i) * entry(&tab, r, col))).collect();
|
||||
// An artificial still basic after phase one stands for a row the others make
|
||||
// redundant, and it is zero there. When this right side moves it, any change alone
|
||||
// leaves no solution, so both ranges are 0.
|
||||
let pinned = (0..tab.t.len()).any(|r| !real_basic(&tab, r) && entry(&tab, r, col).abs() > EPS);
|
||||
let (increase, decrease) =
|
||||
if pinned { (0.0, 0.0) } else { (min_ratio(dir.iter().map(|(v, d)| (*v, -d))), min_ratio(dir.iter().map(|(v, d)| (*v, *d)))) };
|
||||
let b = at.finals.get(i).map_or(f64::NAN, |f| f.1);
|
||||
let mut row = at.row_head(i, r);
|
||||
let price = sens_num(shadow.get(i).copied().unwrap_or(0.0), noise(terms.get(i).copied().unwrap_or(0.0)));
|
||||
row.extend([price, num(b), sens_num(increase, EPS), sens_num(decrease, EPS)]);
|
||||
row
|
||||
})
|
||||
.collect();
|
||||
let nonbasic: Vec<usize> = (0..tab.first_art).filter(|l| !tab.basis.contains(l)).collect();
|
||||
let var_rows = m
|
||||
.vars
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(j, c)| {
|
||||
let coef = cf.get(j).copied().unwrap_or(0.0);
|
||||
let used: f64 = rows.iter().zip(&shadow).map(|((a, _, _), p)| a.get(j).copied().unwrap_or(0.0) * p).sum();
|
||||
// The reduced cost cf_j − Σᵢ a_ij·yᵢ is a sum of cf_j and the prices' terms times
|
||||
// a_ij: the prices' terms rather than the prices, since a price that is a small
|
||||
// difference of large terms carries their rounding error into every reduced cost
|
||||
// it is used in.
|
||||
let used_terms: f64 = rows.iter().zip(&terms).map(|((a, _, _), t)| a.get(j).copied().unwrap_or(0.0).abs() * t).sum();
|
||||
// A change Δ in the coefficient changes each nonbasic reduced cost by Δ·e.
|
||||
let own = |col: usize| tab.cols.get(col).filter(|mp| mp.var == j).map_or(0.0, |mp| mp.sign);
|
||||
let moves: Vec<(f64, f64)> = nonbasic
|
||||
.iter()
|
||||
.map(|&l| {
|
||||
let through_basis: f64 = (0..tab.t.len()).map(|r| tab.basis.get(r).map_or(0.0, |b| own(*b)) * entry(&tab, r, l)).sum();
|
||||
(tab.obj.get(l).copied().unwrap_or(0.0), sign * (own(l) - through_basis))
|
||||
})
|
||||
.collect();
|
||||
let increase = min_ratio(moves.iter().map(|(d, e)| (*d, -e)));
|
||||
let decrease = min_ratio(moves.iter().map(|(d, e)| (*d, *e)));
|
||||
let mut row = at.var_head(j, *c);
|
||||
row.extend([sens_num(coef - used, noise(coef.abs() + used_terms)), num(coef), sens_num(increase, EPS), sens_num(decrease, EPS)]);
|
||||
row
|
||||
})
|
||||
.collect();
|
||||
Ok((
|
||||
var_rows,
|
||||
con_rows,
|
||||
vec!["Cell", "Name", "Final Value", "Reduced Cost", "Objective Coefficient", "Allowable Increase", "Allowable Decrease"],
|
||||
vec!["Cell", "Name", "Final Value", "Shadow Price", "Constraint R.H. Side", "Allowable Increase", "Allowable Decrease"],
|
||||
))
|
||||
}
|
||||
|
||||
/// Whether `x` satisfies the LP (rows and bounds) and reaches its optimum `best` (minimising
|
||||
/// `c·x`), within the constraint precision.
|
||||
fn optimal_point(c: &[f64], rows: &[LpRow], (lo, hi): (&[f64], &[f64]), x: &[f64], best: f64) -> bool {
|
||||
let near = |v: f64, size: f64| v <= super::PRECISION * (1.0 + size);
|
||||
let dot = |a: &[f64]| a.iter().zip(x).fold((0.0, 0.0), |(s, m), (p, q)| (s + p * q, m + (p * q).abs()));
|
||||
let rows_hold = rows.iter().all(|(a, sense, b)| {
|
||||
let (lhs, size) = dot(a);
|
||||
let size = size + b.abs();
|
||||
match sense {
|
||||
Sense::Le => near(lhs - b, size),
|
||||
Sense::Ge => near(b - lhs, size),
|
||||
Sense::Eq => near((lhs - b).abs(), size),
|
||||
}
|
||||
});
|
||||
let bounds_hold = x.iter().enumerate().all(|(j, v)| {
|
||||
let (l, h) = (lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), hi.get(j).copied().unwrap_or(f64::INFINITY));
|
||||
v.is_finite() && near(l - v, v.abs()) && near(v - h, v.abs())
|
||||
});
|
||||
let (value, size) = dot(c);
|
||||
// Two-sided: a point just outside a bound can beat the optimum by far more than the
|
||||
// precision when its costs are large.
|
||||
rows_hold && bounds_hold && near((value - best).abs(), size + best.abs())
|
||||
}
|
||||
|
||||
fn entry(tab: &FinalTableau, r: usize, col: usize) -> f64 {
|
||||
tab.t.get(r).and_then(|row| row.get(col)).copied().unwrap_or(0.0)
|
||||
}
|
||||
|
||||
fn real_basic(tab: &FinalTableau, r: usize) -> bool {
|
||||
tab.basis.get(r).is_some_and(|b| *b < tab.first_art)
|
||||
}
|
||||
|
||||
/// -1 for a row the tableau negated, else 1.
|
||||
fn flip(tab: &FinalTableau, i: usize) -> f64 {
|
||||
if tab.flipped.get(i).copied().unwrap_or(false) { -1.0 } else { 1.0 }
|
||||
}
|
||||
|
||||
/// Central-difference gradients, inside the bounds, of the objective and of each row's
|
||||
/// `lhs − rhs` at the solution.
|
||||
fn gradients(ev: &mut Evaluator, inp: &ReportInput) -> Result<(Vec<f64>, Vec<Vec<f64>>)> {
|
||||
let x = &inp.solution;
|
||||
let n = x.len();
|
||||
let mut df = vec![0.0; n];
|
||||
let mut dg = vec![vec![0.0; n]; ev.rows.len()];
|
||||
for j in 0..n {
|
||||
let xj = x.get(j).copied().unwrap_or(0.0);
|
||||
let (l, u) = (inp.prob.lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), inp.prob.hi.get(j).copied().unwrap_or(f64::INFINITY));
|
||||
let h = 1e-6 * xj.abs().max(1.0);
|
||||
let (a, b) = ((xj - h).max(l), (xj + h).min(u));
|
||||
if b <= a {
|
||||
continue;
|
||||
}
|
||||
let mut point = |v: f64| {
|
||||
let mut p = x.clone();
|
||||
if let Some(e) = p.get_mut(j) {
|
||||
*e = v;
|
||||
}
|
||||
ev.eval(&p).ok_or_else(out_of_evaluations)
|
||||
};
|
||||
let (pa, pb) = (point(a)?, point(b)?);
|
||||
if let Some(e) = df.get_mut(j) {
|
||||
*e = (pb.f - pa.f) / (b - a);
|
||||
}
|
||||
for (row, (ga, gb)) in dg.iter_mut().zip(pa.g.iter().zip(&pb.g)) {
|
||||
if let Some(e) = row.get_mut(j) {
|
||||
*e = (gb - ga) / (b - a);
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok((df, dg))
|
||||
}
|
||||
|
||||
/// The rounding noise in the Lagrange multiplier of row `i`, one of the `active` rows whose
|
||||
/// multipliers `lambda` solve the normal equations of ∇f = Σ λₖ ∇gₖ over the `free` variables.
|
||||
/// Row i's equation gives λᵢ = (Σⱼ gᵢⱼ·dfⱼ − Σ_{k≠i} λₖ·Σⱼ gᵢⱼ·gₖⱼ) / Σⱼ gᵢⱼ²; the difference
|
||||
/// gradients are good to about 1E-6 of each component, and the rounding error of a sum scales
|
||||
/// with the sum of its terms' sizes, so the noise is 1E-6 of the terms' sizes over the row's
|
||||
/// own squared norm. The bound is local to the row: another row enters only through its
|
||||
/// coupling term with this one, so a row whose gradient is orthogonal to the others has a bound
|
||||
/// from its own gradient alone, and a large component of ∇f that the row's gradient doesn't
|
||||
/// touch hides nothing. Nearly parallel active gradients make the normal equations
|
||||
/// ill-conditioned and amplify every multiplier's error beyond any bound local to a row; such
|
||||
/// multipliers are ill-determined whatever the bound, and they print as computed.
|
||||
fn multiplier_noise(i: usize, active: &[usize], lambda: &[f64], free: &[usize], (df, dg): (&[f64], &[Vec<f64>])) -> f64 {
|
||||
let grad = |i: usize, j: usize| dg.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
let own: f64 = free.iter().map(|&j| grad(i, j) * grad(i, j)).sum();
|
||||
let objective: f64 = free.iter().map(|&j| (grad(i, j) * df.get(j).copied().unwrap_or(0.0)).abs()).sum();
|
||||
let coupling: f64 = active
|
||||
.iter()
|
||||
.zip(lambda)
|
||||
.filter(|(k, _)| **k != i)
|
||||
.map(|(&k, l)| free.iter().map(|&j| (grad(i, j) * grad(k, j) * l).abs()).sum::<f64>())
|
||||
.sum();
|
||||
1e-6 * (objective + coupling) / own.max(f64::MIN_POSITIVE)
|
||||
}
|
||||
|
||||
/// Least squares on the free variables: ∇f = Σ λᵢ ∇gᵢ over the `active` rows. An inequality's
|
||||
/// multiplier has the sign that tightening it costs (λ ≥ 0 for `lhs ≤ rhs` when maximising); a
|
||||
/// row whose multiplier has the other sign beyond its own noise ([`multiplier_noise`], the
|
||||
/// criterion the report rounds it with) isn't binding, so it leaves and the rest are solved
|
||||
/// again; the objective's value plays no part. Returns each row's multiplier and its noise (0
|
||||
/// and 0 for a row that doesn't bind). The work is bounded by `budget` (see
|
||||
/// [`MULTIPLIER_WORK`]); past it the report is refused.
|
||||
fn multipliers(
|
||||
rows: &[Row],
|
||||
goal: Goal,
|
||||
mut active: Vec<usize>,
|
||||
free: &[usize],
|
||||
(df, dg): (&[f64], &[Vec<f64>]),
|
||||
budget: u64,
|
||||
) -> Result<(Vec<f64>, Vec<f64>)> {
|
||||
let grad = |i: usize, j: usize| dg.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
let goal_sign = if goal == Goal::Max { 1.0 } else { -1.0 };
|
||||
let mut lambda = vec![0.0; rows.len()];
|
||||
let mut noise = vec![0.0; rows.len()];
|
||||
let mut work = 0u64;
|
||||
while !active.is_empty() {
|
||||
let (k, f) = (active.len() as u64, (free.len() as u64).max(1));
|
||||
work = work.saturating_add(3u64.saturating_mul(k).saturating_mul(k).saturating_mul(f));
|
||||
if work > budget {
|
||||
return Err(EngineError::Other(
|
||||
"The Sensitivity report of this nonlinear model has too many binding constraints to work out in reasonable time.".into(),
|
||||
));
|
||||
}
|
||||
let normal: Vec<Vec<f64>> =
|
||||
active.iter().map(|&i| active.iter().map(|&k| free.iter().map(|&j| grad(i, j) * grad(k, j)).sum()).collect()).collect();
|
||||
let rhs: Vec<f64> = active.iter().map(|&i| free.iter().map(|&j| grad(i, j) * df.get(j).copied().unwrap_or(0.0)).sum()).collect();
|
||||
let solved = solve_dense(normal, rhs);
|
||||
let bounds: Vec<f64> = active.iter().map(|&i| multiplier_noise(i, &active, &solved, free, (df, dg))).collect();
|
||||
let wrong = active
|
||||
.iter()
|
||||
.zip(&solved)
|
||||
.zip(&bounds)
|
||||
.enumerate()
|
||||
.filter_map(|(k, ((&i, l), bound))| {
|
||||
let want = match rows.get(i).map_or(Sense::Eq, |r| r.sense) {
|
||||
Sense::Le => goal_sign,
|
||||
Sense::Ge => -goal_sign,
|
||||
Sense::Eq => return None,
|
||||
};
|
||||
(want * l < -bound).then_some((k, want * l))
|
||||
})
|
||||
.min_by(|a, b| a.1.total_cmp(&b.1));
|
||||
match wrong {
|
||||
Some((k, _)) => {
|
||||
active.remove(k);
|
||||
}
|
||||
None => {
|
||||
for ((i, l), bound) in active.iter().zip(solved).zip(bounds) {
|
||||
if let (Some(e), Some(n)) = (lambda.get_mut(*i), noise.get_mut(*i)) {
|
||||
*e = l;
|
||||
*n = bound;
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok((lambda, noise))
|
||||
}
|
||||
|
||||
/// Whether the solution is a Karush-Kuhn-Tucker point of the model, from the gradients the
|
||||
/// report computed: every row and bound holds, and each reduced gradient (the objective's
|
||||
/// gradient less what the binding rows' multipliers explain) is near zero relative to `terms`,
|
||||
/// the sizes of the terms it is made of, or, for a variable at a bound, points only out of that
|
||||
/// bound in the objective's direction. A binding row whose multiplier had the wrong sign was
|
||||
/// dropped by [`multipliers`], so the part of the gradient it would have explained remains in
|
||||
/// the reduced gradients and the point is refused here. Every tolerance is local to the row or
|
||||
/// variable it tests: a large right side or gradient component elsewhere in the model excuses
|
||||
/// nothing here.
|
||||
fn kkt_point(at: &At, rows: &[Row], terms: &[f64], reduced: &[f64]) -> bool {
|
||||
let (inp, x) = (at.inp, &at.inp.solution);
|
||||
// Feasibility: each row within ten times the constraint precision of its own right side.
|
||||
// The nonlinear method accepts a solution within that of the model's largest constant right
|
||||
// side, which lets a small row go by its own size; a row over by more than its own
|
||||
// precision isn't satisfied, whatever the other rows' sizes, so the report refuses it.
|
||||
let rows_hold = rows.iter().enumerate().all(|(i, r)| {
|
||||
let (g, b) = (at.point.g.get(i).copied().unwrap_or(f64::NAN), at.finals.get(i).map_or(f64::NAN, |f| f.1));
|
||||
let over = match r.sense {
|
||||
Sense::Le => g,
|
||||
Sense::Ge => -g,
|
||||
Sense::Eq => g.abs(),
|
||||
};
|
||||
over <= 10.0 * super::PRECISION * (1.0 + b.abs())
|
||||
});
|
||||
if !rows_hold || !at.point.f.is_finite() {
|
||||
return false;
|
||||
}
|
||||
// Stationarity: each reduced gradient relative to the terms it is made of (the objective's
|
||||
// component and the multipliers' parts, at least 1). At optima the nonlinear method solves,
|
||||
// a reduced gradient is below about 1E-6 of those terms; 1E-3 leaves a wide margin for the
|
||||
// method's stopping test and for rounding in the difference gradients of badly scaled
|
||||
// models, and still refuses a point whose objective improves visibly along a direction the
|
||||
// constraints allow. The bound is the variable's own: with a bound from the gradient's
|
||||
// largest component, max 1E9·x + y at (1, 0.5) passed, although raising y gains 1 per unit
|
||||
// against a tolerance of 1E6.
|
||||
// Moving up a variable improves a maximum by +r and a minimum by −r per unit.
|
||||
let sign = if inp.model.goal == Goal::Max { 1.0 } else { -1.0 };
|
||||
x.iter().enumerate().all(|(j, &v)| {
|
||||
let (l, h) = (inp.prob.lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), inp.prob.hi.get(j).copied().unwrap_or(f64::INFINITY));
|
||||
let near = super::PRECISION * (1.0 + v.abs());
|
||||
if !v.is_finite() || v < l - near || v > h + near {
|
||||
return false;
|
||||
}
|
||||
// The same test for "at a bound" as the choice of free variables above.
|
||||
let edge = 1e-6 * v.abs().max(1.0);
|
||||
let (at_lo, at_hi) = (v <= l + edge, v >= h - edge);
|
||||
let gain = reduced.get(j).map_or(f64::NAN, |r| sign * r);
|
||||
let tol = 1e-3 * terms.get(j).copied().unwrap_or(0.0).max(1.0);
|
||||
gain.is_finite() && (at_hi || gain <= tol) && (at_lo || gain >= -tol)
|
||||
})
|
||||
}
|
||||
|
||||
/// GRG Nonlinear: reduced gradients and the binding rows' Lagrange multipliers.
|
||||
fn nonlinear_body(ev: &mut Evaluator, at: &At) -> Result<Body> {
|
||||
let (inp, m) = (at.inp, at.inp.model);
|
||||
let x = &inp.solution;
|
||||
let (df, dg) = gradients(ev, inp)?;
|
||||
// Binding rows (within the nonlinear method's accuracy) and the variables strictly inside
|
||||
// their bounds.
|
||||
let active: Vec<usize> = ev
|
||||
.rows
|
||||
.iter()
|
||||
.enumerate()
|
||||
.filter(|(i, r)| {
|
||||
let (g, b) = (at.point.g.get(*i).copied().unwrap_or(f64::NAN), at.finals.get(*i).map_or(0.0, |f| f.1));
|
||||
r.sense == Sense::Eq || g.abs() <= 1e-4 * (1.0 + b.abs())
|
||||
})
|
||||
.map(|(i, _)| i)
|
||||
.collect();
|
||||
let free: Vec<usize> = (0..x.len())
|
||||
.filter(|&j| {
|
||||
let v = x.get(j).copied().unwrap_or(0.0);
|
||||
let tol = 1e-6 * v.abs().max(1.0);
|
||||
v > inp.prob.lo.get(j).copied().unwrap_or(f64::NEG_INFINITY) + tol && v < inp.prob.hi.get(j).copied().unwrap_or(f64::INFINITY) - tol
|
||||
})
|
||||
.collect();
|
||||
let (lambda, lambda_noise) = multipliers(&ev.rows, m.goal, active, &free, (&df, &dg), MULTIPLIER_WORK)?;
|
||||
let grad = |i: usize, j: usize| dg.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
let reduced: Vec<f64> =
|
||||
(0..x.len()).map(|j| df.get(j).copied().unwrap_or(0.0) - lambda.iter().enumerate().map(|(i, l)| l * grad(i, j)).sum::<f64>()).collect();
|
||||
// A reduced gradient dfⱼ − Σᵢ λᵢ·gᵢⱼ is a sum of the gradient's own component and the
|
||||
// multipliers' terms; the sizes of those terms bound its rounding error and scale the test
|
||||
// of whether it is zero.
|
||||
let terms: Vec<f64> = (0..x.len())
|
||||
.map(|j| df.get(j).copied().unwrap_or(0.0).abs() + lambda.iter().enumerate().map(|(i, l)| (l * grad(i, j)).abs()).sum::<f64>())
|
||||
.collect();
|
||||
// The report describes the optimum, so the values it shows must be one.
|
||||
if !kkt_point(at, &ev.rows, &terms, &reduced) {
|
||||
return Err(not_optimal());
|
||||
}
|
||||
// Difference gradients are good to about 1E-6 of each component, so a reduced gradient's
|
||||
// noise is 1E-6 of its terms' sizes, plus each multiplier's own noise
|
||||
// ([`multiplier_noise`]) times |gᵢⱼ|. The bound is the variable's own: a small gradient
|
||||
// keeps its value beside a large component of another variable.
|
||||
let reduced_noise = |j: usize| {
|
||||
let carried: f64 = lambda_noise.iter().enumerate().map(|(i, n)| if grad(i, j) == 0.0 { 0.0 } else { n * grad(i, j).abs() }).sum();
|
||||
1e-6 * terms.get(j).copied().unwrap_or(0.0) + carried
|
||||
};
|
||||
let var_rows = m
|
||||
.vars
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(j, c)| {
|
||||
let mut row = at.var_head(j, *c);
|
||||
row.push(sens_num(reduced.get(j).copied().unwrap_or(f64::NAN), reduced_noise(j)));
|
||||
row
|
||||
})
|
||||
.collect();
|
||||
let con_rows = ev
|
||||
.rows
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, r)| {
|
||||
let mut row = at.row_head(i, r);
|
||||
row.push(sens_num(lambda.get(i).copied().unwrap_or(0.0), lambda_noise.get(i).copied().unwrap_or(0.0)));
|
||||
row
|
||||
})
|
||||
.collect();
|
||||
Ok((var_rows, con_rows, vec!["Cell", "Name", "Final Value", "Reduced Gradient"], vec!["Cell", "Name", "Final Value", "Lagrange Multiplier"]))
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::{Goal, MULTIPLIER_WORK, Row, Sense, multipliers};
|
||||
use gridcraft_core::CellRef;
|
||||
|
||||
#[test]
|
||||
fn multipliers_at_the_constraint_limit_stop_at_the_work_limit() {
|
||||
// 1,000 binding `<=` rows over 200 free variables, each multiplier with the wrong sign:
|
||||
// without a limit every pass would drop one row and solve again, about 4·10⁸ operations
|
||||
// a pass for up to 1,000 passes. The limit (lowered in tests) refuses it at once.
|
||||
let (k, f) = (1000usize, 200usize);
|
||||
let rows: Vec<Row> = (0..k).map(|_| Row { lhs: CellRef::new(0, 0), sense: Sense::Le, rhs: None, rhs_const: 0.0 }).collect();
|
||||
let dg: Vec<Vec<f64>> = (0..k).map(|i| (0..f).map(|j| if j == i % f { 1.0 } else { 0.0 }).collect()).collect();
|
||||
let df = vec![-1.0; f];
|
||||
let free: Vec<usize> = (0..f).collect();
|
||||
let got = multipliers(&rows, Goal::Max, (0..k).collect(), &free, (&df, &dg), MULTIPLIER_WORK);
|
||||
assert!(got.is_err());
|
||||
// A small model within the limit still gets its multipliers: ∇f = (1, 2) on two
|
||||
// binding rows with unit gradients.
|
||||
let dg2 = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
|
||||
let (got, noise) = multipliers(&rows[..2], Goal::Max, vec![0, 1], &free[..2], (&[1.0, 2.0], &dg2), MULTIPLIER_WORK).unwrap();
|
||||
assert_eq!(got, vec![1.0, 2.0]);
|
||||
// Each multiplier's noise is 1E-6 of its own row's term (1 · 1 and 1 · 2 over a unit
|
||||
// norm): the rows are orthogonal, so neither enters the other's bound.
|
||||
assert_eq!(noise, vec![1e-6, 2e-6]);
|
||||
}
|
||||
}
|
||||
@@ -283,6 +283,8 @@ fn refused_solver_reports_preserve_the_document_and_history() {
|
||||
(true, "simplex", json!(["A1 <= 4"]), json!(["answer"]), "workbook is protected"),
|
||||
(false, "simplex", json!(["A1 <= 4", "A1 int"]), json!(["answer", "limits"]), "integer or binary"),
|
||||
(false, "evolutionary", json!(["A1 <= 4"]), json!(["answer", "limits"]), "Simplex LP and GRG"),
|
||||
(false, "simplex", json!(["A1 <= 4", "A1 int"]), json!(["answer", "sensitivity"]), "integer or binary"),
|
||||
(false, "evolutionary", json!(["A1 <= 4"]), json!(["sensitivity"]), "Simplex LP and GRG"),
|
||||
] {
|
||||
for (keep, save) in [(true, true), (false, true), (true, false), (false, false)] {
|
||||
let mut s = s();
|
||||
@@ -400,6 +402,823 @@ fn solver_report_insertion_keeps_variable_edit_guards() {
|
||||
}
|
||||
}
|
||||
|
||||
/// The numbers from column D on in the row below the one holding `text` on report `name`.
|
||||
fn report_row(s: &Session, name: &str, text: &str, cols: &[&str]) -> Vec<f64> {
|
||||
let r = row_of(s, name, text) + 1;
|
||||
cols.iter().map(|c| at(s, name, &format!("{c}{r}")).as_f64().unwrap()).collect()
|
||||
}
|
||||
|
||||
fn close(got: &[f64], want: &[f64], tol: f64) -> bool {
|
||||
got.len() == want.len() && got.iter().zip(want).all(|(g, w)| (g - w).abs() <= tol * (1.0 + w.abs()))
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sensitivity_report_gives_the_wyndor_glass_shadow_prices_and_ranges() {
|
||||
// Hillier & Lieberman's Wyndor Glass sensitivity analysis: shadow prices 0, 1.5 and 1
|
||||
// for the three plants; plant 1's hours may fall by 2 and rise without limit, plants 2's
|
||||
// and 3's by 6 either way; the doors' profit 3 may range over 0..7.5, the windows' 5 over
|
||||
// 2 and up.
|
||||
let mut s = s();
|
||||
s.execute(
|
||||
"range.setValues",
|
||||
json!({"range": "A1", "values": [[null, "Doors", "Windows", "Used", "Available"], ["Batches", 0, 0, null, null], ["Plant 1", 1, 0, null, 4], ["Plant 2", 0, 2, null, 12], ["Plant 3", 3, 2, null, 18], ["Profit", 3, 5, null, null]]}),
|
||||
)
|
||||
.unwrap();
|
||||
for r in 3..=6 {
|
||||
set(&mut s, &[(&format!("D{r}"), &format!("=SUMPRODUCT(B{r}:C{r},$B$2:$C$2)"))]);
|
||||
}
|
||||
let solved = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "D6", "variables": "B2:C2", "constraints": ["D3:D5 <= E3:E5"], "method": "simplex", "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
assert_eq!(solved["reports"], json!(["Sensitivity Report 1"]), "{solved}");
|
||||
let rep = "Sensitivity Report 1";
|
||||
assert_eq!(at(&s, rep, "A1"), Value::text("GridCraft Sensitivity Report"));
|
||||
let r = row_of(&s, rep, "Reduced Cost");
|
||||
assert_eq!(at(&s, rep, &format!("F{}", r + 1)), Value::text("Objective Coefficient"));
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
// Final value, reduced cost, objective coefficient, allowable increase and decrease.
|
||||
assert!(close(&report_row(&s, rep, "$B$2", &cols)[1..], &[0.0, 3.0, 4.5, 3.0], 1e-9));
|
||||
assert!(close(&report_row(&s, rep, "$C$2", &cols), &[6.0, 0.0, 5.0, 1e30, 3.0], 1e-9));
|
||||
// Final value, shadow price, right side, allowable increase and decrease.
|
||||
assert!(close(&report_row(&s, rep, "$D$3", &cols), &[2.0, 0.0, 4.0, 1e30, 2.0], 1e-9));
|
||||
assert!(close(&report_row(&s, rep, "$D$4", &cols), &[12.0, 1.5, 12.0, 6.0, 6.0], 1e-9));
|
||||
assert!(close(&report_row(&s, rep, "$D$5", &cols), &[18.0, 1.0, 18.0, 6.0, 6.0], 1e-9));
|
||||
// A minimisation with `>=` rows: min 2x + 3y, x + y ≥ 4, x + 3y ≥ 6 → (3, 1), cost 9.
|
||||
// The duals 1.5 and 0.5 solve u + v = 2, u + 3v = 3; the basis holds for 2 ≤ b₁ ≤ 6,
|
||||
// 4 ≤ b₂ ≤ 12, 1 ≤ c₁ ≤ 3 and 2 ≤ c₂ ≤ 6.
|
||||
set(&mut s, &[("H1", "0"), ("H2", "0"), ("I1", "=2*H1+3*H2"), ("I2", "=H1+H2"), ("I3", "=H1+3*H2")]);
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "I1", "goal": "min", "variables": "H1:H2", "constraints": ["I2 >= 4", "I3 >= 6"], "method": "simplex", "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$H$1", &cols), &[3.0, 0.0, 2.0, 1.0, 1.0], 1e-9), "{r}");
|
||||
assert!(close(&report_row(&s, &rep, "$H$2", &cols), &[1.0, 0.0, 3.0, 3.0, 1.0], 1e-9));
|
||||
assert!(close(&report_row(&s, &rep, "$I$2", &cols), &[4.0, 1.5, 4.0, 2.0, 2.0], 1e-9));
|
||||
assert!(close(&report_row(&s, &rep, "$I$3", &cols), &[6.0, 0.5, 6.0, 6.0, 2.0], 1e-9));
|
||||
// A variable held at zero has a reduced cost: min 2x + 3y with x + y ≥ 4 alone uses
|
||||
// only x; y would add 3 − 2 = 1 per unit.
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "I1", "goal": "min", "variables": "H1:H2", "constraints": ["I2 >= 4"], "method": "simplex", "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$H$2", &cols), &[0.0, 1.0, 3.0, 1e30, 1.0], 1e-9), "{r}");
|
||||
// Integer models, the Evolutionary method and Value Of goals have no Sensitivity report.
|
||||
s.execute("data.solver", json!({"objective": "D6", "variables": "B2:C2", "constraints": ["D3:D5 <= E3:E5"], "method": "simplex"})).unwrap();
|
||||
for p in [
|
||||
json!({"reports": ["sensitivity"], "constraints": ["D3:D5 <= E3:E5", "B2:C2 int"]}),
|
||||
json!({"reports": ["sensitivity"], "method": "evolutionary", "constraints": ["D3:D5 <= E3:E5", "B2:C2 <= 10"]}),
|
||||
json!({"reports": ["sensitivity"], "goal": "value", "value": 30}),
|
||||
] {
|
||||
assert!(s.execute("data.solverReport", p.clone()).is_err(), "{p}");
|
||||
}
|
||||
// data.solver refuses a Value Of goal's Sensitivity report before it changes anything.
|
||||
set(&mut s, &[("B2", "1"), ("C2", "1")]);
|
||||
let undo = s.doc().unwrap().undo.len();
|
||||
let p = json!({"objective": "D6", "goal": "value", "value": 30, "variables": "B2:C2", "constraints": ["D3:D5 <= E3:E5"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let err = s.execute("data.solver", p).unwrap_err();
|
||||
assert!(err.to_string().contains("maximise or minimise"), "{err}");
|
||||
assert_eq!((at(&s, "Sheet1", "B2"), at(&s, "Sheet1", "C2")), (Value::number(1.0), Value::number(1.0)));
|
||||
assert_eq!(s.doc().unwrap().undo.len(), undo);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_ranges_of_a_redundant_equality_are_zero() {
|
||||
// max x with x = 2 and 2x = 4: the second row repeats the first, so any change to either
|
||||
// right side alone has no solution, and both ranges of each row are 0. With x = 2 alone the
|
||||
// right side may rise without limit and fall by 2 (to x = 0).
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("B1", "=A1"), ("C1", "=A1"), ("C2", "=2*A1")]);
|
||||
let cols = ["D", "F", "G", "H"];
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "B1", "variables": "A1", "constraints": ["C1 = 2", "C2 = 4"], "method": "simplex", "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$C$1", &cols), &[2.0, 2.0, 0.0, 0.0], 1e-9), "{:?}", report_row(&s, &rep, "$C$1", &cols));
|
||||
assert!(close(&report_row(&s, &rep, "$C$2", &cols), &[4.0, 4.0, 0.0, 0.0], 1e-9), "{:?}", report_row(&s, &rep, "$C$2", &cols));
|
||||
// The redundant row doesn't restrict the objective coefficient: x = 2 stays the only
|
||||
// solution whatever its profit.
|
||||
assert!(close(&report_row(&s, &rep, "$A$1", &["D", "F", "G", "H"]), &[2.0, 1.0, 1e30, 1e30], 1e-9));
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "B1", "variables": "A1", "constraints": ["C1 = 2"], "method": "simplex", "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$C$1", &cols), &[2.0, 2.0, 1e30, 2.0], 1e-9), "{:?}", report_row(&s, &rep, "$C$1", &cols));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_report_needs_an_optimal_point() {
|
||||
// Wyndor Glass: (2, 6) is the optimum (36). (0, 0) is feasible but not optimal; (4, 4.8)
|
||||
// reaches 36 but breaks plant 3's hours (21.6 > 18).
|
||||
let mut s = s();
|
||||
s.execute(
|
||||
"range.setValues",
|
||||
json!({"range": "A1", "values": [[null, "Doors", "Windows", "Used", "Available"], ["Batches", 0, 0, null, null], ["Plant 1", 1, 0, null, 4], ["Plant 2", 0, 2, null, 12], ["Plant 3", 3, 2, null, 18], ["Profit", 3, 5, null, null]]}),
|
||||
)
|
||||
.unwrap();
|
||||
for r in 3..=6 {
|
||||
set(&mut s, &[(&format!("D{r}"), &format!("=SUMPRODUCT(B{r}:C{r},$B$2:$C$2)"))]);
|
||||
}
|
||||
let model = json!({"objective": "D6", "variables": "B2:C2", "constraints": ["D3:D5 <= E3:E5"], "method": "simplex"});
|
||||
let report = |b: f64, c: f64| json!({"reports": ["sensitivity"], "model": model, "values": {"B2": b, "C2": c}});
|
||||
for (b, c) in [(0.0, 0.0), (4.0, 4.8)] {
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", report(b, c)).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "({b}, {c}): {err}");
|
||||
assert!(std::sync::Arc::ptr_eq(&s.doc().unwrap().wb, &before.wb));
|
||||
}
|
||||
let r = s.execute("data.solverReport", report(2.0, 6.0)).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_report_describes_the_vertex_the_solver_returns() {
|
||||
// max x + y with x + y <= 4, x <= 3, y <= 3: every point from (3, 1) to (1, 3) is optimal,
|
||||
// but the ranges belong to one vertex. The vertex Simplex LP returns is reported (from
|
||||
// several starting values, as Solver Results passes the solve's result on); the other vertex
|
||||
// and the midpoint are refused.
|
||||
let model = json!({"objective": "D1", "variables": "A1:A2", "constraints": ["C1 <= 4", "A1 <= 3", "A2 <= 3"], "method": "simplex"});
|
||||
let mut found = Vec::new();
|
||||
for (x, y) in [("0", "0"), ("0.1", "0.7"), ("2.9", "0.3"), ("1", "1"), ("-5", "12.25")] {
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", x), ("A2", y), ("C1", "=A1+A2"), ("D1", "=A1+A2")]);
|
||||
let mut p = model.clone();
|
||||
p["keep"] = json!(false);
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": r["model"], "values": r["values"]})).unwrap();
|
||||
assert_eq!(rep["sheets"], json!(["Sensitivity Report 1"]), "from ({x}, {y}): {rep}");
|
||||
found.push((r["values"]["A1"].as_f64().unwrap(), r["values"]["A2"].as_f64().unwrap()));
|
||||
}
|
||||
let (vx, vy) = found[0];
|
||||
assert!(found.iter().all(|v| *v == (vx, vy)), "{found:?}");
|
||||
let other = if vx == 3.0 { (1.0, 3.0) } else { (3.0, 1.0) };
|
||||
assert!([(3.0, 1.0), (1.0, 3.0)].contains(&(vx, vy)), "{found:?}");
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("C1", "=A1+A2"), ("D1", "=A1+A2")]);
|
||||
for (x, y) in [other, (2.0, 2.0)] {
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": x, "A2": y}})).unwrap_err();
|
||||
assert!(err.to_string().contains("another optimal solution"), "({x}, {y}): {err}");
|
||||
assert!(std::sync::Arc::ptr_eq(&s.doc().unwrap().wb, &before.wb));
|
||||
// data.solver's Skip mode (as Solver Results passes) skips only the Sensitivity report.
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solverReport",
|
||||
json!({"reports": ["answer", "sensitivity"], "model": model, "values": {"A1": x, "A2": y}, "skipNotOptimal": true}),
|
||||
)
|
||||
.unwrap();
|
||||
assert_eq!(r["skippedReports"][0]["report"], "sensitivity", "{r}");
|
||||
assert_eq!(r["sheets"].as_array().map(Vec::len), Some(1), "{r}");
|
||||
}
|
||||
let r = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": vx, "A2": vy}})).unwrap();
|
||||
assert!(r["sheets"][0].as_str().unwrap().starts_with("Sensitivity Report"), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_report_repeats_the_solve_from_its_start_values() {
|
||||
// max 0.1x + 0.1y with 0.1x + 0.1y <= 0.4, 0.7x + 0.7y <= 2.8, x <= 3, y <= 3: the optimum is
|
||||
// not unique, and the difference coefficients of these formulas differ in their last bits
|
||||
// with the point they are taken at. From these starts the solve returns (3, 1) or (1, 3);
|
||||
// linearised at that solution instead, the report's solve broke the tie the other way and
|
||||
// refused the solve's own values. Built from the same start values, the program is the
|
||||
// same to the bit, so each solve's own report is written.
|
||||
let model =
|
||||
json!({"objective": "D1", "variables": "A1:A2", "constraints": ["C1 <= 0.4", "E1 <= 2.8", "A1 <= 3", "A2 <= 3"], "method": "simplex"});
|
||||
let formulas = [("C1", "=A1*0.1+A2*0.1"), ("D1", "=A1*0.1+A2*0.1"), ("E1", "=A1*0.7+A2*0.7")];
|
||||
let program = |s: &Session, x: &[f64]| {
|
||||
let m = super::parse_model(s, &model).unwrap();
|
||||
let prob = super::build_problem(&m);
|
||||
let mut ev = super::Evaluator {
|
||||
trial: super::Trial::new(s.doc().unwrap()),
|
||||
sheet: m.sheet,
|
||||
vars: m.vars.clone(),
|
||||
objective: m.objective,
|
||||
rows: prob.rows.clone(),
|
||||
evals: 0,
|
||||
budget: 100,
|
||||
target: None,
|
||||
};
|
||||
let p = super::lp::simplex_program(&mut ev, m.goal, x).unwrap();
|
||||
let bits = |v: &[f64]| v.iter().map(|x| x.to_bits()).collect::<Vec<u64>>();
|
||||
let rows: Vec<(Vec<u64>, u64)> = p.rows.iter().map(|(a, _, b)| (bits(a), b.to_bits())).collect();
|
||||
(bits(&p.cf), bits(&p.c), rows, p.work)
|
||||
};
|
||||
let mut vertices = Vec::new();
|
||||
for (x, y) in [("-3", "-5"), ("-3", "2.42"), ("-3", "6.13"), ("-3", "9.84"), ("-3", "13.55"), ("0.37", "0.3"), ("2.55", "-1.2")] {
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", x), ("A2", y)]);
|
||||
set(&mut s, &formulas);
|
||||
let start: Vec<f64> = [x, y].iter().map(|v| v.parse().unwrap()).collect();
|
||||
let before = program(&s, &start);
|
||||
let r = s.execute("data.solver", model.clone()).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
// The program the report builds after the solve (cells now at the solution) from the
|
||||
// solve's `original`, as Solver Results passes it, equals the solve's.
|
||||
let original: Vec<f64> = ["A1", "A2"].iter().map(|c| r["original"][c].as_f64().unwrap()).collect();
|
||||
assert_eq!(original, start);
|
||||
assert_eq!(program(&s, &original), before, "from ({x}, {y})");
|
||||
let rep = s
|
||||
.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": r["model"], "values": r["values"], "original": r["original"]}))
|
||||
.unwrap();
|
||||
assert_eq!(rep["sheets"], json!(["Sensitivity Report 1"]), "from ({x}, {y}): {rep}");
|
||||
let v = (r["values"]["A1"].as_f64().unwrap(), r["values"]["A2"].as_f64().unwrap());
|
||||
vertices.push(if v.0 > v.1 { "(3, 1)" } else { "(1, 3)" });
|
||||
}
|
||||
// Both vertices occur, so the starts exercise the tie.
|
||||
assert!(vertices.contains(&"(3, 1)") && vertices.contains(&"(1, 3)"), "{vertices:?}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_report_refuses_a_point_that_beats_the_optimum() {
|
||||
// min 1E+9·A1 with A1 ≥ 0: the optimum is 0. A1 = -1E-7 is outside the bound by less than
|
||||
// the precision, yet its objective is -100, far from the optimum.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("B1", "=1E+9*A1")]);
|
||||
let model = json!({"objective": "B1", "goal": "min", "variables": "A1", "constraints": ["A1 <= 5"], "method": "simplex"});
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": -1e-7}})).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "{err}");
|
||||
assert!(std::sync::Arc::ptr_eq(&s.doc().unwrap().wb, &before.wb));
|
||||
let r = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": 0}})).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sensitivity_report_of_a_nonlinear_model_gives_lagrange_multipliers() {
|
||||
// The largest rectangle with perimeter at most 20: the area (P/4)² grows by P/8 = 2.5 per
|
||||
// unit of perimeter, and both reduced gradients are 0.
|
||||
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"], "reports": ["sensitivity"]})).unwrap();
|
||||
let rep = "Sensitivity Report 1";
|
||||
assert_eq!(r["reports"], json!([rep]), "{r}");
|
||||
assert_eq!(at(&s, rep, &format!("E{}", row_of(&s, rep, "$A$1"))), Value::text("Reduced Gradient"));
|
||||
assert!(close(&report_row(&s, rep, "$A$1", &["D", "E"]), &[5.0, 0.0], 1e-4));
|
||||
assert!(close(&report_row(&s, rep, "$B$2", &["D", "E"]), &[20.0, 2.5], 1e-4));
|
||||
// min (x + 1)² + (y − 2)² with x + y ≤ 1 and x, y ≥ 0: (0, 1). The constraint's
|
||||
// multiplier is −2; x sits at its bound, where raising it costs 2(x + 1) + 2 = 4.
|
||||
set(&mut s, &[("C1", "0.5"), ("C2", "0.5"), ("D1", "=(C1+1)^2+(C2-2)^2"), ("D2", "=C1+C2")]);
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "D1", "goal": "min", "variables": "C1:C2", "constraints": ["D2 <= 1"], "reports": ["sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$C$1", &["D", "E"]), &[0.0, 4.0], 1e-4), "{r}");
|
||||
assert!(close(&report_row(&s, &rep, "$C$2", &["D", "E"]), &[1.0, 0.0], 1e-4));
|
||||
assert!(close(&report_row(&s, &rep, "$D$2", &["D", "E"]), &[1.0, -2.0], 1e-4));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_report_needs_an_optimal_point() {
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "1"), ("B1", "=A1^2")]);
|
||||
set(&mut s, &[("C1", "0"), ("D1", "=(C1-2)^2")]);
|
||||
set(&mut s, &[("E1", "1"), ("E2", "1"), ("F1", "=E1*E2"), ("F2", "=2*E1+2*E2")]);
|
||||
set(&mut s, &[("G1", "0"), ("G2", "1"), ("H1", "=(G1+1)^2+(G2-2)^2"), ("H2", "=G1+G2")]);
|
||||
let square = json!({"objective": "B1", "goal": "min", "variables": "A1", "constraints": ["A1 <= 5"]});
|
||||
let shifted = json!({"objective": "D1", "goal": "min", "variables": "C1", "constraints": ["C1 <= 5"]});
|
||||
let rectangle = json!({"objective": "F1", "variables": "E1:E2", "constraints": ["F2 <= 20"]});
|
||||
let bounded = json!({"objective": "H1", "goal": "min", "variables": "G1:G2", "constraints": ["H2 <= 1"]});
|
||||
let report = |model: &serde_json::Value, values: serde_json::Value| json!({"reports": ["sensitivity"], "model": model, "values": values});
|
||||
for (model, values, why) in [
|
||||
// min A1² at A1 = 1: the free variable's reduced gradient is 2, not 0.
|
||||
(&square, json!({"A1": 1}), "free gradient"),
|
||||
// The rectangle at 6 × 6 has a larger area than the optimum but a perimeter of 24 > 20.
|
||||
(&rectangle, json!({"E1": 6, "E2": 6}), "infeasible"),
|
||||
// At 4 × 6 the perimeter binds, but its multiplier can't explain the gradient (6, 4).
|
||||
(&rectangle, json!({"E1": 4, "E2": 6}), "unexplained gradient"),
|
||||
// min (C1 − 2)² at the bound C1 = 0: raising C1 lowers the objective.
|
||||
(&shifted, json!({"C1": 0}), "wrong-way bound"),
|
||||
// (1, 0) in the bounded model: the binding row's multiplier would be +4, the wrong sign
|
||||
// for a `<=` row when minimising, so the row is dropped and the gradient is unexplained.
|
||||
(&bounded, json!({"G1": 1, "G2": 0}), "wrong-sign multiplier"),
|
||||
] {
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", report(model, values)).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "{why}: {err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb), "{why}");
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()), "{why}");
|
||||
}
|
||||
// The optima are still reported.
|
||||
for (k, (model, values)) in
|
||||
[(&square, json!({"A1": 0})), (&shifted, json!({"C1": 2})), (&rectangle, json!({"E1": 5, "E2": 5})), (&bounded, json!({"G1": 0, "G2": 1}))]
|
||||
.into_iter()
|
||||
.enumerate()
|
||||
{
|
||||
let r = s.execute("data.solverReport", report(model, values)).unwrap();
|
||||
assert_eq!(r["sheets"], json!([format!("Sensitivity Report {}", k + 1)]), "{r}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn solver_writes_the_other_reports_when_the_run_stops_before_an_optimum() {
|
||||
// Rosenbrock's function from (-1.2, 1) with 60 evaluations: the run stops at the limit with
|
||||
// a solution that is not the optimum (1, 1). The solve applies and the Answer report is
|
||||
// written; the Sensitivity report is left out with its reason.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "-1.2"), ("A2", "1"), ("B1", "=(1-A1)^2+100*(A2-A1^2)^2")]);
|
||||
let undo = s.doc().unwrap().undo.len();
|
||||
let r = s
|
||||
.execute(
|
||||
"data.solver",
|
||||
json!({"objective": "B1", "goal": "min", "variables": "A1:A2", "nonNegative": false, "maxEvaluations": 60, "reports": ["answer", "sensitivity"]}),
|
||||
)
|
||||
.unwrap();
|
||||
assert_eq!(r["status"], "iterationLimit", "{r}");
|
||||
assert_eq!(r["changed"], true);
|
||||
assert_eq!(n(&s, "A1"), r["values"]["A1"].as_f64().unwrap());
|
||||
assert_eq!(r["reports"], json!(["Answer Report 1"]), "{r}");
|
||||
assert_eq!(r["skippedReports"][0]["report"], "sensitivity", "{r}");
|
||||
assert!(r["skippedReports"][0]["reason"].as_str().unwrap().contains("optimal"), "{r}");
|
||||
assert_eq!(s.doc().unwrap().undo.len(), undo + 1);
|
||||
// A solve that reaches the optimum still writes the Sensitivity report and skips nothing.
|
||||
let r = s
|
||||
.execute("data.solver", json!({"objective": "B1", "goal": "min", "variables": "A1:A2", "nonNegative": false, "reports": ["sensitivity"]}))
|
||||
.unwrap();
|
||||
assert_eq!(r["status"], "converged", "{r}");
|
||||
assert_eq!(r["reports"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
assert!(r.get("skippedReports").is_none(), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn solver_report_can_skip_a_sensitivity_report_the_values_do_not_support() {
|
||||
// min A1² at A1 = 1 is not an optimum. Called directly, data.solverReport refuses it; with
|
||||
// skipNotOptimal (as Solver Results passes) the Answer report is still written.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "1"), ("B1", "=A1^2")]);
|
||||
let model = json!({"objective": "B1", "goal": "min", "variables": "A1", "constraints": ["A1 <= 5"]});
|
||||
let p = |reports: serde_json::Value, skip: bool| json!({"reports": reports, "model": model, "values": {"A1": 1}, "skipNotOptimal": skip});
|
||||
let err = s.execute("data.solverReport", p(json!(["answer", "sensitivity"]), false)).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "{err}");
|
||||
let undo = s.doc().unwrap().undo.len();
|
||||
let r = s.execute("data.solverReport", p(json!(["answer", "sensitivity"]), true)).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Answer Report 1"]), "{r}");
|
||||
assert!(r["skippedReports"][0]["reason"].as_str().unwrap().contains("not optimal"), "{r}");
|
||||
assert_eq!(s.doc().unwrap().undo.len(), undo + 1);
|
||||
// With nothing left to write, nothing changes and no Undo step is added.
|
||||
let before = s.doc().unwrap().clone();
|
||||
let r = s.execute("data.solverReport", p(json!(["sensitivity"]), true)).unwrap();
|
||||
assert_eq!(r["sheets"], json!([]), "{r}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
// A model with no Sensitivity report is still refused, skipNotOptimal or not.
|
||||
let err = s.execute("data.solverReport", json!({"reports": ["answer", "sensitivity"], "model": {"objective": "B1", "variables": "A1", "constraints": ["A1 <= 5", "A1 int"]}, "skipNotOptimal": true})).unwrap_err();
|
||||
assert!(err.to_string().contains("integer or binary"), "{err}");
|
||||
// In a workbook whose structure is protected, a request whose only report is skipped
|
||||
// writes no sheet, so it succeeds; one that would write a sheet is still refused.
|
||||
let d = s.doc_mut().unwrap();
|
||||
std::sync::Arc::make_mut(&mut d.wb).protected_structure = true;
|
||||
let before = s.doc().unwrap().clone();
|
||||
let r = s.execute("data.solverReport", p(json!(["sensitivity"]), true)).unwrap();
|
||||
assert_eq!(r["sheets"], json!([]), "{r}");
|
||||
assert!(r["skippedReports"][0]["reason"].as_str().unwrap().contains("not optimal"), "{r}");
|
||||
let err = s.execute("data.solverReport", p(json!(["answer", "sensitivity"]), true)).unwrap_err();
|
||||
assert!(err.to_string().contains("protected"), "{err}");
|
||||
// A request that must write a sheet is refused before its reports are built, so an
|
||||
// integer model's Limits report gives the protection message, not the report's own refusal.
|
||||
let int_model = json!({"objective": "B1", "variables": "A1", "constraints": ["A1 <= 5", "A1 int"]});
|
||||
let err = s.execute("data.solverReport", json!({"reports": ["limits"], "model": int_model, "values": {"A1": 1}})).unwrap_err();
|
||||
assert!(err.to_string().contains("protected"), "{err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_sensitivity_refusal_of_the_model_comes_before_any_report_is_built() {
|
||||
// The refusals that apply to the model (integer or binary constraints, the Evolutionary
|
||||
// method, a Value Of goal) are checked before any report is built, in both modes and from
|
||||
// both commands: a Limits report listed first would otherwise run its two optimisations per
|
||||
// variable, and data.solver its whole solve, for an error. An integer model asking for
|
||||
// limits and sensitivity gets the Sensitivity refusal, not the Limits report's, which the
|
||||
// report loop reached first; a Value Of model, whose Limits report is not refused, gets the
|
||||
// Sensitivity refusal with no sheet written and the workbook as it was.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "1"), ("B1", "=2*A1")]);
|
||||
let before = s.doc().unwrap().clone();
|
||||
let int_model = json!({"objective": "B1", "variables": "A1", "constraints": ["A1 <= 5", "A1 int"]});
|
||||
let value_model = json!({"objective": "B1", "goal": "value", "value": 4, "variables": "A1", "constraints": ["A1 <= 5"]});
|
||||
for (model, values, why) in
|
||||
[(&int_model, json!({"A1": 1}), "The Sensitivity report doesn't apply"), (&value_model, json!({"A1": 2}), "maximise or minimise")]
|
||||
{
|
||||
for skip in [false, true] {
|
||||
let p = json!({"reports": ["limits", "sensitivity"], "model": model, "values": values, "skipNotOptimal": skip});
|
||||
let err = s.execute("data.solverReport", p).unwrap_err();
|
||||
assert!(err.to_string().contains(why), "{skip}: {err}");
|
||||
}
|
||||
}
|
||||
let p = json!({"objective": "B1", "variables": "A1", "constraints": ["A1 <= 5", "A1 int"], "method": "simplex", "reports": ["limits", "sensitivity"]});
|
||||
let err = s.execute("data.solver", p).unwrap_err();
|
||||
assert!(err.to_string().contains("The Sensitivity report doesn't apply"), "{err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_with_too_many_binding_constraints_is_refused_without_change() {
|
||||
// 80 binding rows on two free variables pass the work limit (lowered in tests); 40 don't.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "5"), ("A2", "5"), ("B1", "=A1*A2")]);
|
||||
for r in 1..=80 {
|
||||
set(&mut s, &[(&format!("D{r}"), "=A1+A2")]);
|
||||
}
|
||||
let report = |rows: u32| json!({"reports": ["sensitivity"], "model": {"objective": "B1", "variables": "A1:A2", "constraints": [format!("D1:D{rows} <= 10")]}});
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", report(80)).unwrap_err();
|
||||
assert!(err.to_string().contains("too many binding constraints"), "{err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
let r = s.execute("data.solverReport", report(40)).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sensitivity_applies_agrees_with_the_report() {
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "1"), ("B1", "=3*A1")]);
|
||||
for (extra, ok) in [
|
||||
(json!({}), true),
|
||||
(json!({"goal": "min"}), true),
|
||||
(json!({"method": "grg"}), true),
|
||||
(json!({"goal": "value", "value": 3}), false),
|
||||
(json!({"constraints": ["A1 <= 4", "A1 int"]}), false),
|
||||
(json!({"constraints": ["A1 <= 4", {"lhs": "A1", "op": "bin"}]}), false),
|
||||
(json!({"method": "evolutionary"}), false),
|
||||
] {
|
||||
let mut p = json!({"objective": "B1", "variables": "A1", "constraints": ["A1 <= 4"], "method": "simplex", "keep": false});
|
||||
for (k, v) in extra.as_object().unwrap() {
|
||||
p[k] = v.clone();
|
||||
}
|
||||
assert_eq!(super::sensitivity_applies(&p), ok, "{p}");
|
||||
let r = s.execute("data.solver", p.clone()).unwrap();
|
||||
assert_eq!(super::sensitivity_applies(&r["model"]), ok, "{p}");
|
||||
let rep = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": r["model"], "values": r["values"]}));
|
||||
assert_eq!(rep.is_ok(), ok, "{p}: {rep:?}");
|
||||
}
|
||||
}
|
||||
|
||||
/// How many cells of sheet `name` hold `text`.
|
||||
fn count_of(s: &Session, name: &str, text: &str) -> usize {
|
||||
let d = s.doc().unwrap();
|
||||
let sh = d.wb.sheet(d.wb.sheet_index(name).unwrap()).unwrap();
|
||||
(0..60).flat_map(|r| (0..8).map(move |c| (r, c))).filter(|(r, c)| sh.value(CellRef::new(*r, *c)) == Value::text(text)).count()
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sensitivity_report_is_the_same_whatever_constant_the_objective_adds() {
|
||||
// A constant in the objective changes no price: max 1E9 + x with x ≤ 10 (a row) has the
|
||||
// same report as max x with x ≤ 10. At x = 10 the row's shadow price is 1 (one more unit
|
||||
// of right side adds one to the objective); x is basic, with reduced cost 0 and
|
||||
// coefficient 1, which may rise without limit and fall by 1 (to 0, where the slack ties);
|
||||
// the right side 10 may rise without limit and fall by 10 (to x = 0). A tolerance taken
|
||||
// from the objective's value (1E-9 · 1E9 = 1) wrote the shadow price as 0.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("B1", "=1000000000+A1"), ("C1", "=A1"), ("E1", "0"), ("F1", "=E1"), ("G1", "=E1")]);
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let mut rows = Vec::new();
|
||||
for (obj, var, con, v, c) in [("B1", "A1", "C1", "$A$1", "$C$1"), ("F1", "E1", "G1", "$E$1", "$G$1")] {
|
||||
let p = json!({"objective": obj, "variables": var, "constraints": [format!("{con} <= 10")], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
rows.push((report_row(&s, &rep, v, &cols), report_row(&s, &rep, c, &cols)));
|
||||
}
|
||||
assert!(close(&rows[0].0, &[10.0, 0.0, 1.0, 1e30, 1.0], 1e-9), "{rows:?}");
|
||||
assert!(close(&rows[0].1, &[10.0, 1.0, 10.0, 1e30, 10.0], 1e-9), "{rows:?}");
|
||||
assert_eq!(rows[0], rows[1]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_report_is_the_same_whatever_constant_the_objective_adds() {
|
||||
// min 1E6 + (x − 3)² with x ≤ 2.75 has the same report as min (x − 3)² with x ≤ 2.75. At
|
||||
// x = 2.75 the objective's gradient is 2(2.75 − 3) = −0.5. As a bound on the variable cell
|
||||
// (no constraint row) it is x's reduced gradient, −0.5; as a row y = x ≤ 2.75 it is the
|
||||
// row's Lagrange multiplier, −0.5 (∇f = λ∇y with ∇y = 1), and x's reduced gradient is 0.
|
||||
// A tolerance taken from the objective's value (1E-6 · 1E6 = 1) wrote both as 0. The
|
||||
// difference gradients of the model with the constant lose about four digits to rounding,
|
||||
// so the figures agree to 1E-3.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("B1", "=1000000+(A1-3)^2"), ("C1", "=A1"), ("E1", "0"), ("F1", "=(E1-3)^2"), ("G1", "=E1")]);
|
||||
let mut figures = Vec::new();
|
||||
for (obj, var, con, v, c) in [("B1", "A1", "C1", "$A$1", "$C$1"), ("F1", "E1", "G1", "$E$1", "$G$1")] {
|
||||
let solve = |s: &mut Session, bound: String| {
|
||||
let p = json!({"objective": obj, "goal": "min", "variables": var, "constraints": [bound], "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "converged", "{r}");
|
||||
r["reports"][0].as_str().unwrap().to_string()
|
||||
};
|
||||
let rep = solve(&mut s, format!("{var} <= 2.75"));
|
||||
let at_bound = report_row(&s, &rep, v, &["D", "E"]);
|
||||
assert_eq!(count_of(&s, &rep, c), 0, "a bound has no constraint row");
|
||||
let rep = solve(&mut s, format!("{con} <= 2.75"));
|
||||
figures.push((at_bound, report_row(&s, &rep, v, &["D", "E"]), report_row(&s, &rep, c, &["D", "E"])));
|
||||
}
|
||||
let (at_bound, free, row) = &figures[0];
|
||||
assert!(close(at_bound, &[2.75, -0.5], 1e-3), "{figures:?}");
|
||||
assert!(close(free, &[2.75, 0.0], 1e-3), "{figures:?}");
|
||||
assert!(close(row, &[2.75, -0.5], 1e-3), "{figures:?}");
|
||||
assert!(close(&figures[1].0, at_bound, 1e-3) && close(&figures[1].1, free, 1e-3) && close(&figures[1].2, row, 1e-3), "{figures:?}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_keeps_a_small_price_beside_a_large_coefficient() {
|
||||
// max 1E9·x + y with x ≤ 1 and y ≤ 1 as rows and x, y ≥ 0: (1, 1). Both rows bind, the
|
||||
// basis is {x, y} with B = I, so the shadow prices are the coefficients, 1E9 and 1, and both
|
||||
// reduced costs are 0. x's coefficient may rise without limit and fall by 1E9 (to 0, where
|
||||
// its row's slack ties), y's without limit and by 1. Each right side may rise without limit
|
||||
// (the other row doesn't involve the variable) and fall by 1 (to 0). The second row's price
|
||||
// is the single term 1 · 1, so the first row's coefficient 1E9 is no rounding error in it; a
|
||||
// noise level taken from the largest coefficient (1E-9 · 1E9 = 1) wrote it as 0.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=1000000000*A1+A2"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let p = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["C1 <= 1", "C2 <= 1"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let row = |c: &str| report_row(&s, &rep, c, &cols);
|
||||
assert!(close(&row("$A$1"), &[1.0, 0.0, 1e9, 1e30, 1e9], 1e-9), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[1.0, 0.0, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$C$1"), &[1.0, 1e9, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$C$1"));
|
||||
assert!(close(&row("$C$2"), &[1.0, 1.0, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$C$2"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_keeps_a_small_reduced_cost_beside_a_large_coefficient() {
|
||||
// max 1E9·x + y − z with x ≤ 1 and y ≤ 1 as rows and x, y, z ≥ 0 (z's bound is the only
|
||||
// constraint on it): (1, 1, 0). z is nonbasic at zero with reduced cost −1 − 0 = −1 (no
|
||||
// row involves it): its coefficient must rise by 1 before z is made, so the Allowable
|
||||
// Increase is 1 and the Allowable Decrease 1E+30, as for a maximum's variable at zero in
|
||||
// `simplex_sensitivity_gives_the_dakota_furniture_report`. The other rows are those of
|
||||
// `simplex_sensitivity_keeps_a_small_price_beside_a_large_coefficient`. The reduced cost
|
||||
// is the coefficient alone, so x's coefficient 1E9 is no rounding error in it; a noise
|
||||
// level taken from the largest coefficient wrote it as 0.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("A3", "0"), ("B1", "=1000000000*A1+A2-A3"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let p = json!({"objective": "B1", "variables": "A1:A3", "constraints": ["C1 <= 1", "C2 <= 1"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let row = |c: &str| report_row(&s, &rep, c, &cols);
|
||||
assert!(close(&row("$A$1"), &[1.0, 0.0, 1e9, 1e30, 1e9], 1e-9), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[1.0, 0.0, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$A$3"), &[0.0, -1.0, -1.0, 1.0, 1e30], 1e-9), "{:?}", row("$A$3"));
|
||||
assert!(close(&row("$C$1"), &[1.0, 1e9, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$C$1"));
|
||||
assert!(close(&row("$C$2"), &[1.0, 1.0, 1.0, 1e30, 1.0], 1e-9), "{:?}", row("$C$2"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_keeps_a_small_multiplier_beside_a_large_gradient() {
|
||||
// The model of `simplex_sensitivity_keeps_a_small_price_beside_a_large_coefficient` under
|
||||
// GRG Nonlinear: at (1, 1) both rows bind; with ∇f = (1E9, 1) and unit row gradients the
|
||||
// Lagrange multipliers are 1E9 and 1, and both reduced gradients are 0. The second row's
|
||||
// multiplier is the single term 1 · 1, so the gradient's large first component is no
|
||||
// rounding error in it; a noise level taken from the gradient's size (1E-6 · 1E9 = 1000)
|
||||
// wrote it as 0. The difference gradient of y is the difference of two values near 1E9,
|
||||
// where the last place is 2⁻²² ≈ 2.4E-7: the step 1E-6 rounds to 4 places, 9.54E-7, so
|
||||
// df_y ≈ 0.954 and the multiplier is checked to 0.1. The nonlinear method stops on this
|
||||
// model from (0, 0) with "unbounded", so the report is asked for at the optimum (1, 1).
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=1000000000*A1+A2"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let model = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["C1 <= 1", "C2 <= 1"]});
|
||||
let r = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": 1, "A2": 1}})).unwrap();
|
||||
let rep = r["sheets"][0].as_str().unwrap().to_string();
|
||||
let row = |c: &str| report_row(&s, &rep, c, &["D", "E"]);
|
||||
assert!(close(&row("$A$1"), &[1.0, 0.0], 1e-3), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[1.0, 0.0], 1e-3), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$C$1"), &[1.0, 1e9], 1e-3), "{:?}", row("$C$1"));
|
||||
assert!(close(&row("$C$2"), &[1.0, 1.0], 0.1), "{:?}", row("$C$2"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_keeps_a_small_reduced_gradient_beside_a_large_gradient() {
|
||||
// The model of `simplex_sensitivity_keeps_a_small_reduced_cost_beside_a_large_coefficient`
|
||||
// under GRG Nonlinear: z sits at its bound 0, where its reduced gradient is its own
|
||||
// gradient −1 (no row involves it); a noise level taken from the gradient's size (1000)
|
||||
// wrote it as 0. Its difference gradient is one-sided at the bound, over the step 1E-6
|
||||
// from a value near 1E9 that rounds as in
|
||||
// `nonlinear_sensitivity_keeps_a_small_multiplier_beside_a_large_gradient`: about −0.954,
|
||||
// checked to 0.1. As there, the report is asked for at the optimum (1, 1, 0).
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("A3", "0"), ("B1", "=1000000000*A1+A2-A3"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let model = json!({"objective": "B1", "variables": "A1:A3", "constraints": ["C1 <= 1", "C2 <= 1"]});
|
||||
let r = s.execute("data.solverReport", json!({"reports": ["sensitivity"], "model": model, "values": {"A1": 1, "A2": 1, "A3": 0}})).unwrap();
|
||||
let rep = r["sheets"][0].as_str().unwrap().to_string();
|
||||
let row = |c: &str| report_row(&s, &rep, c, &["D", "E"]);
|
||||
assert!(close(&row("$A$1"), &[1.0, 0.0], 1e-3), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[1.0, 0.0], 1e-3), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$A$3"), &[0.0, -1.0], 0.1), "{:?}", row("$A$3"));
|
||||
assert!(close(&row("$C$1"), &[1.0, 1e9], 1e-3), "{:?}", row("$C$1"));
|
||||
assert!(close(&row("$C$2"), &[1.0, 1.0], 0.1), "{:?}", row("$C$2"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_refuses_a_point_that_is_not_optimal_in_a_small_variable() {
|
||||
// max 1E9·x + y with x ≤ 1 and y ≤ 1 as rows, at (1, 0.5): the second row doesn't bind, y
|
||||
// is free and its reduced gradient is its own gradient, about 1 (0.954 with the rounding
|
||||
// in `nonlinear_sensitivity_keeps_a_small_multiplier_beside_a_large_gradient`), so raising
|
||||
// y to 1 gains 0.5 and the point is not optimal: the report is refused and the workbook
|
||||
// left as it was. A stationarity tolerance taken from the largest gradient component
|
||||
// (1E-3 · 1E9 = 1E6) accepted the point and wrote the report. The optimum (1, 1) is
|
||||
// reported.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=1000000000*A1+A2"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let model = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["C1 <= 1", "C2 <= 1"]});
|
||||
let report = |values: serde_json::Value| json!({"reports": ["sensitivity"], "model": model, "values": values});
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", report(json!({"A1": 1, "A2": 0.5}))).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "{err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
let r = s.execute("data.solverReport", report(json!({"A1": 1, "A2": 1}))).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_refuses_a_point_that_breaks_a_small_row_beside_a_large_one() {
|
||||
// max y − (x − 2)² with x ≤ 1E9 and y ≤ 1 as rows: (2, 1). At (2, 1.0001) the second row
|
||||
// is over by 1E-4, more than ten times the constraint precision of its own right side
|
||||
// (1E-5 · (1 + 1) = 2E-5), so the point is infeasible and the report is refused; it is
|
||||
// stationary (x's gradient is 0 and the row, within the binding tolerance, explains y's),
|
||||
// so feasibility alone decides. A feasibility tolerance relative to the largest constant
|
||||
// right side (1E-5 · 1E9 = 1E4) accepted the point. The optimum (2, 1) is reported.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=A2-(A1-2)^2"), ("C1", "=A1"), ("C2", "=A2")]);
|
||||
let model = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["C1 <= 1000000000", "C2 <= 1"]});
|
||||
let report = |values: serde_json::Value| json!({"reports": ["sensitivity"], "model": model, "values": values});
|
||||
let before = s.doc().unwrap().clone();
|
||||
let err = s.execute("data.solverReport", report(json!({"A1": 2, "A2": 1.0001}))).unwrap_err();
|
||||
assert!(err.to_string().contains("not optimal"), "{err}");
|
||||
let d = s.doc().unwrap();
|
||||
assert!(std::sync::Arc::ptr_eq(&d.wb, &before.wb));
|
||||
assert_eq!((d.undo.len(), d.redo.len()), (before.undo.len(), before.redo.len()));
|
||||
let r = s.execute("data.solverReport", report(json!({"A1": 2, "A2": 1}))).unwrap();
|
||||
assert_eq!(r["sheets"], json!(["Sensitivity Report 1"]), "{r}");
|
||||
assert!(close(&report_row(&s, "Sensitivity Report 1", "$C$2", &["D", "E"]), &[1.0, 1.0], 1e-3), "{r}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_of_a_variable_at_its_upper_bound() {
|
||||
// max 3x + 2y with x + y ≤ 4 (a row), x ≤ 3 and y ≤ 5 (bounds on the variable cells, which
|
||||
// have no row in the Constraints table): x = 3, y = 1, objective 11. The row's shadow price
|
||||
// is y's coefficient 2 (y is basic in it). x, held at its upper bound, has reduced cost
|
||||
// 3 − 1·2 = 1: one more unit of its bound adds 1. This is the convention of Frontline's
|
||||
// "Excel Solver - Interpreting the Sensitivity Report" (solver.com): a variable's dual
|
||||
// value is nonzero only at a bound, measures the objective's change per unit increase in
|
||||
// the variable, and a bound on a variable cell has no row in the Constraints section, its
|
||||
// reduced cost standing for the bound's shadow price. Its coefficient may rise without limit
|
||||
// and fall by 1 (to 2, a tie with y): for a maximum's variable at its upper bound the
|
||||
// Allowable Increase is 1E+30 and the Allowable Decrease the reduced cost, the mirror
|
||||
// image of a variable at zero (`simplex_sensitivity_gives_the_dakota_furniture_report`).
|
||||
// y's coefficient may rise by 1 (to 3, a tie with x) and fall by 2 (to 0, where the row's
|
||||
// slack ties). The right side 4 may rise by 4 (y reaches its bound 5 at 8) and fall by 1
|
||||
// (y reaches 0 at 3).
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=3*A1+2*A2"), ("C1", "=A1+A2")]);
|
||||
let p = json!({"objective": "B1", "variables": "A1:A2", "constraints": ["C1 <= 4", "A1 <= 3", "A2 <= 5"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let row = |c: &str| report_row(&s, &rep, c, &cols);
|
||||
assert!(close(&row("$A$1"), &[3.0, 1.0, 3.0, 1e30, 1.0], 1e-9), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[1.0, 0.0, 2.0, 1.0, 2.0], 1e-9), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$C$1"), &[4.0, 2.0, 4.0, 4.0, 1.0], 1e-9), "{:?}", row("$C$1"));
|
||||
assert_eq!((count_of(&s, &rep, "$A$1"), count_of(&s, &rep, "$A$2"), count_of(&s, &rep, "$C$1")), (1, 1, 1));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_of_a_minimum_with_a_greater_row_and_a_bound() {
|
||||
// min 2x + 3y with x + y ≥ 5 (a row) and x ≤ 2 (a bound): x = 2, y = 3, cost 13. The row's
|
||||
// shadow price is 3: one more unit of right side costs one more y. x, at its upper bound,
|
||||
// has reduced cost 2 − 1·3 = −1 (one more unit of its bound saves 1). Its coefficient may
|
||||
// rise by 1 (to 3, a tie with y) and fall without limit: for a minimum's variable at its
|
||||
// upper bound the Allowable Increase is the reduced cost's size and the Allowable Decrease
|
||||
// 1E+30, the mirror image of a minimum's variable at zero (the reduced cost 1 with
|
||||
// Allowable Increase 1E+30 and Allowable Decrease 1 in
|
||||
// `sensitivity_report_gives_the_wyndor_glass_shadow_prices_and_ranges`). y's coefficient
|
||||
// may rise without limit (x can take no more) and fall by 1 (to 2, a tie with x). The right
|
||||
// side 5 may rise without limit and fall by 3 (y reaches 0 at 2).
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=2*A1+3*A2"), ("C1", "=A1+A2")]);
|
||||
let p = json!({"objective": "B1", "goal": "min", "variables": "A1:A2", "constraints": ["C1 >= 5", "A1 <= 2"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let row = |c: &str| report_row(&s, &rep, c, &cols);
|
||||
assert!(close(&row("$A$1"), &[2.0, -1.0, 2.0, 1.0, 1e30], 1e-9), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[3.0, 0.0, 3.0, 1e30, 1.0], 1e-9), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$C$1"), &[5.0, 3.0, 5.0, 1e30, 3.0], 1e-9), "{:?}", row("$C$1"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simplex_sensitivity_gives_the_dakota_furniture_report() {
|
||||
// Winston, Operations Research: Applications and Algorithms, the Dakota Furniture problem:
|
||||
// max 60x₁ + 30x₂ + 20x₃ (desks, tables, chairs) with 8x₁ + 6x₂ + x₃ ≤ 48 (lumber),
|
||||
// 4x₁ + 2x₂ + 1.5x₃ ≤ 20 (finishing) and 2x₁ + 1.5x₂ + 0.5x₃ ≤ 8 (carpentry): x = (2, 0, 8),
|
||||
// revenue 280. For the basis {s₁, x₃, x₁} the rows of B⁻¹ are (1, 2, −8), (0, 2, −4) and
|
||||
// (0, −0.5, 1.5), so the basic values are B⁻¹b = (24, 8, 2) and the shadow prices
|
||||
// c_B·B⁻¹ = (0, 20, 60)·B⁻¹ = (0, 10, 10). Tables, at zero, have reduced cost
|
||||
// 30 − (6·0 + 2·10 + 1.5·10) = −5: their price must rise by 5 before they are made, so the
|
||||
// Allowable Increase is 5 and the Allowable Decrease 1E+30, as Excel's report shows a
|
||||
// maximum's variable at zero. A change Δ in a basic variable's coefficient moves each
|
||||
// nonbasic reduced cost by Δ times that variable's row of the final tableau; the tables'
|
||||
// column there is B⁻¹(6, 2, 1.5) = (−2, −2, 1.25) and the slacks' columns are B⁻¹'s. Desks
|
||||
// (row 3): −5 − 1.25Δ ≤ 0, −10 + 0.5Δ ≤ 0 and −10 − 1.5Δ ≤ 0 hold for −4 ≤ Δ ≤ 20. Chairs
|
||||
// (row 2): −5 + 2Δ ≤ 0, −10 − 2Δ ≤ 0 and −10 + 4Δ ≤ 0 hold for −5 ≤ Δ ≤ 2.5. A right side
|
||||
// moves the basic values along its column of B⁻¹: lumber (24 + Δ, 8, 2) ≥ 0 for Δ ≥ −24;
|
||||
// finishing (24 + 2Δ, 8 + 2Δ, 2 − 0.5Δ) ≥ 0 for −4 ≤ Δ ≤ 4; carpentry
|
||||
// (24 − 8Δ, 8 − 4Δ, 2 + 1.5Δ) ≥ 0 for −4/3 ≤ Δ ≤ 2.
|
||||
let mut s = s();
|
||||
set(
|
||||
&mut s,
|
||||
&[
|
||||
("A1", "0"),
|
||||
("A2", "0"),
|
||||
("A3", "0"),
|
||||
("B1", "=60*A1+30*A2+20*A3"),
|
||||
("C1", "=8*A1+6*A2+A3"),
|
||||
("C2", "=4*A1+2*A2+1.5*A3"),
|
||||
("C3", "=2*A1+1.5*A2+0.5*A3"),
|
||||
],
|
||||
);
|
||||
let p = json!({"objective": "B1", "variables": "A1:A3", "constraints": ["C1 <= 48", "C2 <= 20", "C3 <= 8"], "method": "simplex", "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
assert_eq!(r["status"], "optimal", "{r}");
|
||||
assert_eq!(r["objective"], 280.0, "{r}");
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
let cols = ["D", "E", "F", "G", "H"];
|
||||
let row = |c: &str| report_row(&s, &rep, c, &cols);
|
||||
assert!(close(&row("$A$1"), &[2.0, 0.0, 60.0, 20.0, 4.0], 1e-9), "{:?}", row("$A$1"));
|
||||
assert!(close(&row("$A$2"), &[0.0, -5.0, 30.0, 5.0, 1e30], 1e-9), "{:?}", row("$A$2"));
|
||||
assert!(close(&row("$A$3"), &[8.0, 0.0, 20.0, 2.5, 5.0], 1e-9), "{:?}", row("$A$3"));
|
||||
assert!(close(&row("$C$1"), &[24.0, 0.0, 48.0, 1e30, 24.0], 1e-9), "{:?}", row("$C$1"));
|
||||
assert!(close(&row("$C$2"), &[20.0, 10.0, 20.0, 4.0, 4.0], 1e-9), "{:?}", row("$C$2"));
|
||||
assert!(close(&row("$C$3"), &[8.0, 10.0, 8.0, 2.0, 4.0 / 3.0], 1e-9), "{:?}", row("$C$3"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonlinear_sensitivity_signs_follow_the_goal() {
|
||||
// A Lagrange multiplier is the derivative of the optimal objective with respect to the
|
||||
// row's right side, and a reduced gradient that with respect to the variable's bound,
|
||||
// whichever way the goal points. Frontline's "Excel Solver - Interpreting the Sensitivity
|
||||
// Report" (solver.com) names a binding constraint's dual value its Lagrange multiplier and a
|
||||
// variable's its reduced gradient, each the objective's change per unit increase in the
|
||||
// constraint's bound or the variable, valid at the optimal point alone for a nonlinear
|
||||
// model.
|
||||
// min (x − 3)² + (y − 1)² with x + y ≤ b: the point of the line nearest (3, 1) is
|
||||
// (2.5, 0.5) at b = 3, f*(b) = (4 − b)²/2 and df*/db = −(4 − b) = −1, the multiplier;
|
||||
// maximising the negated objective reaches the same point with multiplier +1. Both
|
||||
// variables are interior, so their reduced gradients are 0.
|
||||
let mut s = s();
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0"), ("B1", "=(A1-3)^2+(A2-1)^2"), ("B2", "=0-(A1-3)^2-(A2-1)^2"), ("C1", "=A1+A2")]);
|
||||
for (obj, goal, lambda) in [("B1", "min", -1.0), ("B2", "max", 1.0)] {
|
||||
set(&mut s, &[("A1", "0"), ("A2", "0")]);
|
||||
let p = json!({"objective": obj, "goal": goal, "variables": "A1:A2", "constraints": ["C1 <= 3"], "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$A$1", &["D", "E"]), &[2.5, 0.0], 1e-3), "{goal}: {r}");
|
||||
assert!(close(&report_row(&s, &rep, "$A$2", &["D", "E"]), &[0.5, 0.0], 1e-3), "{goal}: {r}");
|
||||
assert!(close(&report_row(&s, &rep, "$C$1", &["D", "E"]), &[3.0, lambda], 1e-3), "{goal}: {r}");
|
||||
}
|
||||
// max −(x + 1)² − (y − 2)² with x + y ≤ 1 and x, y ≥ 0 mirrors the minimum in
|
||||
// `sensitivity_report_of_a_nonlinear_model_gives_lagrange_multipliers`: at (0, 1) the
|
||||
// multiplier is +2 (−2 there) and x's reduced gradient −4 (4 there), since raising x from
|
||||
// its bound lowers the maximum by 2(x + 1) + 2 = 4 per unit.
|
||||
set(&mut s, &[("D1", "0.5"), ("D2", "0.5"), ("E1", "=0-(D1+1)^2-(D2-2)^2"), ("E2", "=D1+D2")]);
|
||||
let p = json!({"objective": "E1", "goal": "max", "variables": "D1:D2", "constraints": ["E2 <= 1"], "reports": ["sensitivity"]});
|
||||
let r = s.execute("data.solver", p).unwrap();
|
||||
let rep = r["reports"][0].as_str().unwrap().to_string();
|
||||
assert!(close(&report_row(&s, &rep, "$D$1", &["D", "E"]), &[0.0, -4.0], 1e-3), "{r}");
|
||||
assert!(close(&report_row(&s, &rep, "$D$2", &["D", "E"]), &[1.0, 0.0], 1e-3), "{r}");
|
||||
assert!(close(&report_row(&s, &rep, "$E$2", &["D", "E"]), &[1.0, 2.0], 1e-3), "{r}");
|
||||
}
|
||||
|
||||
fn formula(s: &Session, a1: &str) -> String {
|
||||
let d = s.doc().unwrap();
|
||||
d.wb.active().unwrap().cell(CellRef::parse(a1).unwrap()).and_then(|c| c.formula.as_ref()).map(|f| f.text.clone()).unwrap_or_default()
|
||||
|
||||
Reference in New Issue
Block a user