Skip to content

Commit 8b44e87

Browse files
murrayrmmurrayrm
committed
update return values to match numpy inf/nan pattern
1 parent 96e44fe commit 8b44e87

7 files changed

Lines changed: 141 additions & 83 deletions

File tree

‎control/margins.py‎

Lines changed: 1 addition & 2 deletions
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‎

Lines changed: 19 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -725,7 +725,9 @@ def slycot_laub(self, x):
725725
Frequency response
726726
"""
727727
from slycot import tb05ad
728-
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)
729731

730732
# Make sure that we are operating on a simple list
731733
if len(x_arr.shape) > 1:
@@ -790,7 +792,9 @@ def horner(self, x, warn_infinite=True):
790792
except (ImportError, Exception):
791793
# Fall back because either Slycot unavailable or cannot handle
792794
# certain cases.
793-
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)
794798

795799
# Make sure that we are operating on a simple list
796800
if len(x_arr.shape) > 1:
@@ -808,11 +812,18 @@ def horner(self, x, warn_infinite=True):
808812
solve(x_idx * eye(self.nstates) - self.A, self.B)) \
809813
+ self.D
810814
except LinAlgError:
815+
# Issue a warning messsage, for consistency with xferfcn
811816
if warn_infinite:
812-
warn("frequency response is not finite",
817+
warn("singular matrix in frequency response",
813818
RuntimeWarning)
814-
# TODO: check for nan cases
815-
out[:,:,idx] = np.inf
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+
816827
return out
817828

818829
def freqresp(self, omega):
@@ -1193,7 +1204,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
11931204
Ad, Bd, C, D, _ = cont2discrete(sys, Twarp, method, alpha)
11941205
return StateSpace(Ad, Bd, C, D, Ts)
11951206

1196-
def dcgain(self):
1207+
def dcgain(self, warn_infinite=False):
11971208
"""Return the zero-frequency gain
11981209
11991210
The zero-frequency gain of a continuous-time state-space
@@ -1212,16 +1223,8 @@ def dcgain(self):
12121223
be the zero-frequency (or DC) gain, or, if the frequency
12131224
response is singular, the array will be filled with np.inf.
12141225
"""
1215-
try:
1216-
if self.isctime():
1217-
gain = np.asarray(self.D -
1218-
self.C.dot(np.linalg.solve(self.A, self.B)))
1219-
else:
1220-
gain = np.squeeze(self.horner(1))
1221-
except LinAlgError:
1222-
# eigenvalue at DC
1223-
gain = np.tile(np.inf, (self.noutputs, self.ninputs))
1224-
return np.squeeze(gain)
1226+
return self(0, warn_infinite=warn_infinite) if self.isctime() \
1227+
else self(1, warn_infinite=warn_infinite)
12251228

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

‎control/tests/freqresp_test.py‎

Lines changed: 76 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -367,7 +367,7 @@ def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
367367

368368

369369
def test_freqresp_warn_infinite():
370-
"""Test evaluation of transfer functions at the origin"""
370+
"""Test evaluation warnings for transfer functions w/ pole at the origin"""
371371
sys_finite = ctrl.tf([1], [1, 0.01])
372372
sys_infinite = ctrl.tf([1], [1, 0.01, 0])
373373

@@ -378,11 +378,13 @@ def test_freqresp_warn_infinite():
378378

379379
# Transfer function with infinite zero frequency gain
380380
with pytest.warns(RuntimeWarning, match="divide by zero"):
381-
np.testing.assert_almost_equal(sys_infinite(0), np.inf)
381+
np.testing.assert_almost_equal(
382+
sys_infinite(0), complex(np.inf, np.nan))
382383
with pytest.warns(RuntimeWarning, match="divide by zero"):
383384
np.testing.assert_almost_equal(
384-
sys_infinite(0, warn_infinite=True), np.inf)
385-
np.testing.assert_almost_equal(sys_infinite(0, warn_infinite=False), np.inf)
385+
sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
386+
np.testing.assert_almost_equal(
387+
sys_infinite(0, warn_infinite=False), complex(np.inf, np.nan))
386388

387389
# Switch to state space
388390
sys_finite = ctrl.tf2ss(sys_finite)
@@ -394,13 +396,15 @@ def test_freqresp_warn_infinite():
394396
np.testing.assert_almost_equal(sys_finite(0, warn_infinite=True), 100)
395397

396398
# State space system with infinite zero frequency gain
397-
with pytest.warns(RuntimeWarning, match="not finite"):
398-
np.testing.assert_almost_equal(sys_infinite(0), np.inf)
399-
with pytest.warns(RuntimeWarning, match="not finite"):
400-
np.testing.assert_almost_equal(sys_infinite(0), np.inf)
401-
np.testing.assert_almost_equal(sys_infinite(0, warn_infinite=True), np.inf)
402-
np.testing.assert_almost_equal(sys_infinite(0, warn_infinite=False), np.inf)
403-
399+
with pytest.warns(RuntimeWarning, match="singular matrix"):
400+
np.testing.assert_almost_equal(
401+
sys_infinite(0), complex(np.inf, np.nan))
402+
with pytest.warns(RuntimeWarning, match="singular matrix"):
403+
np.testing.assert_almost_equal(
404+
sys_infinite(0, warn_infinite=True), complex(np.inf, np.nan))
405+
np.testing.assert_almost_equal(sys_infinite(
406+
0, warn_infinite=False), complex(np.inf, np.nan))
407+
404408

405409
def test_dcgain_consistency():
406410
"""Test to make sure that DC gain is consistently evaluated"""
@@ -412,25 +416,74 @@ def test_dcgain_consistency():
412416
sys_ss = ctrl.tf2ss(sys_tf)
413417
assert 0 in sys_ss.pole()
414418

415-
# Evaluation
416-
np.testing.assert_equal(sys_tf(0), np.inf + 0j)
417-
np.testing.assert_equal(sys_ss(0), np.inf + 0j)
418-
np.testing.assert_equal(sys_tf.dcgain(), np.inf + 0j)
419-
np.testing.assert_equal(sys_ss.dcgain(), np.inf + 0j)
419+
# Finite (real) numerator over 0 denominator => inf + nanj
420+
np.testing.assert_equal(
421+
sys_tf(0, warn_infinite=False), complex(np.inf, np.nan))
422+
np.testing.assert_equal(
423+
sys_ss(0, warn_infinite=False), complex(np.inf, np.nan))
424+
np.testing.assert_equal(
425+
sys_tf(0j, warn_infinite=False), complex(np.inf, np.nan))
426+
np.testing.assert_equal(
427+
sys_ss(0j, warn_infinite=False), complex(np.inf, np.nan))
428+
np.testing.assert_equal(
429+
sys_tf.dcgain(warn_infinite=False), complex(np.inf, np.nan))
430+
np.testing.assert_equal(
431+
sys_ss.dcgain(warn_infinite=False), complex(np.inf, np.nan))
420432

421433
# Set up transfer function with pole, zero at the origin
422434
sys_tf = ctrl.tf([1, 0], [1, 0])
423435
assert 0 in sys_tf.pole()
424436
assert 0 in sys_tf.zero()
425-
437+
426438
sys_ss = ctrl.tf2ss(ctrl.tf([1, 0], [1, 1])) * \
427439
ctrl.tf2ss(ctrl.tf([1], [1, 0]))
428440
assert 0 in sys_ss.pole()
429441
assert 0 in sys_ss.zero()
430442

431-
# Pole and zero at the origin should give nan for the response
432-
np.testing.assert_equal(sys_tf(0), np.nan)
433-
np.testing.assert_equal(sys_tf.dcgain(), np.nan)
434-
# TODO: state space cases not yet working
435-
# np.testing.assert_equal(sys_ss(0), np.nan)
436-
# np.testing.assert_equal(sys_ss.dcgain(), np.nan)
443+
# Pole and zero at the origin should give nan + nanj for the response
444+
np.testing.assert_equal(
445+
sys_tf(0, warn_infinite=False), complex(np.nan, np.nan))
446+
np.testing.assert_equal(
447+
sys_tf(0j, warn_infinite=False), complex(np.nan, np.nan))
448+
np.testing.assert_equal(
449+
sys_tf.dcgain(warn_infinite=False), complex(np.nan, np.nan))
450+
np.testing.assert_equal(
451+
sys_ss(0, warn_infinite=False), complex(np.nan, np.nan))
452+
np.testing.assert_equal(
453+
sys_ss(0j, warn_infinite=False), complex(np.nan, np.nan))
454+
np.testing.assert_equal(
455+
sys_ss.dcgain(warn_infinite=False), complex(np.nan, np.nan))
456+
457+
# Pole with non-zero, complex numerator => inf + infj
458+
s = ctrl.tf('s')
459+
sys_tf = (s + 1) / (s**2 + 1)
460+
assert 1j in sys_tf.pole()
461+
462+
# Set up state space system with pole on imaginary axis
463+
sys_ss = ctrl.tf2ss(sys_tf)
464+
assert 1j in sys_tf.pole()
465+
466+
# Make sure we get correct response if evaluated at the pole
467+
np.testing.assert_equal(
468+
sys_tf(1j, warn_infinite=False), complex(np.inf, np.inf))
469+
470+
# For state space, numerical errors come into play
471+
resp_ss = sys_ss(1j, warn_infinite=False)
472+
if np.isfinite(resp_ss):
473+
assert abs(resp_ss) > 1e15
474+
else:
475+
if resp_ss != complex(np.inf, np.inf):
476+
pytest.xfail("statesp evaluation at poles not fully implemented")
477+
else:
478+
np.testing.assert_equal(resp_ss, complex(np.inf, np.inf))
479+
480+
# DC gain is finite
481+
np.testing.assert_almost_equal(sys_tf.dcgain(), 1.)
482+
np.testing.assert_almost_equal(sys_ss.dcgain(), 1.)
483+
484+
# Make sure that we get the *signed* DC gain
485+
sys_tf = -1 / (s + 1)
486+
np.testing.assert_almost_equal(sys_tf.dcgain(), -1)
487+
488+
sys_ss = ctrl.tf2ss(sys_tf)
489+
np.testing.assert_almost_equal(sys_ss.dcgain(), -1)

‎control/tests/input_element_int_test.py‎

Lines changed: 5 additions & 5 deletions
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‎

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
rss, ss, tf2ss, _statesp_defaults)
2222
from control.tests.conftest import ismatarrayout, slycotonly
2323
from control.xferfcn import TransferFunction, ss2tf
24+
from control.exception import ControlSlycot
2425

2526
from .conftest import editsdefaults
2627

@@ -498,7 +499,7 @@ def test_dc_gain_cont(self):
498499
np.testing.assert_allclose(sys2.dcgain(), expected)
499500

500501
sys3 = StateSpace(0., 1., 1., 0.)
501-
np.testing.assert_equal(sys3.dcgain(), np.inf)
502+
np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
502503

503504
def test_dc_gain_discr(self):
504505
"""Test DC gain for discrete-time state-space systems."""
@@ -516,7 +517,7 @@ def test_dc_gain_discr(self):
516517

517518
# summer
518519
sys = StateSpace(1, 1, 1, 0, True)
519-
np.testing.assert_equal(sys.dcgain(), np.inf)
520+
np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
520521

521522
@pytest.mark.parametrize("outputs", range(1, 6))
522523
@pytest.mark.parametrize("inputs", range(1, 6))
@@ -539,8 +540,15 @@ def test_dc_gain_integrator(self, outputs, inputs, dt):
539540
c = np.eye(max(outputs, states))[:outputs, :states]
540541
d = np.zeros((outputs, inputs))
541542
sys = StateSpace(a, b, c, d, dt)
542-
dc = np.squeeze(np.full_like(d, np.inf))
543-
np.testing.assert_array_equal(dc, sys.dcgain())
543+
dc = np.full_like(d, complex(np.inf, np.nan), dtype=complex)
544+
if sys.issiso():
545+
dc = dc.squeeze()
546+
547+
try:
548+
np.testing.assert_array_equal(dc, sys.dcgain())
549+
except NotImplementedError:
550+
# Skip MIMO tests if there is no slycot
551+
pytest.skip("slycot required for MIMO dcgain")
544552

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

‎control/tests/xferfcn_test.py‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -807,7 +807,7 @@ def test_dcgain_cont(self):
807807
np.testing.assert_equal(sys2.dcgain(), 2)
808808

809809
sys3 = TransferFunction(6, [1, 0])
810-
np.testing.assert_equal(sys3.dcgain(), np.inf)
810+
np.testing.assert_equal(sys3.dcgain(), complex(np.inf, np.nan))
811811

812812
num = [[[15], [21], [33]], [[10], [14], [22]]]
813813
den = [[[1, 3], [2, 3], [3, 3]], [[1, 5], [2, 7], [3, 11]]]
@@ -827,8 +827,13 @@ def test_dcgain_discr(self):
827827

828828
# differencer
829829
sys = TransferFunction(1, [1, -1], True)
830+
np.testing.assert_equal(sys.dcgain(), complex(np.inf, np.nan))
831+
832+
# differencer, with warning
833+
sys = TransferFunction(1, [1, -1], True)
830834
with pytest.warns(RuntimeWarning, match="divide by zero"):
831-
np.testing.assert_equal(sys.dcgain(), np.inf)
835+
np.testing.assert_equal(
836+
sys.dcgain(warn_infinite=True), complex(np.inf, np.nan))
832837

833838
# summer
834839
sys = TransferFunction([1, -1], [1], True)

0 commit comments

Comments
 (0)