|
| 1 | +""" pyplots.ai |
| 2 | +nyquist-basic: Nyquist Plot for Control Systems |
| 3 | +Library: plotnine 0.15.3 | Python 3.14.3 |
| 4 | +Quality: 90/100 | Created: 2026-03-20 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pandas as pd |
| 9 | +from plotnine import ( |
| 10 | + aes, |
| 11 | + annotate, |
| 12 | + arrow, |
| 13 | + coord_fixed, |
| 14 | + element_blank, |
| 15 | + element_line, |
| 16 | + element_rect, |
| 17 | + element_text, |
| 18 | + geom_path, |
| 19 | + geom_point, |
| 20 | + geom_ribbon, |
| 21 | + geom_segment, |
| 22 | + ggplot, |
| 23 | + guide_legend, |
| 24 | + guides, |
| 25 | + labs, |
| 26 | + scale_alpha_manual, |
| 27 | + scale_color_manual, |
| 28 | + scale_linetype_manual, |
| 29 | + scale_size_manual, |
| 30 | + scale_x_continuous, |
| 31 | + scale_y_continuous, |
| 32 | + theme, |
| 33 | + theme_minimal, |
| 34 | +) |
| 35 | + |
| 36 | + |
| 37 | +# Data - Transfer function G(s) = 5 / [(s+1)(0.5s+1)(0.2s+1)] |
| 38 | +# Poles at s = -1, -2, -5 (stable minimum-phase system) |
| 39 | +omega = np.concatenate( |
| 40 | + [np.logspace(-2, -0.5, 150), np.logspace(-0.5, 0.5, 300), np.logspace(0.5, 1.5, 200), np.logspace(1.5, 3, 100)] |
| 41 | +) |
| 42 | + |
| 43 | +K = 5.0 |
| 44 | +jw = 1j * omega |
| 45 | +G = K / ((jw + 1) * (0.5 * jw + 1) * (0.2 * jw + 1)) |
| 46 | + |
| 47 | +real_part = G.real |
| 48 | +imag_part = G.imag |
| 49 | + |
| 50 | +# Build unified DataFrame for both curves |
| 51 | +df_pos = pd.DataFrame({"real": real_part, "imaginary": imag_part, "curve": "Positive freq (ω > 0)"}) |
| 52 | +df_neg = pd.DataFrame({"real": real_part[::-1], "imaginary": -imag_part[::-1], "curve": "Negative freq (ω < 0)"}) |
| 53 | +df_curves = pd.concat([df_pos, df_neg], ignore_index=True) |
| 54 | + |
| 55 | +# Unit circle reference |
| 56 | +theta = np.linspace(0, 2 * np.pi, 200) |
| 57 | +unit_circle_df = pd.DataFrame({"real": np.cos(theta), "imaginary": np.sin(theta)}) |
| 58 | + |
| 59 | +# Stability region shading |
| 60 | +theta_fill = np.linspace(0, 2 * np.pi, 100) |
| 61 | +cos_vals = np.cos(theta_fill) |
| 62 | +sin_vals = np.sin(np.arccos(np.clip(cos_vals, -1, 1))) |
| 63 | +stability_df = pd.DataFrame({"x": cos_vals, "ymin": -sin_vals, "ymax": sin_vals}) |
| 64 | + |
| 65 | +# Critical point (-1, 0) |
| 66 | +critical_df = pd.DataFrame({"real": [-1.0], "imaginary": [0.0]}) |
| 67 | + |
| 68 | +# Frequency annotations - positioned away from crossover markers near origin |
| 69 | +# Use explicit offsets to guarantee no overlap with critical point / crossover labels |
| 70 | +freq_annotations = [ |
| 71 | + (0.1, 0.4, -0.5), # upper right region - offset down-right |
| 72 | + (0.3, 0.6, -0.55), # right side - offset down-right |
| 73 | + (0.7, 0.5, -0.55), # bottom region - offset down |
| 74 | + (1.5, -0.7, 0.5), # lower-left - offset up-left (away from gain crossover label) |
| 75 | + (3.0, 0.75, 0.8), # near origin - offset up-right, well clear of phase crossover marker |
| 76 | +] |
| 77 | +annot_rows = [] |
| 78 | +for wf, ox, oy in freq_annotations: |
| 79 | + idx = np.argmin(np.abs(omega - wf)) |
| 80 | + rx, ry = real_part[idx], imag_part[idx] |
| 81 | + annot_rows.append({"real": rx, "imaginary": ry, "label": f"ω = {wf:g}", "lx": rx + ox, "ly": ry + oy}) |
| 82 | + |
| 83 | +annot_df = pd.DataFrame(annot_rows) |
| 84 | + |
| 85 | +# Direction arrows along positive frequency curve |
| 86 | +arrow_data = [] |
| 87 | +for frac in [0.06, 0.3, 0.65]: |
| 88 | + idx = int(frac * len(df_pos)) |
| 89 | + step = max(2, int(0.015 * len(df_pos))) |
| 90 | + if 0 < idx < len(df_pos) - step: |
| 91 | + arrow_data.append( |
| 92 | + { |
| 93 | + "x": df_pos.iloc[idx]["real"], |
| 94 | + "y": df_pos.iloc[idx]["imaginary"], |
| 95 | + "xend": df_pos.iloc[idx + step]["real"], |
| 96 | + "yend": df_pos.iloc[idx + step]["imaginary"], |
| 97 | + } |
| 98 | + ) |
| 99 | + |
| 100 | +# Direction arrows along negative frequency curve |
| 101 | +for frac in [0.35, 0.75]: |
| 102 | + idx = int(frac * len(df_neg)) |
| 103 | + step = max(2, int(0.015 * len(df_neg))) |
| 104 | + if 0 < idx < len(df_neg) - step: |
| 105 | + arrow_data.append( |
| 106 | + { |
| 107 | + "x": df_neg.iloc[idx]["real"], |
| 108 | + "y": df_neg.iloc[idx]["imaginary"], |
| 109 | + "xend": df_neg.iloc[idx + step]["real"], |
| 110 | + "yend": df_neg.iloc[idx + step]["imaginary"], |
| 111 | + } |
| 112 | + ) |
| 113 | + |
| 114 | +arrow_df = pd.DataFrame(arrow_data) |
| 115 | + |
| 116 | +# Gain and phase margin calculations |
| 117 | +mag = np.abs(G) |
| 118 | +phase = np.degrees(np.angle(G)) |
| 119 | + |
| 120 | +gc_idx = np.argmin(np.abs(mag - 1.0)) |
| 121 | +gc_omega = omega[gc_idx] |
| 122 | +phase_margin = 180 + phase[gc_idx] |
| 123 | + |
| 124 | +pc_idx = np.argmin(np.abs(phase + 180)) |
| 125 | +pc_omega = omega[pc_idx] |
| 126 | +gain_margin_db = -20 * np.log10(mag[pc_idx]) |
| 127 | + |
| 128 | +# Segment from origin to gain crossover point (shows phase margin angle) |
| 129 | +gc_seg_df = pd.DataFrame({"x": [0.0], "y": [0.0], "xend": [G[gc_idx].real], "yend": [G[gc_idx].imag]}) |
| 130 | + |
| 131 | +# Plot |
| 132 | +plot = ( |
| 133 | + ggplot() |
| 134 | + # Stability region shading |
| 135 | + + geom_ribbon(stability_df, aes(x="x", ymin="ymin", ymax="ymax"), fill="#E8EDF2", alpha=0.3) |
| 136 | + # Unit circle |
| 137 | + + geom_path(unit_circle_df, aes(x="real", y="imaginary"), color="#BBBBBB", size=0.8, linetype="dashed") |
| 138 | + # Nyquist curves with unified legend via aes-mapped scales |
| 139 | + + geom_path(df_curves, aes(x="real", y="imaginary", color="curve", linetype="curve", size="curve", alpha="curve")) |
| 140 | + + scale_color_manual( |
| 141 | + values={"Positive freq (ω > 0)": "#306998", "Negative freq (ω < 0)": "#7BA4C7"}, name="Frequency Response" |
| 142 | + ) |
| 143 | + + scale_linetype_manual( |
| 144 | + values={"Positive freq (ω > 0)": "solid", "Negative freq (ω < 0)": "dashed"}, name="Frequency Response" |
| 145 | + ) |
| 146 | + + scale_size_manual(values={"Positive freq (ω > 0)": 1.5, "Negative freq (ω < 0)": 1.0}, name="Frequency Response") |
| 147 | + + scale_alpha_manual( |
| 148 | + values={"Positive freq (ω > 0)": 0.95, "Negative freq (ω < 0)": 0.55}, name="Frequency Response" |
| 149 | + ) |
| 150 | + # Direction arrows |
| 151 | + + geom_segment( |
| 152 | + arrow_df, aes(x="x", y="y", xend="xend", yend="yend"), color="#333333", size=1.0, arrow=arrow(length=0.15) |
| 153 | + ) |
| 154 | + # Phase margin line from origin to gain crossover |
| 155 | + + geom_segment( |
| 156 | + gc_seg_df, aes(x="x", y="y", xend="xend", yend="yend"), color="#9467BD", size=0.7, linetype="dotted", alpha=0.7 |
| 157 | + ) |
| 158 | + # Critical point marker |
| 159 | + + geom_point(critical_df, aes(x="real", y="imaginary"), color="#D62728", size=7, shape="x", stroke=2.5) |
| 160 | + + annotate( |
| 161 | + "text", |
| 162 | + x=-1.95, |
| 163 | + y=0.85, |
| 164 | + label="Critical point\n(−1, 0)", |
| 165 | + color="#D62728", |
| 166 | + size=12, |
| 167 | + fontweight="bold", |
| 168 | + ha="center", |
| 169 | + ) |
| 170 | + + annotate("segment", x=-1.55, y=0.6, xend=-1.05, yend=0.1, color="#D62728", size=0.5, alpha=0.5) |
| 171 | + # Gain crossover marker |
| 172 | + + geom_point( |
| 173 | + pd.DataFrame({"real": [G[gc_idx].real], "imaginary": [G[gc_idx].imag]}), |
| 174 | + aes(x="real", y="imaginary"), |
| 175 | + color="#9467BD", |
| 176 | + size=5.5, |
| 177 | + shape="o", |
| 178 | + stroke=1.8, |
| 179 | + ) |
| 180 | + # Phase crossover marker |
| 181 | + + geom_point( |
| 182 | + pd.DataFrame({"real": [G[pc_idx].real], "imaginary": [G[pc_idx].imag]}), |
| 183 | + aes(x="real", y="imaginary"), |
| 184 | + color="#E8833A", |
| 185 | + size=5.5, |
| 186 | + shape="s", |
| 187 | + stroke=1.8, |
| 188 | + fill="#E8833A", |
| 189 | + ) |
| 190 | + # Gain crossover annotation - positioned in lower-left, well clear of origin and freq labels |
| 191 | + + annotate( |
| 192 | + "text", |
| 193 | + x=-3.2, |
| 194 | + y=-2.6, |
| 195 | + label=f"Gain crossover\nω = {gc_omega:.1f} rad/s · PM = {phase_margin:.0f}°", |
| 196 | + color="#9467BD", |
| 197 | + size=12, |
| 198 | + ha="left", |
| 199 | + fontweight="bold", |
| 200 | + ) |
| 201 | + + annotate( |
| 202 | + "segment", |
| 203 | + x=-2.1, |
| 204 | + y=-2.1, |
| 205 | + xend=G[gc_idx].real - 0.1, |
| 206 | + yend=G[gc_idx].imag - 0.1, |
| 207 | + color="#9467BD", |
| 208 | + size=0.6, |
| 209 | + alpha=0.7, |
| 210 | + ) |
| 211 | + # Phase crossover annotation - positioned in upper-left, well separated |
| 212 | + + annotate( |
| 213 | + "text", |
| 214 | + x=-2.8, |
| 215 | + y=2.6, |
| 216 | + label=f"Phase crossover\nω = {pc_omega:.1f} rad/s · GM = {gain_margin_db:.1f} dB", |
| 217 | + color="#E8833A", |
| 218 | + size=12, |
| 219 | + ha="left", |
| 220 | + fontweight="bold", |
| 221 | + ) |
| 222 | + + annotate( |
| 223 | + "segment", |
| 224 | + x=-1.9, |
| 225 | + y=2.2, |
| 226 | + xend=G[pc_idx].real + 0.05, |
| 227 | + yend=G[pc_idx].imag + 0.05, |
| 228 | + color="#E8833A", |
| 229 | + size=0.6, |
| 230 | + alpha=0.7, |
| 231 | + ) |
| 232 | + # Frequency annotation dots and labels |
| 233 | + + geom_point(annot_df, aes(x="real", y="imaginary"), color="#306998", size=3.5, shape="o", fill="#306998") |
| 234 | +) |
| 235 | + |
| 236 | +# Frequency labels |
| 237 | +for _, row in annot_df.iterrows(): |
| 238 | + plot = plot + annotate("text", x=row["lx"], y=row["ly"], label=row["label"], color="#444444", size=12) |
| 239 | + |
| 240 | +# Axes, scales, and styling |
| 241 | +plot = ( |
| 242 | + plot |
| 243 | + + scale_x_continuous(breaks=[-3, -2, -1, 0, 1, 2, 3, 4, 5]) |
| 244 | + + scale_y_continuous(breaks=[-3, -2, -1, 0, 1, 2, 3]) |
| 245 | + + coord_fixed(ratio=1, xlim=(-3.8, 5.8), ylim=(-3.8, 3.8)) |
| 246 | + + labs( |
| 247 | + title="nyquist-basic · plotnine · pyplots.ai", x="Real Axis [dimensionless]", y="Imaginary Axis [dimensionless]" |
| 248 | + ) |
| 249 | + + guides(color=guide_legend(title="Frequency Response", override_aes={"size": 2})) |
| 250 | + + theme_minimal() |
| 251 | + + theme( |
| 252 | + figure_size=(12, 12), |
| 253 | + plot_title=element_text(size=24, weight="bold", ha="center"), |
| 254 | + axis_title=element_text(size=20), |
| 255 | + axis_text=element_text(size=16, color="#555555"), |
| 256 | + panel_grid_major=element_line(color="#EFEFEF", size=0.3), |
| 257 | + panel_grid_minor=element_blank(), |
| 258 | + plot_background=element_rect(fill="white", color="white"), |
| 259 | + panel_background=element_rect(fill="#FAFBFC", color="white"), |
| 260 | + legend_position="bottom", |
| 261 | + legend_title=element_text(size=14, weight="bold"), |
| 262 | + legend_text=element_text(size=12), |
| 263 | + legend_background=element_rect(fill="white", alpha=0.9), |
| 264 | + ) |
| 265 | +) |
| 266 | + |
| 267 | +plot.save("plot.png", dpi=300, verbose=False) |
0 commit comments