Diagnose full-duration WaveStar NMPC numerical instability

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`

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
+  }
+}

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
+      }
+    }
+  ]
+}

tests/test_wavestar_nmpc_stability_diagnostic.py

@@ -0,0 +1,49 @@
+"""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
+
+
+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)
+    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():
+    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.)