Fix matched c2d/sample returning NaN for a pole or zero at the origin - #1241

coverage: 94.743% (+0.003%) from 94.74% — gaoflow:fix-matched-c2d-origin-pole-zero into python-control:main

Thanks for the PR!

This function has gotten almost no attention and has really atrophied, good find that it fails for systems with poles or zeros at the origin, I agree it shouldn’t do that. So thanks for taking this on.

Something to consider: In the intervening years, scipy has introduced an impulse-matched conversion in cont2discrete called “impulse”. https://github.com/scipy/scipy/blob/v1.18.0/scipy/signal/_lti_conversion.py#L384-L583

Its implementation goes through a state-space conversion and back, but it looks clean and simple (only about 4 lines). Do you think it’s better to stick with your algorithm, that works directly with the poles and zeros of the TF, or switch to calling scipy, which is probably more tested?

(My take - just call the scipy code, but make sure it passes your unit tests. Edit: Unless you know that this direct-from-the-transfer function approach is better!)

I should be clear - I don’t know whether it is numerically more accurate to compute the impulse-matched conversion as it is done in this PR, from the transfer function directly, or as it is done in scipy using a state-space conversion.

Do you know ?

I think the question has a wrinkle in it: matched and scipy's impulse aren't two implementations of the same transform, so it isn't really a question of which is numerically better.

_c2d_matched is matched pole-zero: it maps z = exp(s*T) and then sets the gain from sysC.dcgain(), so DC gain is preserved exactly. scipy's impulse is impulse invariance: it makes the discrete impulse response equal T*h(kT), and DC gain is whatever falls out. Both on T = 0.1:

G(s) = 1/((s+1)(s+3)), continuous DC gain 0.333333

G(s) = 10/(s+2), continuous DC gain 5.0

Each is exact on its own invariant and wrong on the other — 10% off in DC gain for the first-order case. There's also a sample of delay: impulse invariance gives a strictly proper discrete system, matched pole-zero doesn't.

So switching method='matched' to call cont2discrete(..., 'impulse') would change what existing users get, not just how accurately they get it. If you want impulse invariance available I'd add it as its own method name and leave matched alone.

On the narrower numerical question I honestly don't know, and I didn't want to guess: matched needs tf2zpk root-finding, which is badly conditioned for high-order polynomials, while the state-space route builds a companion matrix and takes expm, which has its own conditioning problems. I'd expect both to degrade for high order, differently, but I haven't measured it and wouldn't want to claim a winner without doing so.

Either way this PR is the origin-pole/zero NaN, which is orthogonal — happy to keep it scoped to that.

I measured it. For the same transform neither route wins everywhere, but this PR's algorithm sidesteps the failure mode of both.

Method: a 120-digit mpmath reference per transform (verified against closed forms, including repeated poles, and against an independent residue oracle), both implementations fed identical float64 coefficients, error = max relative error of H(e^jw) on the unit circle. The middle and right columns are the same transform, impulse invariance, computed two ways, so they isolate the route.

The state-space route loses in ss2tf, not in expm: the denominator comes out right either way (coefficients agree to ~1e-15), but the numerator is poly(A - B@C) + (D-1)*poly(A), a difference of two characteristic polynomials. I measured that cancellation at 1e10-1e13 for n=8-10, and 2e12 at T=1e-6 even for n=2, which is where the last row comes from. The TF-direct route has the mirror problem — partial-fraction residues blow up for clustered poles — so it loses badly on Butterworth 10 and on poles -1..-10.

_c2d_matched hits neither, because it maps poles and zeros and multiplies out: no residues, no difference of characteristic polynomials. Against a floor of "compute exactly, then round to float64" it sits at or near that floor on everything I ran (worst case 60x, 8th-order Chebyshev I), including elliptic and cheby2 filters with a full set of finite zeros, so it isn't just riding a trivial numerator.

Independent of accuracy: scipy's impulse requires a strictly proper system, so it raises ValueError for a PI controller, PID, a differentiator, any biproper plant, and elliptic/cheby2 filters — most of what this PR unbreaks.

Where neither helps: at high order the TF representation is the limit, not the algorithm. That rounding floor is 3e-5 for Butterworth 10 and ~1.0 for poles -0.01..-1000 at T=1e-3, so nothing returning a float64 TransferFunction does better there; zpk or state space is the answer. Not tested: MIMO (matched is SISO-only) and anything other than float64.

Also, a correction to my earlier comment: impulse invariance is only strictly proper when the relative degree is at least 2 — for 10/(s+2) scipy returns D = T*C*B != 0.

FYI: CI failures have been addressed in #1243. Once that has been merged into main, you should be able to rebase and clear the errors.

Rebased onto main now that #1243 has merged. The branch head is f80d990; control/tests/discrete_test.py passes 109/109 and ruff is clean on the two touched files.

Sorry, just got back to my computer after being away for awhile. My mistake about the two different conversions - you're right, scipy's matched is matching impulse responses, not poles and zeros.

This looks good, thanks for the pr and interesting to see the result of numerical testing.