GeezerVison is showing github.com — Added hsvd and markov functionality and test script TestModelsimp.py. · python-control/python-control@d1b6b80
File treeExpand file treeCollapse file tree Open diff view settings Expand file treeCollapse file tree Open diff view settings Collapse file Copy file name to clipboard +28Lines changed: 28 additions & 0 deletions | Original file line number | Diff line number | Diff line change |
|---|
| | 1 | +#!/usr/bin/env python | | 2 | + | | 3 | +from modelsimp import * | | 4 | +from matlab import * | | 5 | +import numpy as np | | 6 | +import unittest | | 7 | + | | 8 | +class TestModelsimp(unittest.TestCase): | | 9 | + def testHSVD(self): | | 10 | + A = np.matrix("1. -2.; 3. -4.") | | 11 | + B = np.matrix("5.; 7.") | | 12 | + C = np.matrix("6. 8.") | | 13 | + D = np.matrix("9.") | | 14 | + sys = ss(A,B,C,D) | | 15 | + hsv = hsvd(sys) | | 16 | + hsvtrue = np.matrix("24.42686 0.5731395") | | 17 | + np.testing.assert_array_almost_equal(hsv, hsvtrue) | | 18 | + | | 19 | + def testMarkov(self): | | 20 | + U = np.matrix("1.; 1.; 1.; 1.; 1.") | | 21 | + Y = U | | 22 | + M = 3 | | 23 | + H = markov(Y,U,M) | | 24 | + Htrue = np.matrix("1.; 0.; 0.") | | 25 | + np.testing.assert_array_almost_equal( H, Htrue ) | | 26 | + | | 27 | +if __name__ == '__main__': | | 28 | + unittest.main() |
Collapse file Copy file name to clipboardExpand all lines: src/TestStatefbk.py +1-3Lines changed: 1 addition & 3 deletions | Original file line number | Diff line number | Diff line change |
|---|
@@ -48,8 +48,7 @@ def testGramWc(self): | 48 | 48 | B = np.matrix("5. 6.; 7. 8.") | 49 | 49 | C = np.matrix("4. 5.; 6. 7.") | 50 | 50 | D = np.matrix("13. 14.; 15. 16.") | 51 | | - # sys = ss(A, B, C, D) | 52 | | - sys = 1. | | 51 | + sys = ss(A, B, C, D) | 53 | 52 | Wctrue = np.matrix("18.5 24.5; 24.5 32.5") | 54 | 53 | Wc = gram(sys,'c') | 55 | 54 | np.testing.assert_array_almost_equal(Wc, Wctrue) | @@ -60,7 +59,6 @@ def testGramWo(self): | 60 | 59 | C = np.matrix("4. 5.; 6. 7.") | 61 | 60 | D = np.matrix("13. 14.; 15. 16.") | 62 | 61 | sys = ss(A, B, C, D) | 63 | | - sys = 1. | 64 | 62 | Wotrue = np.matrix("257.5 -94.5; -94.5 56.5") | 65 | 63 | Wo = gram(sys,'o') | 66 | 64 | np.testing.assert_array_almost_equal(Wo, Wotrue) | |
Collapse file Copy file name to clipboardExpand all lines: src/modelsimp.py +32-8Lines changed: 32 additions & 8 deletions | Original file line number | Diff line number | Diff line change |
|---|
@@ -53,6 +53,8 @@ def hsvd(sys): | 53 | 53 | ===== | 54 | 54 | H = hsvd(sys) | 55 | 55 |
| | 56 | + The Hankel singular values are the singular values of the Hankel operator. In practice, we compute the square root of the eigenvalues of the matrix formed by taking the product of the observability and controllability gramians. There are other (more efficient) methods based on solving the Lyapunov equation in a particular way (more details soon). | | 57 | +
| 56 | 58 | Inputs | 57 | 59 | ------ | 58 | 60 | sys : a state space system | @@ -63,15 +65,16 @@ def hsvd(sys): | 63 | 65 |
| 64 | 66 | """ | 65 | 67 |
| 66 | | - # Make sure that SLICOT is installed | 67 | | - try: | 68 | | - from slycot import sb01bd | 69 | | - except ImportError: | 70 | | - raise ControlSlycot("can't find slycot module 'sb01bd'") | 71 | | - | 72 | | - H = 1. | | 68 | + Wc = gram(sys,'c') | | 69 | + Wo = gram(sys,'o') | | 70 | + | | 71 | + WoWc = np.dot(Wo, Wc) | | 72 | + w, v = LA.eig(WoWc) | | 73 | + | | 74 | + hsv = np.sqrt(w) | | 75 | + hsv = np.matrix(hsv) | 73 | 76 | # Return the Hankel singular values | 74 | | - return H | | 77 | + return hsv | 75 | 78 |
| 76 | 79 | def era(YY,m,n,nin,nout,r): | 77 | 80 | """Calculate an ERA model of order r based on the impulse-response data YY | @@ -100,6 +103,7 @@ def markov(Y,U,M): | 100 | 103 | Usage | 101 | 104 | ===== | 102 | 105 | H = markov(Y,U,M) | | 106 | + Currently only works for SISO | 103 | 107 |
| 104 | 108 | Inputs | 105 | 109 | ------ | @@ -113,3 +117,23 @@ def markov(Y,U,M): | 113 | 117 |
| 114 | 118 | """ | 115 | 119 |
| | 120 | + # Convert input parameters to matrices (if they aren't already) | | 121 | + Ymat = np.mat(Y) | | 122 | + Umat = np.mat(U) | | 123 | + n = np.size(U) | | 124 | + | | 125 | + # Construct a matrix of control inputs to invert | | 126 | + UU = Umat | | 127 | + for i in range(1, M-1): | | 128 | + newCol = np.vstack((0, UU[0:n-1,i-2])) | | 129 | + UU = np.hstack((UU, newCol)) | | 130 | + Ulast = np.vstack((0, UU[0:n-1,M-2])) | | 131 | + for i in range(n-1,0,-1): | | 132 | + Ulast[i] = np.sum(Ulast[0:i-1]) | | 133 | + UU = np.hstack((UU, Ulast)) | | 134 | + | | 135 | + # Invert and solve for Markov parameters | | 136 | + H = UU.I | | 137 | + H = np.dot(H, Y) | | 138 | + | | 139 | + return H |
|