|
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