|
81 | 81 | T = feedback(L, 1); |
82 | 82 |
|
83 | 83 | # Compute stability margins |
84 | | -#! Not yet implemented |
85 | | -# (gm, pm, wgc, wpc) = margin(L); |
| 84 | +(gm, pm, wgc, wpc) = margin(L); |
| 85 | +print("Gain margin: %g at %g" % (gm, wgc)) |
| 86 | +print("Phase margin: %g at %g" % (pm, wpc)) |
86 | 87 |
|
87 | | -#! TODO: this figure has something wrong; axis limits mismatch |
88 | 88 | figure(6); clf; |
89 | | -bode(L); |
| 89 | +bode(L, logspace(-4, 3)); |
90 | 90 |
|
91 | | -# Add crossover line |
92 | | -subplot(211); hold(True); |
93 | | -loglog([1e-4, 1e3], [1, 1], 'k-') |
| 91 | +# Add crossover line to the magnitude plot |
| 92 | +# |
| 93 | +# Note: in matplotlib before v2.1, the following code worked: |
| 94 | +# |
| 95 | +# subplot(211); hold(True); |
| 96 | +# loglog([1e-4, 1e3], [1, 1], 'k-') |
| 97 | +# |
| 98 | +# In later versions of matplotlib the call to subplot will clear the |
| 99 | +# axes and so we have to extract the axes that we want to use by hand. |
| 100 | +# In addition, hold() is deprecated so we no longer require it. |
| 101 | +# |
| 102 | +for ax in gcf().axes: |
| 103 | + if ax.get_label() == 'control-bode-magnitude': |
| 104 | + break |
| 105 | +ax.semilogx([1e-4, 1e3], 20 * np.log10([1, 1]), 'k-') |
94 | 106 |
|
| 107 | +# |
95 | 108 | # Replot phase starting at -90 degrees |
96 | | -bode(L, logspace(-4, 3)); |
| 109 | +# |
| 110 | +# Get the phase plot axes |
| 111 | +for ax in gcf().axes: |
| 112 | + if ax.get_label() == 'control-bode-phase': |
| 113 | + break |
| 114 | + |
| 115 | +# Recreate the frequency response and shift the phase |
97 | 116 | (mag, phase, w) = freqresp(L, logspace(-4, 3)); |
98 | 117 | phase = phase - 360; |
99 | | -subplot(212); |
100 | | -semilogx([1e-4, 1e3], [-180, -180], 'k-') |
101 | | -hold(True); |
102 | | -semilogx(w, np.squeeze(phase), 'b-') |
103 | | -axis([1e-4, 1e3, -360, 0]); |
| 118 | + |
| 119 | +# Replot the phase by hand |
| 120 | +ax.semilogx([1e-4, 1e3], [-180, -180], 'k-') |
| 121 | +ax.semilogx(w, np.squeeze(phase), 'b-') |
| 122 | +ax.axis([1e-4, 1e3, -360, 0]); |
104 | 123 | xlabel('Frequency [deg]'); ylabel('Phase [deg]'); |
105 | 124 | # set(gca, 'YTick', [-360, -270, -180, -90, 0]); |
106 | 125 | # set(gca, 'XTick', [10^-4, 10^-2, 1, 100]); |
|
109 | 128 | # Nyquist plot for complete design |
110 | 129 | # |
111 | 130 | figure(7); clf; |
112 | | -axis([-700, 5300, -3000, 3000]); hold(True); |
113 | 131 | nyquist(L, (0.0001, 1000)); |
114 | 132 | axis([-700, 5300, -3000, 3000]); |
115 | 133 |
|
|
118 | 136 |
|
119 | 137 | # Expanded region |
120 | 138 | figure(8); clf; subplot(231); |
121 | | -axis([-10, 5, -20, 20]); hold(True); |
122 | 139 | nyquist(L); |
123 | 140 | axis([-10, 5, -20, 20]); |
124 | 141 |
|
|
136 | 153 |
|
137 | 154 | figure(9); |
138 | 155 | (Yvec, Tvec) = step(T, linspace(0, 20)); |
139 | | -plot(Tvec.T, Yvec.T); hold(True); |
| 156 | +plot(Tvec.T, Yvec.T); |
140 | 157 |
|
141 | 158 | (Yvec, Tvec) = step(Co*S, linspace(0, 20)); |
142 | 159 | plot(Tvec.T, Yvec.T); |
|
0 commit comments