From 11a9f1a91cf07b40fc206bc5b48468b07f13e8eb Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 10:04:38 -0400 Subject: [PATCH 1/2] Document full-duration WaveStar NMPC numerical stability blocker --- .../README.md | 43 +++++ .../autonomous_step_sensitivity_mac_arm.json | 121 ++++++++++++++ .../qp_optimality_mac_arm.json | 114 +++++++++++++ ...test_wavestar_nmpc_stability_diagnostic.py | 39 +++++ ...vestar_nmpc_autonomous_step_sensitivity.py | 156 ++++++++++++++++++ tools/wavestar_nmpc_qp_optimality.py | 155 +++++++++++++++++ 6 files changed, 628 insertions(+) create mode 100644 tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/autonomous_step_sensitivity_mac_arm.json create mode 100644 tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/qp_optimality_mac_arm.json create mode 100644 tests/test_wavestar_nmpc_stability_diagnostic.py create mode 100644 tools/wavestar_nmpc_autonomous_step_sensitivity.py create mode 100644 tools/wavestar_nmpc_qp_optimality.py diff --git a/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/README.md b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/README.md index 90dfacef..9e4e6f3a 100644 --- a/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/README.md +++ b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/README.md @@ -319,6 +319,49 @@ whereas Python uses L-BFGS-B. This bound slack explains that sample's This result identifies the remaining Python controller numerical defect more precisely, but the autonomous full-duration physical gate still fails. +## Full-duration stability blocker, 2026-10-24 + +The checked [225 s source and Python step diagnostic](autonomous_step_sensitivity_mac_arm.json) +recomputes both source self-change and independent physical Python step +sensitivity against the **unchanged** `0.0005` rad pitch, `0.005` rad/s pitch +speed, `0.0001` m stroke, `1` N axial-force, and `0.1` N m request limits. +The two pinned Apple silicon MATLAB traces at 0.125 and 0.0625 ms have +identical 500-bin wave components. They differ by `0.010800` rad pitch and +`0.721063` N m request, first crossing the pitch limit at 17.95 s. The +Python physical runs at 1 and 0.5 ms each use only the pinned HDF5, wave +phases, and their own controller/plant states in feedback. They differ by +`0.005127` rad pitch, `0.024146` rad/s pitch speed, `1.022` mm stroke, +`2.014` N axial force, and `0.414344` N m request; pitch first crosses its +limit at 43.5 s. Both step comparisons fail all five original pointwise +limits. The run used macOS ARM64, NumPy 2.5.4, and SciPy 1.18.1. Reproduce +it with `python tools/wavestar_nmpc_autonomous_step_sensitivity.py --output +/tmp/wavestar-step-sensitivity.json` from the Python repository after +installing its dependencies. The established 225 s aggregate and strict +expected-failure gates are unchanged. + +The checked [fixed-input QP optimality diagnostic](qp_optimality_mac_arm.json) +uses the same logged state, forecast, and preceding MATLAB plan at 17.5, +30, 97.25, and 202.25 s. The Hessian minimum eigenvalues are `0.886–0.899` +and condition numbers are `1.33–1.41`. Python's bound violations are zero; +at the three non-bound examples, its projected KKT residual is at most +`7.1e-9` and its first request is within `3.1e-9` N m of the direct source +QP. At 97.25 s, MATLAB's interior-point result is `1.76e-7` N m above the +lower torque bound while Python returns that bound exactly. Resetting the +preceding horizon plan at 30 s changes the first request from `5.92177` to +`0.681924` N m on the same state and forecast: it changes the actual RTI +linearization and QP, so it cannot be treated as a tighter solve of the +published controller. Reproduce this bounded diagnostic with `python +tools/wavestar_nmpc_qp_optimality.py --output /tmp/wavestar-qp.json`. + +The remaining correction must establish a 225 s MATLAB controller/plant +reference whose own numerical step and supported-platform variation satisfy +the original pointwise limits, then demonstrate independent physical Python +motion, request, force, and work against that reference on supported runners. +The checked QP residuals do not support an ill-conditioned one-step QP as the +sole cause. No such correction or passing autonomous gate has been +established; issue [#623](https://github.com/cmudrc/wec-sim-python/issues/623) +and the qualified case status remain open. + SHA-256: - `source_1ms_repeat.npz`: `42f3264022177833634731711086cc90a1a39eb313ccdc05110123eb9c3035f6` diff --git a/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/autonomous_step_sensitivity_mac_arm.json b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/autonomous_step_sensitivity_mac_arm.json new file mode 100644 index 00000000..5b004b04 --- /dev/null +++ b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/autonomous_step_sensitivity_mac_arm.json @@ -0,0 +1,121 @@ +{ + "purpose": "negative full-duration autonomous pointwise diagnostic", + "duration_s": 225.0, + "controller_step_s": 0.05, + "python_feedback_inputs": "pinned HDF5, wave phases, and own physical states only", + "original_pointwise_limits": { + "pitch_rad": 0.0005, + "pitch_speed_rad_s": 0.005, + "stroke_m": 0.0001, + "axial_force_N": 1.0, + "requested_torque_Nm": 0.1 + }, + "pinned_sha256": { + "source_0p125ms_macos_arm_m1.npz": "ee791a11bb3c199e510dcc33e9f1b718313a34473694258ce5f75fd44ae55df0", + "source_0p0625ms_macos_arm_m1.npz": "b3fe52fe66256a218c924c1c039f60c4b3d9d2192f36b1d9d4efe56831fac26b", + "wavestar.h5": "b2541ff01c319987d9157d1be9fbc9ec39a7660d0572d8e316c428b8aa3359d1" + }, + "platform": "macOS-26.5.1-arm64-arm-64bit-Mach-O", + "numpy": "2.5.4", + "scipy": "1.18.1", + "matlab_0p125ms_vs_0p0625ms": { + "pointwise": { + "pitch_rad": { + "max_abs": 0.010800047009331015, + "first_limit_crossing_s": 17.95 + }, + "pitch_speed_rad_s": { + "max_abs": 0.057609864539682004, + "first_limit_crossing_s": 21.5 + }, + "stroke_m": { + "max_abs": 0.0021581679653852014, + "first_limit_crossing_s": 17.95 + }, + "axial_force_N": { + "max_abs": 3.4746706199945097, + "first_limit_crossing_s": 26.700000000000003 + }, + "requested_torque_Nm": { + "max_abs": 0.7210627885203973, + "first_limit_crossing_s": 21.6 + } + }, + "within_original_pointwise_limits": false + }, + "python_1ms_vs_0p5ms": { + "pointwise": { + "pitch_rad": { + "max_abs": 0.005126777656919185, + "first_limit_crossing_s": 43.5 + }, + "pitch_speed_rad_s": { + "max_abs": 0.02414629731403123, + "first_limit_crossing_s": 43.800000000000004 + }, + "stroke_m": { + "max_abs": 0.0010220081477542164, + "first_limit_crossing_s": 43.5 + }, + "axial_force_N": { + "max_abs": 2.014343434375106, + "first_limit_crossing_s": 44.7 + }, + "requested_torque_Nm": { + "max_abs": 0.41434360671301085, + "first_limit_crossing_s": 41.85 + } + }, + "within_original_pointwise_limits": false + }, + "python_1ms_vs_matlab_0p125ms": { + "pointwise": { + "pitch_rad": { + "max_abs": 0.007225585232541097, + "first_limit_crossing_s": 17.95 + }, + "pitch_speed_rad_s": { + "max_abs": 0.049707098319223625, + "first_limit_crossing_s": 18.1 + }, + "stroke_m": { + "max_abs": 0.001443955853264338, + "first_limit_crossing_s": 17.95 + }, + "axial_force_N": { + "max_abs": 3.055945666091551, + "first_limit_crossing_s": 26.650000000000002 + }, + "requested_torque_Nm": { + "max_abs": 0.6205056457635307, + "first_limit_crossing_s": 17.35 + } + }, + "within_original_pointwise_limits": false + }, + "python_0p5ms_vs_matlab_0p125ms": { + "pointwise": { + "pitch_rad": { + "max_abs": 0.006288137147178071, + "first_limit_crossing_s": 17.95 + }, + "pitch_speed_rad_s": { + "max_abs": 0.046827784730315083, + "first_limit_crossing_s": 18.1 + }, + "stroke_m": { + "max_abs": 0.0012549441223213484, + "first_limit_crossing_s": 17.95 + }, + "axial_force_N": { + "max_abs": 2.62931595366868, + "first_limit_crossing_s": 26.650000000000002 + }, + "requested_torque_Nm": { + "max_abs": 0.5360698159169655, + "first_limit_crossing_s": 17.35 + } + }, + "within_original_pointwise_limits": false + } +} diff --git a/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/qp_optimality_mac_arm.json b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/qp_optimality_mac_arm.json new file mode 100644 index 00000000..1d6b806f --- /dev/null +++ b/tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s/qp_optimality_mac_arm.json @@ -0,0 +1,114 @@ +{ + "purpose": "fixed-input RTI QP optimality and plan-policy diagnostic", + "pinned_sha256": { + "controller_0p25ms_compatible_amd7763.npz": "06d55c1972b6e9130fbeeea52d29cc6691fddfb7408c687c950ce9522c2498e3", + "checked_qp_replay_compatible_amd7763.mat": "cbc5a2acfe651f3570db0405d040ee3d40791b07f30ac4b3d282ae289ec5af7d", + "direct_controller_inputs.npz": "3ff8f4ba8654ef66bf9058a9fe717cfa6c60868300d19e64a089fcb06781dbf0" + }, + "platform": "macOS-26.5.1-arm64-arm-64bit-Mach-O", + "numpy": "2.5.4", + "scipy": "1.18.1", + "active_bound_tolerance_Nm": 1e-09, + "samples": [ + { + "time_s": 17.5, + "prior_plan": "source", + "hessian_min_eigenvalue": 0.8875814894975004, + "hessian_condition_number": 1.412452927392197, + "python_first_request_Nm": 2.2383394523585394, + "matlab_first_request_Nm": 2.238339450716518, + "python_solution": { + "qp_objective": -0.1341191825630914, + "max_projected_kkt_residual": 4.086149217596979e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + }, + "matlab_plan_in_same_qp": { + "qp_objective": -0.13411918256309152, + "max_projected_kkt_residual": 2.48131191635248e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + } + }, + { + "time_s": 30.0, + "prior_plan": "source", + "hessian_min_eigenvalue": 0.8939600873488631, + "hessian_condition_number": 1.3719119443603895, + "python_first_request_Nm": 5.921774300548368, + "matlab_first_request_Nm": 5.921774297497542, + "python_solution": { + "qp_objective": -0.5208669352970017, + "max_projected_kkt_residual": 7.089540218774282e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + }, + "matlab_plan_in_same_qp": { + "qp_objective": -0.5208669352970017, + "max_projected_kkt_residual": 1.0465008359883399e-08, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + } + }, + { + "time_s": 97.25, + "prior_plan": "source", + "hessian_min_eigenvalue": 0.8859142039982983, + "hessian_condition_number": 1.3761506359866307, + "python_first_request_Nm": -12.0, + "matlab_first_request_Nm": -11.999999823862918, + "python_solution": { + "qp_objective": -0.6008552906995042, + "max_projected_kkt_residual": 1.9168391179880118e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 39 + }, + "matlab_plan_in_same_qp": { + "qp_objective": -0.6008552823088922, + "max_projected_kkt_residual": 0.0476369098919201, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + } + }, + { + "time_s": 202.25, + "prior_plan": "source", + "hessian_min_eigenvalue": 0.8988196216308708, + "hessian_condition_number": 1.3308009988604879, + "python_first_request_Nm": 2.792998566072698, + "matlab_first_request_Nm": 2.7929985652597757, + "python_solution": { + "qp_objective": -0.26546019648751623, + "max_projected_kkt_residual": 3.980322182961515e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + }, + "matlab_plan_in_same_qp": { + "qp_objective": -0.2654601964875163, + "max_projected_kkt_residual": 3.848165786024538e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + } + }, + { + "time_s": 30.0, + "prior_plan": "zero", + "hessian_min_eigenvalue": 0.9003463162491234, + "hessian_condition_number": 1.4030200946517661, + "python_first_request_Nm": 0.681924023796555, + "matlab_first_request_Nm": 5.921774297497542, + "python_solution": { + "qp_objective": -2.7823668671492157, + "max_projected_kkt_residual": 4.230424566209479e-09, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + }, + "matlab_plan_in_same_qp": { + "qp_objective": 112.27721456196204, + "max_projected_kkt_residual": 5.23770035051823, + "max_bound_violation_Nm": 0.0, + "free_variables_at_1e-9_Nm_tolerance": 40 + } + } + ] +} diff --git a/tests/test_wavestar_nmpc_stability_diagnostic.py b/tests/test_wavestar_nmpc_stability_diagnostic.py new file mode 100644 index 00000000..8256edb3 --- /dev/null +++ b/tests/test_wavestar_nmpc_stability_diagnostic.py @@ -0,0 +1,39 @@ +"""Preserve the pinned evidence for WaveStar's open pointwise gate.""" + +import numpy as np +import pytest + +from tools.wavestar_nmpc_autonomous_step_sensitivity import ( + LIMITS, _compare, _load_source, +) +from tools.wavestar_nmpc_qp_optimality import audit as qp_audit + + +@pytest.mark.xfail( + strict=True, + reason="published RTI closed loop changes beyond the original pointwise gate", +) +def test_pinned_matlab_source_step_stability_under_original_limits(): + half = _load_source("source_0p125ms_macos_arm_m1.npz") + quarter = _load_source("source_0p0625ms_macos_arm_m1.npz") + np.testing.assert_array_equal(half["components"], quarter["components"]) + comparison = _compare(half, quarter) + for label, limit in LIMITS.items(): + assert comparison["pointwise"][label]["max_abs"] < limit, label + + +def test_pinned_qp_is_well_conditioned_but_plan_policy_changes_objective(): + report = qp_audit() + source_prior = [sample for sample in report["samples"] + if sample["prior_plan"] == "source"] + assert len(source_prior) == 4 + for sample in source_prior: + assert sample["hessian_min_eigenvalue"] > .8 + assert sample["hessian_condition_number"] < 2 + assert sample["python_solution"]["max_bound_violation_Nm"] == 0 + assert (abs(sample["python_first_request_Nm"] - + sample["matlab_first_request_Nm"]) < 2e-7) + reset = report["samples"][-1] + assert reset["time_s"] == 30. and reset["prior_plan"] == "zero" + assert (abs(reset["python_first_request_Nm"] - + reset["matlab_first_request_Nm"]) > 1.) diff --git a/tools/wavestar_nmpc_autonomous_step_sensitivity.py b/tools/wavestar_nmpc_autonomous_step_sensitivity.py new file mode 100644 index 00000000..3a686f19 --- /dev/null +++ b/tools/wavestar_nmpc_autonomous_step_sensitivity.py @@ -0,0 +1,156 @@ +"""Audit 225 s WaveStar closed-loop sensitivity without source-driven feedback. + +This is a negative diagnostic. It measures two pinned MATLAB source plant +steps and two independent Python physical plant steps against the unchanged +pointwise limits. It never supplies MATLAB states, strokes, forces, or torque +requests to the Python controller. +""" + +import argparse +import json +import platform +import sys +from hashlib import sha256 +from pathlib import Path + +import numpy as np +import scipy + +from wecsim.dynamics.wavestar import run_wavestar_published +from wecsim.environment.irregularWave import jonswap_equal_energy_components + + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT / "tests")) +from test_wavestar_split_step_portable import _DelayedSampledPTO # noqa: E402 + + +SOURCE = ROOT / "tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s" +HYDRO = ROOT / "tests/fixtures/matlab/wavestar_nmpc_split_step_225s/wavestar.h5" +FILES = { + "source_0p125ms_macos_arm_m1.npz": + "ee791a11bb3c199e510dcc33e9f1b718313a34473694258ce5f75fd44ae55df0", + "source_0p0625ms_macos_arm_m1.npz": + "b3fe52fe66256a218c924c1c039f60c4b3d9d2192f36b1d9d4efe56831fac26b", + "wavestar.h5": + "b2541ff01c319987d9157d1be9fbc9ec39a7660d0572d8e316c428b8aa3359d1", +} +LIMITS = { + "pitch_rad": .0005, + "pitch_speed_rad_s": .005, + "stroke_m": .0001, + "axial_force_N": 1., + "requested_torque_Nm": .1, +} +CHANNELS = { + "pitch_rad": "pitch", + "pitch_speed_rad_s": "pitch_speed", + "stroke_m": "stroke", + "axial_force_N": "axial_force", + "requested_torque_Nm": "command", +} + + +def _checked(path, name): + actual = sha256(path.read_bytes()).hexdigest() + if actual != FILES[name]: + raise ValueError(f"Pinned WaveStar input {path} changed: {actual}") + return actual + + +def _load_source(name): + path = SOURCE / name + _checked(path, name) + with np.load(path, allow_pickle=False) as saved: + trace = {key: saved[key].copy() for key in saved.files} + if trace["components"].shape != (500, 4): + raise ValueError("The source wave components changed shape") + time = trace["time"] + np.testing.assert_allclose(time, np.arange(4501) * .05, + rtol=0, atol=1e-9) + for key in (*CHANNELS.values(), "pto_speed"): + if trace[key].shape != time.shape or not np.isfinite(trace[key]).all(): + raise ValueError(f"The source {key} trace is invalid") + return trace + + +def _compare(left, right): + time = left["time"] + np.testing.assert_array_equal(time, right["time"]) + pointwise = {} + for label, key in CHANNELS.items(): + difference = np.abs(left[key] - right[key]) + first = np.flatnonzero(difference >= LIMITS[label]) + pointwise[label] = { + "max_abs": float(np.max(difference)), + "first_limit_crossing_s": float(time[first[0]]) if len(first) else None, + } + return { + "pointwise": pointwise, + "within_original_pointwise_limits": all( + pointwise[label]["max_abs"] < limit for label, limit in LIMITS.items() + ), + } + + +def _independent_python(hydro, sea, plant_step): + controller = _DelayedSampledPTO(plant_dt=plant_step) + response = run_wavestar_published( + hydro, sea, dt=plant_step, end_time=225., ramp_time=25., + g=9.80665, pto_controller=controller, + output_stride=round(.05 / plant_step), + ) + trace = { + "time": response.time.copy(), + "pitch": response.angle.copy(), + "pitch_speed": response.angular_speed.copy(), + "stroke": response.pto_stroke.copy(), + "axial_force": response.pto_force.copy(), + "command": np.asarray(controller.commands), + } + if any(value.shape != (4501,) or not np.isfinite(value).all() + for value in trace.values()): + raise ValueError("The independent Python trace is incomplete") + return trace + + +def audit(): + _checked(HYDRO, "wavestar.h5") + half = _load_source("source_0p125ms_macos_arm_m1.npz") + quarter = _load_source("source_0p0625ms_macos_arm_m1.npz") + np.testing.assert_array_equal(half["components"], quarter["components"]) + sea = jonswap_equal_energy_components( + HYDRO, significant_height=.1042, peak_period=1.836, + directions=np.array([0.]), spreading=np.array([1.]), gamma=3.3, + phase=half["components"][:, 3, None], + ) + python_1ms = _independent_python(HYDRO, sea, .001) + python_0p5ms = _independent_python(HYDRO, sea, .0005) + return { + "purpose": "negative full-duration autonomous pointwise diagnostic", + "duration_s": 225., + "controller_step_s": .05, + "python_feedback_inputs": "pinned HDF5, wave phases, and own physical states only", + "original_pointwise_limits": LIMITS, + "pinned_sha256": FILES, + "platform": platform.platform(), + "numpy": np.__version__, + "scipy": scipy.__version__, + "matlab_0p125ms_vs_0p0625ms": _compare(half, quarter), + "python_1ms_vs_0p5ms": _compare(python_1ms, python_0p5ms), + "python_1ms_vs_matlab_0p125ms": _compare(python_1ms, half), + "python_0p5ms_vs_matlab_0p125ms": _compare(python_0p5ms, half), + } + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, required=True) + args = parser.parse_args() + result = audit() + args.output.write_text(json.dumps(result, indent=2) + "\n") + print(json.dumps(result, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tools/wavestar_nmpc_qp_optimality.py b/tools/wavestar_nmpc_qp_optimality.py new file mode 100644 index 00000000..af85685f --- /dev/null +++ b/tools/wavestar_nmpc_qp_optimality.py @@ -0,0 +1,155 @@ +"""Check the pinned WaveStar RTI QP at selected sensitive source samples. + +The source state, forecast, and preceding plan are fixed inputs to this +component diagnostic. No result here constitutes autonomous motion parity. +""" + +import argparse +import json +import platform +from hashlib import sha256 +from pathlib import Path + +import numpy as np +import scipy +from scipy.io import loadmat + +import wecsim.control.wavestarNmpc as nmpc + + +ROOT = Path(__file__).resolve().parents[1] +REF = ROOT / "tests/fixtures/matlab/wavestar_nmpc_reference_convergence_225s" +MATRICES = ROOT / "tests/fixtures/matlab/wavestar_nmpc_full_step/direct_controller_inputs.npz" +HASHES = { + "controller_0p25ms_compatible_amd7763.npz": + "06d55c1972b6e9130fbeeea52d29cc6691fddfb7408c687c950ce9522c2498e3", + "checked_qp_replay_compatible_amd7763.mat": + "cbc5a2acfe651f3570db0405d040ee3d40791b07f30ac4b3d282ae289ec5af7d", + "direct_controller_inputs.npz": + "3ff8f4ba8654ef66bf9058a9fe717cfa6c60868300d19e64a089fcb06781dbf0", +} + + +def _checked(path): + digest = sha256(path.read_bytes()).hexdigest() + if digest != HASHES[path.name]: + raise ValueError(f"Pinned QP input changed: {path}: {digest}") + + +def _qp_metrics(hessian, linear, bounds, delta): + gradient = hessian @ delta + linear + lower, upper = bounds[:, 0], bounds[:, 1] + at_lower = delta <= lower + 1e-9 + at_upper = delta >= upper - 1e-9 + free = ~(at_lower | at_upper) + residual = np.zeros_like(delta) + residual[at_lower] = np.minimum(gradient[at_lower], 0.) + residual[at_upper] = np.maximum(gradient[at_upper], 0.) + residual[free] = gradient[free] + return { + "qp_objective": float(.5 * delta @ hessian @ delta + linear @ delta), + "max_projected_kkt_residual": float(np.max(np.abs(residual))), + "max_bound_violation_Nm": float(max( + 0., np.max(lower - delta), np.max(delta - upper), + )), + "free_variables_at_1e-9_Nm_tolerance": int(np.count_nonzero(free)), + } + + +def _sample(time_s, source_prior, states, forecasts, plans, matrices): + index = round(time_s / .05) + active = index - 300 + controller = nmpc.WaveStarNmpcController() + controller._transition = matrices["A"].copy() + controller._input = matrices["Bu"][:, 0].copy() + controller._wave_input = matrices["Bw"][:, 0].copy() + if source_prior and active > 0: + controller._nominal = plans[active - 1].copy() + shifted = np.r_[controller._nominal[1:], controller._nominal[-1]] + captured = {} + original = nmpc.minimize + + def capture(objective, x0, *, jac, bounds, method, options): + if method == "L-BFGS-B": + zero = np.zeros(40) + linear = jac(zero) + basis = np.eye(40) + hessian = np.column_stack([ + jac(basis[i]) - linear for i in range(40) + ]) + captured.update(hessian=hessian, linear=linear, + bounds=np.asarray(bounds)) + return original(objective, x0, jac=jac, bounds=bounds, + method=method, options=options) + + nmpc.minimize = capture + try: + controller.step(states[index], forecasts[index], time_s) + finally: + nmpc.minimize = original + + hessian = captured["hessian"] + linear = captured["linear"] + bounds = captured["bounds"] + python_delta = controller._nominal - shifted + source_delta = plans[active] - shifted + return { + "time_s": time_s, + "prior_plan": "source" if source_prior else "zero", + "hessian_min_eigenvalue": float(np.linalg.eigvalsh(hessian).min()), + "hessian_condition_number": float(np.linalg.cond(hessian)), + "python_first_request_Nm": float(controller._nominal[0]), + "matlab_first_request_Nm": float(plans[active, 0]), + "python_solution": _qp_metrics(hessian, linear, bounds, python_delta), + "matlab_plan_in_same_qp": _qp_metrics( + hessian, linear, bounds, source_delta, + ), + } + + +def audit(): + controller_file = REF / "controller_0p25ms_compatible_amd7763.npz" + replay_file = REF / "checked_qp_replay_compatible_amd7763.mat" + for path in (controller_file, replay_file, MATRICES): + _checked(path) + with np.load(controller_file, allow_pickle=False) as saved: + time = saved["time"].copy() + states = saved["estimated_state"].copy() + forecasts = saved["forecast"].copy() + with np.load(MATRICES, allow_pickle=False) as saved: + matrices = {key: saved[key].copy() for key in ("A", "Bu", "Bw")} + replay = loadmat(replay_file) + plans = replay["plans"] + np.testing.assert_allclose(time, np.arange(4501) * .05, + rtol=0, atol=1e-9) + np.testing.assert_array_equal(replay["time"].ravel(), time[300:]) + if (states.shape != (4501, 5) or forecasts.shape != (4501, 40) or + plans.shape != (4201, 40)): + raise ValueError("Pinned QP arrays changed shape") + cases = [ + _sample(t, True, states, forecasts, plans, matrices) + for t in (17.5, 30., 97.25, 202.25) + ] + cases.append(_sample(30., False, states, forecasts, plans, matrices)) + return { + "purpose": "fixed-input RTI QP optimality and plan-policy diagnostic", + "pinned_sha256": HASHES, + "platform": platform.platform(), + "numpy": np.__version__, + "scipy": scipy.__version__, + "active_bound_tolerance_Nm": 1e-9, + "samples": cases, + } + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, required=True) + args = parser.parse_args() + report = audit() + args.output.write_text(json.dumps(report, indent=2) + "\n") + print(json.dumps(report, indent=2)) + + +if __name__ == "__main__": + main() From 430ba77f2dfae1d4868f77cf2ecd941b63c8e784 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 10:05:50 -0400 Subject: [PATCH 2/2] Limit WaveStar expected failure to checked step divergence --- tests/test_wavestar_nmpc_stability_diagnostic.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/tests/test_wavestar_nmpc_stability_diagnostic.py b/tests/test_wavestar_nmpc_stability_diagnostic.py index 8256edb3..55637dc6 100644 --- a/tests/test_wavestar_nmpc_stability_diagnostic.py +++ b/tests/test_wavestar_nmpc_stability_diagnostic.py @@ -9,17 +9,27 @@ from tools.wavestar_nmpc_qp_optimality import audit as qp_audit +class PublishedClosedLoopStepGap(AssertionError): + """The checked source pair exceeds an unchanged pointwise limit.""" + + @pytest.mark.xfail( strict=True, reason="published RTI closed loop changes beyond the original pointwise gate", + raises=PublishedClosedLoopStepGap, ) def test_pinned_matlab_source_step_stability_under_original_limits(): half = _load_source("source_0p125ms_macos_arm_m1.npz") quarter = _load_source("source_0p0625ms_macos_arm_m1.npz") np.testing.assert_array_equal(half["components"], quarter["components"]) comparison = _compare(half, quarter) - for label, limit in LIMITS.items(): - assert comparison["pointwise"][label]["max_abs"] < limit, label + failed = { + label: comparison["pointwise"][label]["max_abs"] + for label, limit in LIMITS.items() + if comparison["pointwise"][label]["max_abs"] >= limit + } + if failed: + raise PublishedClosedLoopStepGap(str(failed)) def test_pinned_qp_is_well_conditioned_but_plan_policy_changes_objective():