MinimalRealization over-reduces: 3-state CL system with hidden integrator drops to order 0

[jamestjat]

Summary

(*System).MinimalRealization() reduces a 3-state closed-loop system (PI controller + first-order plant, series-composed with the plant) down to order 0 (a zero static gain), discarding the real dynamics. The original system has one hidden integrator mode from the PI controller that should be removed, leaving a 2-state minimal realization — not an empty one.

Discovered while tracking a step-response instability in a PACE→MDL bridge: the closed-loop Feedback(ctrl, plant, -1) retains the PI integrator as an uncontrollable mode, which IsStable() then flags. MinimalRealization() was the natural fix, but its current output is unusable on this case.

Environment

• controlsys master as of 2026-04-15
• gonum v0.x
• Go 1.22, Windows

Reproduction

package main

import (
	"fmt"

	"github.com/jamestjsp/controlsys"
	"gonum.org/v1/gonum/mat"
)

func main() {
	// PI controller: G_c(s) = K(βs+1)/(τ₁s), K=3.333, τ₁=β=11998.8 s
	// CCF: A=[0], B=[1], C=[K/τ₁], D=[Kβ/τ₁] = D=K (since β=τ₁)
	K, tau1, beta := 3.333, 11998.8, 11998.8
	ctrl, _ := controlsys.New(
		mat.NewDense(1, 1, []float64{0}),
		mat.NewDense(1, 1, []float64{1}),
		mat.NewDense(1, 1, []float64{K / tau1}),
		mat.NewDense(1, 1, []float64{K * beta / tau1}),
		0,
	)

	// First-order plant: G_p(s) = 0.56 / (66240 s + 1)
	Kp, tauP := 0.56, 66240.0
	plant, _ := controlsys.New(
		mat.NewDense(1, 1, []float64{-1.0 / tauP}),
		mat.NewDense(1, 1, []float64{1.0 / tauP}),
		mat.NewDense(1, 1, []float64{Kp}),
		mat.NewDense(1, 1, []float64{0}),
		0,
	)

	// Closed loop SP→OP, then series with plant to get SP→PV.
	cl, _ := controlsys.Feedback(ctrl, plant, -1)
	clToPv, _ := controlsys.Series(cl, plant)

	r, c := clToPv.A.Dims()
	fmt.Printf("clToPv A dims=%dx%d\n", r, c)
	stable, _ := clToPv.IsStable()
	fmt.Printf("clToPv IsStable=%v  (expected: true for SP→PV of stable PI loop)\n", stable)

	mr, err := clToPv.MinimalRealization()
	if err != nil {
		fmt.Printf("MinimalRealization err: %v\n", err)
		return
	}
	fmt.Printf("MinimalRealization: Order=%d\n", mr.Order)
	if mr.Sys != nil {
		rr, cc := mr.Sys.A.Dims()
		fmt.Printf("  min A dims=%dx%d\n", rr, cc)
		fmt.Printf("  min A=%v\n", mat.Formatted(mr.Sys.A))
		mrStable, _ := mr.Sys.IsStable()
		fmt.Printf("  min IsStable=%v\n", mrStable)
	}
}

Actual output

clToPv A dims=3x3
clToPv IsStable=false
MinimalRealization: Order=0
  min A dims=1x1
  min A=[0]
  min IsStable=false

Expected output

MinimalRealization should return a 2-state stable system equivalent in I/O to the closed-loop G_c·G_p/(1+G_c·G_p):

• Minimal order: 2
• IsStable=true
• DC gain ≈ 1 (unity tracking after PI settles)

Analytically the closed-loop characteristic polynomial is
`τ_p·τ_1·s² + (τ_1 + K·K_p·β)·s + K·K_p = 0`
which has two complex-conjugate stable roots near -2.16e-5 ± 4.34e-5 j. The third eigenvalue in the pre-reduced A (exact 0) is an uncontrollable integrator mode that should be discarded.

Why it matters

This blocks using MinimalRealization as a sanity pass after any closed-loop feedback composition where the controller has an integrator — which is the common case for PI/PID control. Downstream consumers (IsStable, Step, FIR export) propagate the hidden-mode artifact as apparent instability.

Notes / hypothesis

Possibly the Hankel / gramian-based reduction is treating the entire system as unobservable because the integrator pole at exactly s=0 causes a singular gramian, and the tolerance collapses all modes rather than peeling off just the hidden one. A polynomial (TF round-trip) minimal realization path would likely avoid this.

Happy to test any proposed fix on the reproducer above.

[jamestjsp]

Fixed in v0.8.0 via a6cf478.