Compute the default axis limits of nyquistplot per subplot from the curves and the critical point

[baggepinnen]

Part of #1081

Problem

The default axis limits of the Plots recipe nyquistplot were computed as

mask = @. (-20 ≤ imdata ≤ 20) & (-20 ≤ redata ≤ 20)
ylims --> (min(minimum(imdata[mask]),-1.05), max(maximum(imdata[mask]),1.05))
xlims --> (min(minimum(redata[mask]),-1.05), max(maximum(redata[mask]),1.05))

outside of the @series block, with the following consequences:

1. The window is fixed to ±20 independently of the gain: nyquistplot(1000/(s+1)^3) shows Re ∈ (-19.8, 1.05), Im ∈ (-1.05, 15.9), although the curve spans Re ∈ [-250, 1000], Im ∈ [-730, 20].
2. A curve located entirely outside of ±20 throws: nyquistplot(ss(100) + tf(1,[1,1])) gives ArgumentError: reducing over an empty collection is not allowed.
3. Since --> does not overwrite, the limits of channel (1,1) of the first system are used for all subplots and all systems: nyquistplot([tf(1,[1,1]) tf(10,[1,1]); tf(0.1,[1,1]) tf(10,[1,2,1])]) gives (-1.05, 1.05) × (-1.05, 1.05) in every subplot, and nyquistplot([tf(1,[1,1])^2, tf(10,[1,1])^2]) gives (-1.05, 1.05) × (-1.05, 1.05).
4. The limits contain the point -1 also if another critical_point is given. The Makie recipe sets no limits; for 1/(s(s+1)) the axis extends to Im = -100, and the region around -1 cannot be read.

Change

A new function ControlSystemsBase.nyquist_limits(systems, w; critical_point, circles, max_factor = 5, pad = 0.1, responses) returns a matrix with the limits (xlims, ylims) of each channel. The limits are the smallest box containing, for all systems,

• the critical point c and the square of half side |c|/2 centered at c,
• the origin,
• the requested Ms, Mt, disk-margin and unit circles of radius at most 2 (the circles are defined relative to -1; Mt circles with Mt close to 1 have a very large radius and are not included),
• the Nyquist curve at the frequencies selected by nyquist_limit_mask(sys, w; max_factor).

The box is enlarged by the fraction pad of its extent along each axis.

nyquist_limit_mask keeps all frequencies if the system has no poles on the imaginary axis (the unit circle in discrete time). Otherwise, it excludes the frequencies at which the magnitude of the factor contributed by these poles exceeds max_factor:

• n poles at ±iωp: factor (ωp²/|ωp² - ω²|)^n, equal to one at ω = 0,
• n poles in the origin: factor (ωl/ω)^n, where ωl is the smallest magnitude of the poles and zeros not located in the origin. If there is no such pole or zero (e.g. 1/s^2), the curve has no intrinsic scale and is not used.

A pole is classified as located on the imaginary axis if count_eigval_multiplicity counts it among the poles at its projection onto the imaginary axis (the unit circle), i.e., with the tolerance (100ϵ)^(1/n)·max|p| used for count_integrators. Discrete-time poles and zeros are mapped by s = log(z)/Ts. The result is invariant under time_scale (tested). For a DelayLtiSystem, the poles and zeros of the delay-free system are used; for other non-rational systems, the entire curve is used.

The magnitude-based alternative (excluding points whose magnitude exceeds a multiple of the median magnitude of the curve) was not chosen, since the median depends on the frequency grid: for a user-provided grid extending far beyond the bandwidth, the median is small and parts of a finite curve would be excluded.

The Plots recipe sets xlims --> and ylims --> in the last series of each subplot (where title and yguide are already set, since Plots applies the subplot attributes of the last series of each subplot), so that each subplot obtains its own limits and limits provided by the user take precedence. The Makie recipe calls limits! on each Cartesian axis unless the keyword argument limits is provided; polar axes are unchanged (rlimits).

Systems with other number types: the default limits require the poles and zeros, so they are computed only if float(numeric_type(sys)) is a LAPACK type (BlasFloat) for all systems. For other number types (e.g. MonteCarloMeasurements.Particles, dual numbers, BigFloat without GenericSchur), the limits are left to the plotting backend, as for any plot without explicit limits. The Plots recipe also skips the computation if both xlims and ylims are provided, and the Makie recipe if limits is provided. This keeps nyquistplot working for systems with uncertain coefficients, which a downstream test of RobustAndOptimalControl uses (nyquistplot(delay(Particles(...))*tf(1, [0.01, 1])*pid(2, 1), w; xlims, ylims) errored in count_eigval_multiplicity with the first commit of this PR).

Effect on users

The default limits of all Nyquist plots change:

• Plots: for systems without poles on the imaginary axis, the whole curve is shown instead of the part inside ±20. The neighborhood of the critical point is (-1.5, -0.5) × (-0.5, 0.5) plus padding instead of (-1.05, 1.05)²; the right half-plane is only included where the curve or a requested circle is located. For 1/(s(s+1)) the window is Re ∈ (-1.65, 0.15), Im ∈ (-5.3, 1.0) instead of (-1.05, 1.05) × (-19.9, 1.05).
• Makie: Cartesian axes now have default limits as in the Plots version.

Remaining differences between the backends (not addressed here): the Makie recipe has no disk_margin_circles keyword argument and no hz hover information.

Visual comparison

Rendered with Plots/GR and CairoMakie (before/after, not committed):

• 1000/(s+1)^3, Plots: before, only the part with Re ∈ (-20, 1) of the high-frequency end of the curve is visible; after, the whole spiral from 1000 down to the origin is visible.
• 2×2 MIMO example, Plots: before, all subplots (-1.05, 1.05)², the curves of gain 10 leave the axes; after, each subplot is fitted to its curve and contains -1.
• 1/(s(s+1)), Makie: before, a thin axis with Im down to -100 and Re ticks overlapping; after, Im ∈ (-5.3, 1.0), the asymptote Re = -1 and the critical point are readable.
• 1/((s²+1)(s+1)), Plots: before, Re, Im ∈ (-15, 15); after, approximately (-2.6, 3) × (-2.8, 2.8) with the critical point in the center.

Tests

New testset nyquistplot limits in test_plots.jl (reproducers above, limits contain the curve of the high-gain system, per-channel limits of the MIMO example and the union over two systems, invariance under time_scale, exact masks for integrators, double integrators, undamped modes and a discrete-time integrator, finite limits containing -1 for 1/(s(s+1)), 1/s^2 and a delayed integrator, critical_point, circles, per-subplot limits in the Plots recipe, user-provided xlims, a BigFloat system with and without user-provided limits), and Makie tests in test_makie_plots.jl (axis limits equal nyquist_limits for each axis, user-provided limits respected, no default limits for a BigFloat system). The RobustAndOptimalControl reproducer with Particles was verified manually, with and without user-provided limits.

Run locally (Julia 1.12.7): test_plots, test_makie_plots, test_freqresp, test_delayed_systems: 363 passed, 0 failed.

---

Host: demeter6, Claude Code session session_01HigNnun4aJKMoWmY7wyULn, working directory ~/.julia/dev/ControlSystems/.claude/worktrees/agent-a8a6bc7868b047e82

🤖 Generated with Claude Code

https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn

[JuliaControlBot]

This is an automated message.
Plots were compared to references. 4/11 images have changed, see differences below.
After pulling this PR, please update the reference images by creating a PR to ControlExamplePlots.jl here.

Difference	Reference Image	New Image
✔️ 0.0
✔️ 0.0
✔️ 0.0
✔️ 0.013

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.09%. Comparing base (f8c5db4) to head (4da51d7).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1084      +/-   ##
==========================================
+ Coverage   91.98%   92.09%   +0.10%     
==========================================
  Files          42       42              
  Lines        5976     6056      +80     
==========================================
+ Hits         5497     5577      +80     
  Misses        479      479              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[JuliaControlBot]

This is an automated message.
Plots were compared to references. 4/11 images have changed, see differences below.
After pulling this PR, please update the reference images by creating a PR to ControlExamplePlots.jl here.

Difference	Reference Image	New Image
✔️ 0.0
✔️ 0.0
✔️ 0.0
✔️ 0.013