Repository navigation
Compute signed phase margins and determine closed-loop stability in delaymargin - #1094
baggepinnen wants to merge 3 commits into
Conversation
|
This is an automated message.
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1094 +/- ##
==========================================
+ Coverage 92.29% 92.49% +0.20%
==========================================
Files 42 42
Lines 6177 6372 +195
==========================================
+ Hits 5701 5894 +193
- Misses 476 478 +2 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Review: signed phase margins and the new
|
| Loop | Parent | This PR | Ground truth (winding count) |
|---|---|---|---|
2.7(s+1)^2/((s-ε)^2(s+ε)), ε = 1e-4, Float64 and Float32 |
pm 413°, dm 2.403 | pm 53.1°, dm 0.3091 | dm 0.3091 ✔️ |
4(s+1)^2/((s-1)^2(0.1s+1)) (2 unstable OL poles) |
pm 458.6°, dm 0.2066 | pm 98.6°, dm 0.0444 | 0 RHP poles at τ = 0.0444, 2 at τ = 0.045 → 0.0444 ✔️ |
(s+0.1)^2/(s^3(0.01s+1)) |
pm 78.1° or 438.1° depending on adjust_phase_start |
78.1° for both | ✔️ |
The second row is a real-world hazard that the parent hides: it reports 5× more delay margin than the loop has. Margins no longer depend on how integrator_excess / count_integrators classify poles near the origin, which removes this function from that blast radius. 👍
1. delaymargin returns a positive, finite margin for loops that are already unstable
| Loop (closed loop unstable without delay) | pm | Parent dm | This PR dm |
|---|---|---|---|
27/(s+1)^3 |
−31.6° | −0.195 | 2.025 |
0.999/(s(s^2+0.002s+1)) |
−89.8° | −1.182 | 3.558 |
2(s^2+0.01s+1.0201)/(s(s^2+0.01s+1)(s+1)) |
[−41.7, 41.4, 39.4] | −0.723 | 0.5645 |
mod(ϕ, 360) turns −31.6° into 328.4° of lag, so the result looks like a healthy margin. The docstring now states the stability assumption, but people tend to call delaymargin precisely to learn whether a loop is OK. The parent's negative value was at least a usable signal. Suggestion: for rational G, check isstable(feedback(G)) (cheap compared with the frequency sweep) and return 0 (or keep a negative value) when it fails. Or at least @warn.
2. delaymargin is discontinuous at pm = 0
L = 1/s · e^{-sπ/2} has pm = −0.033° (barely unstable), and the result stays the same for gains 0.999 to 1.001:
- parent: dm = −0.00058 (≈ 0, correct)
- this PR: dm = 6.28 = 2π/ωc, a whole period of delay
Any loop whose computed pm lands slightly below zero through rounding jumps from about 0 to about 2π/ωc. That's the same root cause as (1), and the same stability check fixes it. The ω = 0 branch has the same edge case: iszero(lag) fails when lag comes out of mod as 360 - eps.
3. allMargins=false can report a positive pm on an unstable loop
For the dipole loop above, the three crossovers are [−41.7°, +41.4°, +39.4°]. The closed loop is unstable, and margin(L).pm is now +39.35 (the parent returned −41.7). "Closest to zero" is defensible geometrically, but margin(L).pm[] > 0 is a common "is this OK" check, and it now passes on an unstable loop. Please at least document this, or break ties toward a negative margin when one exists within a few degrees.
4. Smaller points
- Plot recipes:
plotting.jland the Makie extension still choose the "5 smallest" margins withsortperm(pm), whilemarginnow usesabs. Usesortperm(pm; by=abs)so the plot shows the margin thatmargin()returns. - Semver: the signed convention changes values users may have stored or asserted on (−320 → 40, 458 → 98). That looks like a minor version bump plus a NEWS entry, not a patch.
- Test oracle: the new tests check delay margins with
pade(1.05τ, 6). On the 2-unstable-pole loop above, an 8th-order Padé still reports a stable loop at τ = 0.21 (true limit 0.0444), and on1e7/((s+1e-4)(s+1e4))it reports stable at 1.05τ even though the winding count shows 2 RHP poles. A small winding-number helper (sweep1 + L(jω)e^{-jωτ}and count encirclements, plus open-loop RHP poles) is a more trustworthy oracle for tests like these.
Suggested extra tests (all pass on this PR unless noted)
s = tf("s")
# OL-unstable: delay margin must match the exact limit, not Padé
L = 4(s+1)^2/((s-1)^2*(0.1s+1)); @test delaymargin(L) ≈ 0.0444 rtol=1e-2
# adjust_phase_start must not change margins (parent: 78 vs 438)
L = (s+0.1)^2/(s^3*(0.01s+1))
@test margin(L; adjust_phase_start=true).pm[] == margin(L; adjust_phase_start=false).pm[]
# Float32 agrees with Float64
@test margin(tf(Float32[2.7,5.4,2.7], Float32[1,-1f-4,-1f-8,1f-12])).pm[] ≈ 53.138 atol=1e-2
# Wide scaling
@test delaymargin(1e7/((s+1e-4)*(s+1e4))) ≈ 0.0014786 rtol=1e-3
# Near-cancelling integrator/zero
@test margin((s+1e-9)/(s^2*(s+1))).pm[] ≈ 51.823 atol=1e-2
# Discrete pole near z = 1
@test margin(tf([0.01], [1, -(1-1e-9)], 0.01)).pm[] ≈ 89.713 atol=1e-2
# FAIL today, per (1)/(2): unstable or marginal loops must not report a positive delay margin
@test delaymargin(27/(s+1)^3) <= 0
@test delaymargin(tf(1,[1,0])*delay(π/2)) < 0.01🤖 Posted by an AI agent · harness: Claude Code 2.1.316 · model: claude-opus-5-5
Review: unreviewed
Session: local Claude Code session 796e0e3d-3934-4f65-be25-83ae79ce8750
…y_intervals
The delay margin was computed from the phase margins under the assumption
of a nominally stable closed loop, and returned a positive finite value for
unstable closed loops.
The closed-loop poles of 1 + G(s)e^{-sθ} reach the imaginary axis only at
the gain crossovers of G, for total delays θ_k = (lag + 2πk)/ω, and the
sign of the real part of ds/dθ is -sign(ω d|G|/dω). Starting from the
number of unstable poles of feedback(G), the number of unstable poles for
every delay follows by counting the crossings. This applies to rational
loops and to loops G(s)e^{-sT} with a single input or output delay.
- delaymargin returns NaN for an unstable closed loop, 0 for a closed loop
with poles on the stability boundary, and the first loss of stability
otherwise. Other delay systems retain the previous computation with a
warning.
- delay_stability_intervals returns all intervals of additional delay with
an asymptotically stable closed loop.
- Gain crossover frequencies are refined by Newton iteration, which
corrects the delay margin of the four-crossover test loop from 0.15874
to 0.15644.
- Marginplots select the displayed phase margins by absolute value.
- Tests use an argument-principle count of the exact delayed loop as the
reference instead of Padé approximations.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
baggepinnen
changed the title
|
This is an automated message.
|
- Locate gain crossovers by bisection within the brackets of the default frequency grid, with an analytic derivative and a residual check, instead of Newton iteration, which left the crossover of narrow resonances. - Discrete time: count unstable poles for integer delays exactly. In addition to the interior crossings, each sample of delay adds one unstable pole if |G(-1)| > 1. delay_stability_intervals returns ranges of stable integer delays for discrete-time systems. - Mirror the frequency grid to negative frequencies for systems with complex coefficients. - Evaluate G(0) (G(1) in discrete time) directly to detect a closed-loop pole at the origin for every delay. - Scale the pole tolerance by the norm of the closed-loop system matrix, and snap small positive lags to zero. - Define stability as exponential stability, so that loops with |L(j∞)| = 1 have delay margin 0. - Harden the argument-principle reference in the tests, add a direct discrete-time reference, test all review counterexamples, and check that no stable interval is omitted. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
This is an automated message.
|
Phase margins depend on the chosen unwrapped Bode phase. For
2.7(s+1)^2/((s-ε)^2*(s+ε)), nonzeroεcan produce a phase margin near 413° and a delay margin near 2.403 s, although the relevant values are approximately 53° and 0.309 s.Phase margins
Signed phase margins in
[-180°, 180°)are computed directly from the frequency response at each gain crossover. WithallMargins=false, the margin closest to zero is returned. Integrator-based phase adjustment applies only to the unwrapped phase returned for plotting, so the integrator classification and the initial phase branch do not affect the margins.marginplotselects the displayed margins by absolute value, consistently withmargin. The docstring states that the sign of a phase margin does not determine closed-loop stability.The signed convention preserves negative margins and changes an existing unwrapped result near −320° to its equivalent signed value near +40°. This changes returned values, so the release should be a minor version.
Delay margins and closed-loop stability
A delay margin computed from the phase margins alone requires the assumption of a stable closed loop without additional delay, and returned a positive finite delay margin for unstable closed loops, e.g., 2.03 for
27/(s+1)^3. This cannot be resolved from the frequency response alone, since the number of open-loop unstable poles is required.The closed-loop poles of
1 + G(s)e^{-sθ}are located on the imaginary axis only at gain crossoversωofG, for total delaysθ_k = (mod(arg G(jω) + π, 2π) + 2πk)/ω, and implicit differentiation of the characteristic equation givessign Re(ds/dθ) = -sign(ω d|G(jω)|/dω). Starting from the number of unstable poles offeedback(G), the number of unstable closed-loop poles for every delay follows by counting crossings. This applies to rational loops and to loopsG(s)e^{-sT}with a single input or output delay, for which the nominal delayTis a point on the same delay axis.delaymarginreturnsNaNif the closed loop without additional delay is unstable,0if it has poles on the stability boundary (e.g.,1/s*delay(π/2)), and otherwise the first additional delay at which stability is lost. Stability is defined as exponential stability, which is equivalent to asymptotic stability for|L(j∞)| < 1; loops with|L(j∞)| ≥ 1have delay margin0. For other delay systems (internal delays), the previous computation is retained under the stated assumption, with a warning.delay_stability_intervalsreturns all intervals of additional delay for which the closed loop is stable, including stability switches at larger delays and delays that stabilize an unstable closed loop. For discrete-time systems, it returns ranges of stable integer sample delays.|G(-1)| > 1, e.g.,-1.2/(z+0.5)is unstable with one sample of delay although its phase-based margin is 3.86 samples.marginfrequency grid, which is extended to negative frequencies for systems with complex coefficients. The delay margin of the four-crossover test loop was 0.15874 from the grid-interpolated crossover; the correct value is 0.15644.Fixes ControlSystems.jl#1093. Addresses the review comment and a subsequent review by a Codex agent (counterexamples for near-zero lags, slow unstable poles, discrete Nyquist gains, complex coefficients, zero-frequency roots, narrow resonances, and the test reference), all of which are included as tests.
Validation
The tests verify delay margins and stability intervals against an argument-principle count of the unstable roots of the exact delayed characteristic equation, instead of Padé approximations, which reported stability beyond the delay margin for open-loop unstable and widely scaled loops. The discrete-time results are compared with directly computed closed-loop poles for every integer delay. Covered cases: transfer-function, zero-pole-gain and state-space representations, open-loop unstable plants, unstable closed loops, a delay-stabilized closed loop, stability switches, input and output delays, marginal loops, a neutral-type loop, discrete-time loops, Float32, and internal delays. The full ControlSystemsBase test suite passes (7486 passed, 6 broken as on master) with Julia 1.13.2 on Linux.
Host:
demeter2; Codex session:01a15eb9-1234-7a03-a4be-7f8814d35ca2(initial commit); Claude Code session:285e5228-c064-4641-add1-6fae71325fa0(delay margin); working directory:~/.julia/dev/ControlSystems/.claude/worktrees/crystalline-crunching-gem.🤖 Generated with Claude Code