Add actionable algebraic-loop diagnostics

algebraic_loop.go

@@ -0,0 +1,152 @@
+package controlsys
+
+import (
+	"errors"
+	"fmt"
+	"math"
+	"strings"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+// AlgebraicLoopError describes a singular or numerically singular direct-feedthrough loop.
+type AlgebraicLoopError struct {
+	// Signals contains the named input and output channels active in the
+	// singular feedthrough mode. It is empty when the operation has no names.
+	Signals []string
+	// Condition is the estimated condition number of the feedthrough equation.
+	// An exactly singular equation has an infinite condition number.
+	Condition float64
+
+	loop        *mat.Dense
+	feedthrough *mat.Dense
+	cause       error
+}
+
+func (e *AlgebraicLoopError) Error() string {
+	message := ErrAlgebraicLoop.Error()
+	if e.cause != nil {
+		message = e.cause.Error()
+	} else {
+		message = fmt.Sprintf("%s (condition number %g)", message, e.Condition)
+	}
+	if len(e.Signals) != 0 {
+		message += "; implicated signals: " + strings.Join(e.Signals, ", ")
+	}
+	return message
+}
+
+func (e *AlgebraicLoopError) Unwrap() error {
+	if e.cause != nil {
+		return e.cause
+	}
+	return ErrAlgebraicLoop
+}
+
+func newAlgebraicLoopError(loop, feedthrough *mat.Dense, condition float64) *AlgebraicLoopError {
+	return &AlgebraicLoopError{
+		Condition:   condition,
+		loop:        loop,
+		feedthrough: mat.DenseCopyOf(feedthrough),
+	}
+}
+
+func withAlgebraicLoopSignals(err error, inputNames, outputNames []string) error {
+	var diagnostic *AlgebraicLoopError
+	if !errors.As(err, &diagnostic) {
+		return err
+	}
+	signals := diagnostic.signalNames(inputNames, outputNames)
+	if len(signals) == 0 {
+		return err
+	}
+	enriched := *diagnostic
+	enriched.Signals = signals
+	enriched.loop = nil
+	enriched.feedthrough = nil
+	enriched.cause = err
+	return &enriched
+}
+
+func (e *AlgebraicLoopError) signalNames(inputNames, outputNames []string) []string {
+	if e.loop == nil || e.feedthrough == nil {
+		return nil
+	}
+	n, c := e.loop.Dims()
+	p, m := e.feedthrough.Dims()
+	if n == 0 || c != n || m != n || len(inputNames) != n || len(outputNames) != p {
+		return nil
+	}
+
+	var svd mat.SVD
+	if !svd.Factorize(e.loop, mat.SVDFull) {
+		return nil
+	}
+	values := svd.Values(nil)
+	if len(values) == 0 {
+		return nil
+	}
+	var rightVectors mat.Dense
+	svd.VTo(&rightVectors)
+
+	cutoff := values[0] * eps() * float64(n)
+	firstMode := len(values)
+	for firstMode > 0 && values[firstMode-1] <= cutoff {
+		firstMode--
+	}
+	if firstMode == len(values) {
+		firstMode--
+	}
+
+	activeInputs := make([]bool, n)
+	activeOutputs := make([]bool, p)
+	mode := make([]float64, n)
+	response := make([]float64, p)
+	for k := firstMode; k < len(values); k++ {
+		for i := range n {
+			mode[i] = rightVectors.At(i, k)
+		}
+		markActive(activeInputs, mode)
+
+		raw := e.feedthrough.RawMatrix()
+		for i := range p {
+			sum := 0.0
+			for j := range n {
+				sum += raw.Data[i*raw.Stride+j] * mode[j]
+			}
+			response[i] = sum
+		}
+		markActive(activeOutputs, response)
+	}
+
+	signals := make([]string, 0, n+p)
+	seen := make(map[string]struct{}, n+p)
+	appendActive := func(names []string, active []bool) {
+		for i, name := range names {
+			if !active[i] || name == "" {
+				continue
+			}
+			if _, ok := seen[name]; ok {
+				continue
+			}
+			seen[name] = struct{}{}
+			signals = append(signals, name)
+		}
+	}
+	appendActive(inputNames, activeInputs)
+	appendActive(outputNames, activeOutputs)
+	return signals
+}
+
+func markActive(active []bool, values []float64) {
+	scale := 0.0
+	for _, value := range values {
+		scale = math.Max(scale, math.Abs(value))
+	}
+	tolerance := scale * eps() * float64(len(values))
+	for i, value := range values {
+		if math.Abs(value) > tolerance {
+			active[i] = true
+		}
+	}
+}

connect_test.go

@@ -1121,6 +1121,9 @@ func TestFeedbackApprox_WellPosednessCheck(t *testing.T) {
 	if err == nil {
 		t.Fatal("expected singular algebraic loop error for odd-order Pade with unit gains and negative feedback")
 	}
+	if !errors.Is(err, ErrAlgebraicLoop) {
+		t.Fatalf("err = %v, want ErrAlgebraicLoop", err)
+	}
 
 	// Even-order Pade has D=1: I - (-1)*1*1 = 2 → not singular.
 	cl, err := Feedback(plant, controller, -1, WithPadeOrder(2))
@@ -2309,6 +2312,16 @@ func TestConnect_AlgebraicLoop(t *testing.T) {
 	if !errors.Is(err, ErrAlgebraicLoop) {
 		t.Errorf("got %v, want ErrAlgebraicLoop", err)
 	}
+	var diagnostic *AlgebraicLoopError
+	if !errors.As(err, &diagnostic) {
+		t.Fatalf("err = %T, want *AlgebraicLoopError", err)
+	}
+	if len(diagnostic.Signals) != 0 {
+		t.Errorf("signals = %v, want none for indexed connection", diagnostic.Signals)
+	}
+	if !math.IsInf(diagnostic.Condition, 1) {
+		t.Errorf("condition = %g, want +Inf", diagnostic.Condition)
+	}
 }
 
 func TestConnect_ZeroQ(t *testing.T) {

delay.go

@@ -2151,8 +2151,12 @@ func (sys *System) ZeroDelayApprox() (*System, error) {
 
 	var lu mat.LU
 	lu.Factorize(ImD22)
-	if luNearSingular(&lu) {
-		return nil, ErrAlgebraicLoop
+	condition := lu.Cond()
+	if nearSingularCondition(condition) {
+		return nil, fmt.Errorf(
+			"zero delay approximation: %w",
+			newAlgebraicLoopError(ImD22, sys.LFT.D22, condition),
+		)
 	}
 
 	eye := mat.NewDense(N, N, nil)
@@ -2162,7 +2166,10 @@ func (sys *System) ZeroDelayApprox() (*System, error) {
 	}
 	E := mat.NewDense(N, N, nil)
 	if err := lu.SolveTo(E, false, eye); err != nil {
-		return nil, ErrAlgebraicLoop
+		return nil, fmt.Errorf(
+			"zero delay approximation: %w",
+			newAlgebraicLoopError(ImD22, sys.LFT.D22, condition),
+		)
 	}
 
 	EC2 := mat.NewDense(N, n, nil)

delay_test.go

@@ -3298,6 +3298,13 @@ func TestZeroDelayApproxSingular(t *testing.T) {
 	if !errors.Is(err, ErrAlgebraicLoop) {
 		t.Errorf("expected ErrAlgebraicLoop, got %v", err)
 	}
+	var diagnostic *AlgebraicLoopError
+	if !errors.As(err, &diagnostic) {
+		t.Fatalf("err = %T, want *AlgebraicLoopError", err)
+	}
+	if !math.IsInf(diagnostic.Condition, 1) {
+		t.Errorf("condition = %g, want +Inf", diagnostic.Condition)
+	}
 }
 
 func TestZeroDelayApproxNoInternal(t *testing.T) {

feedback_delay.go

@@ -125,7 +125,7 @@ func (s feedbackDelayStrategy) requireWellPosedPade(plant, controller *System, s
 		return nil
 	}
 	if _, err := solveFeedbackFeedthrough(plant.D, controller.D, sign, pPlant, "Feedback", ErrAlgebraicLoop); err != nil {
-		return fmt.Errorf("Feedback: Pade approximation creates singular algebraic loop; try a different padeOrder (even vs odd) to flip feedthrough sign")
+		return fmt.Errorf("Feedback: Pade approximation creates singular algebraic loop; try a different padeOrder (even vs odd) to flip feedthrough sign: %w", err)
 	}
 	return nil
 }

matutil.go

@@ -79,7 +79,10 @@ func eps() float64 {
 }
 
 func luNearSingular(lu *mat.LU) bool {
-	cond := lu.Cond()
+	return nearSingularCondition(lu.Cond())
+}
+
+func nearSingularCondition(cond float64) bool {
 	return math.IsNaN(cond) || math.IsInf(cond, 1) || cond*eps() >= 1
 }
 

names.go

@@ -372,5 +372,9 @@ func ConnectByName(systems []*System, connections []Connection, inputs, outputs
 		Q.Set(toIdx, fromIdx, gain)
 	}
 
-	return Connect(aug, Q, inIdx, outIdx)
+	result, err := Connect(aug, Q, inIdx, outIdx)
+	if err != nil {
+		return nil, withAlgebraicLoopSignals(err, aug.InputName, aug.OutputName)
+	}
+	return result, nil
 }

names_test.go

@@ -2,6 +2,7 @@ package controlsys
 
 import (
 	"errors"
+	"math"
 	"math/cmplx"
 	"reflect"
 	"strings"
@@ -497,6 +498,114 @@ func TestConnectByName_Feedback(t *testing.T) {
 	}
 }
 
+func TestConnectByName_AlgebraicLoopDiagnostic(t *testing.T) {
+	first, _ := NewGain(mat.NewDense(1, 1, []float64{1}), 0)
+	first.InputName = []string{"first.input"}
+	first.OutputName = []string{"first.output"}
+	second, _ := NewGain(mat.NewDense(1, 1, []float64{1}), 0)
+	second.InputName = []string{"second.input"}
+	second.OutputName = []string{"second.output"}
+
+	_, err := ConnectByName(
+		[]*System{first, second},
+		[]Connection{
+			{From: "first.output", To: "second.input"},
+			{From: "second.output", To: "first.input"},
+		},
+		[]string{"first.input"},
+		[]string{"first.output"},
+	)
+	if !errors.Is(err, ErrAlgebraicLoop) {
+		t.Fatalf("err = %v, want ErrAlgebraicLoop", err)
+	}
+
+	var diagnostic *AlgebraicLoopError
+	if !errors.As(err, &diagnostic) {
+		t.Fatalf("err = %T, want *AlgebraicLoopError", err)
+	}
+	wantSignals := []string{
+		"first.input",
+		"second.input",
+		"first.output",
+		"second.output",
+	}
+	if !reflect.DeepEqual(diagnostic.Signals, wantSignals) {
+		t.Errorf("signals = %v, want %v", diagnostic.Signals, wantSignals)
+	}
+	if !math.IsInf(diagnostic.Condition, 1) {
+		t.Errorf("condition = %g, want +Inf", diagnostic.Condition)
+	}
+}
+
+func TestConnectByName_NearSingularDiagnostic(t *testing.T) {
+	model, _ := NewGain(mat.NewDense(3, 3, []float64{
+		0, 0, 0,
+		0, 0, -1,
+		0, -1, -2 * eps(),
+	}), 0)
+	model.InputName = []string{"stable.input", "critical1.input", "critical2.input"}
+	model.OutputName = []string{"stable.output", "critical1.output", "critical2.output"}
+
+	_, err := ConnectByName(
+		[]*System{model},
+		[]Connection{
+			{From: "stable.output", To: "stable.input"},
+			{From: "critical1.output", To: "critical1.input"},
+			{From: "critical2.output", To: "critical2.input"},
+		},
+		[]string{"stable.input"},
+		[]string{"stable.output"},
+	)
+	if !errors.Is(err, ErrAlgebraicLoop) {
+		t.Fatalf("err = %v, want ErrAlgebraicLoop", err)
+	}
+
+	var diagnostic *AlgebraicLoopError
+	if !errors.As(err, &diagnostic) {
+		t.Fatalf("err = %T, want *AlgebraicLoopError", err)
+	}
+	wantSignals := []string{
+		"critical1.input",
+		"critical2.input",
+		"critical1.output",
+		"critical2.output",
+	}
+	if !reflect.DeepEqual(diagnostic.Signals, wantSignals) {
+		t.Errorf("signals = %v, want %v", diagnostic.Signals, wantSignals)
+	}
+	if math.IsInf(diagnostic.Condition, 0) || diagnostic.Condition*eps() < 1 {
+		t.Errorf("condition = %g, want finite rejected condition", diagnostic.Condition)
+	}
+}
+
+func TestConnectByName_DynamicFeedbackIsWellPosed(t *testing.T) {
+	dynamic, _ := New(
+		mat.NewDense(1, 1, []float64{-1}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	dynamic.InputName = []string{"dynamic.input"}
+	dynamic.OutputName = []string{"dynamic.output"}
+	static, _ := NewGain(mat.NewDense(1, 1, []float64{1}), 0)
+	static.InputName = []string{"static.input"}
+	static.OutputName = []string{"static.output"}
+
+	_, err := ConnectByName(
+		[]*System{dynamic, static},
+		[]Connection{
+			{From: "dynamic.output", To: "static.input"},
+			{From: "static.output", To: "dynamic.input"},
+		},
+		[]string{"dynamic.input"},
+		[]string{"dynamic.output"},
+	)
+	if err != nil {
+		t.Fatalf("dynamic feedback: %v", err)
+	}
+}
+
 func TestConnectByName_WithSumBlk(t *testing.T) {
 	P := makeSISO(-2, 1, 3, 0)
 	P.InputName = []string{"u"}