From 598a91aa0d662467734ddfa2f0042f3b27145de2 Mon Sep 17 00:00:00 2001 From: ipa-lth Date: Wed, 18 Apr 2018 10:33:56 +0200 Subject: [PATCH 1/2] Fix canonical form 'observable' --- control/canonical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/canonical.py b/control/canonical.py index 5eacce372..b5a34ef94 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -133,7 +133,7 @@ def observable_form(xsys): Wrz = obsv(zsys.A, zsys.C) # Transformation from one form to another - Tzx = inv(Wrz) * Wrx + Tzx = solve(Wrz, Wrx) # matrix right division, Tzx = inv(Wrz) * Wrx # Finally, compute the output matrix zsys.B = Tzx * xsys.B From 5009e9293e28b3cc806647c63422ab81acdf45dc Mon Sep 17 00:00:00 2001 From: ipa-lth Date: Mon, 23 Apr 2018 10:54:53 +0200 Subject: [PATCH 2/2] add: test for singularity; unittests --- control/canonical.py | 5 +++- control/tests/canonical_test.py | 44 +++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/control/canonical.py b/control/canonical.py index b5a34ef94..4acd78ff8 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -133,7 +133,10 @@ def observable_form(xsys): Wrz = obsv(zsys.A, zsys.C) # Transformation from one form to another - Tzx = solve(Wrz, Wrx) # matrix right division, Tzx = inv(Wrz) * Wrx + Tzx = solve(Wrz, Wrx) # matrix left division, Tzx = inv(Wrz) * Wrx + + if matrix_rank(Tzx) != xsys.states: + raise ValueError("Transformation matrix singular to working precision.") # Finally, compute the output matrix zsys.B = Tzx * xsys.B diff --git a/control/tests/canonical_test.py b/control/tests/canonical_test.py index f5908a8f4..a450453c1 100644 --- a/control/tests/canonical_test.py +++ b/control/tests/canonical_test.py @@ -52,3 +52,47 @@ def test_unreachable_system(self): # Check if an exception is raised np.testing.assert_raises(ValueError, canonical_form, sys, "reachable") + + def test_observable_form(self): + """Test the observable canonical form""" + + # Create a system in the observable canonical form + coeffs = [1.0, 2.0, 3.0, 4.0, 1.0] + A_true = np.polynomial.polynomial.polycompanion(coeffs) + A_true = np.fliplr(np.flipud(A_true)) + B_true = np.matrix("1.0 1.0 1.0 1.0").T + C_true = np.matrix("1.0 0.0 0.0 0.0") + D_true = 42.0 + + # Perform a coordinate transform with a random invertible matrix + T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + [-0.74855725, -0.39136285, -0.18142339, -0.50356997], + [-0.40688007, 0.81416369, 0.38002113, -0.16483334], + [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) + A = np.linalg.solve(T_true, A_true)*T_true + B = np.linalg.solve(T_true, B_true) + C = C_true*T_true + D = D_true + + # Create a state space system and convert it to the observable canonical form + sys_check, T_check = canonical_form(ss(A, B, C, D), "observable") + + # Check against the true values + np.testing.assert_array_almost_equal(sys_check.A, A_true) + np.testing.assert_array_almost_equal(sys_check.B, B_true) + np.testing.assert_array_almost_equal(sys_check.C, C_true) + np.testing.assert_array_almost_equal(sys_check.D, D_true) + np.testing.assert_array_almost_equal(T_check, T_true) + + def test_unobservable_system(self): + """Test observable canonical form with an unobservable system""" + + # Create an unobservable system + A = np.matrix("1.0 2.0 2.0; 4.0 5.0 5.0; 7.0 8.0 8.0") + B = np.matrix("1.0 1.0 1.0").T + C = np.matrix("1.0 1.0 1.0") + D = 42.0 + sys = ss(A, B, C, D) + + # Check if an exception is raised + np.testing.assert_raises(ValueError, canonical_form, sys, "observable")