Root locus: relative step-size control, improper systems, default gain and plot limits

[baggepinnen]

Problem

rlocus on the following system did not terminate in reasonable time and exhausted the memory:

sys = ss([-0.5 4000; -66.67 -875], [0; 6000;;], [9.5493 0], 0)
rlocus(sys)

The adaptive step-size control in getpoles accepted a step if the sum of the absolute distances moved by the poles was below tol (1e-2 per pole). The number of steps was thus proportional to the distance the poles travel. For this system the open-loop poles have magnitude ≈ 500 and at the default K = 500 the closed-loop poles are located at ≈ -438 ± 3.4e5im, which requires on the order of 3e7 steps (1.4e6 steps already for K = 1). Each step allocated a vector, and the result was finally assembled with hcat(poleout_list...), splatting all steps as arguments.

Change

• The displacement of each pole in a step is divided by the magnitude of the pole before the step, bounded from below by 1e-3 times the smallest nonzero magnitude among the open-loop poles, the zeros and the closed-loop poles at the final gain (magnitudes below sqrt(eps) times the largest are regarded as zero, since a root at the origin is computed with a rounding error). The number of steps then grows logarithmically with the distance travelled and is independent of the time unit of the system. This applies to both the scalar-gain and the matrix-gain variants of getpoles.
• The pole matrices are formed with reduce(hcat, ...) instead of splatting.
• The returned pole matrices are formed with transpose instead of '. The adjoint conjugated all returned poles; this was not visible in plots since the locus is symmetric about the real axis.
• The meaning of tol and initial_stepsize is documented in the docstring of rlocus.

Effect

System	steps before	steps after
system above, K = 500	≈ 3e7 (not completed)	661 (7 ms)
tf([1/5,2],[1,1,1]), K = 200	2954	393
zpk example of #740, K = 200	453	338
tf(1, [1,0,0]), K = 200	1366	784
tf(1, [1,1001,1000]), K = 200	108	105
double_mass_model, state feedback	86	242

The double-mass model requires more steps than before since it has a pole at the origin, which is now resolved on a relative scale as it departs from the origin.

Tests are added for the system above, including a check that the number of steps and the (scaled) poles are unchanged when the time unit is changed from s to ms.

Improper transfer functions

For an improper transfer function P/Q, the closed-loop system has degree(P) poles for k > 0 while the open-loop system has degree(Q) poles. Both SISO methods of getpoles sized the pole matrices from the denominator degree and failed with a DimensionMismatch. They are now sized from the larger of the two degrees.

The adaptive method starts the locus at a gain k0 > 0 (K[1] = k0) instead of approximating k = 0 by k = eps(), which placed the additional poles at a magnitude of ≈ 1e16. k0 follows from Rouché's theorem: if k|P(s)| < |Q(s)| on a closed contour, Q + kP has as many roots inside the contour as Q. The contour consists of a circle of radius ten times the largest magnitude among the open-loop poles and zeros and of circles of radius 0.01|p| around each open-loop pole p, so that at k0 the additional poles are outside the large circle and the remaining poles are close to the open-loop poles. The method with a vector of gains retains the approximation k = eps() for k = 0.

The Plots and Makie recipes mark the open-loop poles by poles(sys) rather than by the first row of the root matrix, which for an improper system also contains the additional poles.

Default gain and plot limits

The default maximum gain K = 500 and the plot limits of ±50 in the Plots recipe depended on the time unit and on the gain of the system. For the system above, the locus was entirely outside the plotted region.

• Default gain (rlocus(P), rlocusplot(P), CSMakie.rlocusplot(P), K = nothing): with ω0 the largest magnitude among the open-loop poles and zeros (at least 1 in discrete time), the default gain is the maximum of |Q(s)/P(s)| over a circle of radius 10ω0 (if degree(Q) > degree(P)) and over circles of radius 0.01max(|z|, 0.01ω0) around each zero z. By Rouché's theorem, at this gain the poles that tend to infinity are located outside the large circle and the remaining poles inside the small circles. The default gain is invariant to a change of the time unit and inversely proportional to the gain of the system. This changes the output of rlocus(P) without a gain argument.
• Default plot limits (ControlSystemsBase.rlocus_limits, used by both the Plots and the Makie recipe): the limits contain the open-loop poles and zeros, the origin, the closed-loop poles at the final gain if the gain is a matrix, and all points of the locus with a magnitude of at most 2ω0. They are symmetric about the real axis and enlarged by 10 %. For discrete-time systems, the limits contain the unit circle, which is drawn as a dashed line.

Automatic limits are not a substitute for these limits: since the default gain places the poles that tend to infinity at a magnitude of at least 10ω0, automatic limits are determined by these branches rather than by the open-loop poles and zeros.

System	default K	max |p| at default K	max |p| at K = 500	ω0
system above	0.137	5640	3.4e5	517
tf([1/5,2],[1,1,1])	4650	920	90	10
1/(s(s+1)(s+2))	9240	22	9.0	2
(-s+1)/(s²+s+1)	303	301	498	1
zpk example of #740	2.5e5	66	11.5	6
zpk([-0.5], [1, 0.6], 0.1, 0.1) (discrete)	3330	330	48	1

---

Host: demeter6, session 6ab4068b-a603-44bd-8c0f-82dee93c0177, working directory ~/.julia/dev/ControlSystems/.claude/worktrees/generic-pondering-babbage

🤖 Generated with Claude Code

[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.0

[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.0

[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.0

[codecov]

Codecov Report

❌ Patch coverage is 94.28571% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.73%. Comparing base (94c867a) to head (ef110c5).

Files with missing lines	Patch %	Lines
lib/ControlSystemsBase/src/root_locus.jl	95.87%	4 Missing ⚠️
...ntrolSystemsBase/ext/ControlSystemsBaseMakieExt.jl	75.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1079      +/-   ##
==========================================
+ Coverage   91.70%   91.73%   +0.02%     
==========================================
  Files          42       42              
  Lines        5796     5865      +69     
==========================================
+ Hits         5315     5380      +65     
- Misses        481      485       +4     

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