From a7e05e83e2aa41448e01334f4a46843384b80499 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Fri, 31 Jul 2026 14:38:33 -0500 Subject: [PATCH] Add algebraic loop diagnostics --- algebraic_loop.go | 152 ++++++++++++++++++++++++++++++++++++++++++++++ connect_test.go | 13 ++++ delay.go | 13 +++- delay_test.go | 7 +++ feedback_delay.go | 2 +- matutil.go | 5 +- names.go | 6 +- names_test.go | 109 +++++++++++++++++++++++++++++++++ wellposed.go | 20 +++++- 9 files changed, 318 insertions(+), 9 deletions(-) create mode 100644 algebraic_loop.go diff --git a/algebraic_loop.go b/algebraic_loop.go new file mode 100644 index 0000000..df551a2 --- /dev/null +++ b/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 + } + } +} diff --git a/connect_test.go b/connect_test.go index 6d98019..a42760d 100644 --- a/connect_test.go +++ b/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) { diff --git a/delay.go b/delay.go index 5b52261..ee7ce54 100644 --- a/delay.go +++ b/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) diff --git a/delay_test.go b/delay_test.go index 89b6616..f23d3fe 100644 --- a/delay_test.go +++ b/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) { diff --git a/feedback_delay.go b/feedback_delay.go index 6197602..7d10f03 100644 --- a/feedback_delay.go +++ b/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 } diff --git a/matutil.go b/matutil.go index 7273b91..99bb306 100644 --- a/matutil.go +++ b/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 } diff --git a/names.go b/names.go index d49ea4a..70d6e06 100644 --- a/names.go +++ b/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 } diff --git a/names_test.go b/names_test.go index 88381d5..fa93f5f 100644 --- a/names_test.go +++ b/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"} diff --git a/wellposed.go b/wellposed.go index e704e0d..fb0785d 100644 --- a/wellposed.go +++ b/wellposed.go @@ -1,6 +1,7 @@ package controlsys import ( + "errors" "fmt" "gonum.org/v1/gonum/mat" @@ -21,17 +22,30 @@ func solveIdentityMinusScaledProduct(left, right *mat.Dense, scale float64, size var lu mat.LU lu.Factorize(loop) - if luNearSingular(&lu) { - return nil, fmt.Errorf("%s: direct feedthrough loop is singular: %w", context, singular) + condition := lu.Cond() + if nearSingularCondition(condition) { + return nil, directFeedthroughSolveError( + context, "is singular", singular, loop, right, condition, + ) } result := mat.NewDense(size, size, nil) if err := lu.SolveTo(result, false, eye); err != nil { - return nil, fmt.Errorf("%s: direct feedthrough loop solve failed: %w", context, singular) + return nil, directFeedthroughSolveError( + context, "solve failed", singular, loop, right, condition, + ) } return result, nil } +func directFeedthroughSolveError(context, failure string, singular error, loop, feedthrough *mat.Dense, condition float64) error { + cause := singular + if errors.Is(singular, ErrAlgebraicLoop) { + cause = newAlgebraicLoopError(loop, feedthrough, condition) + } + return fmt.Errorf("%s: direct feedthrough loop %s: %w", context, failure, cause) +} + func solveFeedbackFeedthrough(plantD, controllerD *mat.Dense, sign float64, size int, context string, singular error) (*mat.Dense, error) { return solveIdentityMinusScaledProduct(plantD, controllerD, sign, size, context, singular) }