TST: Added tests from Byers and Nash

harold/_system_props.py

@@ -155,7 +155,7 @@ def system_norm(G, p=np.inf, hinf_tol=1e-6, eig_tol=1e-8):
         return (gamma_lb + gamma_ub)/2
 
 
-def controllability_indices(A, B):
+def controllability_indices(A, B, tol=None):
     """Computes the controllability indices for a controllable pair (A, B)
 
     Controllability indices are defined as the maximum number of independent
@@ -196,6 +196,9 @@ def controllability_indices(A, B):
         2D (n, n) real-valued array
     B : ndarray
         2D (n, m) real-valued array
+    tol : float
+        Tolerance value for the Arnoldi iteration to decide linear dependence.
+        By default it is sqrt(eps)*n**2
 
     Returns
     -------
@@ -220,7 +223,7 @@ def controllability_indices(A, B):
         raise ValueError("A and B should have the same number of rows")
 
     # FIXME: Tolerance is a bit arbitrary for now!!
-    tol = np.sqrt(np.spacing(1.))*n**2
+    tol = np.sqrt(np.spacing(1.))*n**2 if tol is None else tol
 
     # Will be populated below
     remaining_cols = np.arange(m)

harold/tests/test_static_ctrl_design.py

@@ -1,6 +1,7 @@
 from numpy import eye, array, sort
 from scipy.linalg import block_diag, eigvals
 from numpy.testing import assert_almost_equal, assert_array_almost_equal
+
 from pytest import raises as assert_raises
 from harold import lqr, ackermann, State, Transfer, haroldcompanion
 
@@ -82,3 +83,118 @@ def test_ackermann_uncontrollable():
     B = eye(4)[:, [-2]]
     p = [-10, -9, -8, -7]
     assert_raises(ValueError, ackermann, (A, B), p)
+
+
+def byersnash_A_B_test_pairs():
+    ABs = [
+           # Chemical Reactor (Munro 1979)
+           (array([[1.38, -0.2077, 6.715, -5.676],
+                   [-0.5814, -4.29, 0, 0.675],
+                   [1.067, 4.273, -6.654, 5.893],
+                   [0.048, 4.273, 1.343, -2.104]]),
+            array([[0, 0],
+                   [5.679, 0],
+                   [1.136, -3.146],
+                   [1.136, 0]])),
+           # Distillation Column (Klein, Moore 1977)
+           (array([[-0.1094, 0.0628, 0, 0, 0],
+                   [1.306, -2.132, 0.9807, 0, 0],
+                   [0, 1.595, -3.149, 1.547, 0],
+                   [0, 0.0355, 2.632, -4.257, 1.855],
+                   [0, 0.0023, 0, 0.1636, -0.1625]]),
+            array([[0, 0],
+                   [0.638, 0],
+                   [0.0838, -0.1396],
+                   [0.1004, -0.206],
+                   [0.0063, -0.0128]])),
+           # Nuclear rocket engine (Davison, Chow 1974)
+           (array([[-65.0, 65, -19.5, 19.5],
+                   [0.1, -0.1, 0, 0],
+                   [1, 0, -0.5, -1],
+                   [0, 0, 0.4, 0]]),
+            array([[65., 0],
+                   [0, 0],
+                   [0, 0],
+                   [0, 0.4]])),
+           # MIMO system (Atkinson, 1985)
+           (array([[0, 1, 0],
+                   [0, 0, 1],
+                   [-6, -11, -6]]),
+            array([[1, 1],
+                   [0, 1],
+                   [1, 1]])),
+           # Drum boiler (Bengtsson 1973)
+           (array([[-0.129, 0, 0.396, 0.25, 0.0019],
+                   [0.0329, 0, -0.0078, 0.0122, -0.621],
+                   [0.0072, 0, -0.1, 0.0009, -0.0385],
+                   [0.0041, 0, 0, -0.0822, 0],
+                   [0.0035, 0, 0.0035, 0.0043, -0.0743]]),
+            array([[0, 0.1390],
+                   [0, 0.0359],
+                   [0, -0.0989],
+                   [0.0249, 0],
+                   [0, -0.0053]])),
+           # Miminis random example #1
+           (array([[5.8765, 9.3456, 4.5634, 9.3520],
+                   [6.6526, 0.5867, 3.5829, 0.6534],
+                   [0.0000, 9.6738, 7.4876, 4.7654],
+                   [0.0000, 0.0000, 6.6784, 2.5678]]),
+            array([[3.9878, 0.5432],
+                   [0.0000, 2.7650],
+                   [0.0000, 0.0000],
+                   [0.0000, 0.0000]])),
+           # Miminis random example #2
+           (array([[.5257, .8544, .5596, .5901, .0259, .6213, .7227, .5617],
+                   [.9931, .0643, .1249, .3096, .5174, .3455, .8977, .4682],
+                   [.6489, .8279, .7279, .2552, .3917, .7065, .2428, .7795],
+                   [.9923, .9262, .2678, .6252, .2414, .5211, .4338, .9677],
+                   [.0000, .5667, .5465, .1157, .5064, .2870, .7901, .9809],
+                   [.0000, .0000, .8672, .6117, .4236, .6503, .5069, .8187],
+                   [.0000, .0000, .0000, .0000, .2894, .0881, .5233, .4257],
+                   [.0000, .0000, .0000, .0000, .0000, .4499, .5597, .2462]]),
+            array([[0.9230, 0.3950, 0.8325],
+                   [0.0000, 0.0366, 0.6105],
+                   [0.0000, 0.0000, 0.1871],
+                   [0.0000, 0.0000, 0.0000],
+                   [0.0000, 0.0000, 0.0000],
+                   [0.0000, 0.0000, 0.0000],
+                   [0.0000, 0.0000, 0.0000],
+                   [0.0000, 0.0000, 0.0000]])),
+           # Aircraft control example I (Kautsky and Nichols 1983)
+           (array([[0, 1, 0, 0],
+                   [1.40e-4, -2.04, -1.95, -1.33e-2],
+                   [-2.51e-4, 1, -1.32, -2.38e-2],
+                   [-5.61e-1, 0, 0.358, -2.79e-1]]),
+            array([[0, 0, 0],
+                   [-5.33, 6.45e-3, -2.67e-1],
+                   [-1.60e-1, -1.16e-2, -2.51e-1],
+                   [0, 1.06e-1, 8.62e-2]])),
+           # Aircraft control example II (Kautsky and Nichols 1983)
+           (array([[0, 1, 0, 0],
+                   [5.32e-7, -4.18e-1, -0.12, -2.32e-2],
+                   [-4.62e-9, 1, -0.752, -2.39e-2],
+                   [-5.61e-1, 0, 0.3, -1.74e-2]]),
+            array([[0, 0],
+                   [-1.72e-1, 7.45e-6],
+                   [-2.82e-2, -7.78e-5],
+                   [0, 3.69e-3]])),
+           # Symmetric example (Kautsky and Nichols 1983)
+           (array([[-3.624, 4.9567e-2, -2.4564e-1, 1.3853e-2],
+                   [3.3486e-1, -1.8875, -8.1251e-1, -2.8102e-1],
+                   [-1.9958e-1, -1.1335, -2.2039, -4.5523e-1],
+                   [1.3784e-1, -4.7140e-1, -3.3229e-1, -4.0605]]),
+            array([[2.3122e-1, 3.0761e-1, 3.6164e-1, 3.3217e-1],
+                   [8.8339e-1, 2.1460e-1, 5.6642e-1, 5.0153e-1]]).T),
+           # Ad-hoc ill-conditioned example (Byers and Nash 1989)
+           (array([[0, 0, 0, 0],
+                   [1, 10, 100, 1000],
+                   [0, 1, 10, 100],
+                   [0, 0, 1, 10]]),
+            array([[1, 0],
+                   [0, 1],
+                   [0, 0],
+                   [0, 0]]))
+        ]
+
+    # Return a generator
+    return (x for x in ABs)

harold/tests/test_system_props.py

@@ -22,9 +22,11 @@
 THE SOFTWARE.
 """
 
-from harold import Transfer, system_norm, transfer_to_state
+from harold import (Transfer, system_norm, transfer_to_state,
+                    controllability_indices)
+from test_static_ctrl_design import byersnash_A_B_test_pairs
 
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_equal
 from pytest import raises as assert_raises
 
 
@@ -43,3 +45,12 @@ def test_system_norm_simple():
     F = transfer_to_state(G)
     assert_almost_equal(system_norm(F), 1.1547, decimal=5)
     assert_almost_equal(system_norm(F, p=2), 2.2360679)
+
+
+def test_controllability_indices():
+    test_cont_ind = [[2, 2], [3, 2], [2, 2], [1, 2], [2, 3], [1, 3],
+                     [2, 3, 3], [2, 1, 1], [2, 2], [2, 2], [1, 3]]
+    example_gen = byersnash_A_B_test_pairs()
+    for ind, (A, B) in enumerate(example_gen):
+        indices = controllability_indices(A, B)
+        assert_equal(indices, test_cont_ind[ind])