|
| 1 | +""" pyplots.ai |
| 2 | +root-locus-basic: Root Locus Plot for Control Systems |
| 3 | +Library: pygal 3.1.0 | Python 3.14.3 |
| 4 | +Quality: 82/100 | Created: 2026-03-20 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pygal |
| 9 | +from pygal.style import Style |
| 10 | + |
| 11 | + |
| 12 | +# Data — root locus for G(s) = (s+3) / [s(s+1)(s+2)(s+4)] |
| 13 | +# Open-loop poles: 0, -1, -2, -4 | Open-loop zero: -3 |
| 14 | +num = np.array([1, 3]) |
| 15 | +den = np.polymul(np.polymul([1, 0], [1, 1]), np.polymul([1, 2], [1, 4])) |
| 16 | + |
| 17 | +ol_poles = np.sort(np.roots(den).real) |
| 18 | +ol_zeros = np.sort(np.roots(num).real) |
| 19 | +n_branches = len(den) - 1 |
| 20 | +num_padded = np.zeros(len(den)) |
| 21 | +num_padded[-len(num) :] = num |
| 22 | + |
| 23 | +# Gain sweep with variable density: finer near breakaway and jω crossing |
| 24 | +gains = np.concatenate( |
| 25 | + [np.linspace(0, 2, 200), np.linspace(2, 15, 300), np.linspace(15, 80, 200), np.linspace(80, 500, 150)] |
| 26 | +) |
| 27 | + |
| 28 | +# Compute closed-loop poles for each gain: roots of den(s) + K·num(s) = 0 |
| 29 | +loci = np.zeros((len(gains), n_branches), dtype=complex) |
| 30 | +for i, K in enumerate(gains): |
| 31 | + roots = np.roots(den + K * num_padded) |
| 32 | + if i == 0: |
| 33 | + loci[i] = roots[np.argsort(roots.real)] |
| 34 | + else: |
| 35 | + prev = loci[i - 1] |
| 36 | + available = list(range(n_branches)) |
| 37 | + for j in range(n_branches): |
| 38 | + dists = [abs(roots[k] - prev[j]) if k in available else np.inf for k in range(n_branches)] |
| 39 | + best = int(np.argmin(dists)) |
| 40 | + loci[i, j] = roots[best] |
| 41 | + available.remove(best) |
| 42 | + |
| 43 | +# Find imaginary axis crossings (stability boundary) |
| 44 | +jw_crossings = [] |
| 45 | +for b in range(n_branches): |
| 46 | + reals = loci[:, b].real |
| 47 | + for i in range(len(reals) - 1): |
| 48 | + if reals[i] * reals[i + 1] < 0 and abs(loci[i, b].imag) > 0.1: |
| 49 | + frac = abs(reals[i]) / (abs(reals[i]) + abs(reals[i + 1])) |
| 50 | + im = float(loci[i, b].imag + frac * (loci[i + 1, b].imag - loci[i, b].imag)) |
| 51 | + K_cross = float(gains[i] + frac * (gains[i + 1] - gains[i])) |
| 52 | + jw_crossings.append((round(im, 3), round(K_cross, 2))) |
| 53 | + |
| 54 | +# Find breakaway point: on real axis between poles at -1 and -2 |
| 55 | +# d/ds[1 + K·N(s)/D(s)] = 0 → d/ds[D(s)/N(s)] = 0 |
| 56 | +# Numerically search between -1 and -2 |
| 57 | +s_test = np.linspace(-1.01, -1.99, 500) |
| 58 | +ratio = np.polyval(den, s_test) / np.polyval(num, s_test) |
| 59 | +breakaway_idx = np.argmin(np.abs(np.gradient(ratio, s_test))) |
| 60 | +breakaway_s = round(float(s_test[breakaway_idx]), 3) |
| 61 | +breakaway_K = round(float(-np.polyval(den, breakaway_s) / np.polyval(num, breakaway_s)), 2) |
| 62 | + |
| 63 | +# Real-axis locus segments (to the left of an odd number of real poles+zeros) |
| 64 | +real_segments = [(0, -1), (-2, -3), (-4, -6)] |
| 65 | + |
| 66 | +# Constant damping ratio guide lines (ζ = 0.2, 0.4, 0.6, 0.8) |
| 67 | +zeta_values = [0.2, 0.4, 0.6, 0.8] |
| 68 | +guide_extent = 5.5 |
| 69 | + |
| 70 | +# Constant natural frequency semicircles (ωn = 1, 2, 3, 4, 5) |
| 71 | +wn_values = [1, 2, 3, 4, 5] |
| 72 | + |
| 73 | +# Style — refined palette with warm background tint for polish |
| 74 | +font = "DejaVu Sans, Helvetica, Arial, sans-serif" |
| 75 | +custom_style = Style( |
| 76 | + background="#fafafa", |
| 77 | + plot_background="#fefefe", |
| 78 | + foreground="#1a1a2e", |
| 79 | + foreground_strong="#1a1a2e", |
| 80 | + foreground_subtle="#d8d8d8", |
| 81 | + guide_stroke_color="#eaeaea", |
| 82 | + guide_stroke_dasharray="2, 4", |
| 83 | + colors=( |
| 84 | + "#c0c0c0", # 0: ζ guide lines (soft gray, background) |
| 85 | + "#c0c0c0", # 1: ωn guide semicircles (soft gray, background) |
| 86 | + "#2563EB", # 2: Root locus branches (vivid blue) |
| 87 | + "#EA580C", # 3: Real-axis locus (burnt orange) |
| 88 | + "#DC2626", # 4: Open-loop poles (red) |
| 89 | + "#7C3AED", # 5: Open-loop zero (violet — colorblind-safe vs red) |
| 90 | + "#16A34A", # 6: Breakaway point (green — distinct emphasis) |
| 91 | + "#D97706", # 7: Stability boundary (amber) |
| 92 | + "#1D4ED8", # 8: Gain direction markers (dark blue) |
| 93 | + ), |
| 94 | + font_family=font, |
| 95 | + title_font_family=font, |
| 96 | + title_font_size=52, |
| 97 | + label_font_size=40, |
| 98 | + major_label_font_size=36, |
| 99 | + legend_font_size=26, |
| 100 | + legend_font_family=font, |
| 101 | + value_font_size=24, |
| 102 | + tooltip_font_size=26, |
| 103 | + tooltip_font_family=font, |
| 104 | + opacity=0.92, |
| 105 | + opacity_hover=1.0, |
| 106 | + stroke_width=5, |
| 107 | +) |
| 108 | + |
| 109 | +chart = pygal.XY( |
| 110 | + width=3600, |
| 111 | + height=3600, |
| 112 | + style=custom_style, |
| 113 | + title="root-locus-basic · pygal · pyplots.ai", |
| 114 | + x_title="Real Axis (σ)", |
| 115 | + y_title="Imaginary Axis (jω)", |
| 116 | + show_legend=True, |
| 117 | + legend_at_bottom=True, |
| 118 | + legend_at_bottom_columns=4, |
| 119 | + legend_box_size=22, |
| 120 | + stroke=True, |
| 121 | + dots_size=0, |
| 122 | + show_x_guides=True, |
| 123 | + show_y_guides=True, |
| 124 | + x_value_formatter=lambda v: f"{v:.1f}", |
| 125 | + value_formatter=lambda v: f"{v:.1f}", |
| 126 | + margin_bottom=90, |
| 127 | + margin_left=80, |
| 128 | + margin_right=50, |
| 129 | + margin_top=55, |
| 130 | + xrange=(-6, 4), |
| 131 | + range=(-5, 5), |
| 132 | + print_values=False, |
| 133 | + print_zeroes=False, |
| 134 | + js=[], |
| 135 | + truncate_legend=-1, |
| 136 | + include_x_axis=True, |
| 137 | + allow_interruptions=True, |
| 138 | + spacing=18, |
| 139 | +) |
| 140 | + |
| 141 | +# Constant damping ratio guide lines (ζ rays from origin) — rendered first as background |
| 142 | +zeta_pts = [] |
| 143 | +for zeta in zeta_values: |
| 144 | + theta = np.arccos(zeta) |
| 145 | + for t in np.linspace(0, guide_extent, 25): |
| 146 | + zeta_pts.append((round(-t * np.cos(theta), 3), round(t * np.sin(theta), 3))) |
| 147 | + zeta_pts.append(None) |
| 148 | + for t in np.linspace(0, guide_extent, 25): |
| 149 | + zeta_pts.append((round(-t * np.cos(theta), 3), round(-t * np.sin(theta), 3))) |
| 150 | + zeta_pts.append(None) |
| 151 | +chart.add( |
| 152 | + "ζ guides", zeta_pts, stroke_style={"width": 1.2, "dasharray": "6, 5"}, show_dots=False, allow_interruptions=True |
| 153 | +) |
| 154 | + |
| 155 | +# Constant natural frequency semicircles (ωn arcs in left half-plane) |
| 156 | +wn_pts = [] |
| 157 | +for wn in wn_values: |
| 158 | + angles = np.linspace(np.pi / 2, 3 * np.pi / 2, 50) |
| 159 | + for a in angles: |
| 160 | + wn_pts.append((round(wn * np.cos(a), 3), round(wn * np.sin(a), 3))) |
| 161 | + wn_pts.append(None) |
| 162 | +chart.add( |
| 163 | + "ωn guides", wn_pts, stroke_style={"width": 1.2, "dasharray": "6, 5"}, show_dots=False, allow_interruptions=True |
| 164 | +) |
| 165 | + |
| 166 | +# Root locus branches — skip points near zero at s=-3 to avoid obscuring it |
| 167 | +zero_exclusion_radius = 0.35 |
| 168 | +locus_pts = [] |
| 169 | +for b in range(n_branches): |
| 170 | + branch_data = [] |
| 171 | + for i in range(len(gains)): |
| 172 | + r, im = float(loci[i, b].real), float(loci[i, b].imag) |
| 173 | + if -6 <= r <= 4 and -5 <= im <= 5: |
| 174 | + near_zero = any( |
| 175 | + abs(r - float(z)) < zero_exclusion_radius and abs(im) < zero_exclusion_radius for z in ol_zeros |
| 176 | + ) |
| 177 | + if near_zero: |
| 178 | + if branch_data: |
| 179 | + locus_pts.extend(branch_data) |
| 180 | + locus_pts.append(None) |
| 181 | + branch_data = [] |
| 182 | + else: |
| 183 | + branch_data.append({"value": (round(r, 4), round(im, 4)), "label": f"K = {gains[i]:.2f}"}) |
| 184 | + if branch_data: |
| 185 | + locus_pts.extend(branch_data) |
| 186 | + locus_pts.append(None) |
| 187 | + |
| 188 | +chart.add( |
| 189 | + "Root Locus", locus_pts, stroke_style={"width": 7, "linecap": "round"}, show_dots=False, allow_interruptions=True |
| 190 | +) |
| 191 | + |
| 192 | +# Real-axis locus segments — thick orange line |
| 193 | +real_pts = [] |
| 194 | +for seg_start, seg_end in real_segments: |
| 195 | + for x in np.linspace(seg_start, seg_end, 60): |
| 196 | + real_pts.append((round(float(x), 3), 0.0)) |
| 197 | + real_pts.append(None) |
| 198 | +chart.add( |
| 199 | + "Real-Axis Locus", |
| 200 | + real_pts, |
| 201 | + stroke_style={"width": 12, "linecap": "round"}, |
| 202 | + show_dots=False, |
| 203 | + allow_interruptions=True, |
| 204 | +) |
| 205 | + |
| 206 | +# Open-loop poles (marked with ×) |
| 207 | +pole_pts = [{"value": (round(float(p), 2), 0.0), "label": f"Pole at s = {p:.0f}"} for p in ol_poles] |
| 208 | +chart.add("Poles (×)", pole_pts, stroke=False, dots_size=22) |
| 209 | + |
| 210 | +# Open-loop zero (marked with ○) — large dot, rendered after locus to stay on top |
| 211 | +zero_pts = [{"value": (round(float(z), 2), 0.0), "label": f"Zero at s = {z:.0f}"} for z in ol_zeros] |
| 212 | +chart.add("Zero (○)", zero_pts, stroke=False, dots_size=28) |
| 213 | + |
| 214 | +# Breakaway point — emphasized with distinct green marker |
| 215 | +breakaway_pts = [{"value": (breakaway_s, 0.0), "label": f"Breakaway: s = {breakaway_s}, K = {breakaway_K}"}] |
| 216 | +chart.add("Breakaway Point", breakaway_pts, stroke=False, dots_size=30) |
| 217 | + |
| 218 | +# Imaginary axis crossings — stability boundary markers |
| 219 | +jw_pts = [{"value": (0.0, im), "label": f"jω crossing: s = {im:+.3f}j, K = {K:.2f}"} for im, K in jw_crossings] |
| 220 | +chart.add("Stability Boundary", jw_pts, stroke=False, dots_size=24) |
| 221 | + |
| 222 | +# Gain direction markers along locus at selected gains |
| 223 | +arrow_pts = [] |
| 224 | +arrow_gains = [5, 15, 40, 100, 250] |
| 225 | +for ag in arrow_gains: |
| 226 | + idx = np.argmin(np.abs(gains - ag)) |
| 227 | + for b in range(n_branches): |
| 228 | + r, im = float(loci[idx, b].real), float(loci[idx, b].imag) |
| 229 | + if -6 <= r <= 4 and -5 <= im <= 5: |
| 230 | + arrow_pts.append({"value": (round(r, 3), round(im, 3)), "label": f"K = {ag} →"}) |
| 231 | +chart.add("Gain K →", arrow_pts, stroke=False, dots_size=20) |
| 232 | + |
| 233 | +# Save |
| 234 | +chart.render_to_png("plot.png") |
| 235 | +chart.render_to_file("plot.html") |
0 commit comments