10000 ENH: add `linform` to compute linear system L-infinity norm · python-control/python-control@02b3451 · GitHub
[go: up one dir, main page]
Skip to content

Commit 02b3451

Browse files
roryyorkeroryyorke
committed
ENH: add linform to compute linear system L-infinity norm
linfnorm requires Slycot routine ab13dd, which does the actual work. Added tests to check that wrapping of ab13dd is correct, that static systems are handled, and that discrete-time system frequencies are scaled correctly.
1 parent 983726c commit 02b3451

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
@@ -63,7 +64,12 @@
6364
from . import config
6465
from copy import deepcopy
6566

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

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

@@ -1097,3 +1099,52 @@ def test_latex_repr_testsize(editsdefaults):
10971099

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

0 commit comments

Comments
 (0)
0