Skip to content

Commit ea67935

Browse files
martinupmartinup
committed
Added example as ipython notebook.
Optimizing bode-plot.
1 parent f2c40e6 commit ea67935

4 files changed

Lines changed: 853 additions & 171 deletions

File tree

‎.gitignore‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,3 +12,8 @@ build.log
1212
.coverage
1313
doc/_build
1414
doc/generated
15+
16+
examples/.ipynb_checkpoints/bode-plot-checkpoint.ipynb
17+
.settings/org.eclipse.core.resources.prefs
18+
.pydevproject
19+
.project

‎control/freqplot.py‎

Lines changed: 80 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -44,10 +44,8 @@
4444
import matplotlib.pyplot as plt
4545
import scipy as sp
4646
import numpy as np
47-
#from warnings import warn
4847
from .ctrlutil import unwrap
4948
from .bdalg import feedback
50-
#from .lti import isdtime, timebaseEqual
5149

5250
__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
5351
'bode', 'nyquist', 'gangof4']
@@ -132,14 +130,15 @@ def bode_plot(syslist, omega=None, omega_num=None, dB=None, Hz=None, deg=None,
132130
mags, phases, omegas, nyquistfrqs = [], [], [], []
133131
for sys in syslist:
134132
if (sys.inputs > 1 or sys.outputs > 1):
135-
#TODO: Add MIMO bode plots.
133+
# TODO: Add MIMO bode plots.
136134
raise NotImplementedError("Bode is currently only implemented for SISO systems.")
137135
else:
138136
omega_sys = np.array(omega)
139137
if sys.isdtime(True):
140-
nyquistfrq = 2. * np.pi * 1./sys.dt / 2.
141-
nyquistfrq_limit = nyquistfrq - (omega_sys[-1] - omega_sys[-2]) /2. # floating point!
142-
omega_sys = omega_sys[omega_sys < nyquistfrq_limit]
138+
nyquistfrq = 2. * np.pi * 1. / sys.dt / 2.
139+
omega_sys = omega_sys[omega_sys < nyquistfrq]
140+
nyquistfrq_limit = nyquistfrq - (omega_sys[-1] - omega_sys[-2]) / 2. # floating point!
141+
omega_sys = omega_sys[omega_sys < nyquistfrq_limit] # in two steps to get the right increment
143142
else:
144143
nyquistfrq = None
145144
# Get the magnitude and phase of the system
@@ -162,41 +161,57 @@ def bode_plot(syslist, omega=None, omega_num=None, dB=None, Hz=None, deg=None,
162161
omegas.append(omega_sys)
163162
nyquistfrqs.append(nyquistfrq)
164163
# Get the dimensions of the current axis, which we will divide up
165-
#! TODO: Not current implemented; just use subplot for now
164+
# ! TODO: Not current implemented; just use subplot for now
166165

167166
if (Plot):
168167
# Magnitude plot
169168
ax_mag = plt.subplot(211);
170169
if dB:
171-
pltline = plt.semilogx(omega_plot, 20 * np.log10(mag), *args, **kwargs)
170+
pltline = ax_mag.semilogx(omega_plot, 20 * np.log10(mag), *args, **kwargs)
172171
else:
173-
pltline = plt.loglog(omega_plot, mag, *args, **kwargs)
172+
pltline = ax_mag.loglog(omega_plot, mag, *args, **kwargs)
174173
plt.hold(True);
175174
if nyquistfrq_plot:
176-
plt.axvline(nyquistfrq_plot, color=pltline[0].get_color())
175+
ax_mag.axvline(nyquistfrq_plot, color=pltline[0].get_color())
176+
177177
# Add a grid to the plot + labeling
178-
plt.grid(True)
179-
plt.grid(True, which='minor')
180-
plt.ylabel("Magnitude (dB)" if dB else "Magnitude")
178+
ax_mag.grid(True, which='both')
179+
ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
181180

182181
# Phase plot
183-
plt.subplot(212, sharex=ax_mag);
182+
ax_phase = plt.subplot(212, sharex=ax_mag);
184183
if deg:
185184
phase_plot = phase * 180. / np.pi
186185
else:
187186
phase_plot = phase
188-
plt.semilogx(omega_plot, phase_plot, *args, **kwargs)
189-
plt.hold(True);
187+
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)
188+
ax_phase.hold(True);
190189
if nyquistfrq_plot:
191-
plt.axvline(nyquistfrq_plot, color=pltline[0].get_color())
190+
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())
191+
192192
# Add a grid to the plot + labeling
193-
plt.grid(True)
194-
plt.grid(True, which='minor')
195-
plt.ylabel("Phase (deg)" if deg else "Phase (rad)")
196-
193+
ax_phase.set_ylabel("Phase (deg)" if deg else "Phase (rad)")
194+
def genZeroCenteredSeries(val_min, val_max, period):
195+
v1 = np.ceil(val_min / period - 0.2)
196+
v2 = np.floor(val_max / period + 0.2)
197+
return np.arange(v1, v2 + 1) * period
198+
if deg:
199+
ylim = ax_phase.get_ylim()
200+
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 45.))
201+
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 15.), minor=True)
202+
else:
203+
ylim = ax_phase.get_ylim()
204+
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 4.))
205+
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 12.), minor=True)
206+
ax_phase.grid(True, which='both')
207+
# ax_mag.grid(which='minor', alpha=0.3)
208+
# ax_mag.grid(which='major', alpha=0.9)
209+
# ax_phase.grid(which='minor', alpha=0.3)
210+
# ax_phase.grid(which='major', alpha=0.9)
211+
197212
# Label the frequency axis
198-
plt.xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)")
199-
213+
ax_phase.set_xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)")
214+
200215
if len(syslist) == 1:
201216
return mags[0], phases[0], omegas[0]
202217
else:
@@ -242,19 +257,19 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
242257

243258
# Select a default range if none is provided
244259
if omega is None:
245-
#! TODO: think about doing something smarter for discrete
260+
# ! TODO: think about doing something smarter for discrete
246261
omega = default_frequency_range(syslist)
247262

248263
# Interpolate between wmin and wmax if a tuple or list are provided
249-
elif (isinstance(omega,list) | isinstance(omega,tuple)):
264+
elif (isinstance(omega, list) | isinstance(omega, tuple)):
250265
# Only accept tuple or list of length 2
251266
if (len(omega) != 2):
252267
raise ValueError("Supported frequency arguments are (wmin,wmax) tuple or list, or frequency vector. ")
253268
omega = np.logspace(np.log10(omega[0]), np.log10(omega[1]),
254269
num=50, endpoint=True, base=10.0)
255270
for sys in syslist:
256271
if (sys.inputs > 1 or sys.outputs > 1):
257-
#TODO: Add MIMO nyquist plots.
272+
# TODO: Add MIMO nyquist plots.
258273
raise NotImplementedError("Nyquist is currently only implemented for SISO systems.")
259274
else:
260275
# Get the magnitude and phase of the system
@@ -278,11 +293,11 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
278293
ind = slice(None, None, labelFreq)
279294
for xpt, ypt, omegapt in zip(x[ind], y[ind], omega[ind]):
280295
# Convert to Hz
281-
f = omegapt/(2*sp.pi)
296+
f = omegapt / (2 * sp.pi)
282297

283298
# Factor out multiples of 1000 and limit the
284299
# result to the range [-8, 8].
285-
pow1000 = max(min(get_pow1000(f),8),-8)
300+
pow1000 = max(min(get_pow1000(f), 8), -8)
286301

287302
# Get the SI prefix.
288303
prefix = gen_prefix(pow1000)
@@ -294,12 +309,12 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
294309
# instead of 1.0, and this would otherwise be
295310
# truncated to 0.
296311
plt.text(xpt, ypt,
297-
' ' + str(int(np.round(f/1000**pow1000, 0))) +
312+
' ' + str(int(np.round(f / 1000 ** pow1000, 0))) +
298313
' ' + prefix + 'Hz')
299314
return x, y, omega
300315

301316
# Gang of Four
302-
#! TODO: think about how (and whether) to handle lists of systems
317+
# ! TODO: think about how (and whether) to handle lists of systems
303318
def gangof4_plot(P, C, omega=None):
304319
"""Plot the "Gang of 4" transfer functions for a system
305320
@@ -317,34 +332,34 @@ def gangof4_plot(P, C, omega=None):
317332
-------
318333
None
319334
"""
320-
if (P.inputs > 1 or P.outputs > 1 or C.inputs > 1 or C.outputs >1):
321-
#TODO: Add MIMO go4 plots.
335+
if (P.inputs > 1 or P.outputs > 1 or C.inputs > 1 or C.outputs > 1):
336+
# TODO: Add MIMO go4 plots.
322337
raise NotImplementedError("Gang of four is currently only implemented for SISO systems.")
323338
else:
324339

325340
# Select a default range if none is provided
326-
#! TODO: This needs to be made more intelligent
341+
# ! TODO: This needs to be made more intelligent
327342
if omega is None:
328-
omega = default_frequency_range((P,C))
343+
omega = default_frequency_range((P, C))
329344

330345
# Compute the senstivity functions
331-
L = P*C;
346+
L = P * C;
332347
S = feedback(1, L);
333348
T = L * S;
334349

335350
# Plot the four sensitivity functions
336-
#! TODO: Need to add in the mag = 1 lines
351+
# ! TODO: Need to add in the mag = 1 lines
337352
mag_tmp, phase_tmp, omega = T.freqresp(omega);
338353
mag = np.squeeze(mag_tmp)
339354
phase = np.squeeze(phase_tmp)
340355
plt.subplot(221); plt.loglog(omega, mag);
341356

342-
mag_tmp, phase_tmp, omega = (P*S).freqresp(omega);
357+
mag_tmp, phase_tmp, omega = (P * S).freqresp(omega);
343358
mag = np.squeeze(mag_tmp)
344359
phase = np.squeeze(phase_tmp)
345360
plt.subplot(222); plt.loglog(omega, mag);
346361

347-
mag_tmp, phase_tmp, omega = (C*S).freqresp(omega);
362+
mag_tmp, phase_tmp, omega = (C * S).freqresp(omega);
348363
mag = np.squeeze(mag_tmp)
349364
phase = np.squeeze(phase_tmp)
350365
plt.subplot(223); plt.loglog(omega, mag);
@@ -418,14 +433,13 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe
418433
features_ = features_[features_ != 0.0];
419434
features = np.concatenate((features, features_))
420435
elif sys.isdtime(strict=True):
421-
freq_interesting.append(np.pi*1./sys.dt)
436+
freq_interesting.append(np.pi * 1. / sys.dt)
422437
p = sys.pole()
423438
p = p[p != -1.]
424-
features = np.concatenate((features, np.abs(np.log(p)/sys.dt)))
439+
features = np.concatenate((features, np.abs(np.log(p) / sys.dt)))
425440
z = sys.zero()
426-
z = z[(z.imag != 0.0)]
427-
#z = z[(z.imag != 0.0) | (z.real > 0.0)]
428-
features = np.concatenate((features, np.abs(np.log(z)/sys.dt)))
441+
z = z[(z.imag != 0.0)] # Get rid of poles and zeros at the origin and real <= 0 & imag==0
442+
features = np.concatenate((features, np.abs(np.log(z) / sys.dt)))
429443
else:
430444
pass
431445
# TODO:
@@ -440,20 +454,19 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe
440454
if Hz:
441455
features /= 2.*np.pi
442456
features = np.log10(features)
443-
lsp_min = np.floor(np.min(features)) - feature_periphery_decade
444-
lsp_max = np.ceil(np.max(features)) + feature_periphery_decade
457+
lsp_min = np.floor(np.min(features) - feature_periphery_decade)
458+
lsp_max = np.ceil(np.max(features) + feature_periphery_decade)
445459
lsp_min += np.log10(2.*np.pi)
446460
lsp_max += np.log10(2.*np.pi)
447461
else:
448-
# Take the log of the features
449462
features = np.log10(features)
450-
lsp_min = np.floor(np.min(features)) - feature_periphery_decade
451-
lsp_max = np.ceil(np.max(features)) + feature_periphery_decade
463+
lsp_min = np.floor(np.min(features) - feature_periphery_decade)
464+
lsp_max = np.ceil(np.max(features) + feature_periphery_decade)
452465
if freq_interesting:
453466
lsp_min = min(lsp_min, np.log10(min(freq_interesting)))
454467
lsp_max = max(lsp_max, np.log10(max(freq_interesting)))
455468

456-
#! TODO: Add a check in discrete case to make sure we don't get aliasing (Attention: there is a list of system but only one omega vector)
469+
# ! TODO: Add a check in discrete case to make sure we don't get aliasing (Attention: there is a list of system but only one omega vector)
457470

458471
# Set the range to be an order of magnitude beyond any features
459472
if number_of_samples:
@@ -478,7 +491,7 @@ def get_pow1000(num):
478491
return 0
479492
elif dnum < 0:
480493
dnum = -dnum
481-
return int(floor(dnum.log10()/3))
494+
return int(floor(dnum.log10() / 3))
482495

483496
def gen_prefix(pow1000):
484497
'''Return the SI prefix for a power of 1000.
@@ -487,23 +500,23 @@ def gen_prefix(pow1000):
487500
# deca, deci, and centi).
488501
if pow1000 < -8 or pow1000 > 8:
489502
raise ValueError("Value is out of the range covered by the SI prefixes.")
490-
return ['Y', # yotta (10^24)
491-
'Z', # zetta (10^21)
492-
'E', # exa (10^18)
493-
'P', # peta (10^15)
494-
'T', # tera (10^12)
495-
'G', # giga (10^9)
496-
'M', # mega (10^6)
497-
'k', # kilo (10^3)
498-
'', # (10^0)
499-
'm', # milli (10^-3)
500-
r'$\mu$', # micro (10^-6)
501-
'n', # nano (10^-9)
502-
'p', # pico (10^-12)
503-
'f', # femto (10^-15)
504-
'a', # atto (10^-18)
505-
'z', # zepto (10^-21)
506-
'y'][8 - pow1000] # yocto (10^-24)
503+
return ['Y', # yotta (10^24)
504+
'Z', # zetta (10^21)
505+
'E', # exa (10^18)
506+
'P', # peta (10^15)
507+
'T', # tera (10^12)
508+
'G', # giga (10^9)
509+
'M', # mega (10^6)
510+
'k', # kilo (10^3)
511+
'', # (10^0)
512+
'm', # milli (10^-3)
513+
r'$\mu$', # micro (10^-6)
514+
'n', # nano (10^-9)
515+
'p', # pico (10^-12)
516+
'f', # femto (10^-15)
517+
'a', # atto (10^-18)
518+
'z', # zepto (10^-21)
519+
'y'][8 - pow1000] # yocto (10^-24)
507520

508521
# Function aliases
509522
bode = bode_plot

0 commit comments

Comments
 (0)