Skip to content

Commit 2b6e8e6

Browse files
author
Kevin Chen
committed
Added modelsimp.py with interfaces for hsvd, era, and markov functions. Also, modified TestStatefbk.py so that we import numpy as np, as in other files.
Steven Brunton <sbrunton@princeton.edu>
1 parent 49bdd94 commit 2b6e8e6

2 files changed

Lines changed: 150 additions & 35 deletions

File tree

‎src/TestStatefbk.py‎

Lines changed: 35 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -2,68 +2,68 @@
22

33
from statefbk import ctrb, obsv, place, lqr, gram
44
from matlab import *
5-
import numpy as N
5+
import numpy as np
66
import unittest
77

88
class TestStatefbk(unittest.TestCase):
99
def testCtrbSISO(self):
10-
A = N.matrix("1. 2.; 3. 4.")
11-
B = N.matrix("5.; 7.")
12-
Wctrue = N.matrix("5. 19.; 7. 43.")
10+
A = np.matrix("1. 2.; 3. 4.")
11+
B = np.matrix("5.; 7.")
12+
Wctrue = np.matrix("5. 19.; 7. 43.")
1313
Wc = ctrb(A,B)
14-
N.testing.assert_array_almost_equal(Wc, Wctrue)
14+
np.testing.assert_array_almost_equal(Wc, Wctrue)
1515

1616
def testCtrbMIMO(self):
17-
A = N.matrix("1. 2.; 3. 4.")
18-
B = N.matrix("5. 6.; 7. 8.")
19-
Wctrue = N.matrix("5. 6. 19. 22.; 7. 8. 43. 50.")
17+
A = np.matrix("1. 2.; 3. 4.")
18+
B = np.matrix("5. 6.; 7. 8.")
19+
Wctrue = np.matrix("5. 6. 19. 22.; 7. 8. 43. 50.")
2020
Wc = ctrb(A,B)
21-
N.testing.assert_array_almost_equal(Wc, Wctrue)
21+
np.testing.assert_array_almost_equal(Wc, Wctrue)
2222

2323
def testObsvSISO(self):
24-
A = N.matrix("1. 2.; 3. 4.")
25-
C = N.matrix("5. 7.")
26-
Wotrue = N.matrix("5. 7.; 26. 38.")
24+
A = np.matrix("1. 2.; 3. 4.")
25+
C = np.matrix("5. 7.")
26+
Wotrue = np.matrix("5. 7.; 26. 38.")
2727
Wo = obsv(A,C)
28-
N.testing.assert_array_almost_equal(Wo, Wotrue)
28+
np.testing.assert_array_almost_equal(Wo, Wotrue)
2929

3030
def testObsvMIMO(self):
31-
A = N.matrix("1. 2.; 3. 4.")
32-
C = N.matrix("5. 6.; 7. 8.")
33-
Wotrue = N.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.")
31+
A = np.matrix("1. 2.; 3. 4.")
32+
C = np.matrix("5. 6.; 7. 8.")
33+
Wotrue = np.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.")
3434
Wo = obsv(A,C)
35-
N.testing.assert_array_almost_equal(Wo, Wotrue)
35+
np.testing.assert_array_almost_equal(Wo, Wotrue)
3636

3737
def testCtrbObsvDuality(self):
38-
A = N.matrix("1.2 -2.3; 3.4 -4.5")
39-
B = N.matrix("5.8 6.9; 8. 9.1")
38+
A = np.matrix("1.2 -2.3; 3.4 -4.5")
39+
B = np.matrix("5.8 6.9; 8. 9.1")
4040
Wc = ctrb(A,B);
41-
A = N.transpose(A)
42-
C = N.transpose(B)
43-
Wo = N.transpose(obsv(A,C));
44-
N.testing.assert_array_almost_equal(Wc,Wo)
41+
A = np.transpose(A)
42+
C = np.transpose(B)
43+
Wo = np.transpose(obsv(A,C));
44+
np.testing.assert_array_almost_equal(Wc,Wo)
4545

4646
def testGramWc(self):
47-
A = N.matrix("1. -2.; 3. -4.")
48-
B = N.matrix("5. 6.; 7. 8.")
49-
C = N.matrix("4. 5.; 6. 7.")
50-
D = N.matrix("13. 14.; 15. 16.")
47+
A = np.matrix("1. -2.; 3. -4.")
48+
B = np.matrix("5. 6.; 7. 8.")
49+
C = np.matrix("4. 5.; 6. 7.")
50+
D = np.matrix("13. 14.; 15. 16.")
5151
# sys = ss(A, B, C, D)
5252
sys = 1.
53-
Wctrue = N.matrix("18.5 24.5; 24.5 32.5")
53+
Wctrue = np.matrix("18.5 24.5; 24.5 32.5")
5454
Wc = gram(sys,'c')
55-
N.testing.assert_array_almost_equal(Wc, Wctrue)
55+
np.testing.assert_array_almost_equal(Wc, Wctrue)
5656

5757
def testGramWo(self):
58-
A = N.matrix("1. -2.; 3. -4.")
59-
B = N.matrix("5. 6.; 7. 8.")
60-
C = N.matrix("4. 5.; 6. 7.")
61-
D = N.matrix("13. 14.; 15. 16.")
58+
A = np.matrix("1. -2.; 3. -4.")
59+
B = np.matrix("5. 6.; 7. 8.")
60+
C = np.matrix("4. 5.; 6. 7.")
61+
D = np.matrix("13. 14.; 15. 16.")
6262
sys = ss(A, B, C, D)
6363
sys = 1.
64-
Wotrue = N.matrix("257.5 -94.5; -94.5 56.5")
64+
Wotrue = np.matrix("257.5 -94.5; -94.5 56.5")
6565
Wo = gram(sys,'o')
66-
N.testing.assert_array_almost_equal(Wo, Wotrue)
66+
np.testing.assert_array_almost_equal(Wo, Wotrue)
6767

6868
if __name__ == '__main__':
6969
unittest.main()

‎src/modelsimp.py‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
# modelsimp.py - tools for model simplification
2+
#
3+
# Author: Steve Brunton, Kevin Chen, Lauren Padilla
4+
# Date: 30 Nov 2010
5+
#
6+
# This file contains routines for obtaining reduced order models
7+
#
8+
# Copyright (c) 2010 by California Institute of Technology
9+
# All rights reserved.
10+
#
11+
# Redistribution and use in source and binary forms, with or without
12+
# modification, are permitted provided that the following conditions
13+
# are met:
14+
#
15+
# 1. Redistributions of source code must retain the above copyright
16+
# notice, this list of conditions and the following disclaimer.
17+
#
18+
# 2. Redistributions in binary form must reproduce the above copyright
19+
# notice, this list of conditions and the following disclaimer in the
20+
# documentation and/or other materials provided with the distribution.
21+
#
22+
# 3. Neither the name of the California Institute of Technology nor
23+
# the names of its contributors may be used to endorse or promote
24+
# products derived from this software without specific prior
25+
# written permission.
26+
#
27+
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
28+
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
29+
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
30+
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
31+
# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
32+
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
33+
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
34+
# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
35+
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
36+
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
37+
# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
38+
# SUCH DAMAGE.
39+
#
40+
# $Id$
41+
42+
# External packages and modules
43+
import numpy as np
44+
import ctrlutil
45+
from control.exception import *
46+
47+
# Hankel Singular Value Decomposition
48+
# The following returns the Hankel singular values, which are singular values of the matrix formed by multiplying the controllability and observability grammians
49+
def hsvd(sys):
50+
"""Calculate the Hankel singular values
51+
52+
Usage
53+
=====
54+
H = hsvd(sys)
55+
56+
Inputs
57+
------
58+
sys : a state space system
59+
60+
Outputs
61+
-------
62+
H : a list of Hankel singular values
63+
64+
"""
65+
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.
73+
# Return the Hankel singular values
74+
return H
75+
76+
def era(YY,m,n,nin,nout,r):
77+
"""Calculate an ERA model of order r based on the impulse-response data YY
78+
79+
Usage
80+
=====
81+
sys = era(YY,m,n,nin,nout,r)
82+
83+
Inputs
84+
------
85+
YY : nout x nin dimensional impulse-response data
86+
m : number of rows in Hankel matrix
87+
n : number of columns in Hankel matrix
88+
nin : number of input variables
89+
nout : number of output variables
90+
r : order of model
91+
92+
Outputs
93+
-------
94+
sys : a reduced order model sys=ss(Ar,Br,Cr,Dr)
95+
96+
"""
97+
def markov(Y,U,M):
98+
"""Calculate the first M Markov parameters [D CB CAB ...] from input U, output Y
99+
100+
Usage
101+
=====
102+
H = markov(Y,U,M)
103+
104+
Inputs
105+
------
106+
Y : output data
107+
U : input data
108+
M : number of Markov parameters to output
109+
110+
Outputs
111+
-------
112+
H : first M Markov parameters
113+
114+
"""
115+

0 commit comments

Comments
 (0)