DARE convergence test is absolute-only, and `Lqr` still aborts when the iteration budget is short

[gabrielsantosphilips]

Summary

DiscreteAlgebraicRiccatiEquation defaults to MaxIterations = 300. Lqr instantiates it with that default and converts a non-converged result into really_assert(r.converged) — an abort. IntegralStateFeedbackLqi builds an Lqr internally, so it inherits the same failure.

The solver documentation already anticipates this case:

Slow convergence. Systems with eigenvalues near the unit circle converge slowly. Increase the maximum iteration bound if needed.

There is no way to act on that advice through Lqr or IntegralStateFeedbackLqi — MaxIterations is not forwarded, and the failure surfaces as a hard abort rather than a recoverable error.

Impact

Finely sampled plants have discrete poles near 1, which is the normal situation for an embedded control loop. A first-order mechanical plant with an integral augmentation, sampled at 1 kHz, needs 9480 sweeps. At the default 300 the solver returns converged = false with P off by a factor of 10, and the Lqr constructor aborts on target.

doc/solvers/DiscreteAlgebraicRiccatiEquation.md states iteration counts are "typically 10–30", which understates this class of problem by more than two orders of magnitude.

Reproduction

#include "numerical/solvers/DiscreteAlgebraicRiccatiEquation.hpp"
#include <cstdio>

// Speed plant w[k+1] = 0.999 w[k] + 0.5 u[k], augmented with an integral state
// scaled by Ts = 1 ms, as IntegralStateFeedbackLqi builds it.
int main()
{
    const math::SquareMatrix<float, 2> a{ { 0.999f, 0.0f }, { -0.001f, 1.0f } };
    const math::Matrix<float, 2, 1>    b{ { 0.5f }, { 0.0f } };
    const math::SquareMatrix<float, 2> q{ { 1.0f, 0.0f }, { 0.0f, 0.1f } };
    const math::SquareMatrix<float, 1> r{ 0.01f };

    solvers::DiscreteAlgebraicRiccatiEquation<float, 2, 1, 300>   atCap;
    solvers::DiscreteAlgebraicRiccatiEquation<float, 2, 1, 30000> raised;

    auto low  = atCap.Solve(a, b, q, r);
    auto high = raised.Solve(a, b, q, r);

    std::printf("cap=300   converged=%d  P22=%.4g\n", low.converged,  low.value.at(1, 1));
    std::printf("cap=30000 converged=%d  P22=%.4g\n", high.converged, high.value.at(1, 1));
}

Output:

cap=300   converged=0  P22=30.01
cap=30000 converged=1  P22=314.8

Residual trace of the same iteration:

sweeps	max abs(dP)	P22
300	9.9e-2	30.0
1000	9.1e-2	96.9
5000	1.6e-2	290.7
20000	3.0e-8	316.4

Sweeps to reach the absolute tolerance 1e-3: 9480. Sweeps to reach a relative tolerance of 1e-6: 11271.

Going through Lqr instead aborts:

Assertion `r.converged' failed.
numerical/controllers/implementations/Lqr.hpp:54

Root cause

Value iteration on the DARE converges linearly at a rate set by the closed-loop spectral radius. As Ts shrinks, that radius approaches 1 and the sweep count grows roughly as 1/Ts. 300 is only adequate for coarsely sampled or aggressively weighted problems.

Secondary issue — absolute convergence threshold

if (math::Abs(pij - math::ToFloat(Pprev.at(i, j))) > tolerance)   // tolerance = math::Tolerance<T>() = 1e-3f

The test is absolute, so the effective stopping precision depends on the magnitude of P, which depends on the units the user happens to work in. The same design applies equally to a P of order 1e-2 and one of order 1e3. This is the same class of problem as #259 (Cholesky absolute pivot threshold).

Suggested direction

1. Forward MaxIterations from Lqr and IntegralStateFeedbackLqi, or raise the default to something realistic for Ts-scaled plants.
2. Replace the absolute test with a mixed absolute/relative criterion, e.g. abs(dP) <= tol * (1 + abs(P)).
3. Give Lqr a non-aborting path — a factory returning std::optional, or an explicit converged() accessor — so a caller can fall back instead of taking down the device. really_assert in a member-initializer lambda is a hard abort in an embedded runtime.

Found while building an LQI speed controller on top of IntegralStateFeedbackLqi<float, 1, 1, 1>. Worked around by expressing the design in per-unit input with unit-sample integration, which moves the closed-loop poles away from 1 and converges in well under 300 sweeps — but that is a change of problem formulation, not a fix.

[gabrielsantosphilips]

Re-verified against v3.1.0. The MaxIterations forwarding from #279 works — IntegralStateFeedbackLqi<float, 1, 1, 1, 30000> constructs the Ts-scaled design that used to abort here.

Points 2 (absolute convergence threshold) and 3 (really_assert on non-convergence) are still open, and while re-checking I found the tolerance issue is worse than described above: with small Q/R the solver reports converged = true on a P that has barely moved from Q, yielding gains that are off by 15 % / 67 % with no error signalled. Carried over to #288 with a reproduction, since this issue is closed by a merged PR and cannot be reopened.