|
| 1 | +""" pyplots.ai |
| 2 | +nyquist-basic: Nyquist Plot for Control Systems |
| 3 | +Library: letsplot 4.9.0 | Python 3.14.3 |
| 4 | +Quality: 91/100 | Created: 2026-03-20 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pandas as pd |
| 9 | +from lets_plot import * # noqa: F403 |
| 10 | +from lets_plot.export import ggsave as export_ggsave |
| 11 | +from scipy import signal |
| 12 | + |
| 13 | + |
| 14 | +LetsPlot.setup_html() # noqa: F405 |
| 15 | + |
| 16 | +# Curated color palette |
| 17 | +COLOR_PRIMARY = "#1B4F72" |
| 18 | +COLOR_MIRROR = "#5DADE2" |
| 19 | +COLOR_CRITICAL = "#C0392B" |
| 20 | +COLOR_STABILITY = "#E74C3C" |
| 21 | +COLOR_UNIT_CIRCLE = "#AEB6BF" |
| 22 | +COLOR_GAIN_MARGIN = "#E74C3C" |
| 23 | +COLOR_ANNOTATION = "#2C3E50" |
| 24 | +COLOR_ARROW = "#1A5276" |
| 25 | + |
| 26 | +# Data - Third-order system: G(s) = 2 / ((s+1)(0.5s+1)(0.2s+1)) |
| 27 | +num = [2.0] |
| 28 | +den = np.polymul(np.polymul([1, 1], [0.5, 1]), [0.2, 1]) |
| 29 | +system = signal.TransferFunction(num, den) |
| 30 | + |
| 31 | +omega = np.logspace(-2, 2, 500) |
| 32 | +_, H = signal.freqresp(system, omega) |
| 33 | + |
| 34 | +real_part = H.real |
| 35 | +imag_part = H.imag |
| 36 | + |
| 37 | +df = pd.DataFrame({"real": real_part, "imaginary": imag_part, "frequency": omega}) |
| 38 | + |
| 39 | +# Mirror (negative frequencies) for complete Nyquist contour |
| 40 | +df_mirror = pd.DataFrame({"real": real_part, "imaginary": -imag_part, "frequency": -omega}) |
| 41 | + |
| 42 | +# Unit circle |
| 43 | +theta = np.linspace(0, 2 * np.pi, 200) |
| 44 | +df_circle = pd.DataFrame({"real": np.cos(theta), "imaginary": np.sin(theta)}) |
| 45 | + |
| 46 | +# Critical point (-1, 0) |
| 47 | +df_critical = pd.DataFrame({"real": [-1.0], "imaginary": [0.0]}) |
| 48 | + |
| 49 | +# Gain margin line: from critical point to nearest point on curve where it crosses negative real axis |
| 50 | +real_cross_mask = (imag_part[1:] * imag_part[:-1] < 0) & (real_part[:-1] < 0) |
| 51 | +cross_indices = np.where(real_cross_mask)[0] |
| 52 | +cross_idx = cross_indices[0] if len(cross_indices) > 0 else np.argmin(np.abs(imag_part[200:])) + 200 |
| 53 | +df_gain_margin = pd.DataFrame({"x": [-1.0], "y": [0.0], "xend": [real_part[cross_idx]], "yend": [0.0]}) |
| 54 | + |
| 55 | +# Arrow indicators along the curve showing direction of increasing frequency |
| 56 | +arrow_indices = [50, 150, 300] |
| 57 | +df_arrows = df.iloc[arrow_indices].copy() |
| 58 | +df_arrows_next = df.iloc[[i + 5 for i in arrow_indices]].copy() |
| 59 | +df_segments = pd.DataFrame( |
| 60 | + { |
| 61 | + "x": df_arrows["real"].values, |
| 62 | + "y": df_arrows["imaginary"].values, |
| 63 | + "xend": df_arrows_next["real"].values, |
| 64 | + "yend": df_arrows_next["imaginary"].values, |
| 65 | + } |
| 66 | +) |
| 67 | + |
| 68 | +# Frequency labels at key points |
| 69 | +label_indices = [0, 80, 200, 350] |
| 70 | +df_labels = df.iloc[label_indices].copy() |
| 71 | +df_labels["label"] = [f"ω={omega[i]:.2f}" if omega[i] < 10 else f"ω={omega[i]:.0f}" for i in label_indices] |
| 72 | +df_labels["nudge_x"] = [0.12, -0.18, -0.12, 0.12] |
| 73 | +df_labels["nudge_y"] = [0.10, -0.14, 0.12, 0.10] |
| 74 | + |
| 75 | +# Stability region highlight - shaded area around critical point |
| 76 | +stability_theta = np.linspace(0, 2 * np.pi, 100) |
| 77 | +stability_r = 0.15 |
| 78 | +df_stability_zone = pd.DataFrame( |
| 79 | + {"real": -1.0 + stability_r * np.cos(stability_theta), "imaginary": stability_r * np.sin(stability_theta)} |
| 80 | +) |
| 81 | + |
| 82 | +# Gain margin annotation |
| 83 | +gain_margin_db = 20 * np.log10(1.0 / abs(real_part[cross_idx])) |
| 84 | +df_gm_label = pd.DataFrame( |
| 85 | + {"real": [(-1.0 + real_part[cross_idx]) / 2], "imaginary": [0.12], "label": [f"GM = {gain_margin_db:.1f} dB"]} |
| 86 | +) |
| 87 | + |
| 88 | +# Plot |
| 89 | +plot = ( |
| 90 | + ggplot() |
| 91 | + # Stability zone highlight around critical point |
| 92 | + + geom_polygon(aes(x="real", y="imaginary"), data=df_stability_zone, fill=COLOR_STABILITY, alpha=0.07) |
| 93 | + # Unit circle |
| 94 | + + geom_path(aes(x="real", y="imaginary"), data=df_circle, color=COLOR_UNIT_CIRCLE, size=0.6, linetype="dashed") |
| 95 | + # Gain margin indicator line |
| 96 | + + geom_segment( |
| 97 | + aes(x="x", y="y", xend="xend", yend="yend"), |
| 98 | + data=df_gain_margin, |
| 99 | + color=COLOR_GAIN_MARGIN, |
| 100 | + size=0.7, |
| 101 | + linetype="dotted", |
| 102 | + ) |
| 103 | + # Gain margin label |
| 104 | + + geom_label( |
| 105 | + aes(x="real", y="imaginary", label="label"), |
| 106 | + data=df_gm_label, |
| 107 | + size=16, |
| 108 | + color=COLOR_CRITICAL, |
| 109 | + fill="#FDEDEC", |
| 110 | + alpha=0.9, |
| 111 | + label_padding=0.3, |
| 112 | + label_r=0.15, |
| 113 | + label_size=0.5, |
| 114 | + fontface="bold", |
| 115 | + ) |
| 116 | + # Mirror curve (negative frequencies) |
| 117 | + + geom_path( |
| 118 | + aes(x="real", y="imaginary"), data=df_mirror, color=COLOR_MIRROR, size=1.0, alpha=0.45, linetype="dashed" |
| 119 | + ) |
| 120 | + # Main Nyquist curve with tooltips |
| 121 | + + geom_path( |
| 122 | + aes(x="real", y="imaginary"), |
| 123 | + data=df, |
| 124 | + color=COLOR_PRIMARY, |
| 125 | + size=2.2, |
| 126 | + tooltips=layer_tooltips() |
| 127 | + .format("real", ".3f") |
| 128 | + .format("imaginary", ".3f") |
| 129 | + .format("frequency", ".3f") |
| 130 | + .line("Re: @real") |
| 131 | + .line("Im: @imaginary") |
| 132 | + .line("ω: @frequency rad/s"), |
| 133 | + ) |
| 134 | + # Direction arrows |
| 135 | + + geom_segment( |
| 136 | + aes(x="x", y="y", xend="xend", yend="yend"), |
| 137 | + data=df_segments, |
| 138 | + color=COLOR_ARROW, |
| 139 | + size=1.5, |
| 140 | + arrow=arrow(length=14, type="closed"), |
| 141 | + ) |
| 142 | + # Critical point (-1, 0) - prominent marker |
| 143 | + + geom_point(aes(x="real", y="imaginary"), data=df_critical, color=COLOR_CRITICAL, size=11, shape=4, stroke=3.5) |
| 144 | + # Critical point label |
| 145 | + + geom_text( |
| 146 | + aes(x="real", y="imaginary"), |
| 147 | + data=pd.DataFrame({"real": [-0.72], "imaginary": [0.28]}), |
| 148 | + label="(-1, 0)", |
| 149 | + size=18, |
| 150 | + color=COLOR_CRITICAL, |
| 151 | + fontface="bold", |
| 152 | + ) |
| 153 | + # Origin marker |
| 154 | + + geom_point( |
| 155 | + aes(x="real", y="imaginary"), |
| 156 | + data=pd.DataFrame({"real": [0.0], "imaginary": [0.0]}), |
| 157 | + color="#7F8C8D", |
| 158 | + size=4, |
| 159 | + shape=3, |
| 160 | + stroke=2.0, |
| 161 | + ) |
| 162 | + # Frequency annotations along curve |
| 163 | + + geom_label( |
| 164 | + aes(x="real", y="imaginary", label="label"), |
| 165 | + data=pd.DataFrame( |
| 166 | + { |
| 167 | + "real": df_labels["real"].values + df_labels["nudge_x"].values, |
| 168 | + "imaginary": df_labels["imaginary"].values + df_labels["nudge_y"].values, |
| 169 | + "label": df_labels["label"].values, |
| 170 | + } |
| 171 | + ), |
| 172 | + size=17, |
| 173 | + color=COLOR_ANNOTATION, |
| 174 | + fill="#F8F9F9", |
| 175 | + alpha=0.9, |
| 176 | + label_padding=0.35, |
| 177 | + label_r=0.2, |
| 178 | + label_size=0.3, |
| 179 | + ) |
| 180 | + # Styling |
| 181 | + + labs( |
| 182 | + x="Re{G(jω)}", |
| 183 | + y="Im{G(jω)}", |
| 184 | + title="nyquist-basic · letsplot · pyplots.ai", |
| 185 | + subtitle="G(s) = 2 / [(s+1)(0.5s+1)(0.2s+1)] — Stable: curve does not encircle (-1, 0)", |
| 186 | + ) |
| 187 | + + coord_fixed(ratio=1) |
| 188 | + + ggsize(1200, 1200) |
| 189 | + + theme_minimal() |
| 190 | + + theme( |
| 191 | + axis_text=element_text(size=16, color="#566573"), |
| 192 | + axis_title=element_text(size=22, color="#2C3E50", face="bold"), |
| 193 | + plot_title=element_text(size=26, color="#1B2631", face="bold"), |
| 194 | + plot_subtitle=element_text(size=16, color="#5D6D7E", face="italic"), |
| 195 | + panel_grid_major=element_line(color="#EAECEE", size=0.2), |
| 196 | + panel_grid_minor=element_blank(), |
| 197 | + plot_background=element_rect(fill="#FDFEFE", color="#FDFEFE"), |
| 198 | + panel_background=element_rect(fill="#FFFFFF", color="#D5D8DC", size=0.3), |
| 199 | + axis_ticks=element_blank(), |
| 200 | + axis_ticks_length=0, |
| 201 | + plot_margin=[40, 40, 25, 25], |
| 202 | + ) |
| 203 | +) |
| 204 | + |
| 205 | +# Save |
| 206 | +export_ggsave(plot, filename="plot.png", path=".", scale=3) |
| 207 | +export_ggsave(plot, filename="plot.html", path=".") |
0 commit comments