Fix default margin frequency ranges and missing phase guides

[dongxuelian2]

Closes #771.

The original system, tf([1], [1, 13, 40, 0]), has a unity-gain crossover near 0.025 rad/s, below the original default grid's 0.1 rad/s lower bound. Select a margin-specific default range using pole/zero features and limiting gains while preserving the original Bode samples. The example now reports a phase margin of 89.53446° at 0.02500054 rad/s in both plotting backends.

The limiting-gain checks use type-aware isapprox comparisons against unity, so rounding error in unit DC or high-frequency gain does not force extension into noise crossings. Near-origin classification calls count_eigval_multiplicity and adds a relative tolerance locally to margin-grid selection; the shared utility and its defaults are unchanged. Use that classification consistently for features, limiting gains, and endpoint evaluations. This keeps the double-mass model's default range at 0.01–1000 rad/s with 300 points while retaining the 1e-6 resonance. Nonrational delay/interconnected models retain their existing range selection because their augmented channel gains do not define the external limiting gain.

Plots and Makie share a backend-independent helper for critical-phase guides, including a fallback when no finite phase margin exists. The helper does not add a margin or a crossover frequency. Keep fullPhase aligned with pm and wpm when Plots retains only five phase margins. Backend tests use stable smoke/regression checks, with guide selection tested directly. Reuse pole/zero features to avoid repeated zpkdata calls and document the default range in marginplot.

Explicit frequency vectors are used without extending their range; ordinary Bode grids and public APIs remain unchanged. Discrete default margin grids are capped at Nyquist.

The newly detected discrete-integrator crossover also exposed a pre-existing phase-correction defect: unstable discrete poles were counted using the continuous-time right-half-plane criterion. Use the unit circle for discrete systems, with regressions for an integrator and additional stable/unstable positive poles. This corrects the phase margin from approximately −270° to 90° for the stable integrator case; the same erroneous explicit-grid result was reproduced on untouched master.

Validation with Julia 1.10.12 on Windows (Plots 1.41.7, Makie 0.24.15/CairoMakie 0.15.15):

• test_analysis.jl: 226/226 passed.
• test_plots.jl: 31/31 passed.
• test_makie_plots.jl: 103/103 passed.
• test_freqresp.jl: 101 passed, 1 allocation-limit failure at line 69. Untouched master has the identical failure: 123,504 bytes against a 68,006.4-byte limit.
• All reviewer reproducers completed, including the explicit-grid Tss truncation case.
• 180 additional adversarial checks passed for tolerance boundaries, Float32/Float64/BigFloat, unchanged master Bode grids, Nyquist bounds, explicit vectors, discrete phase correction, and guide-only behavior. Both backends were rendered and visually inspected.
• Scoped JuliaFormatter, parsing, git diff --check, and final diff review passed. The full package suite was not run.

No package dependencies were added.

[baggepinnen]

Thank you for the PR. The original example from #771 is resolved by margin and by the Plots recipe for marginplot, and the modified test files pass locally (Linux, Julia 1.13; the allocation test in test_freqresp.jl also passes there). The following issues were found.

1. marginplot throws for systems whose DC gain equals one up to rounding error

The low-frequency extension compares dc != 1 exactly. When the DC gain differs from one only by rounding error, the loop lowers w1 until the magnitude lies on the same side of one as dc, and the extended range then contains crossings of |G| = 1 that are caused by rounding error alone.

Tss = feedback(ss(tf(1, [1, 1, 0])) * ss(pid(1.0, 1.0, 0.1; Tf=0.01)))
margin(Tss; allMargins=true).wpm[] # 19 crossings, 18 of them near 1e-8 rad/s with pm = 180°; master: 1 crossing
marginplot(Tss)                    # DimensionMismatch at plotting.jl:876

L = ss(tf([1, 1], [1, 10])) * ss(10.0) # unit DC gain
margin(L; allMargins=true).wpm[]       # 15 crossings between 1e-8 and 2e-8 rad/s

The exception originates from a defect that is already present on master: when more than five margins are found, the recipe truncates pm and wpm but not fullPhase, so fullPhase .- pm cannot be broadcast. On master the defect is reachable only with an explicit frequency vector (marginplot(Tss, exp10.(range(-8, 4, length=1000))) also throws on master), whereas this PR makes it reachable with the default arguments. Comparing dc and hf against one with a tolerance, and truncating fullPhase together with pm, would address both. The high-frequency check hf != 1 has the same structure.

2. Poles that are located at the origin up to rounding error determine the lower frequency bound

For Val{:margin} the fzp .> -4 filter is removed, so poles that are zero up to rounding error now determine the lower bound of the frequency grid. DemoSystems.double_mass_model() has a pole at -1.1e-14 (a numerical integrator); its margin grid spans 1e-16 to 1e3 rad/s with 1138 points, compared to 1e-2 to 1e3 for the Bode grid, and marginplot displays 20 decades of frequency with a phase transition from 0° to -90° near 1e-14 rad/s.

Outside the range of the poles and zeros the magnitude is monotone, so the gain-based extension already locates crossings below that range; removing the filter is required only for the test case with a resonance at 1e-6 rad/s and gain 1e-13. A criterion that distinguishes eigenvalues located at the origin up to rounding error from slow dynamics would retain both behaviors. count_eigval_multiplicity (used by count_integrators and integrator_excess_with_tol) is the existing utility for this classification, but with its default tolerance eps(maximum(abs, p)) it does not classify the pole at -1.1e-14 as an integrator (count_integrators(double_mass_model()) == 0). A threshold relative to the largest pole magnitude, such as sqrt(eps())*maximum(abs, p) (the tolerance that utility applies for multiplicity two), removes this pole while keeping the resonance at 1e-6 rad/s in the test case. The same classification could replace the exact iszero test that determines the number of integrators and differentiators in the computation of dc.

3. The Makie recipe is not updated

ext/ControlSystemsBaseMakieExt.jl still calls _processfreqplot(Val{:bode}(), ...) (lines 386-387) and draws a phase guide only for each phase margin (lines 491-494). CSMakie.marginplot(tf(1, [1, 13, 40, 0])) therefore still uses a frequency grid that starts at 0.1 rad/s, has the title Pm: [] ωpm: [], and contains no phase guide. Both changes should be applied to the Makie recipe as well, with a corresponding test in test_makie_plots.jl.

4. The new plot tests depend on internals of Plots.jl

The recipe itself uses only RecipesBase, but the tests access p[2].series_list and p[1].series_list[1][:x] (an internal field of Plots.Subplot), as well as the representation into which Plots expands an :hline series (series[:y]) and the stored value of series[:linestyle]. The existing tests in test_plots.jl are smoke tests only, so a change in a Plots release could cause these tests to fail without any change in this package. Extracting the selection of the phase guides (from fullPhase, pm and phasedata) into a helper function, testing that function directly, and calling it from both the Plots and the Makie recipe would remove this dependency and keep the two recipes consistent.

Minor

• _default_freq_vector(systems, Val{:margin}()) calls itself with Val{:bode}(), so zpkdata is computed twice for each system.
• The docstring of marginplot is not updated to describe the new default frequency range.

[baggepinnen]

I'll look into adjusting the default tolerance of count_eigval_multiplicity, I have to do a lot of downstream testing before changing such a fundamental utility

[baggepinnen]

Commit 2811b2f addresses the points of the previous review. On top of it, two commits were pushed to this branch:

• A merge of master, which contains Make the tolerance of the integrator count independent of the time unit #1077. That PR changed the default radius of count_eigval_multiplicity to (100ϵ)^(1/i) * scale and added the keyword scale.
• 1ff5e00, which changes _margin_nonintegrators as follows:
  • The relative floor sqrt(eps)*scale is removed. It was suggested in the previous review because the former tolerance did not classify the pole of DemoSystems.double_mass_model() at -1.1e-14 as an integrator. Since Make the tolerance of the integrator count independent of the time unit #1077, count_eigval_multiplicity classifies this pole without the floor. The floor made the classification for the frequency grid differ from the one used for the phase adjustment in sisomargin.
  • The zeros are classified with scale = max(maximum(abs, p), maximum(abs, z)), as in integrator_excess_with_tol. Previously, a zero in the origin that is not computed exactly defined the scale itself and determined the lower bound of the grid. For balreal(ss(s^2/((s+1)*(s^2 + 0.2s + 4))))[1], the zeros are computed as ±1.4e-8. The margin grid started at 9.6e-11 rad/s, and margin(...; allMargins=true) reported 61 gain margins between 9.8e-11 and 2.8e-8 rad/s, at which rounding error determines the phase. This system and balreal(ss(s/(s+1)^2))[1] are added as tests.

The full test suite of ControlSystemsBase passes (Linux, Julia 1.13).

The remaining gain margin at ω = 0 with gm ≈ 1e16 for systems with a zero in the origin is present on master as well. It will be addressed in a separate PR.

[codecov]

Codecov Report

❌ Patch coverage is 98.93617% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 91.91%. Comparing base (59c065a) to head (a7cba57).
⚠️ Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
lib/ControlSystemsBase/src/plotting.jl	93.33%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1076      +/-   ##
==========================================
+ Coverage   91.73%   91.91%   +0.18%     
==========================================
  Files          42       42              
  Lines        5865     5935      +70     
==========================================
+ Hits         5380     5455      +75     
+ Misses        485      480       -5     

☔ 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.

[baggepinnen]

Two further commits were pushed to this branch:

• A merge of master, which contains Do not report a gain margin at zero frequency for systems with a zero in the origin #1078 (removal of the gain crossover at ω = 0 for systems with a zero in the origin) and Root locus: relative step-size control, improper systems, default gain and plot limits #1079.
• a7cba57, which addresses two problems of MIMO systems.

One default frequency vector per channel. margin(sys) and marginplot (Plots and Makie) used one default frequency vector for all channels, and marginplot also for all systems. The extension of the frequency range required by a slow pole in one channel was therefore applied to every channel, including channels for which _bounds_and_features deliberately does not extend the range because their gain at low frequencies equals one. Example, with a non-diagonal realization of the first block:

s = tf("s")
sys = append(ss(1/((s + a)*(s + 1))), ss(1e4/(s + 1e4)))

Number of phase crossovers of channel (2, 2), at which rounding error determines whether the gain exceeds one:

	a = 1e-12	a = 1e-6
frequency vector shared by all channels (previous commit)	600	241
default frequency vector of master	1	1
frequency vector of the channel (this commit)	1	1

The remaining crossover is located at ω = 0, where the gain equals one exactly. Each channel now uses its own default frequency vector, and margin(sys) agrees with margin(sys[i, j]) for every channel. In marginplot, the curves of different channels and systems may therefore cover different frequency ranges. The magnitude is now computed per channel; previously, the Plots recipe applied the conversion to dB to the magnitudes of all channels once per channel.

Structurally zero channels. This problem is also present on master. In the example above, the off-diagonal channels are structurally zero, and for a = 1e-12, margin(sys; allMargins=true) reported 334 gain crossovers for channel (2, 1), both with the default frequency vector of master and of this PR. Two changes in sisomargin address this:

• The state variables that are structurally uncontrollable or unobservable are removed (_sminreal), as zpkdata does when the default frequency vector is computed. Otherwise, the Hessenberg reduction of A mixes the state variables, and the computed frequency response of a structurally zero channel is rounding error. Of 300 random state-space systems with a structurally zero transfer function, 29 obtained gain crossovers without this reduction, with gain margins such as 3.9e16, 1.7e17 and 2.3e18. After the reduction, the state dimension of such a channel is zero and its frequency response is exactly zero. As a consequence, the poles and zeros in the origin are classified on the same realization in sisomargin and in _margin_nonintegrators.
• Gain crossovers at which the computed frequency response is exactly zero are removed, since the Nyquist curve passes through the origin there and the computed phase is determined by the signs of the real and imaginary parts of zero. For example, channel (2, 1) of -ss([1/(s+1) 0; 0 1/(s+2)]) has the gain -0.0, its computed phase is 180°, and every frequency of the grid was reported as a gain crossover.

The test set Margins of MIMO systems covers both problems. Without the per-channel frequency vectors, 6 of its 28 assertions fail; without the structural reduction, 3 fail; without the removal of crossovers at which the frequency response is zero, 4 fail. marginplot of the example system, of a vector of systems and with the dB scale is tested for both backends.

The full test suite of ControlSystemsBase passes (Linux, Julia 1.13): 5773 passed, 6 broken. In a first run, the gang-of-four test at test_analysis.jl:609, which uses unseeded ssrand systems, failed with 3.8e-9 < 1e-9; this test does not call the changed functions and fails in 1 of 5000 random draws independently of this PR.