margin and delaymargin are 360° off when rounding splits a multiple integrator into poles off the origin #1093

Expected. margin gives the phase margin in [0°, 360°). delaymargin gives the delay that makes the loop unstable. For L = 2.7(s+1)²/s³, these are 53.130° at 3 rad/s and 0.30910 s. The result is the same when rounding moves the three integrator poles by ε = 1e-6 or 1e-4. Linearized models often have such poles.

Actual. For ε = 1e-6 and ε = 1e-4, margin gives 413.14°. This is 360° too large. delaymargin gives 2.403 s. This is almost 8 times the correct value. A sixth-order Padé check shows that the loop becomes unstable between 0.95τ and 1.05τ with τ = 0.309 s, so the true delay margin is unchanged.

Likely cause (code reading, not tested): sisomargin corrects the phase only when integrator_excess is not zero. Since #1090, a cluster counts as integrators only when its polynomial is within about 100·eps of zⁿ, so these clusters are deliberately not counted. The phase then starts 360° off. delaymargin (analysis.jl:723) divides that phase margin by the crossover frequency without wrapping it to [0°, 360°). A fix that wraps the phase margin in margin and delaymargin would not depend on the integrator count. Related: #1081.

Reproducer (ControlSystems.jl master 2737b53; ControlSystemsBase 1.25.1):

using ControlSystemsBase
println("ControlSystemsBase ", pkgversion(ControlSystemsBase))
s = tf("s")
# L0 = 2.7 (s+1)^2 / s^3: |L0(3i)| = 1, ∠L0(3i) = -270° + 2 atan(3), so
# PM = 2 atan(3) - 90° = 53.130°, delay margin = deg2rad(PM) / 3 = 0.30910 s.
pm0 = 2rad2deg(atan(3)) - 90
println("closed form: PM = ", pm0, "°, delay margin = ", deg2rad(pm0) / 3, " s")
for ε in (0.0, 1e-6, 1e-4)
    # The triple integrator moved by ε, as rounding moves it in a linearized model
    L = ss(2.7(s + 1)^2 / ((s - ε)^2 * (s + ε)))
    m = margin(L)
    println("ε = $ε: closed loop stable: ", isstable(feedback(L)), ", poles(L) = ", round.(poles(L), sigdigits = 3))
    println("  margin:      PM = ", m.pm[1], "° at ", m.wpm[1], " rad/s")
    println("  delaymargin: ", delaymargin(L), " s")
    τ = deg2rad(mod(m.pm[1], 360)) / m.wpm[1] # wrapped phase margin over crossover frequency
    println("  Padé(6) check, closed loop stable with delay 0.95τ: ", isstable(feedback(L * pade(0.95τ, 6))),
        ", with 1.05τ: ", isstable(feedback(L * pade(1.05τ, 6))), " (τ = $τ s)")
end

Output:

ControlSystemsBase 1.25.1
closed form: PM = 53.13010235415598°, delay margin = 0.3090984060005374 s
ε = 0.0: closed loop stable: true, poles(L) = ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]
  margin:      PM = 53.1306178712631° at 3.000044987961516 rad/s
  delaymargin: 0.30909676994629687 s
  Padé(6) check, closed loop stable with delay 0.95τ: true, with 1.05τ: false (τ = 0.30909676994629687 s)
ε = 1.0e-6: closed loop stable: true, poles(L) = ComplexF64[-1.0e-6 + 0.0im, 1.0e-6 + 2.1e-14im, 1.0e-6 - 2.1e-14im]
  margin:      PM = 413.140879627855° at 3.000942427073555 rad/s
  delaymargin: 2.4028013863375914 s
  Padé(6) check, closed loop stable with delay 0.95τ: true, with 1.05τ: false (τ = 0.30906401553878105 s)
ε = 0.0001: closed loop stable: true, poles(L) = ComplexF64[-0.0001 + 0.0im, 0.0001 + 0.0im, 0.0001 + 0.0im]
  margin:      PM = 413.13755576698037° at 3.000817254571634 rad/s
  delaymargin: 2.402882281748546 s
  Padé(6) check, closed loop stable with delay 0.95τ: true, with 1.05τ: false (τ = 0.30905757529320055 s)

Workaround. Wrap the result: pm = mod(m.pm[1], 360) and τ = deg2rad(pm) / m.wpm[1]. This gives the correct values above.

Environment: Julia 1.12.7 (macOS arm64).

🤖 Posted by an AI agent — harness: Claude Code 2.1.295 (desktop app) · model: claude-opus-5-5
Session: local Claude Code session 98c0c801-f493-4ae7-bb02-47ce9ed80bd5