Skip to content

Commit b8b828e

Browse files
murrayrmmurrayrm
authored
Merge pull request #542 from murrayrm/infinite_gain
Fix inconsistencies in frequency response at poles
2 parents 9c6c619 + 598058c commit b8b828e

8 files changed

+229
-72
lines changed

‎control/lti.py

+4-2
Original file line numberDiff line numberDiff line change
@@ -644,8 +644,10 @@ def dcgain(sys):
644644
Returns
645645
-------
646646
gain : ndarray
647-
The zero-frequency gain, or np.nan if the system has a pole
648-
at the origin
647+
The zero-frequency gain, or (inf + nanj) if the system has a pole at
648+
the origin, (nan + nanj) if there is a pole/zero cancellation at the
649+
origin.
650+
649651
"""
650652
return sys.dcgain()
651653

‎control/margins.py

+1-2
Original file line numberDiff line numberDiff line change
@@ -294,8 +294,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
294294
num_iw, den_iw = _poly_iw(sys)
295295
# frequency for gain margin: phase crosses -180 degrees
296296
w_180 = _poly_iw_real_crossing(num_iw, den_iw, epsw)
297-
with np.errstate(all='ignore'): # den=0 is okay
298-
w180_resp = sys(1J * w_180)
297+
w180_resp = sys(1J * w_180, warn_infinite=False) # den=0 is okay
299298

300299
# frequency for phase margin : gain crosses magnitude 1
301300
wc = _poly_iw_mag1_crossing(num_iw, den_iw, epsw)

‎control/statesp.py

+37-23
Original file line numberDiff line numberDiff line change
@@ -668,7 +668,7 @@ def __rdiv__(self, other):
668668
raise NotImplementedError(
669669
"StateSpace.__rdiv__ is not implemented yet.")
670670

671-
def __call__(self, x, squeeze=None):
671+
def __call__(self, x, squeeze=None, warn_infinite=True):
672672
"""Evaluate system's transfer function at complex frequency.
673673
674674
Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -689,6 +689,8 @@ def __call__(self, x, squeeze=None):
689689
keep all indices (output, input and, if omega is array_like,
690690
frequency) even if the system is SISO. The default value can be
691691
set using config.defaults['control.squeeze_frequency_response'].
692+
warn_infinite : bool, optional
693+
If set to `False`, don't warn if frequency response is infinite.
692694
693695
Returns
694696
-------
@@ -702,7 +704,7 @@ def __call__(self, x, squeeze=None):
702704
703705
"""
704706
# Use Slycot if available
705-
out = self.horner(x)
707+
out = self.horner(x, warn_infinite=warn_infinite)
706708
return _process_frequency_response(self, x, out, squeeze=squeeze)
707709

708710
def slycot_laub(self, x):
@@ -723,7 +725,9 @@ def slycot_laub(self, x):
723725
Frequency response
724726
"""
725727
from slycot import tb05ad
726-
x_arr = np.atleast_1d(x) # array-like version of x
728+
729+
# Make sure the argument is a 1D array of complex numbers
730+
x_arr = np.atleast_1d(x).astype(complex, copy=False)
727731

728732
# Make sure that we are operating on a simple list
729733
if len(x_arr.shape) > 1:
@@ -758,7 +762,7 @@ def slycot_laub(self, x):
758762
out[:, :, kk+1] = result[0] + self.D
759763
return out
760764

761-
def horner(self, x):
765+
def horner(self, x, warn_infinite=True):
762766
"""Evaluate system's transfer function at complex frequency
763767
using Laub's or Horner's method.
764768
@@ -788,21 +792,38 @@ def horner(self, x):
788792
except (ImportError, Exception):
789793
# Fall back because either Slycot unavailable or cannot handle
790794
# certain cases.
791-
x_arr = np.atleast_1d(x) # force to be an array
795+
796+
# Make sure the argument is a 1D array of complex numbers
797+
x_arr = np.atleast_1d(x).astype(complex, copy=False)
792798

793799
# Make sure that we are operating on a simple list
794800
if len(x_arr.shape) > 1:
795801
raise ValueError("input list must be 1D")
796802

797803
# Preallocate
798-
out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex)
804+
out = empty((self.noutputs, self.ninputs, len(x_arr)),
805+
dtype=complex)
799806

800807
#TODO: can this be vectorized?
801808
for idx, x_idx in enumerate(x_arr):
802-
out[:,:,idx] = \
803-
np.dot(self.C,
809+
try:
810+
out[:,:,idx] = np.dot(
811+
self.C,
804812
solve(x_idx * eye(self.nstates) - self.A, self.B)) \
805-
+ self.D
813+
+ self.D
814+
except LinAlgError:
815+
# Issue a warning messsage, for consistency with xferfcn
816+
if warn_infinite:
817+
warn("singular matrix in frequency response",
818+
RuntimeWarning)
819+
820+
# Evaluating at a pole. Return value depends if there
821+
# is a zero at the same point or not.
822+
if x_idx in self.zero():
823+
out[:,:,idx] = complex(np.nan, np.nan)
824+
else:
825+
out[:,:,idx] = complex(np.inf, np.nan)
826+
806827
return out
807828

808829
def freqresp(self, omega):
@@ -1183,7 +1204,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
11831204
Ad, Bd, C, D, _ = cont2discrete(sys, Twarp, method, alpha)
11841205
return StateSpace(Ad, Bd, C, D, Ts)
11851206

1186-
def dcgain(self):
1207+
def dcgain(self, warn_infinite=False):
11871208
"""Return the zero-frequency gain
11881209
11891210
The zero-frequency gain of a continuous-time state-space
@@ -1198,20 +1219,13 @@ def dcgain(self):
11981219
Returns
11991220
-------
12001221
gain : ndarray
1201-
An array of shape (outputs,inputs); the array will either
1202-
be the zero-frequency (or DC) gain, or, if the frequency
1203-
response is singular, the array will be filled with np.nan.
1222+
An array of shape (outputs,inputs); the array will either be the
1223+
zero-frequency (or DC) gain, or, if the frequency response is
1224+
singular, the array will be filled with (inf + nanj).
1225+
12041226
"""
1205-
try:
1206-
if self.isctime():
1207-
gain = np.asarray(self.D -
1208-
self.C.dot(np.linalg.solve(self.A, self.B)))
1209-
else:
1210-
gain = np.squeeze(self.horner(1))
1211-
except LinAlgError:
1212-
# eigenvalue at DC
1213-
gain = np.tile(np.nan, (self.noutputs, self.ninputs))
1214-
return np.squeeze(gain)
1227+
return self(0, warn_infinite=warn_infinite) if self.isctime() \
1228+
else self(1, warn_infinite=warn_infinite)
12151229

12161230
def _isstatic(self):
12171231
"""True if and only if the system has no dynamics, that is,

‎control/tests/freqresp_test.py

+137
Original file line numberDiff line numberDiff line change
@@ -371,3 +371,140 @@ def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
371371
mag, phase, omega = ctrl.bode(TF, wrap_phase=wrap_phase)
372372
assert(min(phase) >= min_phase)
373373
assert(max(phase) <= max_phase)
374+
375+
376+
def test_freqresp_warn_infinite():
377+
"""Test evaluation warnings for transfer functions w/ pole at the origin"""
378+
sys_finite = ctrl.tf([1], [1, 0.01])
379+
sys_infinite = ctrl.tf([1], [1, 0.01, 0])
380+
381+
# Transfer function with finite zero frequency gain
382+
np.testing.assert_almost_equal(sys_finite(0), 100)
383+
np.testing.assert_almost_equal(sys_finite(0, warn_infinite=False), 100)
384+
np.testing.assert_almost_equal(sys_finite(0, warn_infinite=True), 100)
385+
386+
# Transfer function with infinite zero frequency gain
387+
with pytest.warns(RuntimeWarning, match="divide by zero"):
388+
np.testing.assert_almost_equal(
389+
sys_infinite(0), complex(np.inf, np.nan))
390+
with pytest.warns(RuntimeWarning, match="divide by zero"):
391+
np.testing.assert_almost_equal(
392+
sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
393+
np.testing.assert_almost_equal(
394+
sys_infinite(0, warn_infinite=False), complex(np.inf, np.nan))
395+
396+
# Switch to state space
397+
sys_finite = ctrl.tf2ss(sys_finite)
398+
sys_infinite = ctrl.tf2ss(sys_infinite)
399+
400+
# State space system with finite zero frequency gain
401+
np.testing.assert_almost_equal(sys_finite(0), 100)
402+
np.testing.assert_almost_equal(sys_finite(0, warn_infinite=False), 100)
403+
np.testing.assert_almost_equal(sys_finite(0, warn_infinite=True), 100)
404+
405+
# State space system with infinite zero frequency gain
406+
with pytest.warns(RuntimeWarning, match="singular matrix"):
407+
np.testing.assert_almost_equal(
408+
sys_infinite(0), complex(np.inf, np.nan))
409+
with pytest.warns(RuntimeWarning, match="singular matrix"):
410+
np.testing.assert_almost_equal(
411+
sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
412+
np.testing.assert_almost_equal(sys_infinite(
413+
0, warn_infinite=False), complex(np.inf, np.nan))
414+
415+
416+
def test_dcgain_consistency():
417+
"""Test to make sure that DC gain is consistently evaluated"""
418+
# Set up transfer function with pole at the origin
419+
sys_tf = ctrl.tf([1], [1, 0])
420+
assert 0 in sys_tf.pole()
421+
422+
# Set up state space system with pole at the origin
423+
sys_ss = ctrl.tf2ss(sys_tf)
424+
assert 0 in sys_ss.pole()
425+
426+
# Finite (real) numerator over 0 denominator => inf + nanj
427+
np.testing.assert_equal(
428+
sys_tf(0, warn_infinite=False), complex(np.inf, np.nan))
429+
np.testing.assert_equal(
430+
sys_ss(0, warn_infinite=False), complex(np.inf, np.nan))
431+
np.testing.assert_equal(
432+
sys_tf(0j, warn_infinite=False), complex(np.inf, np.nan))
433+
np.testing.assert_equal(
434+
sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
435+
np.testing.assert_equal(
436+
sys_tf.dcgain(warn_infinite=False), complex(np.inf, np.nan))
437+
np.testing.assert_equal(
438+
sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
439+
440+
# Set up transfer function with pole, zero at the origin
441+
sys_tf = ctrl.tf([1, 0], [1, 0])
442+
assert 0 in sys_tf.pole()
443+
assert 0 in sys_tf.zero()
444+
445+
# Pole and zero at the origin should give nan + nanj for the response
446+
np.testing.assert_equal(
447+
sys_tf(0, warn_infinite=False), complex(np.nan, np.nan))
448+
np.testing.assert_equal(
449+
sys_tf(0j, warn_infinite=False), complex(np.nan, np.nan))
450+
np.testing.assert_equal(
451+
sys_tf.dcgain(warn_infinite=False), complex(np.nan, np.nan))
452+
453+
# Set up state space version
454+
sys_ss = ctrl.tf2ss(ctrl.tf([1, 0], [1, 1])) * \
455+
ctrl.tf2ss(ctrl.tf([1], [1, 0]))
456+
457+
# Different systems give different representations => test accordingly
458+
if 0 in sys_ss.pole() and 0 in sys_ss.zero():
459+
# Pole and zero at the origin => should get (nan + nanj)
460+
np.testing.assert_equal(
461+
sys_ss(0, warn_infinite=False), complex(np.nan, np.nan))
462+
np.testing.assert_equal(
463+
sys_ss(0j, warn_infinite=False), complex(np.nan, np.nan))
464+
np.testing.assert_equal(
465+
sys_ss.dcgain(warn_infinite=False), complex(np.nan, np.nan))
466+
elif 0 in sys_ss.pole():
467+
# Pole at the origin, but zero elsewhere => should get (inf + nanj)
468+
np.testing.assert_equal(
469+
sys_ss(0, warn_infinite=False), complex(np.inf, np.nan))
470+
np.testing.assert_equal(
471+
sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
472+
np.testing.assert_equal(
473+
sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
474+
else:
475+
# Near pole/zero cancellation => nothing sensible to check
476+
pass
477+
478+
# Pole with non-zero, complex numerator => inf + infj
479+
s = ctrl.tf('s')
480+
sys_tf = (s + 1) / (s**2 + 1)
481+
assert 1j in sys_tf.pole()
482+
483+
# Set up state space system with pole on imaginary axis
484+
sys_ss = ctrl.tf2ss(sys_tf)
485+
assert 1j in sys_tf.pole()
486+
487+
# Make sure we get correct response if evaluated at the pole
488+
np.testing.assert_equal(
489+
sys_tf(1j, warn_infinite=False), complex(np.inf, np.inf))
490+
491+
# For state space, numerical errors come into play
492+
resp_ss = sys_ss(1j, warn_infinite=False)
493+
if np.isfinite(resp_ss):
494+
assert abs(resp_ss) > 1e15
495+
else:
496+
if resp_ss != complex(np.inf, np.inf):
497+
pytest.xfail("statesp evaluation at poles not fully implemented")
498+
else:
499+
np.testing.assert_equal(resp_ss, complex(np.inf, np.inf))
500+
501+
# DC gain is finite
502+
np.testing.assert_almost_equal(sys_tf.dcgain(), 1.)
503+
np.testing.assert_almost_equal(sys_ss.dcgain(), 1.)
504+
505+
# Make sure that we get the *signed* DC gain
506+
sys_tf = -1 / (s + 1)
507+
np.testing.assert_almost_equal(sys_tf.dcgain(), -1)
508+
509+
sys_ss = ctrl.tf2ss(sys_tf)
510+
np.testing.assert_almost_equal(sys_ss.dcgain(), -1)

‎control/tests/input_element_int_test.py

+5-5
Original file line numberDiff line numberDiff line change
@@ -23,15 +23,15 @@ def test_tf_den_with_numpy_int_element(self):
2323

2424
sys = tf(num, den)
2525

26-
np.testing.assert_array_max_ulp(1., dcgain(sys))
26+
np.testing.assert_almost_equal(1., dcgain(sys))
2727

2828
def test_tf_num_with_numpy_int_element(self):
2929
num = np.convolve([1], [1, 1])
3030
den = np.convolve([1, 2, 1], [1, 1, 1])
3131

3232
sys = tf(num, den)
3333

34-
np.testing.assert_array_max_ulp(1., dcgain(sys))
34+
np.testing.assert_almost_equal(1., dcgain(sys))
3535

3636
# currently these pass
3737
def test_tf_input_with_int_element(self):
@@ -40,7 +40,7 @@ def test_tf_input_with_int_element(self):
4040

4141
sys = tf(num, den)
4242

43-
np.testing.assert_array_max_ulp(1., dcgain(sys))
43+
np.testing.assert_almost_equal(1., dcgain(sys))
4444

4545
def test_ss_input_with_int_element(self):
4646
a = np.array([[0, 1],
@@ -52,7 +52,7 @@ def test_ss_input_with_int_element(self):
5252

5353
sys = ss(a, b, c, d)
5454
sys2 = tf(sys)
55-
np.testing.assert_array_max_ulp(dcgain(sys), dcgain(sys2))
55+
np.testing.assert_almost_equal(dcgain(sys), dcgain(sys2))
5656

5757
def test_ss_input_with_0int_dcgain(self):
5858
a = np.array([[0, 1],
@@ -63,4 +63,4 @@ def test_ss_input_with_0int_dcgain(self):
6363
d = 0
6464
sys = ss(a, b, c, d)
6565
np.testing.assert_allclose(dcgain(sys), 0,
66-
atol=np.finfo(np.float).epsneg)
66+
atol=np.finfo(np.float).epsneg)

‎control/tests/statesp_test.py

+12-5
Original file line numberDiff line numberDiff line change
@@ -498,7 +498,7 @@ def test_dc_gain_cont(self):
498498
np.testing.assert_allclose(sys2.dcgain(), expected)
499499

500500
sys3 = StateSpace(0., 1., 1., 0.)
501-
np.testing.assert_equal(sys3.dcgain(), np.nan)
501+
np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
502502

503503
def test_dc_gain_discr(self):
504504
"""Test DC gain for discrete-time state-space systems."""
@@ -516,14 +516,14 @@ def test_dc_gain_discr(self):
516516

517517
# summer
518518
sys = StateSpace(1, 1, 1, 0, True)
519-
np.testing.assert_equal(sys.dcgain(), np.nan)
519+
np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
520520

521521
@pytest.mark.parametrize("outputs", range(1, 6))
522522
@pytest.mark.parametrize("inputs", range(1, 6))
523523
@pytest.mark.parametrize("dt", [None, 0, 1, True],
524524
ids=["dtNone", "c", "dt1", "dtTrue"])
525525
def test_dc_gain_integrator(self, outputs, inputs, dt):
526-
"""DC gain when eigenvalue at DC returns appropriately sized array of nan.
526+
"""DC gain w/ pole at origin returns appropriately sized array of inf.
527527
528528
the SISO case is also tested in test_dc_gain_{cont,discr}
529529
time systems (dt=0)
@@ -539,8 +539,15 @@ def test_dc_gain_integrator(self, outputs, inputs, dt):
539539
c = np.eye(max(outputs, states))[:outputs, :states]
540540
d = np.zeros((outputs, inputs))
541541
sys = StateSpace(a, b, c, d, dt)
542-
dc = np.squeeze(np.full_like(d, np.nan))
543-
np.testing.assert_array_equal(dc, sys.dcgain())
542+
dc = np.full_like(d, complex(np.inf, np.nan), dtype=complex)
543+
if sys.issiso():
544+
dc = dc.squeeze()
545+
546+
try:
547+
np.testing.assert_array_equal(dc, sys.dcgain())
548+
except NotImplementedError:
549+
# Skip MIMO tests if there is no slycot
550+
pytest.skip("slycot required for MIMO dcgain")
544551

545552
def test_scalar_static_gain(self):
546553
"""Regression: can we create a scalar static gain?

0 commit comments

Comments
 (0)