Skip to content

Commit f201a91

Browse files
author
Kevin Chen
committed
Added rss-balred.py to /examples and bode now returns mag, phase, freq, with the option to supress plotting.
Lauren Padilla <lpadilla@princeton.edu>
1 parent c95353a commit f201a91

5 files changed

Lines changed: 96 additions & 46 deletions

File tree

‎examples/rss-balred.py‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
#!/usr/bin/env python
2+
3+
import numpy as np
4+
import control.modelsimp as msimp
5+
import control.matlab as mt
6+
from control.statesp import StateSpace
7+
import matplotlib.pyplot as plt
8+
9+
plt.close('all')
10+
11+
#controlable canonical realization computed in matlab for the transfer function:
12+
# num = [1 11 45 32], den = [1 15 60 200 60]
13+
A = np.matrix('-15., -7.5, -6.25, -1.875; \
14+
8., 0., 0., 0.; \
15+
0., 4., 0., 0.; \
16+
0., 0., 1., 0.')
17+
B = np.matrix('2.; 0.; 0.; 0.')
18+
C = np.matrix('0.5, 0.6875, 0.7031, 0.5')
19+
D = np.matrix('0.')
20+
21+
# The full system
22+
fsys = StateSpace(A,B,C,D)
23+
# The reduced system, truncating the order by 1
24+
ord = 3
25+
rsys = msimp.balred(fsys,ord, method = 'truncate')
26+
27+
# Comparison of the step responses of the full and reduced systems
28+
plt.figure(1)
29+
t, y = mt.step(fsys)
30+
tr, yr = mt.step(rsys)
31+
plt.plot(t.T, y.T)
32+
plt.hold(True)
33+
plt.plot(tr.T, yr.T)
34+
35+
# Repeat balanced reduction, now with 100-dimensional random state space
36+
sysrand = mt.rss(100, 1, 1)
37+
rsysrand = msimp.balred(sysrand,10,method ='truncate')
38+
39+
# Comparison of the impulse responses of the full and reduced random systems
40+
plt.figure(2)
41+
trand, yrand = mt.impulse(sysrand)
42+
trandr, yrandr = mt.impulse(rsysrand)
43+
plt.plot(trand.T, yrand.T,trandr.T, yrandr.T)
44+

‎examples/secord-matlab.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222

2323
# Bode plot for the system
2424
figure(2)
25-
bode(sys, logspace(-2, 2))
25+
mag,phase,om = bode(sys, logspace(-2, 2),Plot=True)
2626

2727
# Nyquist plot for the system
2828
figure(3)

‎src/freqplot.py‎

Lines changed: 47 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -54,12 +54,12 @@
5454
#
5555

5656
# Bode plot
57-
def bode(syslist, omega=None, dB=False, Hz=False, color=None):
57+
def bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True):
5858
"""Bode plot for a system
5959
6060
Usage
6161
=====
62-
(magh, phaseh) = bode(syslist, omega=None, dB=False, Hz=False)
62+
(magh, phaseh, omega) = bode(syslist, omega=None, dB=False, Hz=False, color=None, Plot=True)
6363
6464
Plots a Bode plot for the system over a (optional) frequency range.
6565
@@ -73,16 +73,18 @@ def bode(syslist, omega=None, dB=False, Hz=False, color=None):
7373
If True, plot result in dB
7474
Hz : boolean
7575
If True, plot frequency in Hz (omega must be provided in rad/sec)
76+
Plot : boolean
77+
If True, plot magnitude and phase
7678
7779
Return values
7880
-------------
79-
magh : graphics handle to magnitude plot (for rescaling, etc)
80-
phaseh : graphics handle to phase plot
81+
magh : magnitude array
82+
phaseh : phase array
83+
omega : frequency array
8184
8285
Notes
8386
-----
84-
1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the
85-
frequency response for a system.
87+
1. Alternatively, may use (mag, phase, freq) = sys.freqresp(freq) to generate the frequency response for a system.
8688
"""
8789
# If argument was a singleton, turn it into a list
8890
if (not getattr(syslist, '__iter__', False)):
@@ -108,46 +110,47 @@ def bode(syslist, omega=None, dB=False, Hz=False, color=None):
108110
# Get the dimensions of the current axis, which we will divide up
109111
#! TODO: Not current implemented; just use subplot for now
110112

111-
# Magnitude plot
112-
plt.subplot(211);
113-
if dB:
114-
if color==None:
115-
plt.semilogx(omega, mag)
113+
if (Plot):
114+
# Magnitude plot
115+
plt.subplot(211);
116+
if dB:
117+
if color==None:
118+
plt.semilogx(omega, mag)
119+
else:
120+
plt.semilogx(omega, mag, color=color)
121+
plt.ylabel("Magnitude (dB)")
116122
else:
117-
plt.semilogx(omega, mag, color=color)
118-
plt.ylabel("Magnitude (dB)")
119-
else:
123+
if color==None:
124+
plt.loglog(omega, mag)
125+
else:
126+
plt.loglog(omega, mag, color=color)
127+
plt.ylabel("Magnitude")
128+
129+
# Add a grid to the plot
130+
plt.grid(True)
131+
plt.grid(True, which='minor')
132+
plt.hold(True);
133+
134+
# Phase plot
135+
plt.subplot(212);
120136
if color==None:
121-
plt.loglog(omega, mag)
122-
else:
123-
plt.loglog(omega, mag, color=color)
124-
plt.ylabel("Magnitude")
125-
126-
# Add a grid to the plot
127-
plt.grid(True)
128-
plt.grid(True, which='minor')
129-
plt.hold(True);
130-
131-
# Phase plot
132-
plt.subplot(212);
133-
if color==None:
134-
plt.semilogx(omega, phase)
135-
else:
136-
plt.semilogx(omega, phase, color=color)
137-
plt.hold(True)
138-
139-
# Add a grid to the plot
140-
plt.grid(True)
141-
plt.grid(True, which='minor')
142-
plt.ylabel("Phase (deg)")
143-
144-
# Label the frequency axis
145-
if Hz:
146-
plt.xlabel("Frequency (Hz)")
147-
else:
148-
plt.xlabel("Frequency (rad/sec)")
149-
150-
return (211, 212)
137+
plt.semilogx(omega, phase)
138+
else:
139+
plt.semilogx(omega, phase, color=color)
140+
plt.hold(True)
141+
142+
# Add a grid to the plot
143+
plt.grid(True)
144+
plt.grid(True, which='minor')
145+
plt.ylabel("Phase (deg)")
146+
147+
# Label the frequency axis
148+
if Hz:
149+
plt.xlabel("Frequency (Hz)")
150+
else:
151+
plt.xlabel("Frequency (rad/sec)")
152+
153+
return mag, phase, omega
151154

152155
# Nyquist plot
153156
def nyquist(syslist, omega=None):

‎src/modelsimp.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,8 @@ def balred(sys,orders,method='truncate'):
207207

208208
#Check system is stable
209209
D,V = np.linalg.eig(sys.A)
210+
print D.shape
211+
print D
210212
for e in D:
211213
if e.real >= 0:
212214
raise ValueError, "Oops, the system is unstable!"

‎src/xferfcn.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -710,7 +710,8 @@ def _convertToTransferFunction(sys, **kw):
710710
num[i][j] = list(tfout[6][i, j, :])
711711
# Each transfer function matrix row has a common denominator.
712712
den[i][j] = list(tfout[5][i, :])
713-
713+
print num
714+
print den
714715
return TransferFunction(num, den)
715716
elif isinstance(sys, (int, long, float, complex)):
716717
if "inputs" in kw:

0 commit comments

Comments
 (0)