|
| 1 | +# =================================================== |
| 2 | +# ab08n* tests |
| 3 | + |
| 4 | +import unittest |
| 5 | +from slycot import analysis |
| 6 | +import numpy as np |
| 7 | + |
| 8 | +from scipy import linalg |
| 9 | +from scipy import signal |
| 10 | +from numpy.testing import assert_equal, assert_allclose, assert_array_equal |
| 11 | + |
| 12 | +class test_ab13bd(unittest.TestCase): |
| 13 | + """ Test regular pencil construction ab08nX with input parameters |
| 14 | + according to example in documentation """ |
| 15 | + |
| 16 | + A = np.array([[0.0, 1.0],[-0.5, -0.1]]) |
| 17 | + B = np.array([[0.],[1.]]) |
| 18 | + C = np.eye(2) |
| 19 | + D = np.zeros((2,1)) |
| 20 | + |
| 21 | + Ad, Bd, Cd, Dd, dt = signal.cont2discrete((A, B, C, D), 0.1, method='zoh') |
| 22 | + |
| 23 | + def test_no_change_args_ccase(self): |
| 24 | + """ ab13md must not change its arguments. continuous system case. |
| 25 | + """ |
| 26 | + |
| 27 | + acopy = self.A.copy() |
| 28 | + bcopy = self.B.copy() |
| 29 | + ccopy = self.C.copy() |
| 30 | + dcopy = self.D.copy() |
| 31 | + |
| 32 | + dico = 'C' |
| 33 | + jobn = 'H' |
| 34 | + |
| 35 | + n, m = self.B.shape |
| 36 | + p = self.C.shape[0] |
| 37 | + |
| 38 | + analysis.ab13bd(dico, jobn, n, m, p, self.A, self.B, self.C, self.D) |
| 39 | + assert_array_equal(self.A, acopy) |
| 40 | + assert_array_equal(self.B, bcopy) |
| 41 | + assert_array_equal(self.C, ccopy) |
| 42 | + assert_array_equal(self.D, dcopy) |
| 43 | + |
| 44 | + def test_no_change_args_dcase(self): |
| 45 | + """ ab13md must not change its arguments. discrete system case. |
| 46 | + """ |
| 47 | + |
| 48 | + acopy = self.Ad.copy() |
| 49 | + bcopy = self.Bd.copy() |
| 50 | + ccopy = self.Cd.copy() |
| 51 | + dcopy = self.Dd.copy() |
| 52 | + |
| 53 | + dico = 'D' |
| 54 | + jobn = 'H' |
| 55 | + |
| 56 | + n, m = self.Bd.shape |
| 57 | + p = self.Cd.shape[0] |
| 58 | + |
| 59 | + analysis.ab13bd(dico, jobn, n, m, p, self.Ad, self.Bd, self.Cd, self.Dd) |
| 60 | + assert_array_equal(self.Ad, acopy) |
| 61 | + assert_array_equal(self.Bd, bcopy) |
| 62 | + assert_array_equal(self.Cd, ccopy) |
| 63 | + assert_array_equal(self.Dd, dcopy) |
| 64 | + |
| 65 | + def test_ab13bd_2norm_ccase(self): |
| 66 | + """ Compare ab13md with scipy solution (Lyapunov Equation). |
| 67 | + continuous system case. |
| 68 | + """ |
| 69 | + |
| 70 | + A = self.A |
| 71 | + B = self.B |
| 72 | + C = self.C |
| 73 | + D = self.D |
| 74 | + |
| 75 | + n, m = self.B.shape |
| 76 | + p = self.C.shape[0] |
| 77 | + |
| 78 | + dico = 'C' |
| 79 | + jobn = 'H' |
| 80 | + |
| 81 | + h2norm = analysis.ab13bd(dico, jobn, n, m, p, A, B, C, D) |
| 82 | + |
| 83 | + Lc = linalg.solve_continuous_lyapunov(A, -B@B.T) |
| 84 | + h2norm_Lc = np.sqrt(np.trace(C@Lc@C.T)) |
| 85 | + print(h2norm_Lc, h2norm) |
| 86 | + assert_allclose(h2norm_Lc, h2norm, atol=1e-5) |
| 87 | + |
| 88 | + Lo = linalg.solve_continuous_lyapunov(A.T, -C.T@C) |
| 89 | + h2norm_Lo = np.sqrt(np.trace(B.T@Lo@B)) |
| 90 | + print(h2norm_Lo, h2norm) |
| 91 | + assert_allclose(h2norm_Lo, h2norm, atol=1e-5) |
| 92 | + |
| 93 | + def test_ab13bd_2norm_dcase(self): |
| 94 | + """ Compare ab13md with scipy solution (Lyapunov Equation). |
| 95 | + discrete system case. |
| 96 | + """ |
| 97 | + |
| 98 | + Ad = self.Ad |
| 99 | + Bd = self.Bd |
| 100 | + Cd = self.Cd |
| 101 | + Dd = self.Dd |
| 102 | + |
| 103 | + n, m = Bd.shape |
| 104 | + p = Cd.shape[0] |
| 105 | + |
| 106 | + dico = 'D' |
| 107 | + jobn = 'H' |
| 108 | + |
| 109 | + h2norm = analysis.ab13bd(dico, jobn, n, m, p, Ad, Bd, Cd, Dd) |
| 110 | + |
| 111 | + Lc = linalg.solve_discrete_lyapunov(Ad, Bd@Bd.T) |
| 112 | + h2norm_Lc = np.sqrt(np.trace(Cd@Lc@Cd.T + Dd@Dd.T)) |
| 113 | + print(h2norm, h2norm_Lc) |
| 114 | + assert_allclose(h2norm_Lc, h2norm, atol=1e-5) |
| 115 | + |
| 116 | + Lo = linalg.solve_discrete_lyapunov(Ad.T, Cd.T@Cd) |
| 117 | + h2norm_Lo = np.sqrt(np.trace(Bd.T@Lo@Bd + Dd.T@Dd)) |
| 118 | + print(h2norm, h2norm_Lo) |
| 119 | + assert_allclose(h2norm_Lo, h2norm, atol=1e-5) |
| 120 | + |
| 121 | + |
| 122 | +if __name__ == "__main__": |
| 123 | + unittest.main() |
0 commit comments