Skip to content

Commit dd95f35

Browse files
murrayrmmurrayrm
authored
Merge pull request #729 from roryyorke/rory/linfnorm
ENH: add `linform` to compute linear system L-infinity norm
2 parents 46b3c53 + 02b3451 commit dd95f35

File tree

2 files changed

+117
-2
lines changed

2 files changed

+117
-2
lines changed

‎control/statesp.py

Lines changed: 65 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,7 @@
5555
from numpy.linalg import solve, eigvals, matrix_rank
5656
from numpy.linalg.linalg import LinAlgError
5757
import scipy as sp
58+
import scipy.linalg
5859
from scipy.signal import cont2discrete
5960
from scipy.signal import StateSpace as signalStateSpace
6061
from warnings import warn
@@ -65,7 +66,12 @@
6566
from . import config
6667
from copy import deepcopy
6768

68-
__all__ = ['StateSpace', 'tf2ss', 'ssdata']
69+
try:
70+
from slycot import ab13dd
71+
except ImportError:
72+
ab13dd = None
73+
74+
__all__ = ['StateSpace', 'tf2ss', 'ssdata', 'linfnorm']
6975

7076
# Define module default parameter values
7177
_statesp_defaults = {
@@ -1895,3 +1901,61 @@ def ssdata(sys):
18951901
"""
18961902
ss = _convert_to_statespace(sys)
18971903
return ss.A, ss.B, ss.C, ss.D
1904+
1905+
1906+
def linfnorm(sys, tol=1e-10):
1907+
"""L-infinity norm of a linear system
1908+
1909+
Parameters
1910+
----------
1911+
sys : LTI (StateSpace or TransferFunction)
1912+
system to evalute L-infinity norm of
1913+
tol : real scalar
1914+
tolerance on norm estimate
1915+
1916+
Returns
1917+
-------
1918+
gpeak : non-negative scalar
1919+
L-infinity norm
1920+
fpeak : non-negative scalar
1921+
Frequency, in rad/s, at which gpeak occurs
1922+
1923+
For stable systems, the L-infinity and H-infinity norms are equal;
1924+
for unstable systems, the H-infinity norm is infinite, while the
1925+
L-infinity norm is finite if the system has no poles on the
1926+
imaginary axis.
1927+
1928+
See also
1929+
--------
1930+
slycot.ab13dd : the Slycot routine linfnorm that does the calculation
1931+
"""
1932+
1933+
if ab13dd is None:
1934+
raise ControlSlycot("Can't find slycot module 'ab13dd'")
1935+
1936+
a, b, c, d = ssdata(_convert_to_statespace(sys))
1937+
e = np.eye(a.shape[0])
1938+
1939+
n = a.shape[0]
1940+
m = b.shape[1]
1941+
p = c.shape[0]
1942+
1943+
if n == 0:
1944+
# ab13dd doesn't accept empty A, B, C, D;
1945+
# static gain case is easy enough to compute
1946+
gpeak = scipy.linalg.svdvals(d)[0]
1947+
# max svd is constant with freq; arbitrarily choose 0 as peak
1948+
fpeak = 0
1949+
return gpeak, fpeak
1950+
1951+
dico = 'C' if sys.isctime() else 'D'
1952+
jobe = 'I'
1953+
equil = 'S'
1954+
jobd = 'Z' if all(0 == d.flat) else 'D'
1955+
1956+
gpeak, fpeak = ab13dd(dico, jobe, equil, jobd, n, m, p, a, e, b, c, d, tol)
1957+
1958+
if dico=='D':
1959+
fpeak /= sys.dt
1960+
1961+
return gpeak, fpeak

‎control/tests/statesp_test.py

Lines changed: 52 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,13 +19,15 @@
1919
from control.dtime import sample_system
2020
from control.lti import evalfr
2121
from control.statesp import StateSpace, _convert_to_statespace, tf2ss, \
22-
_statesp_defaults, _rss_generate
22+
_statesp_defaults, _rss_generate, linfnorm
2323
from control.iosys import ss, rss, drss
2424
from control.tests.conftest import ismatarrayout, slycotonly
2525
from control.xferfcn import TransferFunction, ss2tf
2626

27+
2728
from .conftest import editsdefaults
2829

30+
2931
class TestStateSpace:
3032
"""Tests for the StateSpace class."""
3133

@@ -1107,3 +1109,52 @@ def test_latex_repr_testsize(editsdefaults):
11071109

11081110
gstatic = ss([], [], [], 1)
11091111
assert gstatic._repr_latex_() is None
1112+
1113+
1114+
class TestLinfnorm:
1115+
# these are simple tests; we assume ab13dd is correct
1116+
# python-control specific behaviour is:
1117+
# - checking for continuous- and discrete-time
1118+
# - scaling fpeak for discrete-time
1119+
# - handling static gains
1120+
1121+
# the underdamped gpeak and fpeak are found from
1122+
# gpeak = 1/(2*zeta*(1-zeta**2)**0.5)
1123+
# fpeak = wn*(1-2*zeta**2)**0.5
1124+
@pytest.fixture(params=[
1125+
('static', ct.tf, ([1.23],[1]), 1.23, 0),
1126+
('underdamped', ct.tf, ([100],[1, 2*0.5*10, 100]), 1.1547005, 7.0710678),
1127+
])
1128+
def ct_siso(self, request):
1129+
name, systype, sysargs, refgpeak, reffpeak = request.param
1130+
return systype(*sysargs), refgpeak, reffpeak
1131+
1132+
@pytest.fixture(params=[
1133+
('underdamped', ct.tf, ([100],[1, 2*0.5*10, 100]), 1e-4, 1.1547005, 7.0710678),
1134+
])
1135+
def dt_siso(self, request):
1136+
name, systype, sysargs, dt, refgpeak, reffpeak = request.param
1137+
return ct.c2d(systype(*sysargs), dt), refgpeak, reffpeak
1138+
1139+
@slycotonly
1140+
def test_linfnorm_ct_siso(self, ct_siso):
1141+
sys, refgpeak, reffpeak = ct_siso
1142+
gpeak, fpeak = linfnorm(sys)
1143+
np.testing.assert_allclose(gpeak, refgpeak)
1144+
np.testing.assert_allclose(fpeak, reffpeak)
1145+
1146+
@slycotonly
1147+
def test_linfnorm_dt_siso(self, dt_siso):
1148+
sys, refgpeak, reffpeak = dt_siso
1149+
gpeak, fpeak = linfnorm(sys)
1150+
# c2d pole-mapping has round-off
1151+
np.testing.assert_allclose(gpeak, refgpeak)
1152+
np.testing.assert_allclose(fpeak, reffpeak)
1153+
1154+
@slycotonly
1155+
def test_linfnorm_ct_mimo(self, ct_siso):
1156+
siso, refgpeak, reffpeak = ct_siso
1157+
sys = ct.append(siso, siso)
1158+
gpeak, fpeak = linfnorm(sys)
1159+
np.testing.assert_allclose(gpeak, refgpeak)
1160+
np.testing.assert_allclose(fpeak, reffpeak)

0 commit comments

Comments
 (0)