|
| 1 | +""" pyplots.ai |
| 2 | +nyquist-basic: Nyquist Plot for Control Systems |
| 3 | +Library: altair 6.0.0 | Python 3.14.3 |
| 4 | +Quality: 91/100 | Created: 2026-03-20 |
| 5 | +""" |
| 6 | + |
| 7 | +import altair as alt |
| 8 | +import numpy as np |
| 9 | +import pandas as pd |
| 10 | + |
| 11 | + |
| 12 | +# Data - Nyquist plot for G(s) = 5 / ((s+1)(0.5s+1)(0.1s+1)) |
| 13 | +# A stable third-order system with DC gain = 5 |
| 14 | +omega = np.concatenate( |
| 15 | + [np.logspace(-2, -0.5, 100), np.logspace(-0.5, 0.5, 250), np.logspace(0.5, 1.5, 200), np.logspace(1.5, 3, 100)] |
| 16 | +) |
| 17 | + |
| 18 | +K = 5.0 |
| 19 | +s = 1j * omega |
| 20 | +G = K / ((s + 1.0) * (0.5 * s + 1.0) * (0.1 * s + 1.0)) |
| 21 | + |
| 22 | +real_part = G.real |
| 23 | +imag_part = G.imag |
| 24 | + |
| 25 | +# Positive frequency branch |
| 26 | +pos_df = pd.DataFrame( |
| 27 | + { |
| 28 | + "real": real_part, |
| 29 | + "imaginary": imag_part, |
| 30 | + "frequency": omega, |
| 31 | + "branch": "ω ≥ 0 (positive)", |
| 32 | + "idx": np.arange(len(omega)), |
| 33 | + } |
| 34 | +) |
| 35 | + |
| 36 | +# Negative frequency branch (mirror about real axis) |
| 37 | +neg_df = pd.DataFrame( |
| 38 | + { |
| 39 | + "real": real_part, |
| 40 | + "imaginary": -imag_part, |
| 41 | + "frequency": -omega[::-1], |
| 42 | + "branch": "ω ≤ 0 (negative)", |
| 43 | + "idx": np.arange(len(omega)), |
| 44 | + } |
| 45 | +) |
| 46 | + |
| 47 | +nyquist_df = pd.concat([pos_df, neg_df], ignore_index=True) |
| 48 | + |
| 49 | +# Unit circle for reference |
| 50 | +theta = np.linspace(0, 2 * np.pi, 200) |
| 51 | +unit_circle_df = pd.DataFrame({"ux": np.cos(theta), "uy": np.sin(theta), "idx": np.arange(len(theta))}) |
| 52 | + |
| 53 | +# Critical point (-1, 0) |
| 54 | +critical_df = pd.DataFrame({"real": [-1.0], "imaginary": [0.0], "label": ["Critical Point (−1, 0)"]}) |
| 55 | + |
| 56 | +# Find gain crossover (|G(jω)| = 1) |
| 57 | +magnitude = np.abs(G) |
| 58 | +gain_cross_idx = np.argmin(np.abs(magnitude - 1.0)) |
| 59 | +gain_cross_omega = omega[gain_cross_idx] |
| 60 | + |
| 61 | +# Phase crossover: where imaginary part crosses zero (excluding near DC) |
| 62 | +sign_changes = np.where(np.diff(np.sign(imag_part[30:])))[0] + 30 |
| 63 | +phase_cross_idx = sign_changes[0] if len(sign_changes) > 0 else None |
| 64 | +phase_cross_omega = omega[phase_cross_idx] if phase_cross_idx is not None else None |
| 65 | + |
| 66 | +# Frequency markers along curve — offset label_y to avoid overlaps |
| 67 | +freq_annotations = [] |
| 68 | +for w_mark, dy_px in [(0.5, -18), (1.0, -18), (5.0, -18)]: |
| 69 | + idx = np.argmin(np.abs(omega - w_mark)) |
| 70 | + freq_annotations.append( |
| 71 | + {"real": real_part[idx], "imaginary": imag_part[idx], "label": f"ω = {w_mark}", "dy_px": dy_px} |
| 72 | + ) |
| 73 | + |
| 74 | +# ω = 2.0 and gain crossover (ω ≈ 2.7) are close — place on opposite sides |
| 75 | +idx_2 = np.argmin(np.abs(omega - 2.0)) |
| 76 | +freq_annotations.append({"real": real_part[idx_2], "imaginary": imag_part[idx_2], "label": "ω = 2.0", "dy_px": -20}) |
| 77 | +freq_annotations.append( |
| 78 | + { |
| 79 | + "real": real_part[gain_cross_idx], |
| 80 | + "imaginary": imag_part[gain_cross_idx], |
| 81 | + "label": f"ω ≈ {gain_cross_omega:.1f} (|G|=1)", |
| 82 | + "dy_px": 20, |
| 83 | + } |
| 84 | +) |
| 85 | + |
| 86 | +# Phase crossover is near the critical point — annotate it separately |
| 87 | +# to avoid overlap with the critical point label |
| 88 | +if phase_cross_idx is not None: |
| 89 | + phase_cross_real = real_part[phase_cross_idx] |
| 90 | + phase_cross_imag = imag_part[phase_cross_idx] |
| 91 | + phase_cross_label = f"ω ≈ {phase_cross_omega:.1f} (phase crossover)" |
| 92 | + |
| 93 | +freq_df = pd.DataFrame(freq_annotations) |
| 94 | + |
| 95 | +# Arrow indicators for direction of increasing frequency |
| 96 | +arrow_rows = [] |
| 97 | +for target_w in [0.8, 3.0]: |
| 98 | + idx = np.argmin(np.abs(omega - target_w)) |
| 99 | + im_val = imag_part[idx] |
| 100 | + if abs(im_val) > 0.05: |
| 101 | + shape = "down" if im_val < 0 else "up" |
| 102 | + arrow_rows.append({"ax": real_part[idx], "ay": imag_part[idx], "branch": "ω ≥ 0 (positive)", "shape": shape}) |
| 103 | + arrow_rows.append( |
| 104 | + { |
| 105 | + "ax": real_part[idx], |
| 106 | + "ay": -imag_part[idx], |
| 107 | + "branch": "ω ≤ 0 (negative)", |
| 108 | + "shape": "up" if shape == "down" else "down", |
| 109 | + } |
| 110 | + ) |
| 111 | +arrow_df = pd.DataFrame(arrow_rows) |
| 112 | + |
| 113 | +# Axis scales - 1:1 aspect ratio, tighter domain around data |
| 114 | +plot_range = 5.5 |
| 115 | +x_scale = alt.Scale(domain=[-plot_range, plot_range], nice=False) |
| 116 | +y_scale = alt.Scale(domain=[-plot_range, plot_range], nice=False) |
| 117 | + |
| 118 | +branch_palette = ["#306998", "#e07b39"] |
| 119 | +branch_domain = ["ω ≥ 0 (positive)", "ω ≤ 0 (negative)"] |
| 120 | + |
| 121 | +# Layer: Nyquist curve with conditional selection highlight |
| 122 | +highlight = alt.selection_point(fields=["branch"], bind="legend") |
| 123 | + |
| 124 | +nyquist_layer = ( |
| 125 | + alt.Chart(nyquist_df) |
| 126 | + .mark_line(strokeWidth=2.8) |
| 127 | + .encode( |
| 128 | + x=alt.X( |
| 129 | + "real:Q", |
| 130 | + scale=x_scale, |
| 131 | + title="Real Part — Re[G(jω)]", |
| 132 | + axis=alt.Axis( |
| 133 | + labelFontSize=16, |
| 134 | + titleFontSize=20, |
| 135 | + titleFontWeight="bold", |
| 136 | + titleColor="#2a2a2a", |
| 137 | + labelColor="#444444", |
| 138 | + grid=False, |
| 139 | + titlePadding=14, |
| 140 | + domainColor="#aaaaaa", |
| 141 | + tickColor="#aaaaaa", |
| 142 | + ), |
| 143 | + ), |
| 144 | + y=alt.Y( |
| 145 | + "imaginary:Q", |
| 146 | + scale=y_scale, |
| 147 | + title="Imaginary Part — Im[G(jω)]", |
| 148 | + axis=alt.Axis( |
| 149 | + labelFontSize=16, |
| 150 | + titleFontSize=20, |
| 151 | + titleFontWeight="bold", |
| 152 | + titleColor="#2a2a2a", |
| 153 | + labelColor="#444444", |
| 154 | + grid=False, |
| 155 | + titlePadding=14, |
| 156 | + domainColor="#aaaaaa", |
| 157 | + tickColor="#aaaaaa", |
| 158 | + ), |
| 159 | + ), |
| 160 | + color=alt.Color( |
| 161 | + "branch:N", |
| 162 | + scale=alt.Scale(domain=branch_domain, range=branch_palette), |
| 163 | + legend=alt.Legend( |
| 164 | + title="Frequency Branch", |
| 165 | + titleFontSize=16, |
| 166 | + labelFontSize=14, |
| 167 | + symbolSize=180, |
| 168 | + symbolStrokeWidth=3, |
| 169 | + orient="top-right", |
| 170 | + offset=5, |
| 171 | + ), |
| 172 | + ), |
| 173 | + opacity=alt.condition(highlight, alt.value(0.9), alt.value(0.2)), |
| 174 | + order="idx:Q", |
| 175 | + tooltip=[ |
| 176 | + alt.Tooltip("branch:N", title="Branch"), |
| 177 | + alt.Tooltip("real:Q", title="Re(G)", format=".3f"), |
| 178 | + alt.Tooltip("imaginary:Q", title="Im(G)", format=".3f"), |
| 179 | + alt.Tooltip("frequency:Q", title="ω (rad/s)", format=".3f"), |
| 180 | + ], |
| 181 | + ) |
| 182 | + .add_params(highlight) |
| 183 | +) |
| 184 | + |
| 185 | +# Layer: Unit circle |
| 186 | +unit_circle_layer = ( |
| 187 | + alt.Chart(unit_circle_df) |
| 188 | + .mark_line(strokeWidth=1.5, strokeDash=[6, 4], color="#cccccc", opacity=0.6) |
| 189 | + .encode(x=alt.X("ux:Q", scale=x_scale), y=alt.Y("uy:Q", scale=y_scale), order="idx:Q") |
| 190 | +) |
| 191 | + |
| 192 | +# Layer: Critical point (-1, 0) with pulsing ring |
| 193 | +critical_ring = ( |
| 194 | + alt.Chart(critical_df) |
| 195 | + .mark_point(shape="circle", size=800, strokeWidth=2, color="#d62728", filled=False, opacity=0.3) |
| 196 | + .encode(x=alt.X("real:Q", scale=x_scale), y=alt.Y("imaginary:Q", scale=y_scale)) |
| 197 | +) |
| 198 | + |
| 199 | +critical_layer = ( |
| 200 | + alt.Chart(critical_df) |
| 201 | + .mark_point(shape="cross", size=500, strokeWidth=4, color="#d62728", filled=False) |
| 202 | + .encode( |
| 203 | + x=alt.X("real:Q", scale=x_scale), |
| 204 | + y=alt.Y("imaginary:Q", scale=y_scale), |
| 205 | + tooltip=[alt.Tooltip("label:N", title="Critical Point")], |
| 206 | + ) |
| 207 | +) |
| 208 | + |
| 209 | +critical_text = ( |
| 210 | + alt.Chart(critical_df) |
| 211 | + .mark_text(fontSize=16, fontWeight="bold", color="#c5211e", align="right", dx=-18, dy=-16) |
| 212 | + .encode(x=alt.X("real:Q", scale=x_scale), y=alt.Y("imaginary:Q", scale=y_scale), text="label:N") |
| 213 | +) |
| 214 | + |
| 215 | +# Layer: Phase crossover annotation (positioned below critical point to avoid overlap) |
| 216 | +phase_cross_layers = [] |
| 217 | +if phase_cross_idx is not None: |
| 218 | + pc_df = pd.DataFrame({"real": [phase_cross_real], "imaginary": [phase_cross_imag], "label": [phase_cross_label]}) |
| 219 | + phase_cross_point = ( |
| 220 | + alt.Chart(pc_df) |
| 221 | + .mark_point(shape="diamond", size=250, color="#8b4513", filled=True, opacity=0.85) |
| 222 | + .encode( |
| 223 | + x=alt.X("real:Q", scale=x_scale), |
| 224 | + y=alt.Y("imaginary:Q", scale=y_scale), |
| 225 | + tooltip=[alt.Tooltip("label:N", title="Phase Crossover")], |
| 226 | + ) |
| 227 | + ) |
| 228 | + phase_cross_text = ( |
| 229 | + alt.Chart(pc_df) |
| 230 | + .mark_text(fontSize=13, color="#8b4513", fontWeight="bold", align="left", dx=14, dy=18) |
| 231 | + .encode(x=alt.X("real:Q", scale=x_scale), y=alt.Y("imaginary:Q", scale=y_scale), text="label:N") |
| 232 | + ) |
| 233 | + phase_cross_layers = [phase_cross_point, phase_cross_text] |
| 234 | + |
| 235 | +# Layer: Frequency annotation points |
| 236 | +freq_points = ( |
| 237 | + alt.Chart(freq_df) |
| 238 | + .mark_point(shape="circle", size=200, color="#306998", filled=True, opacity=0.85) |
| 239 | + .encode( |
| 240 | + x=alt.X("real:Q", scale=x_scale), |
| 241 | + y=alt.Y("imaginary:Q", scale=y_scale), |
| 242 | + tooltip=[alt.Tooltip("label:N", title="Frequency")], |
| 243 | + ) |
| 244 | +) |
| 245 | + |
| 246 | +# Split labels by dy offset to avoid overlap (Altair mark-level dy is static) |
| 247 | +freq_above = freq_df[freq_df["dy_px"] < 0].copy() |
| 248 | +freq_below = freq_df[freq_df["dy_px"] > 0].copy() |
| 249 | + |
| 250 | +freq_labels_above = ( |
| 251 | + alt.Chart(freq_above) |
| 252 | + .mark_text(fontSize=14, color="#333333", fontWeight="bold", align="left", dx=14, dy=-18) |
| 253 | + .encode(x=alt.X("real:Q", scale=x_scale), y=alt.Y("imaginary:Q", scale=y_scale), text="label:N") |
| 254 | +) |
| 255 | + |
| 256 | +freq_labels_below = ( |
| 257 | + alt.Chart(freq_below) |
| 258 | + .mark_text(fontSize=14, color="#333333", fontWeight="bold", align="left", dx=14, dy=20) |
| 259 | + .encode(x=alt.X("real:Q", scale=x_scale), y=alt.Y("imaginary:Q", scale=y_scale), text="label:N") |
| 260 | +) |
| 261 | + |
| 262 | +# Layer: Direction arrows |
| 263 | +arrow_up_df = arrow_df[arrow_df["shape"] == "up"].copy() |
| 264 | +arrow_down_df = arrow_df[arrow_df["shape"] == "down"].copy() |
| 265 | + |
| 266 | +arrow_up_layer = ( |
| 267 | + alt.Chart(arrow_up_df) |
| 268 | + .mark_point(shape="triangle-up", size=280, filled=True, opacity=0.85) |
| 269 | + .encode( |
| 270 | + x=alt.X("ax:Q", scale=x_scale), |
| 271 | + y=alt.Y("ay:Q", scale=y_scale), |
| 272 | + color=alt.Color("branch:N", scale=alt.Scale(domain=branch_domain, range=branch_palette), legend=None), |
| 273 | + ) |
| 274 | +) |
| 275 | + |
| 276 | +arrow_down_layer = ( |
| 277 | + alt.Chart(arrow_down_df) |
| 278 | + .mark_point(shape="triangle-down", size=280, filled=True, opacity=0.85) |
| 279 | + .encode( |
| 280 | + x=alt.X("ax:Q", scale=x_scale), |
| 281 | + y=alt.Y("ay:Q", scale=y_scale), |
| 282 | + color=alt.Color("branch:N", scale=alt.Scale(domain=branch_domain, range=branch_palette), legend=None), |
| 283 | + ) |
| 284 | +) |
| 285 | + |
| 286 | +# Compose all layers |
| 287 | +combined = ( |
| 288 | + unit_circle_layer |
| 289 | + + nyquist_layer |
| 290 | + + critical_ring |
| 291 | + + critical_layer |
| 292 | + + critical_text |
| 293 | + + freq_points |
| 294 | + + freq_labels_above |
| 295 | + + freq_labels_below |
| 296 | + + arrow_up_layer |
| 297 | + + arrow_down_layer |
| 298 | +) |
| 299 | +for layer in phase_cross_layers: |
| 300 | + combined = combined + layer |
| 301 | + |
| 302 | +chart = ( |
| 303 | + combined.properties( |
| 304 | + width=1200, |
| 305 | + height=1200, |
| 306 | + title=alt.Title( |
| 307 | + "nyquist-basic · altair · pyplots.ai", |
| 308 | + fontSize=28, |
| 309 | + fontWeight="bold", |
| 310 | + color="#1a1a1a", |
| 311 | + subtitle="G(s) = 5 / (s+1)(0.5s+1)(0.1s+1) · Open-Loop Frequency Response", |
| 312 | + subtitleFontSize=18, |
| 313 | + subtitleColor="#555555", |
| 314 | + subtitlePadding=10, |
| 315 | + anchor="start", |
| 316 | + offset=10, |
| 317 | + ), |
| 318 | + ) |
| 319 | + .configure_view(strokeWidth=0) |
| 320 | + .interactive() |
| 321 | +) |
| 322 | + |
| 323 | +chart.save("plot.png", scale_factor=3.0) |
| 324 | +chart.save("plot.html") |
0 commit comments