From ed73d670c345ec12d9e845cfaceb2128c111cdef Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:11:08 -0500 Subject: [PATCH 01/10] Correct model array stacking semantics --- README.md | 6 +++-- model_array.go | 54 ++++++++++++++++++++++++++++++++++++++++++--- model_array_test.go | 50 ++++++++++++++++++++++++++++++++++++++--- 3 files changed, 102 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 9095d22..56d5d68 100644 --- a/README.md +++ b/README.md @@ -90,7 +90,8 @@ func main() { | `NewZPKMIMO` | MIMO zero-pole-gain model | | `NewFRD` | Frequency-response data model from sampled complex responses | | `NewModelArray` | Compatible array of state-space models for sweeps or model grids | -| `StackModelArrays` | Concatenate compatible model arrays along a new leading axis | +| `StackModelArrays` | Stack equal-shaped compatible arrays along a new leading axis | +| `ConcatModelArrays` | Flatten and concatenate compatible model arrays | | `NewGeneralizedModel` | Wrap a fixed or tunable block and attach analysis points | | `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with an analysis point | | `NewPhysicalComponent` | Wrap a model with named physical ports | @@ -334,7 +335,8 @@ func main() { | Function/Method | Description | |-----------------|-------------| | `NewModelArray` | Create a shaped array of compatible state-space models | -| `StackModelArrays` | Stack compatible model arrays | +| `StackModelArrays` | Stack equal-shaped compatible arrays along a new leading axis | +| `ConcatModelArrays` | Flatten and concatenate compatible model arrays | | `(*ModelArray).Model` / `ModelFlat` | Retrieve a model by multidimensional or flat index | | `(*ModelArray).SelectFlat` | Select a flat-index subset of models | | `(*ModelArray).FreqResponse` | Frequency response for every model in the array | diff --git a/model_array.go b/model_array.go index 0a9615f..e4d010b 100644 --- a/model_array.go +++ b/model_array.go @@ -80,16 +80,52 @@ func StackModelArrays(arrays ...*ModelArray) (*ModelArray, error) { return nil, fmt.Errorf("StackModelArrays: no arrays: %w", ErrDimensionMismatch) } var ref *ModelArray - total := 0 for _, arr := range arrays { if arr == nil { return nil, fmt.Errorf("StackModelArrays: nil array: %w", ErrDimensionMismatch) } if ref == nil { ref = arr - } else if err := validateModelArrayHeadersCompatible(ref, arr); err != nil { + continue + } + if !sameModelArrayShape(ref.shape, arr.shape) { + return nil, fmt.Errorf("StackModelArrays: shape %v != %v: %w", arr.shape, ref.shape, ErrDimensionMismatch) + } + if err := validateModelArrayHeadersCompatible(ref, arr); err != nil { return nil, fmt.Errorf("StackModelArrays: %w", err) } + } + shape := make([]int, len(ref.shape)+1) + shape[0] = len(arrays) + copy(shape[1:], ref.shape) + models := copyModelArrays(arrays) + return NewModelArray(shape, models) +} + +func ConcatModelArrays(arrays ...*ModelArray) (*ModelArray, error) { + if len(arrays) == 0 { + return nil, fmt.Errorf("ConcatModelArrays: no arrays: %w", ErrDimensionMismatch) + } + var ref *ModelArray + total := 0 + for _, arr := range arrays { + if arr == nil { + return nil, fmt.Errorf("ConcatModelArrays: nil array: %w", ErrDimensionMismatch) + } + if ref == nil { + ref = arr + } else if err := validateModelArrayHeadersCompatible(ref, arr); err != nil { + return nil, fmt.Errorf("ConcatModelArrays: %w", err) + } + total += arr.Len() + } + models := copyModelArrays(arrays) + return NewModelArray([]int{total}, models) +} + +func copyModelArrays(arrays []*ModelArray) []*System { + total := 0 + for _, arr := range arrays { total += arr.Len() } models := make([]*System, 0, total) @@ -102,7 +138,19 @@ func StackModelArrays(arrays ...*ModelArray) (*ModelArray, error) { } } } - return NewModelArray([]int{total}, models) + return models +} + +func sameModelArrayShape(a, b []int) bool { + if len(a) != len(b) { + return false + } + for i := range a { + if a[i] != b[i] { + return false + } + } + return true } func (a *ModelArray) Shape() []int { diff --git a/model_array_test.go b/model_array_test.go index 4a9b725..78ab7bc 100644 --- a/model_array_test.go +++ b/model_array_test.go @@ -65,7 +65,7 @@ func TestModelArraySelectAndStack(t *testing.T) { if err != nil { t.Fatalf("left array: %v", err) } - right, err := NewModelArray([]int{1}, []*System{b}) + right, err := NewModelArray([]int{2}, []*System{b, c}) if err != nil { t.Fatalf("right array: %v", err) } @@ -84,17 +84,61 @@ func TestModelArraySelectAndStack(t *testing.T) { if err != nil { t.Fatalf("StackModelArrays: %v", err) } - if got, want := stacked.Shape(), []int{3}; !sameInts(got, want) { + if got, want := stacked.Shape(), []int{2, 2}; !sameInts(got, want) { t.Fatalf("stacked Shape() = %v, want %v", got, want) } + if got, ok, err := stacked.Model(0, 1); err != nil || ok || got != nil { + t.Fatalf("stacked void = (%v, %v, %v), want nil,false,nil", got, ok, err) + } + if got, ok, err := stacked.Model(1, 0); err != nil || !ok || got.A.At(0, 0) != b.A.At(0, 0) { + t.Fatalf("stacked Model(1,0) = (%v, %v, %v), want copied b", got, ok, err) + } incompatible, err := NewModelArray([]int{1}, []*System{c}) if err != nil { t.Fatal(err) } + if _, err := StackModelArrays(left, incompatible); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("StackModelArrays shape mismatch err = %v, want ErrDimensionMismatch", err) + } + + incompatible, err = NewModelArray([]int{2}, []*System{b, c}) + if err != nil { + t.Fatal(err) + } incompatible.models[0].Dt = 0.2 if _, err := StackModelArrays(left, incompatible); !errors.Is(err, ErrDimensionMismatch) { - t.Fatalf("StackModelArrays incompatible err = %v, want ErrDimensionMismatch", err) + t.Fatalf("StackModelArrays header mismatch err = %v, want ErrDimensionMismatch", err) + } + if _, err := StackModelArrays(); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("StackModelArrays empty err = %v, want ErrDimensionMismatch", err) + } + if _, err := StackModelArrays(left, nil); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("StackModelArrays nil err = %v, want ErrDimensionMismatch", err) + } +} + +func TestConcatModelArraysFlattensCompatibleArrays(t *testing.T) { + a := modelArrayTestSystem(t, 0, []float64{-1, 2, -3, -4}) + b := modelArrayTestSystem(t, 0, []float64{-2, 1, -4, -5}) + left, err := NewModelArray([]int{1, 2}, []*System{a, nil}) + if err != nil { + t.Fatal(err) + } + right, err := NewModelArray([]int{1}, []*System{b}) + if err != nil { + t.Fatal(err) + } + + concatenated, err := ConcatModelArrays(left, right) + if err != nil { + t.Fatalf("ConcatModelArrays: %v", err) + } + if got, want := concatenated.Shape(), []int{3}; !sameInts(got, want) { + t.Fatalf("concatenated Shape() = %v, want %v", got, want) + } + if _, ok, err := concatenated.ModelFlat(2); err != nil || !ok { + t.Fatalf("concatenated final model ok=%v err=%v, want ok", ok, err) } } From b105aa0cdeecd6c3d5c76b6df28ebe3ae6435f32 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:15:28 -0500 Subject: [PATCH 02/10] Clarify sampled passivity guarantees --- README.md | 3 +- passivity.go | 133 +++++++++++++++++++++++++++++++++++++++++----- passivity_test.go | 61 +++++++++++++++++++++ 3 files changed, 183 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 56d5d68..7d759a6 100644 --- a/README.md +++ b/README.md @@ -178,7 +178,8 @@ func main() { | `Bandwidth` | -3 dB bandwidth | | `RootLocus` | Root locus as a function of loop gain | | `Pzmap` | Poles and transmission zeros for plotting/inspection | -| `Passive` / `FRDPassive` | Passivity check from a state-space model or FRD samples | +| `SampledPassive` / `FRDPassive` | Passivity evidence on an explicit frequency grid; passing samples are not an analytic certificate | +| `Passive` | Compatibility alias for `SampledPassive` | | `SpectralFactor` | Spectral factor for supported static-gain models | ### Control Design diff --git a/passivity.go b/passivity.go index 90a3ae0..4e11c0b 100644 --- a/passivity.go +++ b/passivity.go @@ -13,30 +13,55 @@ type PassivityOptions struct { Tol float64 } +// PassivityStatus distinguishes a sampled pass from a conclusive violation or +// an analytic certificate. +type PassivityStatus string + +const ( + PassivityViolated PassivityStatus = "violated" + PassivitySampled PassivityStatus = "sampled-pass" + PassivityCertified PassivityStatus = "certified" +) + type PassivityResult struct { + Status PassivityStatus Passive bool MinHermitianPart float64 Frequency float64 + Omega []float64 + Tolerance float64 } +// Passive reports passivity evidence over a finite frequency grid. +// Deprecated: use SampledPassive to make the sampled guarantee explicit. func Passive(sys *System, opts *PassivityOptions) (*PassivityResult, error) { + return SampledPassive(sys, opts) +} + +// SampledPassive checks the positive-real frequency-domain condition over a +// finite grid. PassivitySampled is evidence, not an analytic certificate; +// PassivityViolated identifies a conclusive sampled counterexample. +func SampledPassive(sys *System, opts *PassivityOptions) (*PassivityResult, error) { if sys == nil { - return nil, fmt.Errorf("Passive: nil system: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("SampledPassive: nil system: %w", ErrDimensionMismatch) } - if err := newDescriptorPolicy(sys).requireStandard("Passive"); err != nil { + if err := newDescriptorPolicy(sys).requireStandard("SampledPassive"); err != nil { return nil, err } if sys.HasDelay() { - return nil, fmt.Errorf("Passive: delayed systems are not supported: %w", ErrDescriptorUnsupported) + return nil, fmt.Errorf("SampledPassive: delayed systems are not supported: %w", ErrDescriptorUnsupported) } stable, err := sys.IsStable() if err != nil { return nil, err } if !stable { - return nil, fmt.Errorf("Passive: %w", ErrUnstable) + return nil, fmt.Errorf("SampledPassive: %w", ErrUnstable) + } + omega, tol, err := passivityGrid(opts, sys.Dt) + if err != nil { + return nil, fmt.Errorf("SampledPassive: %w", err) } - omega, tol := passivityGrid(opts) frd, err := sys.FRD(omega) if err != nil { return nil, err @@ -48,13 +73,44 @@ func FRDPassive(frd *FRD, opts *PassivityOptions) (*PassivityResult, error) { if frd == nil || len(frd.Response) == 0 { return nil, fmt.Errorf("FRDPassive: insufficient data: %w", ErrInsufficientData) } + if len(frd.Omega) != len(frd.Response) { + return nil, fmt.Errorf("FRDPassive: %d responses for %d frequencies: %w", len(frd.Response), len(frd.Omega), ErrDimensionMismatch) + } + if err := validatePassivityGrid(frd.Omega, frd.Dt); err != nil { + return nil, fmt.Errorf("FRDPassive: %w", err) + } + if len(frd.Response[0]) == 0 || len(frd.Response[0][0]) == 0 { + return nil, fmt.Errorf("FRDPassive: empty response matrix: %w", ErrDimensionMismatch) + } p, m := len(frd.Response[0]), len(frd.Response[0][0]) if p != m { return nil, fmt.Errorf("FRDPassive: model must be square, got %dx%d: %w", p, m, ErrDimensionMismatch) } - _, tol := passivityGrid(opts) - result := &PassivityResult{Passive: true, MinHermitianPart: math.Inf(1)} + tol, err := passivityTolerance(opts) + if err != nil { + return nil, fmt.Errorf("FRDPassive: %w", err) + } + result := &PassivityResult{ + Status: PassivitySampled, + Passive: true, + MinHermitianPart: math.Inf(1), + Omega: copyFloatSlice(frd.Omega), + Tolerance: tol, + } for k, h := range frd.Response { + if len(h) != p { + return nil, fmt.Errorf("FRDPassive: response[%d] has %d rows, want %d: %w", k, len(h), p, ErrDimensionMismatch) + } + for i, row := range h { + if len(row) != m { + return nil, fmt.Errorf("FRDPassive: response[%d][%d] has %d columns, want %d: %w", k, i, len(row), m, ErrDimensionMismatch) + } + for j, value := range row { + if !finiteComplex(value) { + return nil, fmt.Errorf("FRDPassive: response[%d][%d][%d] is not finite: %w", k, i, j, ErrInsufficientData) + } + } + } minPart := minHermitianPart(h) if minPart < result.MinHermitianPart { result.MinHermitianPart = minPart @@ -64,6 +120,9 @@ func FRDPassive(frd *FRD, opts *PassivityOptions) (*PassivityResult, error) { } } result.Passive = result.MinHermitianPart >= -tol + if !result.Passive { + result.Status = PassivityViolated + } return result, nil } @@ -102,15 +161,63 @@ func SpectralFactor(sys *System) (*System, error) { return NewGain(D, sys.Dt) } -func passivityGrid(opts *PassivityOptions) ([]float64, float64) { - tol := 1e-9 - if opts != nil && opts.Tol > 0 { - tol = opts.Tol +func passivityGrid(opts *PassivityOptions, dt float64) ([]float64, float64, error) { + tol, err := passivityTolerance(opts) + if err != nil { + return nil, 0, err } if opts != nil && len(opts.Omega) > 0 { - return opts.Omega, tol + omega := copyFloatSlice(opts.Omega) + if err := validatePassivityGrid(omega, dt); err != nil { + return nil, 0, err + } + return omega, tol, nil + } + upper := 1e2 + if dt > 0 { + upper = math.Pi / dt + } + lower := math.Min(1e-2, upper*1e-4) + omega := make([]float64, 121) + copy(omega[1:], logspace(math.Log10(lower), math.Log10(upper), 120)) + return omega, tol, nil +} + +func passivityTolerance(opts *PassivityOptions) (float64, error) { + tol := 1e-9 + if opts == nil || opts.Tol == 0 { + return tol, nil } - return logspace(-2, 2, 120), tol + if math.IsNaN(opts.Tol) || math.IsInf(opts.Tol, 0) || opts.Tol < 0 { + return 0, fmt.Errorf("invalid tolerance %g: %w", opts.Tol, ErrDimensionMismatch) + } + return opts.Tol, nil +} + +func validatePassivityGrid(omega []float64, dt float64) error { + if len(omega) == 0 { + return fmt.Errorf("frequency grid is empty: %w", ErrInsufficientData) + } + upper := math.Inf(1) + if dt > 0 { + upper = math.Pi / dt + } + for i, frequency := range omega { + if math.IsNaN(frequency) || math.IsInf(frequency, 0) || frequency < 0 { + return fmt.Errorf("omega[%d]=%g is invalid: %w", i, frequency, ErrDimensionMismatch) + } + if frequency > upper*(1+1e-12) { + return fmt.Errorf("omega[%d]=%g exceeds Nyquist frequency %g: %w", i, frequency, upper, ErrDimensionMismatch) + } + if i > 0 && frequency < omega[i-1] { + return fmt.Errorf("frequency grid is not sorted at index %d: %w", i, ErrDimensionMismatch) + } + } + return nil +} + +func finiteComplex(value complex128) bool { + return !math.IsNaN(real(value)) && !math.IsNaN(imag(value)) && !math.IsInf(real(value), 0) && !math.IsInf(imag(value), 0) } func minHermitianPart(h [][]complex128) float64 { diff --git a/passivity_test.go b/passivity_test.go index 3388988..e98b8c9 100644 --- a/passivity_test.go +++ b/passivity_test.go @@ -16,6 +16,9 @@ func TestPassivityDenseAndFRD(t *testing.T) { if !result.Passive || result.MinHermitianPart < 0 { t.Fatalf("expected passive result, got %#v", result) } + if result.Status != PassivitySampled || len(result.Omega) != 121 { + t.Fatalf("expected disclosed sampled result, got %#v", result) + } nonpassive := makeSISO(-1, 1, -1, 0) bad, err := Passive(nonpassive, nil) @@ -25,6 +28,9 @@ func TestPassivityDenseAndFRD(t *testing.T) { if bad.Passive { t.Fatalf("expected nonpassive result, got %#v", bad) } + if bad.Status != PassivityViolated { + t.Fatalf("nonpassive status = %q, want %q", bad.Status, PassivityViolated) + } frd, err := passive.FRD([]float64{0.1, 1, 10}) if err != nil { @@ -37,6 +43,49 @@ func TestPassivityDenseAndFRD(t *testing.T) { if !frdResult.Passive { t.Fatalf("expected passive FRD result, got %#v", frdResult) } + if got, want := frdResult.Omega, []float64{0.1, 1, 10}; !sameFloatSlice(got, want) { + t.Fatalf("checked frequencies = %v, want %v", got, want) + } +} + +func TestSampledPassiveDoesNotClaimCertification(t *testing.T) { + sys := makeSISO(-1, 1, -2, 1) + result, err := SampledPassive(sys, &PassivityOptions{Omega: []float64{10, 100}}) + if err != nil { + t.Fatalf("SampledPassive high-frequency grid: %v", err) + } + if !result.Passive || result.Status != PassivitySampled { + t.Fatalf("high-frequency sampled result = %#v, want sampled pass", result) + } + + result, err = SampledPassive(sys, &PassivityOptions{Omega: []float64{0, 10}}) + if err != nil { + t.Fatalf("SampledPassive grid including DC: %v", err) + } + if result.Passive || result.Status != PassivityViolated || result.Frequency != 0 { + t.Fatalf("grid including DC result = %#v, want violation at DC", result) + } +} + +func TestSampledPassiveValidatesGridAndUsesDiscreteNyquist(t *testing.T) { + discrete := makeSISO(0.5, 1, 1, 0) + discrete.Dt = 0.2 + result, err := SampledPassive(discrete, nil) + if err != nil { + t.Fatalf("SampledPassive discrete default: %v", err) + } + if got, want := result.Omega[len(result.Omega)-1], 3.141592653589793/0.2; got != want { + t.Fatalf("discrete grid upper frequency = %g, want Nyquist %g", got, want) + } + if _, err := SampledPassive(discrete, &PassivityOptions{Omega: []float64{0, 20}}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("above-Nyquist grid err = %v, want ErrDimensionMismatch", err) + } + if _, err := SampledPassive(discrete, &PassivityOptions{Omega: []float64{1, 0}}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("unsorted grid err = %v, want ErrDimensionMismatch", err) + } + if _, err := SampledPassive(discrete, &PassivityOptions{Tol: -1}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("negative tolerance err = %v, want ErrDimensionMismatch", err) + } } func TestFRDPassiveUsesMIMOHermitianEigenvalue(t *testing.T) { @@ -84,3 +133,15 @@ func TestSpectralFactorStaticGainAndUnsupportedCases(t *testing.T) { t.Fatalf("delayed spectral factor err = %v, want ErrDescriptorUnsupported", err) } } + +func sameFloatSlice(a, b []float64) bool { + if len(a) != len(b) { + return false + } + for i := range a { + if a[i] != b[i] { + return false + } + } + return true +} From e029f71a069f98e965f0f9778881ac0e184097ad Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:24:24 -0500 Subject: [PATCH 03/10] Implement basis-invariant modal truncation --- README.md | 2 +- modal_reduction.go | 280 ++++++++++++++++++++++++++++++++++------ modal_reduction_test.go | 141 ++++++++++++++++++-- 3 files changed, 374 insertions(+), 49 deletions(-) diff --git a/README.md b/README.md index 7d759a6..c8a2874 100644 --- a/README.md +++ b/README.md @@ -242,7 +242,7 @@ func main() { | `Sminreal` | Minimal realization via staircase reduction | | `Stabsep` | Stable/unstable decomposition | | `Modsep` | Modal decomposition around a cutoff | -| `ModalTruncate` | Modal truncation result with kept-state metadata | +| `ModalTruncate` | Basis-invariant real-Schur modal truncation with retained poles and projection bases | | `Canon` | Modal or companion canonical form | | `SS2SS` | Similarity transform with a user-supplied state basis | | `StateTransform` | Alias-style state-basis transform helper | diff --git a/modal_reduction.go b/modal_reduction.go index 8aeb1d6..95b9e7f 100644 --- a/modal_reduction.go +++ b/modal_reduction.go @@ -1,8 +1,14 @@ package controlsys import ( + "fmt" "math" - "sort" + "math/cmplx" + + "gonum.org/v1/gonum/blas" + "gonum.org/v1/gonum/blas/blas64" + "gonum.org/v1/gonum/lapack" + "gonum.org/v1/gonum/mat" ) type ModalTruncateOptions struct { @@ -11,10 +17,13 @@ type ModalTruncateOptions struct { } type ModalReductionResult struct { - Sys *System - Order int - Method string - Kept []int + Sys *System + Order int + Method string + Kept []int + KeptPoles []complex128 + Basis *mat.Dense + Projection *mat.Dense } func ModalTruncate(sys *System, opts *ModalTruncateOptions) (*ModalReductionResult, error) { @@ -28,63 +37,258 @@ func ModalTruncate(sys *System, opts *ModalTruncateOptions) (*ModalReductionResu if err := policy.requireDelayFree("ModalTruncate"); err != nil { return nil, err } - n, _, _ := sys.Dims() + n, m, p := sys.Dims() if n == 0 { - return &ModalReductionResult{Sys: sys.Copy(), Order: 0, Method: "modal-truncate"}, nil + return &ModalReductionResult{Sys: sys.Copy(), Order: 0, Method: "real-schur-modal-truncate", Basis: &mat.Dense{}, Projection: &mat.Dense{}}, nil } - order := opts.Order - if order == 0 { - order = modalAutoOrder(sys, opts.MaxRealPart) + if opts.Order < 0 || opts.Order > n { + return nil, ErrInvalidOrder + } + + t, z, err := modalSchur(sys) + if err != nil { + return nil, err } - if order < 1 || order > n { + threshold := opts.Order == 0 && opts.MaxRealPart != 0 + if err := orderModalSchur(t, z, n, sys.Dt, opts.MaxRealPart, threshold); err != nil { + return nil, err + } + poles := schurEigenvaluesRaw(t, n) + order := modalReductionOrder(t, poles, n, sys.Dt, opts) + if order < 1 || order > n || splitsSchurBlock(t, n, order) { return nil, ErrInvalidOrder } + if order < unstableModalCount(poles, sys.Dt) { + return nil, fmt.Errorf("ModalTruncate: order %d would discard an unstable mode: %w", order, ErrInvalidOrder) + } + keptPoles := append([]complex128(nil), poles[:order]...) if order == n { - return &ModalReductionResult{Sys: sys.Copy(), Order: n, Method: "modal-truncate", Kept: rangeInts(n)}, nil + identity := eyeDense(n) + return &ModalReductionResult{ + Sys: sys.Copy(), + Order: n, + Method: "real-schur-modal-truncate", + Kept: rangeInts(n), + KeptPoles: keptPoles, + Basis: identity, + Projection: eyeDense(n), + }, nil + } + + x, err := separateSchurBlocks(t, n, order) + if err != nil { + return nil, err + } + bModal, cModal := modalInputOutput(sys, z, n, m, p) + aReduced := extractModalBlock(t, n, 0, order, 0, order) + bReduced := mat.NewDense(order, m, nil) + for i := range order { + for j := range m { + value := bModal[i*m+j] + for k := order; k < n; k++ { + value -= x[i*(n-order)+k-order] * bModal[k*m+j] + } + bReduced.Set(i, j, value) + } } - elim := make([]int, 0, n-order) - for i := order; i < n; i++ { - elim = append(elim, i) + cReduced := mat.NewDense(p, order, nil) + for i := range p { + copy(cReduced.RawMatrix().Data[i*order:(i+1)*order], cModal[i*n:i*n+order]) } - reduced, err := Modred(sys, elim, Truncate) + reduced, err := policy.resultWithOriginalFeedthrough(aReduced, bReduced, cReduced) if err != nil { return nil, err } - return &ModalReductionResult{Sys: reduced, Order: order, Method: "modal-truncate", Kept: rangeInts(order)}, nil + basis, projection := modalProjectionBases(z, x, n, order) + return &ModalReductionResult{ + Sys: reduced, + Order: order, + Method: "real-schur-modal-truncate", + Kept: rangeInts(order), + KeptPoles: keptPoles, + Basis: basis, + Projection: projection, + }, nil +} + +func modalSchur(sys *System) (t, z []float64, err error) { + n, _, _ := sys.Dims() + t = make([]float64, n*n) + aRaw := sys.A.RawMatrix() + copyStrided(t, n, aRaw.Data, aRaw.Stride, n, n) + z = make([]float64, n*n) + wr := make([]float64, n) + wi := make([]float64, n) + bwork := make([]bool, n) + query := make([]float64, 1) + impl.Dgees(lapack.SchurHess, lapack.SortNone, nil, n, t, n, wr, wi, z, n, query, -1, bwork) + work := make([]float64, int(query[0])) + _, ok := impl.Dgees(lapack.SchurHess, lapack.SortNone, nil, n, t, n, wr, wi, z, n, work, len(work), bwork) + if !ok { + return nil, nil, ErrSchurFailed + } + return t, z, nil +} + +func orderModalSchur(t, z []float64, n int, dt, maxRealPart float64, threshold bool) error { + work := make([]float64, n) + placed := 0 + for placed < n { + evals := schurEigenvaluesRaw(t, n) + best := placed + for i := placed; i < n; i += schurBlockSize(t, n, i) { + if modalBlockBefore(evals[i], evals[best], dt, maxRealPart, threshold) { + best = i + } + } + if best != placed { + _, _, ok := impl.Dtrexc(lapack.UpdateSchur, n, t, n, z, n, best, placed, work) + if !ok { + return ErrSchurFailed + } + } + placed += schurBlockSize(t, n, placed) + } + return nil } -func modalAutoOrder(sys *System, maxRealPart float64) int { - poles, err := sys.Poles() - if err != nil || len(poles) == 0 { - return 0 +func modalBlockBefore(candidate, current complex128, dt, maxRealPart float64, threshold bool) bool { + if threshold { + candidateSelected := real(candidate) >= maxRealPart || modalPoleUnstable(candidate, dt) + currentSelected := real(current) >= maxRealPart || modalPoleUnstable(current, dt) + if candidateSelected != currentSelected { + return candidateSelected + } } - if maxRealPart == 0 { - maxRealPart = math.Inf(1) + candidateScore := real(candidate) + currentScore := real(current) + if dt > 0 { + candidateScore = cmplx.Abs(candidate) + currentScore = cmplx.Abs(current) } - type poleScore struct { - index int - real float64 + if candidateScore != currentScore { + return candidateScore > currentScore } - scores := make([]poleScore, len(poles)) - for i, p := range poles { - scores[i] = poleScore{index: i, real: real(p)} + return math.Abs(imag(candidate)) < math.Abs(imag(current)) +} + +func modalReductionOrder(t []float64, poles []complex128, n int, dt float64, opts *ModalTruncateOptions) int { + if opts.Order > 0 { + return opts.Order } - sort.Slice(scores, func(i, j int) bool { return scores[i].real > scores[j].real }) - order := 0 - for _, score := range scores { - if score.real >= maxRealPart { - order++ + if opts.MaxRealPart != 0 { + order := 0 + for order < n && (real(poles[order]) >= opts.MaxRealPart || modalPoleUnstable(poles[order], dt)) { + order += schurBlockSize(t, n, order) + } + if order > 0 { + return order } } + order := n / 2 if order == 0 { - order = len(poles) / 2 - if order == 0 { - order = 1 - } + order = 1 + } + if splitsSchurBlock(t, n, order) { + order++ + } + unstable := unstableModalCount(poles, dt) + if order < unstable { + order = unstable } return order } +func separateSchurBlocks(t []float64, n, order int) ([]float64, error) { + discarded := n - order + x := make([]float64, order*discarded) + for i := range order { + for j := range discarded { + x[i*discarded+j] = -t[i*n+order+j] + } + } + scale, ok := impl.Dtrsyl(blas.NoTrans, blas.NoTrans, -1, order, discarded, t, n, t[order*n+order:], n, x, discarded) + if !ok || scale == 0 || scale < 1e-12 { + return nil, fmt.Errorf("ModalTruncate: retained and discarded modes cannot be separated reliably: %w", ErrSchurFailed) + } + if scale != 1 { + for i := range x { + x[i] /= scale + } + } + return x, nil +} + +func modalInputOutput(sys *System, z []float64, n, m, p int) (bModal, cModal []float64) { + zGeneral := blas64.General{Rows: n, Cols: n, Stride: n, Data: z} + bRaw := sys.B.RawMatrix() + bModal = make([]float64, n*m) + blas64.Gemm(blas.Trans, blas.NoTrans, 1, zGeneral, + blas64.General{Rows: n, Cols: m, Stride: bRaw.Stride, Data: bRaw.Data}, + 0, blas64.General{Rows: n, Cols: m, Stride: m, Data: bModal}) + cRaw := sys.C.RawMatrix() + cModal = make([]float64, p*n) + blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, + blas64.General{Rows: p, Cols: n, Stride: cRaw.Stride, Data: cRaw.Data}, zGeneral, + 0, blas64.General{Rows: p, Cols: n, Stride: n, Data: cModal}) + return bModal, cModal +} + +func modalProjectionBases(z, x []float64, n, order int) (*mat.Dense, *mat.Dense) { + basis := mat.NewDense(n, order, nil) + for i := range n { + copy(basis.RawMatrix().Data[i*order:(i+1)*order], z[i*n:i*n+order]) + } + schurProjection := mat.NewDense(order, n, nil) + for i := range order { + schurProjection.Set(i, i, 1) + for j := order; j < n; j++ { + schurProjection.Set(i, j, -x[i*(n-order)+j-order]) + } + } + projection := mat.NewDense(order, n, nil) + projection.Mul(schurProjection, mat.NewDense(n, n, z).T()) + return basis, projection +} + +func extractModalBlock(data []float64, stride, rowStart, rowEnd, colStart, colEnd int) *mat.Dense { + rows := rowEnd - rowStart + cols := colEnd - colStart + out := mat.NewDense(rows, cols, nil) + for i := range rows { + copy(out.RawMatrix().Data[i*cols:(i+1)*cols], data[(rowStart+i)*stride+colStart:(rowStart+i)*stride+colEnd]) + } + return out +} + +func schurBlockSize(t []float64, n, start int) int { + if start+1 < n && t[(start+1)*n+start] != 0 { + return 2 + } + return 1 +} + +func splitsSchurBlock(t []float64, n, order int) bool { + return order > 0 && order < n && t[order*n+order-1] != 0 +} + +func unstableModalCount(poles []complex128, dt float64) int { + count := 0 + for _, pole := range poles { + if modalPoleUnstable(pole, dt) { + count++ + } + } + return count +} + +func modalPoleUnstable(pole complex128, dt float64) bool { + if dt > 0 { + return cmplx.Abs(pole) >= 1 + } + return real(pole) >= 0 +} + func rangeInts(n int) []int { out := make([]int, n) for i := range out { diff --git a/modal_reduction_test.go b/modal_reduction_test.go index 247885b..3c88162 100644 --- a/modal_reduction_test.go +++ b/modal_reduction_test.go @@ -1,9 +1,16 @@ package controlsys -import "testing" +import ( + "errors" + "math" + "math/cmplx" + "testing" -func TestModalTruncateReducesOrderAndPreservesMetadata(t *testing.T) { - sys := benchSysNonSym(4, 1, 1) + "gonum.org/v1/gonum/mat" +) + +func TestModalTruncateRetainsDominantModesAndMetadata(t *testing.T) { + sys := modalTestSystem(t) sys.InputName = []string{"u"} sys.OutputName = []string{"y"} sys.StateName = []string{"x1", "x2", "x3", "x4"} @@ -12,7 +19,7 @@ func TestModalTruncateReducesOrderAndPreservesMetadata(t *testing.T) { if err != nil { t.Fatalf("ModalTruncate: %v", err) } - if result.Method != "modal-truncate" || result.Order != 2 { + if result.Method != "real-schur-modal-truncate" || result.Order != 2 { t.Fatalf("metadata = %#v", result) } if n, _, _ := result.Sys.Dims(); n != 2 { @@ -21,18 +28,132 @@ func TestModalTruncateReducesOrderAndPreservesMetadata(t *testing.T) { if !sameStrings(result.Sys.InputName, []string{"u"}) || !sameStrings(result.Sys.OutputName, []string{"y"}) { t.Fatalf("names = %v/%v", result.Sys.InputName, result.Sys.OutputName) } + if !sameComplexApprox(result.KeptPoles, []complex128{-0.1, -1}, 1e-12) { + t.Fatalf("kept poles = %v, want [-0.1 -1]", result.KeptPoles) + } + var identity mat.Dense + identity.Mul(result.Projection, result.Basis) + if !mat.EqualApprox(&identity, mat.NewDiagDense(2, []float64{1, 1}), 1e-11) { + t.Fatalf("projection*basis =\n%v", mat.Formatted(&identity)) + } +} + +func TestModalTruncateIsInvariantUnderSimilarityTransform(t *testing.T) { + base := modalTestSystem(t) + transform := mat.NewDense(4, 4, []float64{ + 1, 0.2, -0.1, 0.3, + 0.1, 1.2, 0.4, -0.2, + -0.2, 0.1, 0.9, 0.2, + 0.3, -0.1, 0.2, 1.1, + }) + equivalent, err := SS2SS(base, transform) + if err != nil { + t.Fatalf("SS2SS: %v", err) + } + + left, err := ModalTruncate(base, &ModalTruncateOptions{Order: 2}) + if err != nil { + t.Fatalf("ModalTruncate base: %v", err) + } + right, err := ModalTruncate(equivalent, &ModalTruncateOptions{Order: 2}) + if err != nil { + t.Fatalf("ModalTruncate equivalent: %v", err) + } + if !sameComplexApprox(left.KeptPoles, right.KeptPoles, 1e-10) { + t.Fatalf("kept poles differ: %v vs %v", left.KeptPoles, right.KeptPoles) + } + omega := []float64{0, 0.1, 1, 10} + leftResponse, err := left.Sys.FreqResponse(omega) + if err != nil { + t.Fatal(err) + } + rightResponse, err := right.Sys.FreqResponse(omega) + if err != nil { + t.Fatal(err) + } + for k := range omega { + if delta := cmplx.Abs(leftResponse.At(k, 0, 0) - rightResponse.At(k, 0, 0)); delta > 1e-9 { + t.Fatalf("response difference at omega=%g is %g", omega[k], delta) + } + } +} + +func TestModalTruncatePreservesComplexPairs(t *testing.T) { + sys, err := New( + mat.NewDense(3, 3, []float64{ + -0.1, -2, 0.3, + 2, -0.1, -0.2, + 0, 0, -5, + }), + mat.NewDense(3, 1, []float64{1, 2, 3}), + mat.NewDense(1, 3, []float64{2, -1, 0.5}), + mat.NewDense(1, 1, nil), + 0, + ) + if err != nil { + t.Fatal(err) + } + if _, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 1}); !errors.Is(err, ErrInvalidOrder) { + t.Fatalf("split-pair order err = %v, want ErrInvalidOrder", err) + } + result, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 2}) + if err != nil { + t.Fatalf("pair-preserving reduction: %v", err) + } + if len(result.KeptPoles) != 2 || math.Abs(imag(result.KeptPoles[0])) < 1.9 || cmplx.Conj(result.KeptPoles[0]) != result.KeptPoles[1] { + t.Fatalf("kept poles = %v, want conjugate pair", result.KeptPoles) + } } -func TestModalTruncateAutoOrderAndInvalidOptions(t *testing.T) { - sys := benchSysNonSym(4, 1, 1) - result, err := ModalTruncate(sys, &ModalTruncateOptions{MaxRealPart: -1.0}) +func TestModalTruncatePreservesUnstableModesAndAutoSelects(t *testing.T) { + sys, err := New( + mat.DenseCopyOf(mat.NewDiagDense(4, []float64{0.2, 0.1, -1, -4})), + mat.NewDense(4, 1, []float64{1, 1, 1, 1}), + mat.NewDense(1, 4, []float64{1, 1, 1, 1}), + mat.NewDense(1, 1, nil), + 0, + ) + if err != nil { + t.Fatal(err) + } + if _, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 1}); !errors.Is(err, ErrInvalidOrder) { + t.Fatalf("unstable-discard order err = %v, want ErrInvalidOrder", err) + } + result, err := ModalTruncate(sys, &ModalTruncateOptions{MaxRealPart: -0.5}) if err != nil { t.Fatalf("ModalTruncate auto: %v", err) } - if result.Order == 0 || result.Order >= 4 { - t.Fatalf("auto order = %d, want partial reduction", result.Order) + if result.Order != 2 || !sameComplexApprox(result.KeptPoles, []complex128{0.2, 0.1}, 1e-12) { + t.Fatalf("auto-selected result = order %d poles %v", result.Order, result.KeptPoles) } - if _, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 5}); err != ErrInvalidOrder { + if _, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 5}); !errors.Is(err, ErrInvalidOrder) { t.Fatalf("invalid order err = %v, want ErrInvalidOrder", err) } } + +func modalTestSystem(t *testing.T) *System { + t.Helper() + sys, err := New( + mat.DenseCopyOf(mat.NewDiagDense(4, []float64{-0.1, -1, -5, -10})), + mat.NewDense(4, 1, []float64{1, 2, 3, 4}), + mat.NewDense(1, 4, []float64{4, -1, 2, 0.5}), + mat.NewDense(1, 1, []float64{0.2}), + 0, + ) + if err != nil { + t.Fatal(err) + } + return sys +} + +func sameComplexApprox(a, b []complex128, tolerance float64) bool { + if len(a) != len(b) { + return false + } + for i := range a { + if cmplx.Abs(a[i]-b[i]) > tolerance { + return false + } + } + return true +} From e7552cf47698163a12eb138009db7dd71596b18e Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:27:26 -0500 Subject: [PATCH 04/10] Update Gonum fork to v0.17.7 --- go.mod | 2 +- go.sum | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/go.mod b/go.mod index 7c11a8b..b560322 100644 --- a/go.mod +++ b/go.mod @@ -4,4 +4,4 @@ go 1.26.3 require gonum.org/v1/gonum v0.15.0 -replace gonum.org/v1/gonum => github.com/jamestjsp/gonum v0.17.6-fork +replace gonum.org/v1/gonum => github.com/jamestjsp/gonum v0.17.7-fork diff --git a/go.sum b/go.sum index f6227c7..5b429a6 100644 --- a/go.sum +++ b/go.sum @@ -1,2 +1,2 @@ -github.com/jamestjsp/gonum v0.17.6-fork h1:92RTsQJ3+diY4MSw9wNFFs9LrHnLzAzP9ZAT2C2ZYkw= -github.com/jamestjsp/gonum v0.17.6-fork/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E= +github.com/jamestjsp/gonum v0.17.7-fork h1:evxkama96reFIjJEZiHYhRGtibxK4d076+3vB/Q0t+c= +github.com/jamestjsp/gonum v0.17.7-fork/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E= From 244a2068b0cc2fb04af3247770fb255504d6e64d Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:40:04 -0500 Subject: [PATCH 05/10] Assemble physical port constraint dynamics --- README.md | 8 +- descriptor_test.go | 40 ++- frequency.go | 46 +++- physical_assembly.go | 514 +++++++++++++++++++++++++++++++++++--- physical_assembly_test.go | 149 ++++++++--- 5 files changed, 672 insertions(+), 85 deletions(-) diff --git a/README.md b/README.md index c8a2874..0028cce 100644 --- a/README.md +++ b/README.md @@ -35,7 +35,7 @@ This package is intended to be usable in production control and estimation code, - **State estimation:** Extended Kalman Filter (EKF) for nonlinear systems - **System identification:** Eigensystem Realization Algorithm (ERA) and frequency-response estimation from I/O data - **Nonlinear systems:** Jacobian linearization around operating points; Smith predictor for time-delay plants -- **Model arrays and physical assembly:** compatible model grids for parameter sweeps and port-checked physical component assembly +- **Model arrays and physical assembly:** compatible model grids for parameter sweeps and descriptor assembly of connected physical ports - **Model reduction & decomposition:** controllability/observability staircase, balanced realization, balanced truncation, stable/unstable and modal separation, modal truncation - **System norms & covariance:** H2/H-infinity norms, Hankel singular values, state covariance - **Interconnection:** series, parallel, feedback, safe feedback, append, block diagonal, named/indexed connect, FRD interconnections, sum blocks, and LFT @@ -95,7 +95,7 @@ func main() { | `NewGeneralizedModel` | Wrap a fixed or tunable block and attach analysis points | | `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with an analysis point | | `NewPhysicalComponent` | Wrap a model with named physical ports | -| `AssemblePhysical` | Validate physical port compatibility and append component models | +| `AssemblePhysical` | Assemble connected across/through port equations into a descriptor model | | `NewDescriptor` | Descriptor state-space model with explicit E matrix | | `Rss` | Random stable continuous-time state-space model | | `Drss` | Random stable discrete-time state-space model | @@ -348,8 +348,8 @@ func main() { | Function/Type | Description | |---------------|-------------| | `NewPhysicalComponent` | Wrap a model with named physical ports | -| `AssemblePhysical` | Validate physical port compatibility and append component models | -| `PhysicalPort` / `PhysicalConnection` | Port and connection metadata for physical assembly | +| `AssemblePhysical` | Assemble connected across/through port equations into a descriptor model | +| `PhysicalPort` / `PhysicalConnection` | Bind component input/output channels to physical ports and node connections | ### Transport Delays diff --git a/descriptor_test.go b/descriptor_test.go index bb2a102..e26ef80 100644 --- a/descriptor_test.go +++ b/descriptor_test.go @@ -147,14 +147,6 @@ func TestDescriptorSystem_UnsupportedOperationsRejectExplicitly(t *testing.T) { _, err := sys.TransferFunction(nil) return err }}, - {name: "FreqResponse", run: func() error { - _, err := sys.FreqResponse([]float64{1}) - return err - }}, - {name: "EvalFr", run: func() error { - _, err := sys.EvalFr(1i) - return err - }}, {name: "Simulate", run: func() error { discrete := sys.Copy() discrete.Dt = 0.1 @@ -176,6 +168,38 @@ func TestDescriptorSystem_UnsupportedOperationsRejectExplicitly(t *testing.T) { } } +func TestDescriptorSystem_FrequencyResponseUsesDescriptorPencil(t *testing.T) { + sys, err := NewDescriptor( + mat.NewDense(1, 1, []float64{-1}), + mat.NewDense(1, 1, []float64{3}), + mat.NewDense(1, 1, []float64{4}), + mat.NewDense(1, 1, []float64{0.5}), + mat.NewDense(1, 1, []float64{2}), + 0, + ) + if err != nil { + t.Fatal(err) + } + omega := []float64{0, 1, 5} + response, err := sys.FreqResponse(omega) + if err != nil { + t.Fatalf("FreqResponse: %v", err) + } + for k, frequency := range omega { + want := 12/complex(1, 2*frequency) + 0.5 + if got := response.At(k, 0, 0); cmplx.Abs(got-want) > 1e-12 { + t.Fatalf("response at %g = %v, want %v", frequency, got, want) + } + } + evaluated, err := sys.EvalFr(1i) + if err != nil { + t.Fatalf("EvalFr: %v", err) + } + if want := 12/complex(1, 2) + 0.5; cmplx.Abs(evaluated[0][0]-want) > 1e-12 { + t.Fatalf("EvalFr = %v, want %v", evaluated[0][0], want) + } +} + func TestDescriptorSystem_IdentityDescriptorUsesStandardOperations(t *testing.T) { sys, err := New( mat.NewDense(2, 2, []float64{-1, 2, 0.5, -3}), diff --git a/frequency.go b/frequency.go index 1939a06..22d1eab 100644 --- a/frequency.go +++ b/frequency.go @@ -110,9 +110,6 @@ func unwrapBodePhase(phase []float64, p, m, nw int) { } func (sys *System) FreqResponse(omega []float64) (*FreqResponseMatrix, error) { - if err := newDescriptorPolicy(sys).requireStandard("FreqResponse"); err != nil { - return nil, err - } return newFrequencyEvaluator(sys).response(omega) } @@ -137,9 +134,6 @@ func (sys *System) Bode(omega []float64, nPoints int) (*BodeResult, error) { } func (sys *System) EvalFr(s complex128) ([][]complex128, error) { - if err := newDescriptorPolicy(sys).requireStandard("EvalFr"); err != nil { - return nil, err - } return newFrequencyEvaluator(sys).eval(s) } @@ -160,6 +154,9 @@ func (e frequencyEvaluator) response(omega []float64) (*FreqResponseMatrix, erro return nil, nil } if e.sys.HasInternalDelay() { + if e.sys.IsDescriptor() { + return nil, fmt.Errorf("FreqResponse: descriptor models with internal delays are not supported: %w", ErrDescriptorUnsupported) + } resp, err := freqResponseLFT(e.sys, omega, e.p, e.m) if err != nil { return nil, err @@ -173,6 +170,16 @@ func (e frequencyEvaluator) response(omega []float64) (*FreqResponseMatrix, erro nw := len(omega) pm := e.p * e.m data := make([]complex128, nw*pm) + if e.sys.IsDescriptor() { + for k, w := range omega { + s := e.sAt(w) + if err := e.evalStateSpaceInto(s, data[k*pm:(k+1)*pm]); err != nil { + return nil, err + } + } + applyIODelayPhase(e.sys, omega, data, e.p, e.m, true) + return e.matrix(data, omega), nil + } if nw == 1 { s := e.sAt(omega[0]) if err := e.evalStateSpaceInto(s, data); err != nil { @@ -197,6 +204,9 @@ func (e frequencyEvaluator) eval(s complex128) ([][]complex128, error) { pm := e.p * e.m if e.sys.HasInternalDelay() { + if e.sys.IsDescriptor() { + return nil, fmt.Errorf("EvalFr: descriptor models with internal delays are not supported: %w", ErrDescriptorUnsupported) + } g, err := evalFrLFT(e.sys, s, e.p, e.m) if err != nil { return nil, err @@ -373,7 +383,14 @@ func evalFrSSInto(ws *ssEvalWorkspace, sys *System, s complex128, n, p, m int) e } aRaw := sys.A.RawMatrix() - if err := cResolventInto(ws.resolvent, ws.sIA, ws.invBuf, aRaw.Data, aRaw.Stride, s, n); err != nil { + var err error + if sys.E == nil { + err = cResolventInto(ws.resolvent, ws.sIA, ws.invBuf, aRaw.Data, aRaw.Stride, s, n) + } else { + eRaw := sys.E.RawMatrix() + err = cDescriptorResolventInto(ws.resolvent, ws.sIA, ws.invBuf, aRaw.Data, aRaw.Stride, eRaw.Data, eRaw.Stride, s, n) + } + if err != nil { return err } bRaw := sys.B.RawMatrix() @@ -549,6 +566,21 @@ func cResolventInto(dst, sIA, invBuf []complex128, aData []float64, aStride int, return cInvertInto(dst, invBuf, sIA, n) } +func cDescriptorResolventInto(dst, sEA, invBuf []complex128, aData []float64, aStride int, eData []float64, eStride int, s complex128, n int) error { + if n == 0 { + return nil + } + for i := range n { + row := i * n + aRow := i * aStride + eRow := i * eStride + for j := range n { + sEA[row+j] = s*complex(eData[eRow+j], 0) - complex(aData[aRow+j], 0) + } + } + return cInvertInto(dst, invBuf, sEA, n) +} + func cComputeHInto(dst, temp, resolvent []complex128, cData []float64, cStride int, bData []float64, bStride int, dData []float64, dStride, n, rows, cols int) { for i := range dst[:rows*cols] { dst[i] = 0 diff --git a/physical_assembly.go b/physical_assembly.go index 4afac0e..fe9fa65 100644 --- a/physical_assembly.go +++ b/physical_assembly.go @@ -1,6 +1,11 @@ package controlsys -import "fmt" +import ( + "fmt" + "sort" + + "gonum.org/v1/gonum/mat" +) type PhysicalPortKind int @@ -13,6 +18,8 @@ type PhysicalPort struct { Name string Kind PhysicalPortKind Dimension int + Input []int + Output []int } type PhysicalComponent struct { @@ -30,75 +37,508 @@ type PhysicalConnection struct { } func NewPhysicalComponent(name string, sys *System, ports []PhysicalPort) PhysicalComponent { - return PhysicalComponent{Name: name, System: sys.Copy(), Ports: append([]PhysicalPort(nil), ports...)} + component := PhysicalComponent{Name: name, Ports: copyPhysicalPorts(ports)} + if sys != nil { + component.System = sys.Copy() + } + return component } func AssemblePhysical(name string, components []PhysicalComponent, connections []PhysicalConnection) (*System, error) { + plan, err := newPhysicalAssemblyPlan(name, components, connections) + if err != nil { + return nil, err + } + return plan.assemble() +} + +type physicalAssemblyPlan struct { + name string + components []PhysicalComponent + ports []physicalPortBinding + groups []physicalPortGroup + stateOffsets []int + inputOffsets []int + outputOffsets []int + totalStates int + totalInputs int + totalOutputs int + internalInputs []int + internalOutputs []int + externalInputs []int + externalOutputs []int + internalInputPos map[int]int +} + +type physicalPortBinding struct { + key string + component int + port PhysicalPort + inputs []int + outputs []int +} + +type physicalPortGroup struct { + ports []int + grounded bool +} + +func newPhysicalAssemblyPlan(name string, components []PhysicalComponent, connections []PhysicalConnection) (*physicalAssemblyPlan, error) { + if name == "" { + return nil, fmt.Errorf("AssemblePhysical: empty assembly name: %w", ErrDimensionMismatch) + } if len(components) == 0 { return nil, fmt.Errorf("AssemblePhysical: no components: %w", ErrDimensionMismatch) } - byName := make(map[string]PhysicalComponent, len(components)) - for _, component := range components { + plan := &physicalAssemblyPlan{ + name: name, + components: make([]PhysicalComponent, len(components)), + stateOffsets: make([]int, len(components)), + inputOffsets: make([]int, len(components)), + outputOffsets: make([]int, len(components)), + } + byComponent := make(map[string]int, len(components)) + for i, component := range components { if component.Name == "" || component.System == nil { - return nil, fmt.Errorf("AssemblePhysical: invalid component: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("AssemblePhysical: invalid component at index %d: %w", i, ErrDimensionMismatch) } - if _, exists := byName[component.Name]; exists { + if _, exists := byComponent[component.Name]; exists { return nil, fmt.Errorf("AssemblePhysical: duplicate component %q: %w", component.Name, ErrDimensionMismatch) } - byName[component.Name] = component + if err := component.System.Validate(); err != nil { + return nil, fmt.Errorf("AssemblePhysical: component %q: %w", component.Name, err) + } + if i > 0 && component.System.Dt != components[0].System.Dt { + return nil, fmt.Errorf("AssemblePhysical: component %q has sample time %g, want %g: %w", component.Name, component.System.Dt, components[0].System.Dt, ErrDomainMismatch) + } + plan.stateOffsets[i] = plan.totalStates + plan.inputOffsets[i] = plan.totalInputs + plan.outputOffsets[i] = plan.totalOutputs + n, m, p := component.System.Dims() + plan.totalStates += n + plan.totalInputs += m + plan.totalOutputs += p + plan.components[i] = NewPhysicalComponent(component.Name, component.System, component.Ports) + byComponent[component.Name] = i + } + portIndex, err := plan.bindPorts() + if err != nil { + return nil, err + } + if err := plan.bindConnections(connections, portIndex); err != nil { + return nil, err + } + return plan, nil +} + +func (p *physicalAssemblyPlan) bindPorts() (map[string]int, error) { + portIndex := make(map[string]int) + for componentIndex, component := range p.components { + _, m, outputs := component.System.Dims() + inputCursor := 0 + outputCursor := 0 + usedInputs := make([]bool, m) + usedOutputs := make([]bool, outputs) + for _, port := range component.Ports { + if port.Name == "" || port.Dimension <= 0 { + return nil, fmt.Errorf("AssemblePhysical: invalid port on %q: %w", component.Name, ErrDimensionMismatch) + } + if port.Kind != PhysicalPortDisplacement && port.Kind != PhysicalPortEffort { + return nil, fmt.Errorf("AssemblePhysical: invalid port kind on %s.%s: %w", component.Name, port.Name, ErrDimensionMismatch) + } + key := component.Name + "." + port.Name + if _, exists := portIndex[key]; exists { + return nil, fmt.Errorf("AssemblePhysical: duplicate port %q: %w", key, ErrDimensionMismatch) + } + inputs := copyIntSlice(port.Input) + outputsForPort := copyIntSlice(port.Output) + if len(inputs) == 0 && len(outputsForPort) == 0 { + inputs = integerRange(inputCursor, port.Dimension) + outputsForPort = integerRange(outputCursor, port.Dimension) + inputCursor += port.Dimension + outputCursor += port.Dimension + } + if len(inputs) != port.Dimension || len(outputsForPort) != port.Dimension { + return nil, fmt.Errorf("AssemblePhysical: port %q must bind %d input and output channels: %w", key, port.Dimension, ErrDimensionMismatch) + } + for _, channel := range inputs { + if channel < 0 || channel >= m || usedInputs[channel] { + return nil, fmt.Errorf("AssemblePhysical: invalid or reused input channel %d on %q: %w", channel, key, ErrDimensionMismatch) + } + usedInputs[channel] = true + } + for _, channel := range outputsForPort { + if channel < 0 || channel >= outputs || usedOutputs[channel] { + return nil, fmt.Errorf("AssemblePhysical: invalid or reused output channel %d on %q: %w", channel, key, ErrDimensionMismatch) + } + usedOutputs[channel] = true + } + binding := physicalPortBinding{key: key, component: componentIndex, port: port, inputs: make([]int, port.Dimension), outputs: make([]int, port.Dimension)} + for i := range port.Dimension { + binding.inputs[i] = p.inputOffsets[componentIndex] + inputs[i] + binding.outputs[i] = p.outputOffsets[componentIndex] + outputsForPort[i] + } + portIndex[key] = len(p.ports) + p.ports = append(p.ports, binding) + } } - for _, conn := range connections { - from, err := lookupPhysicalPort(byName, conn.FromComponent, conn.FromPort) + return portIndex, nil +} + +func (p *physicalAssemblyPlan) bindConnections(connections []PhysicalConnection, portIndex map[string]int) error { + parent := make([]int, len(p.ports)) + for i := range parent { + parent[i] = i + } + active := make([]bool, len(p.ports)) + grounded := make([]bool, len(p.ports)) + seen := make(map[string]bool) + for _, connection := range connections { + from, err := physicalPortIndex(portIndex, connection.FromComponent, connection.FromPort) if err != nil { - return nil, err + return err } - if conn.Grounded { + active[from] = true + if connection.Grounded { + key := "ground:" + p.ports[from].key + if seen[key] { + return fmt.Errorf("AssemblePhysical: duplicate grounding of %s: %w", p.ports[from].key, ErrDimensionMismatch) + } + seen[key] = true + grounded[from] = true continue } - to, err := lookupPhysicalPort(byName, conn.ToComponent, conn.ToPort) + to, err := physicalPortIndex(portIndex, connection.ToComponent, connection.ToPort) if err != nil { - return nil, err + return err + } + if from == to { + return fmt.Errorf("AssemblePhysical: port %s cannot connect to itself: %w", p.ports[from].key, ErrDimensionMismatch) + } + fromPort, toPort := p.ports[from].port, p.ports[to].port + if fromPort.Kind != toPort.Kind || fromPort.Dimension != toPort.Dimension { + return fmt.Errorf("AssemblePhysical: incompatible ports %s and %s: %w", p.ports[from].key, p.ports[to].key, ErrDimensionMismatch) + } + edge := normalizedPhysicalEdge(p.ports[from].key, p.ports[to].key) + if seen[edge] { + return fmt.Errorf("AssemblePhysical: duplicate connection %s: %w", edge, ErrDimensionMismatch) + } + seen[edge] = true + active[to] = true + unionPhysicalPorts(parent, from, to) + } + groups := make(map[int]*physicalPortGroup) + for i := range p.ports { + if !active[i] { + continue + } + root := findPhysicalPort(parent, i) + group := groups[root] + if group == nil { + group = &physicalPortGroup{} + groups[root] = group + } + group.ports = append(group.ports, i) + group.grounded = group.grounded || grounded[i] + } + roots := make([]int, 0, len(groups)) + for root := range groups { + roots = append(roots, root) + } + sort.Ints(roots) + for _, root := range roots { + group := groups[root] + if !group.grounded && len(group.ports) < 2 { + return fmt.Errorf("AssemblePhysical: ungrounded node has fewer than two ports: %w", ErrDimensionMismatch) } - if from.Kind != to.Kind || from.Dimension != to.Dimension { - return nil, fmt.Errorf("AssemblePhysical: incompatible ports %s.%s and %s.%s: %w", - conn.FromComponent, conn.FromPort, conn.ToComponent, conn.ToPort, ErrDimensionMismatch) + p.groups = append(p.groups, *group) + for _, portIndex := range group.ports { + p.internalInputs = append(p.internalInputs, p.ports[portIndex].inputs...) + p.internalOutputs = append(p.internalOutputs, p.ports[portIndex].outputs...) } } + p.internalInputPos = make(map[int]int, len(p.internalInputs)) + for position, channel := range p.internalInputs { + p.internalInputPos[channel] = position + } + p.externalInputs = complementChannels(p.totalInputs, p.internalInputs) + p.externalOutputs = complementChannels(p.totalOutputs, p.internalOutputs) + return nil +} - systems := make([]*System, len(components)) - for i, component := range components { - sys := component.System.Copy() - prefixSystemMetadata(sys, component.Name) - systems[i] = sys +func (p *physicalAssemblyPlan) assemble() (*System, error) { + if len(p.groups) == 0 { + return p.appendComponents() + } + for _, component := range p.components { + if component.System.HasDelay() { + return nil, fmt.Errorf("AssemblePhysical: connected delayed components are not supported: %w", ErrDescriptorUnsupported) + } + } + a, b, c, d, e := p.aggregateMatrices() + q := len(p.internalInputs) + n := p.totalStates + nAugmented := n + q + aAugmented := newDense(nAugmented, nAugmented) + eAugmented := newDense(nAugmented, nAugmented) + bAugmented := newDense(nAugmented, len(p.externalInputs)) + setBlock(aAugmented, 0, 0, a) + setBlock(eAugmented, 0, 0, e) + for state := range n { + for position, channel := range p.internalInputs { + aAugmented.Set(state, n+position, b.At(state, channel)) + } + for position, channel := range p.externalInputs { + bAugmented.Set(state, position, b.At(state, channel)) + } + } + constraint := 0 + for _, group := range p.groups { + dimension := p.ports[group.ports[0]].port.Dimension + if group.grounded { + for _, portIndex := range group.ports { + for coordinate := range dimension { + p.addAcrossConstraint(aAugmented, bAugmented, c, d, n+constraint, portIndex, coordinate, 1) + constraint++ + } + } + continue + } + reference := group.ports[0] + for _, portIndex := range group.ports[1:] { + for coordinate := range dimension { + row := n + constraint + p.addAcrossConstraint(aAugmented, bAugmented, c, d, row, portIndex, coordinate, 1) + p.addAcrossConstraint(aAugmented, bAugmented, c, d, row, reference, coordinate, -1) + constraint++ + } + } + for coordinate := range dimension { + row := n + constraint + for _, portIndex := range group.ports { + channel := p.ports[portIndex].inputs[coordinate] + aAugmented.Set(row, n+p.internalInputPos[channel], aAugmented.At(row, n+p.internalInputPos[channel])+1) + } + constraint++ + } + } + if constraint != q { + return nil, fmt.Errorf("AssemblePhysical: generated %d constraints for %d internal variables: %w", constraint, q, ErrDimensionMismatch) + } + cAugmented := newDense(len(p.externalOutputs), nAugmented) + dAugmented := newDense(len(p.externalOutputs), len(p.externalInputs)) + for row, output := range p.externalOutputs { + for state := range n { + cAugmented.Set(row, state, c.At(output, state)) + } + for position, input := range p.internalInputs { + cAugmented.Set(row, n+position, d.At(output, input)) + } + for column, input := range p.externalInputs { + dAugmented.Set(row, column, d.At(output, input)) + } + } + var constructorB, constructorC, constructorD *mat.Dense + if len(p.externalInputs) > 0 { + constructorB = bAugmented + } + if len(p.externalOutputs) > 0 { + constructorC = cAugmented } - if len(systems) == 1 { - return systems[0], nil + if len(p.externalInputs) > 0 && len(p.externalOutputs) > 0 { + constructorD = dAugmented } - assembled := systems[0] - for _, sys := range systems[1:] { - next, err := Append(assembled, sys) + result, err := NewDescriptor(aAugmented, constructorB, constructorC, constructorD, eAugmented, p.components[0].System.Dt) + if err != nil { + return nil, err + } + result.InputName = selectedPhysicalNames(p.aggregateInputNames(), p.externalInputs) + result.OutputName = selectedPhysicalNames(p.aggregateOutputNames(), p.externalOutputs) + result.StateName = append(p.aggregateStateNames(), p.internalStateNames()...) + result.Notes = p.name + return result, nil +} + +func (p *physicalAssemblyPlan) addAcrossConstraint(a, b, c, d *mat.Dense, row, portIndex, coordinate int, factor float64) { + output := p.ports[portIndex].outputs[coordinate] + for state := range p.totalStates { + a.Set(row, state, a.At(row, state)+factor*c.At(output, state)) + } + for position, input := range p.internalInputs { + a.Set(row, p.totalStates+position, a.At(row, p.totalStates+position)+factor*d.At(output, input)) + } + for position, input := range p.externalInputs { + b.Set(row, position, b.At(row, position)+factor*d.At(output, input)) + } +} + +func (p *physicalAssemblyPlan) aggregateMatrices() (a, b, c, d, e *mat.Dense) { + a = newDense(p.totalStates, p.totalStates) + b = newDense(p.totalStates, p.totalInputs) + c = newDense(p.totalOutputs, p.totalStates) + d = newDense(p.totalOutputs, p.totalInputs) + e = newDense(p.totalStates, p.totalStates) + for i, component := range p.components { + n, _, _ := component.System.Dims() + stateOffset := p.stateOffsets[i] + inputOffset := p.inputOffsets[i] + outputOffset := p.outputOffsets[i] + setBlock(a, stateOffset, stateOffset, component.System.A) + setBlock(b, stateOffset, inputOffset, component.System.B) + setBlock(c, outputOffset, stateOffset, component.System.C) + setBlock(d, outputOffset, inputOffset, component.System.D) + if component.System.E == nil { + for state := range n { + e.Set(stateOffset+state, stateOffset+state, 1) + } + } else { + setBlock(e, stateOffset, stateOffset, component.System.E) + } + } + return a, b, c, d, e +} + +func (p *physicalAssemblyPlan) appendComponents() (*System, error) { + assembled := p.components[0].System.Copy() + prefixSystemMetadata(assembled, p.components[0].Name) + for _, component := range p.components[1:] { + next := component.System.Copy() + prefixSystemMetadata(next, component.Name) + var err error + assembled, err = Append(assembled, next) if err != nil { return nil, err } - assembled = next } + assembled.Notes = p.name return assembled, nil } -func lookupPhysicalPort(components map[string]PhysicalComponent, componentName, portName string) (PhysicalPort, error) { - component, ok := components[componentName] - if !ok { - return PhysicalPort{}, fmt.Errorf("AssemblePhysical: component %q not found: %w", componentName, ErrSignalNotFound) +func (p *physicalAssemblyPlan) aggregateInputNames() []string { + names := make([]string, 0, p.totalInputs) + for _, component := range p.components { + _, m, _ := component.System.Dims() + names = append(names, physicalSignalNames(component.System.InputName, m, component.Name, "input")...) } - for _, port := range component.Ports { - if port.Name == portName { - if port.Dimension <= 0 { - return PhysicalPort{}, fmt.Errorf("AssemblePhysical: invalid port dimension: %w", ErrDimensionMismatch) + return names +} + +func (p *physicalAssemblyPlan) aggregateOutputNames() []string { + names := make([]string, 0, p.totalOutputs) + for _, component := range p.components { + _, _, outputs := component.System.Dims() + names = append(names, physicalSignalNames(component.System.OutputName, outputs, component.Name, "output")...) + } + return names +} + +func (p *physicalAssemblyPlan) aggregateStateNames() []string { + names := make([]string, 0, p.totalStates) + for _, component := range p.components { + n, _, _ := component.System.Dims() + names = append(names, physicalSignalNames(component.System.StateName, n, component.Name, "state")...) + } + return names +} + +func (p *physicalAssemblyPlan) internalStateNames() []string { + names := make([]string, 0, len(p.internalInputs)) + for _, group := range p.groups { + for _, portIndex := range group.ports { + binding := p.ports[portIndex] + for coordinate := range binding.port.Dimension { + names = append(names, fmt.Sprintf("%s.%s.%s.through[%d]", p.name, p.components[binding.component].Name, binding.port.Name, coordinate)) } - return port, nil } } - return PhysicalPort{}, fmt.Errorf("AssemblePhysical: port %q not found on %q: %w", portName, componentName, ErrSignalNotFound) + return names +} + +func copyPhysicalPorts(ports []PhysicalPort) []PhysicalPort { + out := make([]PhysicalPort, len(ports)) + for i, port := range ports { + out[i] = port + out[i].Input = copyIntSlice(port.Input) + out[i].Output = copyIntSlice(port.Output) + } + return out +} + +func physicalPortIndex(index map[string]int, component, port string) (int, error) { + key := component + "." + port + value, ok := index[key] + if !ok { + return 0, fmt.Errorf("AssemblePhysical: port %q not found: %w", key, ErrSignalNotFound) + } + return value, nil +} + +func normalizedPhysicalEdge(a, b string) string { + if a > b { + a, b = b, a + } + return a + "<->" + b +} + +func findPhysicalPort(parent []int, value int) int { + for parent[value] != value { + parent[value] = parent[parent[value]] + value = parent[value] + } + return value +} + +func unionPhysicalPorts(parent []int, a, b int) { + aRoot := findPhysicalPort(parent, a) + bRoot := findPhysicalPort(parent, b) + if aRoot != bRoot { + if aRoot > bRoot { + aRoot, bRoot = bRoot, aRoot + } + parent[bRoot] = aRoot + } +} + +func complementChannels(total int, internal []int) []int { + used := make([]bool, total) + for _, channel := range internal { + used[channel] = true + } + external := make([]int, 0, total-len(internal)) + for channel := range total { + if !used[channel] { + external = append(external, channel) + } + } + return external +} + +func selectedPhysicalNames(names []string, channels []int) []string { + selected := make([]string, len(channels)) + for i, channel := range channels { + selected[i] = names[channel] + } + return selected +} + +func physicalSignalNames(names []string, count int, prefix, kind string) []string { + out := make([]string, count) + for i := range count { + name := fmt.Sprintf("%s%d", kind, i+1) + if i < len(names) && names[i] != "" { + name = names[i] + } + out[i] = prefix + "." + name + } + return out +} + +func integerRange(start, count int) []int { + values := make([]int, count) + for i := range count { + values[i] = start + i + } + return values } func prefixSystemMetadata(sys *System, prefix string) { diff --git a/physical_assembly_test.go b/physical_assembly_test.go index dcae691..47449e5 100644 --- a/physical_assembly_test.go +++ b/physical_assembly_test.go @@ -2,69 +2,160 @@ package controlsys import ( "errors" + "math/cmplx" "testing" "gonum.org/v1/gonum/mat" ) -func TestPhysicalAssemblyGroundedDescriptorComponent(t *testing.T) { - component := physicalTestComponent(t, "mass") - asm, err := AssemblePhysical("grounded", []PhysicalComponent{component}, []PhysicalConnection{ +func TestPhysicalAssemblyConnectedPortsChangeTransferBehavior(t *testing.T) { + left := physicalCoupledComponent(t, "left", 1) + right := physicalCoupledComponent(t, "right", 2) + assembled, err := AssemblePhysical("pair", []PhysicalComponent{left, right}, []PhysicalConnection{ + {FromComponent: "left", FromPort: "mount", ToComponent: "right", ToPort: "mount"}, + }) + if err != nil { + t.Fatalf("AssemblePhysical: %v", err) + } + if n, m, p := assembled.Dims(); n != 4 || m != 2 || p != 2 { + t.Fatalf("dims = (%d,%d,%d), want (4,2,2)", n, m, p) + } + if !assembled.IsDescriptor() { + t.Fatal("connected physical assembly must retain algebraic port constraints") + } + response, err := assembled.FreqResponse([]float64{1}) + if err != nil { + t.Fatalf("FreqResponse: %v", err) + } + want := 1 / complex(3, 2) + if got := response.At(0, 0, 1); cmplx.Abs(got-want) > 1e-12 { + t.Fatalf("cross response = %v, want %v", got, want) + } + if cmplx.Abs(response.At(0, 0, 1)) == 0 { + t.Fatal("connected assembly retained block-diagonal transfer behavior") + } + if !sameStrings(assembled.InputName, []string{"left.force", "right.force"}) { + t.Fatalf("input names = %v", assembled.InputName) + } + if assembled.Notes != "pair" { + t.Fatalf("assembly notes = %q, want pair", assembled.Notes) + } +} + +func TestPhysicalAssemblyGroundingIntroducesReactionConstraint(t *testing.T) { + component := physicalCoupledComponent(t, "mass", 1) + assembled, err := AssemblePhysical("grounded", []PhysicalComponent{component}, []PhysicalConnection{ {FromComponent: "mass", FromPort: "mount", Grounded: true}, }) if err != nil { t.Fatalf("AssemblePhysical: %v", err) } - if !asm.IsDescriptor() { - t.Fatal("assembled descriptor component should remain descriptor") + if n, m, p := assembled.Dims(); n != 2 || m != 1 || p != 1 { + t.Fatalf("dims = (%d,%d,%d), want (2,1,1)", n, m, p) } - if !matEqual(asm.A, component.System.A, 1e-12) || !matEqual(asm.E, component.System.E, 1e-12) { - t.Fatalf("assembled matrices do not match component") + response, err := assembled.FreqResponse([]float64{0, 1, 10}) + if err != nil { + t.Fatalf("FreqResponse: %v", err) } - if !sameStrings(asm.StateName, []string{"mass.x1", "mass.x2"}) { - t.Fatalf("state names = %v", asm.StateName) + for k := range response.NFreq { + if got := cmplx.Abs(response.At(k, 0, 0)); got > 1e-12 { + t.Fatalf("grounded response[%d] = %g, want zero", k, got) + } + } + if assembled.E.At(1, 1) != 0 || assembled.A.At(1, 0) != 1 || assembled.A.At(0, 1) != 1 { + t.Fatalf("ground reaction descriptor equations are incorrect: E=%v A=%v", mat.Formatted(assembled.E), mat.Formatted(assembled.A)) } } -func TestPhysicalAssemblyTwoComponentCoupling(t *testing.T) { - left := physicalTestComponent(t, "left") - right := physicalTestComponent(t, "right") - asm, err := AssemblePhysical("pair", []PhysicalComponent{left, right}, []PhysicalConnection{ - {FromComponent: "left", FromPort: "mount", ToComponent: "right", ToPort: "mount"}, +func TestPhysicalAssemblyBuildsMultiPortNode(t *testing.T) { + components := []PhysicalComponent{ + physicalCoupledComponent(t, "one", 1), + physicalCoupledComponent(t, "two", 2), + physicalCoupledComponent(t, "three", 3), + } + assembled, err := AssemblePhysical("node", components, []PhysicalConnection{ + {FromComponent: "one", FromPort: "mount", ToComponent: "two", ToPort: "mount"}, + {FromComponent: "two", FromPort: "mount", ToComponent: "three", ToPort: "mount"}, }) if err != nil { t.Fatalf("AssemblePhysical: %v", err) } - if n, m, p := asm.Dims(); n != 4 || m != 2 || p != 2 { - t.Fatalf("dims = (%d,%d,%d), want (4,2,2)", n, m, p) + if n, m, p := assembled.Dims(); n != 6 || m != 3 || p != 3 { + t.Fatalf("dims = (%d,%d,%d), want (6,3,3)", n, m, p) } - if asm.A.At(0, 1) != left.System.A.At(0, 1) || asm.A.At(2, 3) != right.System.A.At(0, 1) { - t.Fatalf("block diagonal A not preserved: %v", mat.Formatted(asm.A)) + response, err := assembled.FreqResponse([]float64{1}) + if err != nil { + t.Fatal(err) } - if !sameStrings(asm.InputName, []string{"left.force", "right.force"}) { - t.Fatalf("input names = %v", asm.InputName) + want := 1 / complex(6, 3) + if got := response.At(0, 0, 2); cmplx.Abs(got-want) > 1e-12 { + t.Fatalf("three-component cross response = %v, want %v", got, want) } } -func TestPhysicalAssemblyRejectsIncompatiblePorts(t *testing.T) { - left := physicalTestComponent(t, "left") - right := physicalTestComponent(t, "right") +func TestPhysicalAssemblyPreservesUnconnectedDescriptorComponent(t *testing.T) { + component := physicalDescriptorComponent(t, "mass") + assembled, err := AssemblePhysical("unconnected", []PhysicalComponent{component}, nil) + if err != nil { + t.Fatalf("AssemblePhysical: %v", err) + } + if !assembled.IsDescriptor() || !matEqual(assembled.A, component.System.A, 1e-12) || !matEqual(assembled.E, component.System.E, 1e-12) { + t.Fatalf("unconnected descriptor component was not preserved") + } + if !sameStrings(assembled.StateName, []string{"mass.x1", "mass.x2"}) { + t.Fatalf("state names = %v", assembled.StateName) + } +} + +func TestPhysicalAssemblyRejectsInvalidTopologies(t *testing.T) { + left := physicalCoupledComponent(t, "left", 1) + right := physicalCoupledComponent(t, "right", 2) right.Ports[0].Kind = PhysicalPortEffort - _, err := AssemblePhysical("bad", []PhysicalComponent{left, right}, []PhysicalConnection{ - {FromComponent: "left", FromPort: "mount", ToComponent: "right", ToPort: "mount"}, - }) - if !errors.Is(err, ErrDimensionMismatch) { + connection := PhysicalConnection{FromComponent: "left", FromPort: "mount", ToComponent: "right", ToPort: "mount"} + if _, err := AssemblePhysical("bad", []PhysicalComponent{left, right}, []PhysicalConnection{connection}); !errors.Is(err, ErrDimensionMismatch) { t.Fatalf("incompatible port err = %v, want ErrDimensionMismatch", err) } + right.Ports[0].Kind = PhysicalPortDisplacement + if _, err := AssemblePhysical("duplicate", []PhysicalComponent{left, right}, []PhysicalConnection{connection, connection}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("duplicate connection err = %v, want ErrDimensionMismatch", err) + } + badMapping := left + badMapping.Ports[0].Input = []int{4} + if _, err := AssemblePhysical("mapping", []PhysicalComponent{badMapping}, []PhysicalConnection{{FromComponent: "left", FromPort: "mount", Grounded: true}}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("bad channel mapping err = %v, want ErrDimensionMismatch", err) + } + if component := NewPhysicalComponent("nil", nil, nil); component.System != nil { + t.Fatal("nil-system constructor should retain a nil system for assembly validation") + } +} + +func physicalCoupledComponent(t *testing.T, name string, decay float64) PhysicalComponent { + t.Helper() + sys, err := New( + mat.NewDense(1, 1, []float64{-decay}), + mat.NewDense(1, 2, []float64{1, 1}), + mat.NewDense(2, 1, []float64{1, 1}), + mat.NewDense(2, 2, nil), + 0, + ) + if err != nil { + t.Fatal(err) + } + sys.InputName = []string{"force", "mount.force"} + sys.OutputName = []string{"position", "mount.position"} + sys.StateName = []string{"position"} + return NewPhysicalComponent(name, sys, []PhysicalPort{{ + Name: "mount", Kind: PhysicalPortDisplacement, Dimension: 1, Input: []int{1}, Output: []int{1}, + }}) } -func physicalTestComponent(t *testing.T, name string) PhysicalComponent { +func physicalDescriptorComponent(t *testing.T, name string) PhysicalComponent { t.Helper() sys, err := NewDescriptor( mat.NewDense(2, 2, []float64{-1, 0.5, -2, -3}), mat.NewDense(2, 1, []float64{1, -1}), mat.NewDense(1, 2, []float64{2, -0.5}), - mat.NewDense(1, 1, []float64{0}), + mat.NewDense(1, 1, nil), mat.NewDense(2, 2, []float64{2, 0, 0.1, 3}), 0, ) From 9efb609bbe5d07ccc3eab19a74d919ae6caffd4e Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 15:57:43 -0500 Subject: [PATCH 06/10] Bind tuning goals to loop topology --- README.md | 8 +- docs/api-mutation-audit.md | 6 +- docs/codebase-interface-diagram.md | 2 + generalized.go | 65 ++++-- generalized_test.go | 62 ++++++ tuning.go | 71 +++++-- tuning_goal.go | 309 +++++++++++++++++++++++++++-- tuning_goal_test.go | 179 ++++++++++++++++- tuning_test.go | 23 ++- 9 files changed, 663 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index 0028cce..21b7e65 100644 --- a/README.md +++ b/README.md @@ -93,7 +93,7 @@ func main() { | `StackModelArrays` | Stack equal-shaped compatible arrays along a new leading axis | | `ConcatModelArrays` | Flatten and concatenate compatible model arrays | | `NewGeneralizedModel` | Wrap a fixed or tunable block and attach analysis points | -| `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with an analysis point | +| `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with a plant-output analysis point | | `NewPhysicalComponent` | Wrap a model with named physical ports | | `AssemblePhysical` | Assemble connected across/through port equations into a descriptor model | | `NewDescriptor` | Descriptor state-space model with explicit E matrix | @@ -119,7 +119,8 @@ func main() { | Function/Type | Description | |---------------|-------------| | `NewGeneralizedModel` | Wrap a fixed or tunable block and attach analysis points | -| `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with an analysis point | +| `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with a plant-output analysis point | +| `(*GeneralizedClosedLoop).InsertAnalysisPoint` | Bind another named break to the plant input or output | | `TunableReal` | Bounded scalar parameter used by tunable blocks | | `TunableGain` | Tunable static-gain block | | `TunablePID` | Tunable PID controller block | @@ -129,7 +130,8 @@ func main() { | `NewSensitivityGoal` / `NewWeightedGainGoal` | Tuning-goal constructors for gain and sensitivity limits | | `NewLoopShapeGoal` / `NewMarginGoal` | Tuning-goal constructors for loop-shape and robustness constraints | | `NewPoleGoal` / `NewOvershootGoal` | Tuning-goal constructors for pole-location and step-response constraints | -| `Systune` / `Looptune` | Fixed-structure tuning over free tunable parameters | +| `GridTune` | Bounded Cartesian-grid tuning over free parameters and point-specific goals | +| `Systune` / `Looptune` | Compatibility wrappers around `GridTune` | ### Frequency Response & Plotting diff --git a/docs/api-mutation-audit.md b/docs/api-mutation-audit.md index 47f0db3..defd8dc 100644 --- a/docs/api-mutation-audit.md +++ b/docs/api-mutation-audit.md @@ -288,6 +288,7 @@ ownership clarity, and release readiness, not shrinking the toolbox shape. | `(*GeneralizedModel).AnalysisPoint` | `pure` | Returns value object. | | `(*GeneralizedModel).CurrentSystem` | `returns-mutable` | Returns current numeric model. | | `NewGeneralizedClosedLoop` | `copy-in`, `returns-mutable` | Copies plant; stores controller block. | +| `(*GeneralizedClosedLoop).InsertAnalysisPoint` | `mutates` | Binds a named plant-input or plant-output loop break. | | `(*GeneralizedClosedLoop).AnalysisPoint` | `pure` | Returns value object. | | `(*GeneralizedClosedLoop).OpenLoop` | `returns-mutable` | Builds current open-loop model. | | `(*GeneralizedClosedLoop).ClosedLoop` | `returns-mutable` | Alias for complementary sensitivity. | @@ -326,8 +327,9 @@ ownership clarity, and release readiness, not shrinking the toolbox shape. | `(*TunableSS).RandomSample` | `copy-out`, `returns-mutable` | Returns sampled block copy. | | `(*TunableSS).FreeParameters` | `alias-risk` | Returns pointers to internal free parameters. | | `(*TunableSS).SampleBlock` | `copy-out`, `returns-mutable` | Interface wrapper around `Sample`. | -| `Systune` | `view-in`, `returns-mutable` | Returns mutable controller/closed-loop result. | -| `Looptune` | `view-in`, `returns-mutable` | Same implementation as `Systune`. | +| `GridTune` | `view-in`, `returns-mutable` | Returns a mutable result from a bounded Cartesian search. | +| `Systune` | `view-in`, `returns-mutable` | Compatibility wrapper around `GridTune`. | +| `Looptune` | `view-in`, `returns-mutable` | Compatibility wrapper around `GridTune`. | ## Controller, PID, EKF, Physical, and Goal APIs diff --git a/docs/codebase-interface-diagram.md b/docs/codebase-interface-diagram.md index dd4b988..eeb30bf 100644 --- a/docs/codebase-interface-diagram.md +++ b/docs/codebase-interface-diagram.md @@ -229,6 +229,8 @@ classDiagram } class GeneralizedClosedLoop { + +InsertAnalysisPoint() + +AnalysisPoint() +OpenLoop() +ClosedLoop() +Sensitivity() diff --git a/generalized.go b/generalized.go index 24d6c3b..3e594d5 100644 --- a/generalized.go +++ b/generalized.go @@ -31,8 +31,19 @@ type GeneralizedModel struct { analysisPoints map[string]AnalysisPoint } +// AnalysisPointLocation identifies the signal where a feedback loop is broken. +type AnalysisPointLocation uint8 + +const ( + AnalysisPointUnspecified AnalysisPointLocation = iota + AnalysisPointPlantOutput + AnalysisPointPlantInput +) + +// AnalysisPoint binds a name to a loop location. type AnalysisPoint struct { - Name string + Name string + Location AnalysisPointLocation } func NewGeneralizedModel(name string, block any) (*GeneralizedModel, error) { @@ -125,6 +136,9 @@ func NewGeneralizedClosedLoop(name string, plant *System, controller any, analys if plant == nil { return nil, fmt.Errorf("NewGeneralizedClosedLoop: nil plant: %w", ErrDimensionMismatch) } + if analysisPoint == "" { + return nil, fmt.Errorf("NewGeneralizedClosedLoop: analysis point is empty: %w", ErrDimensionMismatch) + } g := &GeneralizedClosedLoop{ name: name, plant: plant.Copy(), @@ -135,10 +149,28 @@ func NewGeneralizedClosedLoop(name string, plant *System, controller any, analys if tunable, ok := controller.(TunableBlock); ok { g.tunableController = tunable } - g.analysisPoints[analysisPoint] = AnalysisPoint{Name: analysisPoint} + g.analysisPoints[analysisPoint] = AnalysisPoint{Name: analysisPoint, Location: AnalysisPointPlantOutput} return g, nil } +// InsertAnalysisPoint binds name to the plant-input or plant-output loop break. +func (g *GeneralizedClosedLoop) InsertAnalysisPoint(name string, location AnalysisPointLocation) error { + if g == nil { + return fmt.Errorf("GeneralizedClosedLoop.InsertAnalysisPoint: nil model: %w", ErrDimensionMismatch) + } + if name == "" { + return fmt.Errorf("GeneralizedClosedLoop.InsertAnalysisPoint: name is empty: %w", ErrDimensionMismatch) + } + if location != AnalysisPointPlantOutput && location != AnalysisPointPlantInput { + return fmt.Errorf("GeneralizedClosedLoop.InsertAnalysisPoint: invalid location %d: %w", location, ErrDimensionMismatch) + } + if g.analysisPoints == nil { + g.analysisPoints = make(map[string]AnalysisPoint) + } + g.analysisPoints[name] = AnalysisPoint{Name: name, Location: location} + return nil +} + func (g *GeneralizedClosedLoop) AnalysisPoint(name string) (AnalysisPoint, error) { if g == nil { return AnalysisPoint{}, fmt.Errorf("GeneralizedClosedLoop.AnalysisPoint: nil model: %w", ErrDimensionMismatch) @@ -151,14 +183,22 @@ func (g *GeneralizedClosedLoop) AnalysisPoint(name string) (AnalysisPoint, error } func (g *GeneralizedClosedLoop) OpenLoop(name string) (*System, error) { - if _, err := g.AnalysisPoint(name); err != nil { + point, err := g.AnalysisPoint(name) + if err != nil { return nil, err } controller, err := g.controller.CurrentSystem() if err != nil { return nil, err } - return Series(controller, g.plant) + switch point.Location { + case AnalysisPointPlantOutput: + return Series(controller, g.plant) + case AnalysisPointPlantInput: + return Series(g.plant, controller) + default: + return nil, fmt.Errorf("GeneralizedClosedLoop.OpenLoop: analysis point %q has no loop location: %w", name, ErrDimensionMismatch) + } } func (g *GeneralizedClosedLoop) ClosedLoop(name string) (*System, error) { @@ -166,26 +206,23 @@ func (g *GeneralizedClosedLoop) ClosedLoop(name string) (*System, error) { } func (g *GeneralizedClosedLoop) ComplementarySensitivity(name string) (*System, error) { - L, err := g.OpenLoop(name) + loop, err := g.OpenLoop(name) if err != nil { return nil, err } - return Feedback(L, nil, -1) + return Feedback(loop, nil, -1) } func (g *GeneralizedClosedLoop) Sensitivity(name string) (*System, error) { - T, err := g.ComplementarySensitivity(name) + loop, err := g.OpenLoop(name) if err != nil { return nil, err } - one, err := NewGain(eyeDense(1), T.Dt) - if err != nil { - return nil, err + _, inputs, outputs := loop.Dims() + if inputs != outputs { + return nil, fmt.Errorf("GeneralizedClosedLoop.Sensitivity: loop at %q is %dx%d: %w", name, outputs, inputs, ErrDimensionMismatch) } - negT := T.Copy() - negT.C.Scale(-1, negT.C) - negT.D.Scale(-1, negT.D) - return Parallel(one, negT) + return Feedback(makeIdentityGain(outputs, loop.Dt), loop, -1) } func (g *GeneralizedClosedLoop) primaryAnalysisPointName() string { diff --git a/generalized_test.go b/generalized_test.go index e9f1f67..085f963 100644 --- a/generalized_test.go +++ b/generalized_test.go @@ -90,3 +90,65 @@ func TestGeneralizedClosedLoopAnalysisHelpers(t *testing.T) { } } } + +func TestGeneralizedClosedLoopAnalysisPointsBindDistinctMIMOBreaks(t *testing.T) { + plant, err := NewGain(mat.NewDense(2, 2, []float64{1, 2, 0, 1}), 0) + if err != nil { + t.Fatal(err) + } + controller, err := NewGain(mat.NewDense(2, 2, []float64{1, 0, 3, 1}), 0) + if err != nil { + t.Fatal(err) + } + loop, err := NewGeneralizedClosedLoop("mimo", plant, controller, "plant_output") + if err != nil { + t.Fatal(err) + } + if err := loop.InsertAnalysisPoint("plant_input", AnalysisPointPlantInput); err != nil { + t.Fatalf("InsertAnalysisPoint: %v", err) + } + + outputLoop, err := loop.OpenLoop("plant_output") + if err != nil { + t.Fatalf("output OpenLoop: %v", err) + } + inputLoop, err := loop.OpenLoop("plant_input") + if err != nil { + t.Fatalf("input OpenLoop: %v", err) + } + wantOutput := mat.NewDense(2, 2, []float64{7, 2, 3, 1}) + wantInput := mat.NewDense(2, 2, []float64{1, 2, 3, 7}) + if !mat.EqualApprox(outputLoop.D, wantOutput, 1e-14) { + t.Fatalf("plant-output loop =\n%v\nwant\n%v", mat.Formatted(outputLoop.D), mat.Formatted(wantOutput)) + } + if !mat.EqualApprox(inputLoop.D, wantInput, 1e-14) { + t.Fatalf("plant-input loop =\n%v\nwant\n%v", mat.Formatted(inputLoop.D), mat.Formatted(wantInput)) + } + + for _, name := range []string{"plant_output", "plant_input"} { + sensitivity, err := loop.Sensitivity(name) + if err != nil { + t.Fatalf("Sensitivity(%q): %v", name, err) + } + complementary, err := loop.ComplementarySensitivity(name) + if err != nil { + t.Fatalf("ComplementarySensitivity(%q): %v", name, err) + } + sum := mat.NewDense(2, 2, nil) + sum.Add(sensitivity.D, complementary.D) + if !mat.EqualApprox(sum, eyeDense(2), 1e-12) { + t.Fatalf("S+T at %q =\n%v\nwant identity", name, mat.Formatted(sum)) + } + } +} + +func TestGeneralizedClosedLoopRejectsInvalidAnalysisPointLocation(t *testing.T) { + plant := makeSISO(-1, 1, 1, 0) + loop, err := NewGeneralizedClosedLoop("loop", plant, plant, "output") + if err != nil { + t.Fatal(err) + } + if err := loop.InsertAnalysisPoint("bad", AnalysisPointUnspecified); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("InsertAnalysisPoint error = %v, want ErrDimensionMismatch", err) + } +} diff --git a/tuning.go b/tuning.go index 6416f12..663d0aa 100644 --- a/tuning.go +++ b/tuning.go @@ -7,10 +7,12 @@ import ( ) type SystuneOptions struct { - GridPoints int + GridPoints int + MaxEvaluations int } type SystuneResult struct { + Method string Pass bool Score float64 Iterations int @@ -21,34 +23,47 @@ type SystuneResult struct { } func Systune(model *GeneralizedClosedLoop, goals []TuningGoal, opts *SystuneOptions) (*SystuneResult, error) { - return tuneFixedStructure(model, goals, opts) + return GridTune(model, goals, opts) } func Looptune(model *GeneralizedClosedLoop, goals []TuningGoal, opts *SystuneOptions) (*SystuneResult, error) { - return tuneFixedStructure(model, goals, opts) + return GridTune(model, goals, opts) } -func tuneFixedStructure(model *GeneralizedClosedLoop, goals []TuningGoal, opts *SystuneOptions) (*SystuneResult, error) { +// GridTune searches a bounded Cartesian grid of the controller's free parameters. +func GridTune(model *GeneralizedClosedLoop, goals []TuningGoal, opts *SystuneOptions) (*SystuneResult, error) { if model == nil { - return nil, fmt.Errorf("Systune: nil model: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("GridTune: nil model: %w", ErrDimensionMismatch) } controller := model.tunableController if controller == nil { - return nil, fmt.Errorf("Systune: controller is not tunable: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("GridTune: controller is not tunable: %w", ErrDimensionMismatch) } if len(goals) == 0 { - return nil, fmt.Errorf("Systune: no goals: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("GridTune: no goals: %w", ErrDimensionMismatch) } gridPoints := 5 + maxEvaluations := 100_000 if opts != nil && opts.GridPoints > 0 { gridPoints = opts.GridPoints } + if opts != nil && opts.MaxEvaluations > 0 { + maxEvaluations = opts.MaxEvaluations + } params := controller.FreeParameters() if len(params) == 0 { - return nil, fmt.Errorf("Systune: no free tunable parameters: %w", ErrDimensionMismatch) + return nil, fmt.Errorf("GridTune: no free tunable parameters: %w", ErrDimensionMismatch) + } + evaluationCount := 1 + for _, param := range params { + count := len(parameterGrid(param, gridPoints)) + if evaluationCount > maxEvaluations/count { + return nil, fmt.Errorf("GridTune: Cartesian grid exceeds %d evaluations: %w", maxEvaluations, ErrDimensionMismatch) + } + evaluationCount *= count } - best := &SystuneResult{Score: math.Inf(1)} + best := &SystuneResult{Method: "cartesian-grid", Score: math.Inf(1)} iterations := 0 values := make(map[string]float64, len(params)) var search func(int) error @@ -66,7 +81,7 @@ func tuneFixedStructure(model *GeneralizedClosedLoop, goals []TuningGoal, opts * if err != nil { return err } - goalResults, score, pass, err := evaluateTuningGoals(closed, goals) + goalResults, score, pass, err := evaluateTuningGoals(&candidate, closed, goals) if err != nil { return err } @@ -76,6 +91,7 @@ func tuneFixedStructure(model *GeneralizedClosedLoop, goals []TuningGoal, opts * return err } best = &SystuneResult{ + Method: "cartesian-grid", Pass: pass, Score: score, Parameters: copyStringFloatMap(values), @@ -103,12 +119,30 @@ func tuneFixedStructure(model *GeneralizedClosedLoop, goals []TuningGoal, opts * return best, nil } -func evaluateTuningGoals(sys *System, goals []TuningGoal) ([]TuningGoalResult, float64, bool, error) { +func evaluateTuningGoals(model *GeneralizedClosedLoop, primaryClosedLoop *System, goals []TuningGoal) ([]TuningGoalResult, float64, bool, error) { results := make([]TuningGoalResult, len(goals)) score := 0.0 pass := true + primary := model.primaryAnalysisPointName() + cache := map[tuningGoalResponseKey]*System{ + {point: primary, response: tuningGoalClosedLoopResponse}: primaryClosedLoop, + } for i, goal := range goals { - result, err := goal.Evaluate(sys) + point := goal.spec.AnalysisPoint + if point == "" { + point = primary + } + key := tuningGoalResponseKey{point: point, response: tuningGoalResponseForType(goal.spec.Type)} + sys := cache[key] + if sys == nil { + var err error + sys, err = tuningGoalSystem(model, goal.spec) + if err != nil { + return nil, 0, false, err + } + cache[key] = sys + } + result, err := goal.evaluateSystem(sys) if err != nil { return nil, 0, false, err } @@ -121,14 +155,13 @@ func evaluateTuningGoals(sys *System, goals []TuningGoal) ([]TuningGoalResult, f return results, score, pass, nil } +type tuningGoalResponseKey struct { + point string + response tuningGoalResponse +} + func goalViolation(result TuningGoalResult) float64 { - if result.Pass { - return 0 - } - if result.Limit == 0 { - return math.Abs(result.Value) - } - return math.Abs((result.Value - result.Limit) / result.Limit) + return result.Violation } func uniqueFreeTunableReals(params [][]*TunableReal) []*TunableReal { diff --git a/tuning_goal.go b/tuning_goal.go index 90ce7ae..2f4e9f2 100644 --- a/tuning_goal.go +++ b/tuning_goal.go @@ -3,7 +3,8 @@ package controlsys import ( "fmt" "math" - "math/cmplx" + + "gonum.org/v1/gonum/mat" ) type TuningGoalType int @@ -20,10 +21,14 @@ const ( ) type TuningGoalSpec struct { - Name string - Type TuningGoalType - Max float64 - Min float64 + Name string + Type TuningGoalType + Max float64 + Min float64 + AnalysisPoint string + Omega []float64 + InputWeight *System + OutputWeight *System } type TuningGoal struct { @@ -35,6 +40,7 @@ type TuningGoalResult struct { Pass bool Value float64 Limit float64 + Violation float64 Diagnostics map[string]float64 } @@ -45,6 +51,28 @@ func NewTuningGoal(spec TuningGoalSpec) (TuningGoal, error) { if spec.Type != TuningGoalPole && (spec.Max < 0 || spec.Min < 0) { return TuningGoal{}, fmt.Errorf("NewTuningGoal: negative bound: %w", ErrDimensionMismatch) } + if spec.Type < TuningGoalTracking || spec.Type > TuningGoalOvershoot { + return TuningGoal{}, fmt.Errorf("NewTuningGoal: unsupported goal type %d: %w", spec.Type, ErrDimensionMismatch) + } + if spec.Type == TuningGoalLoopShape && spec.Min > spec.Max { + return TuningGoal{}, fmt.Errorf("NewTuningGoal: loop-shape minimum %g exceeds maximum %g: %w", spec.Min, spec.Max, ErrDimensionMismatch) + } + if spec.Omega != nil && !tuningGoalUsesFrequencyGrid(spec.Type) { + return TuningGoal{}, fmt.Errorf("NewTuningGoal: goal type %d does not use a frequency grid: %w", spec.Type, ErrDimensionMismatch) + } + if err := validateTuningGoalFrequencyGrid(spec.Omega); err != nil { + return TuningGoal{}, err + } + if (spec.InputWeight != nil || spec.OutputWeight != nil) && spec.Type != TuningGoalWeightedGain { + return TuningGoal{}, fmt.Errorf("NewTuningGoal: weights require a weighted-gain goal: %w", ErrDimensionMismatch) + } + spec.Omega = copyFloatSlice(spec.Omega) + if spec.InputWeight != nil { + spec.InputWeight = spec.InputWeight.Copy() + } + if spec.OutputWeight != nil { + spec.OutputWeight = spec.OutputWeight.Copy() + } return TuningGoal{spec: spec}, nil } @@ -109,10 +137,14 @@ func (g TuningGoal) Name() string { } func (g TuningGoal) Evaluate(model any) (TuningGoalResult, error) { - sys, err := tuningGoalSystem(model) + sys, err := tuningGoalSystem(model, g.spec) if err != nil { return TuningGoalResult{}, err } + return g.evaluateSystem(sys) +} + +func (g TuningGoal) evaluateSystem(sys *System) (TuningGoalResult, error) { switch g.spec.Type { case TuningGoalTracking: return g.evaluateTracking(sys) @@ -131,19 +163,58 @@ func (g TuningGoal) Evaluate(model any) (TuningGoalResult, error) { } } -func tuningGoalSystem(model any) (*System, error) { +func tuningGoalSystem(model any, spec TuningGoalSpec) (*System, error) { switch v := model.(type) { case *System: return v.Copy(), nil case *GeneralizedModel: return v.CurrentSystem() case *GeneralizedClosedLoop: - return v.ClosedLoop(v.primaryAnalysisPointName()) + point := spec.AnalysisPoint + if point == "" { + point = v.primaryAnalysisPointName() + } + switch tuningGoalResponseForType(spec.Type) { + case tuningGoalSensitivityResponse: + return v.Sensitivity(point) + case tuningGoalOpenLoopResponse: + return v.OpenLoop(point) + default: + return v.ClosedLoop(point) + } default: return nil, fmt.Errorf("TuningGoal.Evaluate: unsupported model %T: %w", model, ErrDimensionMismatch) } } +type tuningGoalResponse uint8 + +const ( + tuningGoalClosedLoopResponse tuningGoalResponse = iota + tuningGoalSensitivityResponse + tuningGoalOpenLoopResponse +) + +func tuningGoalResponseForType(goalType TuningGoalType) tuningGoalResponse { + switch goalType { + case TuningGoalRejection, TuningGoalSensitivity: + return tuningGoalSensitivityResponse + case TuningGoalLoopShape, TuningGoalMargin: + return tuningGoalOpenLoopResponse + default: + return tuningGoalClosedLoopResponse + } +} + +func tuningGoalUsesFrequencyGrid(goalType TuningGoalType) bool { + switch goalType { + case TuningGoalRejection, TuningGoalSensitivity, TuningGoalWeightedGain, TuningGoalLoopShape: + return true + default: + return false + } +} + func firstAnalysisPointName(points map[string]AnalysisPoint) string { for name := range points { return name @@ -161,7 +232,7 @@ func (g TuningGoal) evaluateTracking(sys *System) (TuningGoalResult, error) { } func (g TuningGoal) evaluateMaxGain(sys *System) (TuningGoalResult, error) { - value, err := maxFrequencyGain(sys, nil) + value, err := maxFrequencyGain(sys, g.spec.Omega, g.spec.OutputWeight, g.spec.InputWeight) if err != nil { return TuningGoalResult{}, err } @@ -169,12 +240,18 @@ func (g TuningGoal) evaluateMaxGain(sys *System) (TuningGoalResult, error) { } func (g TuningGoal) evaluateLoopShape(sys *System) (TuningGoalResult, error) { - value, err := maxFrequencyGain(sys, nil) + minimum, maximum, err := frequencyGainRange(sys, g.spec.Omega, nil, nil) if err != nil { return TuningGoalResult{}, err } - pass := value >= g.spec.Min && value <= g.spec.Max - return g.scalarResult(value, g.spec.Max, pass, map[string]float64{"max_gain": value, "min_gain": g.spec.Min}), nil + pass := minimum >= g.spec.Min && maximum <= g.spec.Max + result := g.scalarResult(maximum, g.spec.Max, pass, map[string]float64{ + "sampled_min_gain": minimum, + "sampled_max_gain": maximum, + "required_min_gain": g.spec.Min, + }) + result.Violation = math.Max(normalizedLowerViolation(minimum, g.spec.Min), normalizedUpperViolation(maximum, g.spec.Max)) + return result, nil } func (g TuningGoal) evaluateMargin(sys *System) (TuningGoalResult, error) { @@ -184,7 +261,9 @@ func (g TuningGoal) evaluateMargin(sys *System) (TuningGoalResult, error) { } pass := margin.GainMargin >= g.spec.Min && margin.PhaseMargin >= g.spec.Max diag := map[string]float64{"gain_margin_db": margin.GainMargin, "phase_margin_deg": margin.PhaseMargin} - return g.scalarResult(math.Min(margin.GainMargin, margin.PhaseMargin), g.spec.Max, pass, diag), nil + result := g.scalarResult(math.Min(margin.GainMargin, margin.PhaseMargin), g.spec.Max, pass, diag) + result.Violation = math.Max(normalizedLowerViolation(margin.GainMargin, g.spec.Min), normalizedLowerViolation(margin.PhaseMargin, g.spec.Max)) + return result, nil } func (g TuningGoal) evaluatePole(sys *System) (TuningGoalResult, error) { @@ -219,7 +298,34 @@ func (g TuningGoal) evaluateOvershoot(sys *System) (TuningGoalResult, error) { } func (g TuningGoal) scalarResult(value, limit float64, pass bool, diag map[string]float64) TuningGoalResult { - return TuningGoalResult{GoalName: g.spec.Name, Pass: pass, Value: value, Limit: limit, Diagnostics: diag} + return TuningGoalResult{ + GoalName: g.spec.Name, + Pass: pass, + Value: value, + Limit: limit, + Violation: normalizedUpperViolation(value, limit), + Diagnostics: diag, + } +} + +func normalizedUpperViolation(value, limit float64) float64 { + if value <= limit { + return 0 + } + if limit == 0 { + return value - limit + } + return (value - limit) / math.Abs(limit) +} + +func normalizedLowerViolation(value, limit float64) float64 { + if value >= limit { + return 0 + } + if limit == 0 { + return limit - value + } + return (limit - value) / math.Abs(limit) } func maxDCErrorFromOne(dc interface { @@ -242,23 +348,182 @@ func maxDCErrorFromOne(dc interface { return maxErr } -func maxFrequencyGain(sys *System, omega []float64) (float64, error) { +func validateTuningGoalFrequencyGrid(omega []float64) error { + if omega != nil && len(omega) == 0 { + return fmt.Errorf("NewTuningGoal: frequency grid is empty: %w", ErrDimensionMismatch) + } + for i, w := range omega { + if math.IsNaN(w) || math.IsInf(w, 0) || w < 0 { + return fmt.Errorf("NewTuningGoal: invalid frequency omega[%d]=%g: %w", i, w, ErrDimensionMismatch) + } + if i > 0 && w <= omega[i-1] { + return fmt.Errorf("NewTuningGoal: frequencies must be strictly increasing: %w", ErrDimensionMismatch) + } + } + return nil +} + +func maxFrequencyGain(sys *System, omega []float64, outputWeight, inputWeight *System) (float64, error) { + _, maximum, err := frequencyGainRange(sys, omega, outputWeight, inputWeight) + return maximum, err +} + +func frequencyGainRange(sys *System, omega []float64, outputWeight, inputWeight *System) (float64, float64, error) { if omega == nil { omega = logspace(-2, 2, 80) } + if len(omega) == 0 { + return 0, 0, fmt.Errorf("frequency gain: empty grid: %w", ErrDimensionMismatch) + } resp, err := sys.FreqResponse(omega) if err != nil { - return 0, err + return 0, 0, err + } + var outputResponse, inputResponse *FreqResponseMatrix + if outputWeight != nil { + if err := domainMatch(sys, outputWeight); err != nil { + return 0, 0, fmt.Errorf("output weight: %w", err) + } + outputResponse, err = outputWeight.FreqResponse(omega) + if err != nil { + return 0, 0, err + } + } + if inputWeight != nil { + if err := domainMatch(sys, inputWeight); err != nil { + return 0, 0, fmt.Errorf("input weight: %w", err) + } + inputResponse, err = inputWeight.FreqResponse(omega) + if err != nil { + return 0, 0, err + } } maxGain := 0.0 + minGain := math.Inf(1) + responseData := make([]complex128, resp.P*resp.M) + var singularValues complexSingularValueWorkspace for k := range omega { - for i := 0; i < resp.P; i++ { - for j := 0; j < resp.M; j++ { - if gain := cmplx.Abs(resp.At(k, i, j)); gain > maxGain { - maxGain = gain - } + gain := complexResponseAt(resp, k, responseData) + if outputResponse != nil { + gain, err = multiplyComplexMatrices(complexResponseAt(outputResponse, k, nil), gain) + if err != nil { + return 0, 0, fmt.Errorf("output weight: %w", err) + } + } + if inputResponse != nil { + gain, err = multiplyComplexMatrices(gain, complexResponseAt(inputResponse, k, nil)) + if err != nil { + return 0, 0, fmt.Errorf("input weight: %w", err) + } + } + sigma, err := singularValues.maximum(gain) + if err != nil { + return 0, 0, err + } + if sigma > maxGain { + maxGain = sigma + } + if sigma < minGain { + minGain = sigma + } + } + return minGain, maxGain, nil +} + +type complexMatrix struct { + rows int + cols int + data []complex128 +} + +func complexResponseAt(response *FreqResponseMatrix, frequency int, data []complex128) complexMatrix { + if len(data) != response.P*response.M { + data = make([]complex128, response.P*response.M) + } + for i := range response.P { + for j := range response.M { + data[i*response.M+j] = response.At(frequency, i, j) + } + } + return complexMatrix{rows: response.P, cols: response.M, data: data} +} + +func multiplyComplexMatrices(a, b complexMatrix) (complexMatrix, error) { + if a.cols != b.rows { + return complexMatrix{}, fmt.Errorf("matrix dimensions %dx%d and %dx%d: %w", a.rows, a.cols, b.rows, b.cols, ErrDimensionMismatch) + } + result := complexMatrix{rows: a.rows, cols: b.cols, data: make([]complex128, a.rows*b.cols)} + for i := range a.rows { + for k := range a.cols { + aik := a.data[i*a.cols+k] + for j := range b.cols { + result.data[i*result.cols+j] += aik * b.data[k*b.cols+j] } } } - return maxGain, nil + return result, nil +} + +type complexSingularValueWorkspace struct { + realForm *mat.Dense + svd mat.SVD +} + +func (w *complexSingularValueWorkspace) maximum(a complexMatrix) (float64, error) { + if a.rows == 0 || a.cols == 0 { + return 0, nil + } + if a.rows == 1 || a.cols == 1 { + norm := 0.0 + for _, value := range a.data { + norm = math.Hypot(norm, math.Hypot(real(value), imag(value))) + } + return norm, nil + } + if a.rows == 2 && a.cols == 2 { + return maximumComplex2x2SingularValue(a), nil + } + rows, cols := 2*a.rows, 2*a.cols + if w.realForm == nil { + w.realForm = mat.NewDense(rows, cols, nil) + } else if r, c := w.realForm.Dims(); r != rows || c != cols { + w.realForm.Reset() + w.realForm.ReuseAs(rows, cols) + } + for i := range a.rows { + for j := range a.cols { + value := a.data[i*a.cols+j] + w.realForm.Set(i, j, real(value)) + w.realForm.Set(i, j+a.cols, -imag(value)) + w.realForm.Set(i+a.rows, j, imag(value)) + w.realForm.Set(i+a.rows, j+a.cols, real(value)) + } + } + if ok := w.svd.Factorize(w.realForm, mat.SVDNone); !ok { + return 0, fmt.Errorf("maximum singular value: decomposition failed: %w", ErrSingularTransform) + } + values := w.svd.Values(nil) + return values[0], nil +} + +func maximumComplex2x2SingularValue(a complexMatrix) float64 { + scale := 0.0 + for _, value := range a.data { + scale = math.Max(scale, math.Hypot(real(value), imag(value))) + } + if scale == 0 { + return 0 + } + a00 := a.data[0] / complex(scale, 0) + a01 := a.data[1] / complex(scale, 0) + a10 := a.data[2] / complex(scale, 0) + a11 := a.data[3] / complex(scale, 0) + frobeniusSquared := complexMagnitudeSquared(a00) + complexMagnitudeSquared(a01) + complexMagnitudeSquared(a10) + complexMagnitudeSquared(a11) + determinantSquared := complexMagnitudeSquared(a00*a11 - a01*a10) + discriminant := math.Max(0, frobeniusSquared*frobeniusSquared-4*determinantSquared) + return scale * math.Sqrt((frobeniusSquared+math.Sqrt(discriminant))/2) +} + +func complexMagnitudeSquared(value complex128) float64 { + return real(value)*real(value) + imag(value)*imag(value) } diff --git a/tuning_goal_test.go b/tuning_goal_test.go index aa899d8..4af808e 100644 --- a/tuning_goal_test.go +++ b/tuning_goal_test.go @@ -1,6 +1,12 @@ package controlsys -import "testing" +import ( + "errors" + "math" + "testing" + + "gonum.org/v1/gonum/mat" +) func TestTuningGoalsEvaluatePassFailFamilies(t *testing.T) { sys := makeSISO(-2, 2, 1, 0) @@ -77,3 +83,174 @@ func TestTuningGoalValidation(t *testing.T) { t.Fatalf("expected pole goal failure, got %#v", fail) } } + +func TestTuningGoalUsesMaximumSingularValueForMIMO(t *testing.T) { + sys, err := NewGain(mat.NewDense(2, 2, []float64{1, 1, 1, 1}), 0) + if err != nil { + t.Fatal(err) + } + result, err := NewWeightedGainGoal("sigma_max", 1.5).Evaluate(sys) + if err != nil { + t.Fatalf("Evaluate: %v", err) + } + if result.Pass || math.Abs(result.Value-2) > 1e-12 { + t.Fatalf("result = %#v, want sigma_max=2 failure", result) + } +} + +func TestTuningGoalHonorsFrequencyGridAndDynamicWeights(t *testing.T) { + lowpass := makeSISO(-1, 1, 1, 0) + bandGoal, err := NewTuningGoal(TuningGoalSpec{ + Name: "high_frequency", + Type: TuningGoalWeightedGain, + Max: 0.1, + Omega: []float64{100}, + }) + if err != nil { + t.Fatal(err) + } + bandResult, err := bandGoal.Evaluate(lowpass) + if err != nil { + t.Fatalf("band Evaluate: %v", err) + } + if !bandResult.Pass { + t.Fatalf("high-frequency result = %#v, want pass", bandResult) + } + + weight, err := NewGain(mat.NewDense(1, 1, []float64{2}), 0) + if err != nil { + t.Fatal(err) + } + weightedGoal, err := NewTuningGoal(TuningGoalSpec{ + Name: "weighted", + Type: TuningGoalWeightedGain, + Max: 1.5, + Omega: []float64{0}, + OutputWeight: weight, + }) + if err != nil { + t.Fatal(err) + } + weightedResult, err := weightedGoal.Evaluate(lowpass) + if err != nil { + t.Fatalf("weighted Evaluate: %v", err) + } + if weightedResult.Pass || math.Abs(weightedResult.Value-2) > 1e-12 { + t.Fatalf("weighted result = %#v, want gain=2 failure", weightedResult) + } +} + +func TestLoopShapeChecksBothSampledEnvelopeBounds(t *testing.T) { + goal, err := NewTuningGoal(TuningGoalSpec{ + Name: "envelope", + Type: TuningGoalLoopShape, + Min: 0.5, + Max: 2, + Omega: []float64{0, 100}, + }) + if err != nil { + t.Fatal(err) + } + result, err := goal.Evaluate(makeSISO(-1, 1, 1, 0)) + if err != nil { + t.Fatalf("Evaluate: %v", err) + } + if result.Pass || result.Diagnostics["sampled_min_gain"] >= 0.5 || result.Violation <= 0 { + t.Fatalf("result = %#v, want lower-envelope violation", result) + } +} + +func TestTuningGoalRoutesGeneralizedLoopResponses(t *testing.T) { + plant, _ := NewGain(mat.NewDense(1, 1, []float64{2}), 0) + controller, _ := NewGain(mat.NewDense(1, 1, []float64{1}), 0) + loop, err := NewGeneralizedClosedLoop("loop", plant, controller, "output") + if err != nil { + t.Fatal(err) + } + sensitivity, err := NewSensitivityGoal("sensitivity", 0.4).Evaluate(loop) + if err != nil { + t.Fatalf("sensitivity Evaluate: %v", err) + } + closedLoop, err := NewWeightedGainGoal("closed_loop", 0.4).Evaluate(loop) + if err != nil { + t.Fatalf("closed-loop Evaluate: %v", err) + } + if !sensitivity.Pass || closedLoop.Pass { + t.Fatalf("sensitivity=%#v closed-loop=%#v, want distinct routed responses", sensitivity, closedLoop) + } +} + +func TestTuningGoalTargetsNamedRectangularLoopBreak(t *testing.T) { + plant, _ := NewGain(mat.NewDense(1, 2, []float64{1, 2}), 0) + controller, _ := NewGain(mat.NewDense(2, 1, []float64{3, 4}), 0) + loop, err := NewGeneralizedClosedLoop("loop", plant, controller, "plant_output") + if err != nil { + t.Fatal(err) + } + if err := loop.InsertAnalysisPoint("plant_input", AnalysisPointPlantInput); err != nil { + t.Fatal(err) + } + outputGoal, err := NewTuningGoal(TuningGoalSpec{ + Name: "output_sensitivity", + Type: TuningGoalSensitivity, + Max: 0.2, + AnalysisPoint: "plant_output", + Omega: []float64{0}, + }) + if err != nil { + t.Fatal(err) + } + inputGoal, err := NewTuningGoal(TuningGoalSpec{ + Name: "input_sensitivity", + Type: TuningGoalSensitivity, + Max: 0.2, + AnalysisPoint: "plant_input", + Omega: []float64{0}, + }) + if err != nil { + t.Fatal(err) + } + outputResult, err := outputGoal.Evaluate(loop) + if err != nil { + t.Fatalf("output Evaluate: %v", err) + } + inputResult, err := inputGoal.Evaluate(loop) + if err != nil { + t.Fatalf("input Evaluate: %v", err) + } + if !outputResult.Pass || inputResult.Pass { + t.Fatalf("output=%#v input=%#v, want distinct point constraints", outputResult, inputResult) + } +} + +func TestTuningGoalRejectsInvalidFrequencyAndWeightDimensions(t *testing.T) { + if _, err := NewTuningGoal(TuningGoalSpec{ + Name: "grid", + Type: TuningGoalWeightedGain, + Max: 1, + Omega: []float64{1, 1}, + }); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("frequency validation error = %v", err) + } + if _, err := NewTuningGoal(TuningGoalSpec{ + Name: "poles", + Type: TuningGoalPole, + Max: -0.1, + Omega: []float64{1}, + }); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("unused frequency validation error = %v", err) + } + badWeight, _ := NewGain(mat.NewDense(2, 2, nil), 0) + goal, err := NewTuningGoal(TuningGoalSpec{ + Name: "weight", + Type: TuningGoalWeightedGain, + Max: 1, + OutputWeight: badWeight, + }) + if err != nil { + t.Fatal(err) + } + if _, err := goal.Evaluate(makeSISO(-1, 1, 1, 0)); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("weight dimension error = %v, want ErrDimensionMismatch", err) + } +} diff --git a/tuning_test.go b/tuning_test.go index b4a3b64..501c1bf 100644 --- a/tuning_test.go +++ b/tuning_test.go @@ -1,6 +1,9 @@ package controlsys -import "testing" +import ( + "errors" + "testing" +) func TestSystuneTunesSISOTunableGain(t *testing.T) { plant := makeSISO(-2, 1, 1, 0) @@ -18,6 +21,9 @@ func TestSystuneTunesSISOTunableGain(t *testing.T) { if !result.Pass { t.Fatalf("expected tuned result to pass, got %#v", result) } + if result.Method != "cartesian-grid" { + t.Fatalf("method = %q, want cartesian-grid", result.Method) + } if result.Parameters["K"] <= 0.1 { t.Fatalf("K was not increased: %#v", result.Parameters) } @@ -26,6 +32,21 @@ func TestSystuneTunesSISOTunableGain(t *testing.T) { } } +func TestGridTuneLimitsCartesianSearch(t *testing.T) { + plant := benchSysNonSym(2, 2, 2) + k1, _ := NewTunableReal("K1", 0.1, TunableBounds{Lower: 0.1, Upper: 2}) + k2, _ := NewTunableReal("K2", 0.1, TunableBounds{Lower: 0.1, Upper: 2}) + controller := NewTunableGain("Kblock", [][]*TunableReal{{k1, fixedReal(t, "z12_limit", 0)}, {fixedReal(t, "z21_limit", 0), k2}}, 0) + closed, err := NewGeneralizedClosedLoop("loop", plant, controller, "u") + if err != nil { + t.Fatal(err) + } + _, err = GridTune(closed, []TuningGoal{NewWeightedGainGoal("bounded", 10)}, &SystuneOptions{GridPoints: 5, MaxEvaluations: 24}) + if !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("GridTune error = %v, want ErrDimensionMismatch", err) + } +} + func TestSystuneTunesSmallMIMOTunableGain(t *testing.T) { plant := benchSysNonSym(2, 2, 2) k1, _ := NewTunableReal("K1", 0.1, TunableBounds{Lower: 0.1, Upper: 2}) From 623be2aaf14a11d6a6cb6ee9f04cb278f6af5c04 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 16:09:43 -0500 Subject: [PATCH 07/10] plan: complete architectural correctness --- .ergo/plans.jsonl | 48 ++++++++++++++++++++++++++ docs/api-mutation-audit.md | 2 +- docs/codebase-interface-diagram.md | 47 ++++++++++++------------- docs/codebase-internal-seam-map.svg | 2 +- docs/codebase-public-interface-map.svg | 2 +- 5 files changed, 74 insertions(+), 27 deletions(-) create mode 100644 .ergo/plans.jsonl diff --git a/.ergo/plans.jsonl b/.ergo/plans.jsonl new file mode 100644 index 0000000..4db7eae --- /dev/null +++ b/.ergo/plans.jsonl @@ -0,0 +1,48 @@ +{"type":"new_task","ts":"2026-07-18T20:09:01.023735Z","data":{"id":"FDNFWB","uuid":"b29e262f-3f32-4561-9dbc-bd7ee8e2b00d","epic_id":"","state":"todo","title":"Architectural and mathematical correctness","body":"","created_at":"2026-07-18T20:09:01.023735Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.02797Z","data":{"id":"2H6QTC","uuid":"02f5b1b3-2fed-42d4-9793-79a169303da0","epic_id":"FDNFWB","state":"todo","title":"Correct ModelArray stacking semantics","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02797Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027977Z","data":{"id":"XIMOJQ","uuid":"d8ca4cc0-a62e-4010-9046-93c452307a8c","epic_id":"FDNFWB","state":"todo","title":"Make passivity results mathematically honest","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027977Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.02798Z","data":{"id":"KFFJUD","uuid":"9c568d6f-cbb7-425b-9b86-f79988d40d7d","epic_id":"FDNFWB","state":"todo","title":"Implement basis-invariant modal truncation","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02798Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027982Z","data":{"id":"V7TYI7","uuid":"a43a226c-82ba-4bb1-bbde-18543108fc56","epic_id":"FDNFWB","state":"todo","title":"Assemble physical port connections into dynamics","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027982Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027985Z","data":{"id":"JA5SHN","uuid":"0417cefc-9a27-4c1d-8b9f-63a952377f8c","epic_id":"FDNFWB","state":"todo","title":"Deepen generalized loops and tuning goals","body":"## Goal\n- Bind analysis points to actual interconnection channels and extract the requested loop.\n- Make MIMO sensitivity and tuning-goal evaluation mathematically correct.\n- Give grid tuning an honest, explicit policy while preserving useful extension seams.\n\n## Acceptance Criteria\n- Different analysis points can identify different loops in one topology.\n- Sensitivity reuses dimension-aware loop-sensitivity behavior.\n- MIMO gain goals use maximum singular value and honor frequency bands/weights.\n- Systune and Looptune either have distinct documented policies or delegate to an explicitly named shared grid tuner.\n- Tests cover multiple points, MIMO loops, custom TunableBlock implementations, and deterministic tuning limits.\n\n## Validation Gates\n- go test -run 'Generalized|AnalysisPoint|Tuning|Systune|Looptune|Loopsens'\n- go test -race -run 'Generalized|Tuning'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027985Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027987Z","data":{"id":"CMLZFF","uuid":"d1dad137-0c41-411b-8c6f-82c3233ac95c","epic_id":"FDNFWB","state":"todo","title":"Integrate architectural correctness changes","body":"## Goal\n- Align package documentation and examples with the corrected mathematical contracts.\n- Run repository-wide static, race, numerical, and benchmark checks before publication.\n\n## Acceptance Criteria\n- README and package docs no longer overstate tracer capabilities.\n- All exported API docs match implemented guarantees and compatibility notes.\n- Full tests, race tests, go fix, go vet, and relevant benchmarks pass.\n- Commits remain grouped by feature and the branch is ready for a draft PR.\n\n## Validation Gates\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n- go test -bench=. -benchmem ./...","created_at":"2026-07-18T20:09:01.027987Z"}} +{"type":"link","ts":"2026-07-18T20:09:11.615311Z","data":{"from_id":"CMLZFF","to_id":"2H6QTC","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.712916Z","data":{"from_id":"CMLZFF","to_id":"XIMOJQ","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.862506Z","data":{"from_id":"CMLZFF","to_id":"KFFJUD","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.955566Z","data":{"from_id":"CMLZFF","to_id":"V7TYI7","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:12.049981Z","data":{"from_id":"CMLZFF","to_id":"JA5SHN","type":"depends"}} +{"type":"claim","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","agent_id":"gpt-5@codex","ts":"2026-07-18T20:09:21.963819Z"}} +{"type":"state","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","state":"doing","ts":"2026-07-18T20:09:21.963819Z"}} +{"type":"result","ts":"2026-07-18T20:11:27.026516Z","data":{"task_id":"2H6QTC","summary":"model_array.go","path":"model_array.go","sha256_at_attach":"e305638757980cb032cfa577899f61559b360137c69a46db0c92acf575e33796","mtime_at_attach":"2026-07-18T20:10:13.854733865Z","git_commit_at_attach":"ed73d670c345ec12d9e845cfaceb2128c111cdef","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"body","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...\n\n## Completion\n- StackModelArrays now creates a leading axis and rejects unequal shapes.\n- ConcatModelArrays preserves the prior flatten-and-append operation under an honest name.\n- Added shape, ordering, void-entry, incompatibility, and concatenation coverage.\n- Validated with gopls diagnostics, focused tests, the full suite, and go vet.\n","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"state","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","state":"done","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"claim","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","agent_id":"gpt-5@codex","ts":"2026-07-18T20:11:27.119203Z"}} +{"type":"state","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","state":"doing","ts":"2026-07-18T20:11:27.119203Z"}} +{"type":"result","ts":"2026-07-18T20:15:56.288913Z","data":{"task_id":"XIMOJQ","summary":"passivity.go","path":"passivity.go","sha256_at_attach":"fdbcd7cf267a25e8857ba5145f40483f809088cf224c03d818d4c7ae6fe07d68","mtime_at_attach":"2026-07-18T20:14:46.772346488Z","git_commit_at_attach":"b105aa0cdeecd6c3d5c76b6df28ebe3ae6435f32","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"body","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...\n\n## Completion\n- Added explicit violated, sampled-pass, and certified result statuses and retained the exact grid and tolerance.\n- Added SampledPassive while preserving Passive as a compatibility alias.\n- Included DC, bounded discrete-time checks at Nyquist, and rejected malformed grids, tolerances, FRD shapes, and non-finite data.\n- Added a persistent restricted-grid counterexample and completed full, race, diagnostic, and vet validation.\n","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"state","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","state":"done","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"claim","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","agent_id":"gpt-5@codex","ts":"2026-07-18T20:15:56.38886Z"}} +{"type":"state","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","state":"doing","ts":"2026-07-18T20:15:56.38886Z"}} +{"type":"result","ts":"2026-07-18T20:24:42.638055Z","data":{"task_id":"KFFJUD","summary":"modal_reduction.go","path":"modal_reduction.go","sha256_at_attach":"a8eb1f89dabe851fd369880a046b852679a1ec3299414e46770d521f22f35a5f","mtime_at_attach":"2026-07-18T20:20:16.446837181Z","git_commit_at_attach":"e029f71a069f98e965f0f9778881ac0e184097ad","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"body","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...\n\n## Completion\n- Replaced source-coordinate deletion with ordered real Schur decomposition and Sylvester block separation.\n- Preserved unstable modes and conjugate pairs and exposed retained poles plus left/right projection bases.\n- Added persistent similarity-invariance, projection, pair-preservation, unstable-mode, and threshold-selection tests.\n- Full tests, focused race tests, diagnostics, and vet pass.\n- Benchstat records 492.8 us/op and 146.8 KiB/op for the corrected N=50/order=10 algorithm versus 5.75 us/op and 25.5 KiB/op for the previous non-modal state slice.\n","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"state","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","state":"done","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"new_task","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","uuid":"8ed87628-27b4-4194-800d-8de87d884765","epic_id":"FDNFWB","state":"todo","title":"Update Gonum fork dependency","body":"","created_at":"2026-07-18T20:24:42.731138Z"}} +{"type":"claim","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","agent_id":"gpt-5@codex","ts":"2026-07-18T20:24:42.731138Z"}} +{"type":"state","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","state":"doing","ts":"2026-07-18T20:24:42.731138Z"}} +{"type":"body","ts":"2026-07-18T20:25:03.048658Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the current v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n","ts":"2026-07-18T20:25:03.048658Z"}} +{"type":"link","ts":"2026-07-18T20:25:03.144556Z","data":{"from_id":"CMLZFF","to_id":"4NPDOO","type":"depends"}} +{"type":"result","ts":"2026-07-18T20:27:45.387877Z","data":{"task_id":"4NPDOO","summary":"go.mod","path":"go.mod","sha256_at_attach":"4757c07b5ab37cbaea17c69fd4dc06762d2ed0e6c9216b65e76f4de2dfee3811","mtime_at_attach":"2026-07-18T20:25:07.260545573Z","git_commit_at_attach":"e7552cf47698163a12eb138009db7dd71596b18e","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"body","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the previous v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n\n## Completion\n- Updated the replace target and checksums to v0.17.7-fork, the latest live fork tag.\n- go mod tidy changed no unrelated modules.\n- Full tests, race tests, go fix, go vet, and diagnostics pass.\n- Vulnerability analysis found only two Go 1.26.4 standard-library advisories and no Gonum/dependency findings.\n- Canon/modal benchmark comparison showed no significant timing or allocation change.\n","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"state","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","state":"done","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"claim","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","agent_id":"gpt-5@codex","ts":"2026-07-18T20:27:45.479569Z"}} +{"type":"state","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","state":"doing","ts":"2026-07-18T20:27:45.479569Z"}} +{"type":"result","ts":"2026-07-18T20:40:29.598727Z","data":{"task_id":"V7TYI7","summary":"physical_assembly.go","path":"physical_assembly.go","sha256_at_attach":"0cc166244f3a0d904d0df1cf625ff075c91473983efc13aa639f20a71ad7a876","mtime_at_attach":"2026-07-18T20:38:52.841498285Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"body","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...\n\n## Completion\n- Added explicit physical-port input/output channel ownership with compatibility mapping for legacy positional ports.\n- Added an assembly-plan owner for component validation, node unioning, external-channel selection, descriptor constraints, and metadata.\n- Connected nodes enforce across-variable equality and through-variable conservation; grounded nodes preserve reaction variables.\n- Added descriptor frequency response through the regular pencil sE-A and direct scalar-oracle coverage.\n- Added pair, grounded, three-component, descriptor-preservation, and invalid-topology behavior tests.\n- Full and focused race tests, diagnostics, and vet pass.\n- Correct assembly benchmark is 24.9 us/op and 92.2 KiB/op versus 14.9 us/op and 55.1 KiB/op for the prior validation-only append.\n","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"state","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","state":"done","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"claim","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","agent_id":"gpt-5@codex","ts":"2026-07-18T20:40:29.696252Z"}} +{"type":"state","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","state":"doing","ts":"2026-07-18T20:40:29.696252Z"}} +{"type":"result","ts":"2026-07-18T20:57:32.859488Z","data":{"task_id":"JA5SHN","summary":"generalized.go","path":"generalized.go","sha256_at_attach":"04348604c2ce16487fcc6c5c920df195c2246f523037c689ead6fe0efd91157d","mtime_at_attach":"2026-07-18T20:56:28.108894409Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:57:32.859488Z"}} +{"type":"state","ts":"2026-07-18T20:57:32.859488Z","data":{"id":"JA5SHN","state":"done","ts":"2026-07-18T20:57:32.859488Z"}} +{"type":"claim","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","agent_id":"codex@gpt-5","ts":"2026-07-18T20:57:52.859874Z"}} +{"type":"state","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","state":"doing","ts":"2026-07-18T20:57:52.859874Z"}} +{"type":"result","ts":"2026-07-18T21:09:30.134851Z","data":{"task_id":"CMLZFF","summary":"docs/codebase-interface-diagram.md","path":"docs/codebase-interface-diagram.md","sha256_at_attach":"926f3544aa5bf37f9e532cb81a2595c257a2d5f0691118000c4cf9348c4fb00b","mtime_at_attach":"2026-07-18T20:59:46.593472404Z","git_commit_at_attach":"9efb609bbe5d07ccc3eab19a74d919ae6caffd4e","ts":"2026-07-18T21:09:30.134851Z"}} +{"type":"state","ts":"2026-07-18T21:09:30.134851Z","data":{"id":"CMLZFF","state":"done","ts":"2026-07-18T21:09:30.134851Z"}} diff --git a/docs/api-mutation-audit.md b/docs/api-mutation-audit.md index defd8dc..9ba006b 100644 --- a/docs/api-mutation-audit.md +++ b/docs/api-mutation-audit.md @@ -396,7 +396,7 @@ ownership clarity, and release readiness, not shrinking the toolbox shape. | Model containers | `ModelArray`, `GeneralizedModel`, `GeneralizedClosedLoop`, `TunableReal`, `TunableGain`, `TunablePID`, `TunableTF`, `TunableSS` | Mostly private fields with mutating methods; `FreeParameters` exposes parameter pointers. | | Options/workspaces | `C2DOptions`, `TransferFuncOpts`, `StateSpaceOpts`, `FreqRespEstOpts`, `StepInfoOptions`, `SimulateOpts`, `RiccatiOpts`, `RiccatiWorkspace`, `LyapunovOpts`, `LyapunovWorkspace`, `PidtuneOptions`, `SystuneOptions`, `PassivityOptions`, `ReduceOpts`, `ModalTruncateOptions` | Options are caller-owned; workspaces and simulation buffers are mutable and should not be shared concurrently. | | Result structs | `BalrealResult`, `CanonResult`, `GramResult`, `H2SynResult`, `HinfSynResult`, `LqgResult`, `LoopsensResult`, `MarginResult`, `AllMarginResult`, `DiskMarginResult`, `ReduceResult`, `ModalReductionResult`, `ModsepResult`, `PrescaleResult`, `PzmapResult`, `TimeResponse`, `StepInfoResult`, `RiccatiResult`, `RootLocusResult`, `Response`, `SsbalResult`, `StabsepResult`, `StaircaseResult`, `StateSpaceResult`, `TransferFuncResult`, `SystuneResult`, `TuningGoalResult`, `ZerosResult`, `ZPKResult`, `ERAResult`, `FRDPeakGainResult`, `FreqRespEstResult`, `ModelArrayFreqResponse`, `ModelArrayTimeResponse` | Results are mutable data containers; callers should treat them as owned outputs unless workspace-backed docs say otherwise. | -| Value and enum types | `BalredMethod`, `CanonForm`, `AbsorbScope`, `GramType`, `PhysicalPortKind`, `PIDForm`, `PidtuneType`, `C2DMethod`, `C2DDelayModeling`, `FreqRespEstMethod`, `ReduceMode`, `TuningGoalType`, `TuningGoalSpec`, `TuningGoal`, `TunableBounds`, `AnalysisPoint`, `PhysicalPort`, `PhysicalConnection`, `Connection`, `DampInfo`, `StepMetric`, `NonlinearModel`, `EKFModel`, `NumericBlock`, `TunableBlock`, `FRDResponseMapper`, `PIDOption`, `SafeFeedbackOption` | Mostly value types; callback and option function types may retain references through user code. | +| Value and enum types | `BalredMethod`, `CanonForm`, `AbsorbScope`, `GramType`, `PhysicalPortKind`, `PIDForm`, `PidtuneType`, `C2DMethod`, `C2DDelayModeling`, `FreqRespEstMethod`, `ReduceMode`, `TuningGoalType`, `TuningGoalSpec`, `TuningGoal`, `TunableBounds`, `AnalysisPointLocation`, `AnalysisPoint`, `PhysicalPort`, `PhysicalConnection`, `Connection`, `DampInfo`, `StepMetric`, `NonlinearModel`, `EKFModel`, `NumericBlock`, `TunableBlock`, `FRDResponseMapper`, `PIDOption`, `SafeFeedbackOption` | Mostly value types; callback and option function types may retain references through user code. | ## Release Gates diff --git a/docs/codebase-interface-diagram.md b/docs/codebase-interface-diagram.md index eeb30bf..3c38957 100644 --- a/docs/codebase-interface-diagram.md +++ b/docs/codebase-interface-diagram.md @@ -1,6 +1,6 @@ # Controlsys Codebase Interface Diagram -This diagram shows the current module interfaces on `main`. It is a codebase-level view, not a complete call graph: the public model interfaces are centered, and the internal seams show where recurring rules are localized. +This diagram shows the current package interfaces. It is a codebase-level view, not a complete call graph: the public model interfaces are centered, and the internal seams show where recurring rules are localized. ## Public Interface Map @@ -28,12 +28,12 @@ flowchart LR realizeTransferFunc["TransferFunc.StateSpace"] constructZPK["NewZPK / NewZPKMIMO"] constructFRD["NewFRD"] - constructModelArray["NewModelArray / StackModelArrays"] + constructModelArray["NewModelArray / StackModelArrays
ConcatModelArrays"] constructGeneralized["NewGeneralizedModel
NewGeneralizedClosedLoop"] identifyERA["ERA
Markov parameters to state-space model"] estimateFreqResponse["FreqRespEst
sampled input/output to response estimate"] linearizeNonlinear["Linearize / NewEKF
local nonlinear approximation"] - assemblePhysical["NewPhysicalComponent / AssemblePhysical
port-checked component assembly"] + assemblePhysical["NewPhysicalComponent / AssemblePhysical
constraint-based descriptor assembly"] end subgraph interconnection["Interconnection interfaces"] @@ -62,7 +62,7 @@ flowchart LR energyAnalysis["Gram / HSV / Norm
H2Norm / HinfNorm / Covar"] structureAnalysis["Ctrb / Obsv / IsStabilizable / IsDetectable"] loopAnalysis["Loopsens / RootLocus / FRDMargin"] - passivityAnalysis["Passive / FRDPassive / SpectralFactor"] + passivityAnalysis["SampledPassive / FRDPassive
Passive compatibility alias / SpectralFactor"] end subgraph transforms["Transformation and reduction"] @@ -78,7 +78,7 @@ flowchart LR observerDesign["Kalman / Kalmd / Lqe / Lqg
Estim / Reg"] robustSynthesis["H2Syn / HinfSyn"] pidDesign["NewPID / NewPIDStd / Pidtune
PID / PID2 / SmithPredictor"] - fixedStructureTuning["Systune / Looptune
tuning goals"] + fixedStructureTuning["GridTune / Systune / Looptune
point-specific tuning goals"] end caller --> constructStateSpace @@ -322,38 +322,37 @@ classDiagram class generalizedTuningSeam { +numericBlockFromAny() +primaryAnalysisPointName() - +tuneFixedStructure() + +GridTune() +evaluateTuningGoals() } class tuningGoalEvaluator { - +tracking() - +maxGain() - +loopShape() - +margin() - +pole() - +overshoot() + +evaluateSystem() + +tuningGoalSystem() + +frequencyGainRange() + +maximumComplexSingularValue() } class passivitySeam { + +SampledPassive() +Passive() +FRDPassive() +SpectralFactor() } - class physicalAssemblySeam { - +NewPhysicalComponent() - +lookupPhysicalPort() - +prefixSystemMetadata() - +AssemblePhysical() + class physicalAssemblyPlan { + +bindPorts() + +bindConnections() + +assemble() + +aggregateMatrices() } class modelArraySeam { +validateModelArrayCompatible() + +validateModelArrayHeadersCompatible() + +StackModelArrays() + +ConcatModelArrays() +flatIndex() - +ModelFlat() - +FreqResponse() - +Step() } class controllerObserverPolicy { @@ -378,7 +377,7 @@ classDiagram System --> timeResponsePlanner : time response System --> simulationDispatcher : sampled simulation System --> passivitySeam : passivity and spectral factor - System --> physicalAssemblySeam : physical component assembly + System --> physicalAssemblyPlan : physical constraint assembly System --> ModelArray : compatible model arrays TransferFunc --> System : realization @@ -392,7 +391,7 @@ classDiagram generalizedPlantPartition --> System : H2/Hinf controller synthesis generalizedTuningSeam --> tuningGoalEvaluator : fixed-structure tuning - tuningGoalEvaluator --> System : evaluates closed-loop model + tuningGoalEvaluator --> System : evaluates point-specific loop responses controllerObserverPolicy --> System : regulator and estimator assembly matrixEquationProblem --> controllerObserverPolicy : Riccati and Lyapunov validation ``` @@ -404,5 +403,5 @@ classDiagram - Interconnection routines concentrate compatibility checks, direct feedthrough handling, delay movement, and metadata propagation behind a small caller-facing interface. - Delay behavior is intentionally split between topology and conversion seams: topology answers what delay structure exists; conversion decides whether it remains explicit, becomes a delay bank, or moves into LFT form. - Analysis routines share sampled-response layouts so frequency-response data, Bode results, singular-value analysis, and frequency-response estimates use the same output/input/frequency indexing. -- Model-array, physical-assembly, and state-space utility seams support model-grid, port-checked assembly, and signal-selection workflows while keeping compatibility checks and metadata rules localized. -- Synthesis routines route generalized-plant, generalized tuning, and controller/observer rules through policy modules before returning controller or closed-loop state-space models. +- Model-array, physical-assembly, and state-space utility seams support model-grid semantics, across/through descriptor constraints, and signal-selection workflows while keeping compatibility checks and metadata rules localized. +- Synthesis routines route generalized-plant, point-specific tuning, and controller/observer rules through policy modules before returning controller or closed-loop state-space models. diff --git a/docs/codebase-internal-seam-map.svg b/docs/codebase-internal-seam-map.svg index 7988c5e..5a07a6a 100644 --- a/docs/codebase-internal-seam-map.svg +++ b/docs/codebase-internal-seam-map.svg @@ -1 +1 @@ -

descriptor gate

time-domain rules

external delay decomposition

delay conversion

interconnection planning

realization assembly

state-space utilities

frequency response

time response

sampled simulation

passivity and spectral factor

physical component assembly

compatible model arrays

realization

rational-channel conversion

sampled response access

flat response layout

array validation and analysis

analysis-point model wrapper

loop extraction

sampled controller blocks

H2/Hinf controller synthesis

fixed-structure tuning

evaluates closed-loop model

regulator and estimator assembly

Riccati and Lyapunov validation

System

+Dims()

+Validate()

+Poles()

+IsStable()

+DCGain()

+IsContinuous()

+IsDiscrete()

+Simulate()

+FreqResponse()

+TransferFunction()

TransferFunc

+Dims()

+Eval()

+EvalMulti()

+StateSpace()

+ZPK()

ZPK

+Dims()

+Eval()

+FreqResponse()

+TransferFunction()

+StateSpace()

FRD

+Dims()

+NumFrequencies()

+At()

+Abs()

+SelectFrequencies()

+MapResponse()

+PeakGain()

+FreqResponse()

+Bode()

ModelArray

+Shape()

+Model()

+SelectFlat()

+ModelFlat()

+FreqResponse()

+Step()

GeneralizedModel

+InsertAnalysisPoint()

+AnalysisPoint()

+HasAnalysisPoint()

+CurrentSystem()

GeneralizedClosedLoop

+OpenLoop()

+ClosedLoop()

+Sensitivity()

+ComplementarySensitivity()

TunableBlock

+CurrentSystem()

+FreeParameters()

+SampleBlock()

descriptorPolicy

+validate()

+poles()

+requireStandard()

+requireRiccatiStandard()

timeDomain

+validateSampleTime()

+frequencyVariable()

+ensureCompatible()

delayTopology

+totalExternal()

+decomposedExternal()

+decomposableExternal()

delayConversionPolicy

+applyDiscreteDelayFields()

+applyContinuousDelayFields()

+replaceDiscreteExternal()

+replaceContinuousExternal()

interconnectionTopology

+seriesDelayPlan()

+parallelDelayPlan()

+seriesRequiresLFT()

+parallelRequiresLFT()

realizationTransformPolicy

+requireStandard()

+requireDelayFree()

+result()

+resultWithOriginalFeedthrough()

+resultWithZeroFeedthrough()

stateSpaceUtilitySeam

+NewDescriptor()

+ToExplicit()

+EliminateStates()

+FixedInputReduction()

+SelectByName()

+SelectByIndex()

+AugmentInternalDelayOutputs()

frequencyEvaluator

+response()

+eval()

timeResponsePlanner

+auto()

+lsim()

simulationDispatcher

+run()

sampledResponseLayout

+offset()

+blockOffset()

generalizedPlantPartition

+validateControllerChannels()

+newController()

+closedLoopPoles()

generalizedTuningSeam

+numericBlockFromAny()

+primaryAnalysisPointName()

+tuneFixedStructure()

+evaluateTuningGoals()

tuningGoalEvaluator

+tracking()

+maxGain()

+loopShape()

+margin()

+pole()

+overshoot()

passivitySeam

+Passive()

+FRDPassive()

+SpectralFactor()

physicalAssemblySeam

+NewPhysicalComponent()

+lookupPhysicalPort()

+prefixSystemMetadata()

+AssemblePhysical()

modelArraySeam

+validateModelArrayCompatible()

+flatIndex()

+ModelFlat()

+FreqResponse()

+Step()

controllerObserverPolicy

+validateNoise()

+regulator()

+estimator()

matrixEquationProblem

+riccatiProblem

+lyapunovProblem

 \ No newline at end of file +

descriptor gate

time-domain rules

external delay decomposition

delay conversion

interconnection planning

realization assembly

state-space utilities

frequency response

time response

sampled simulation

passivity and spectral factor

physical constraint assembly

compatible model arrays

realization

rational-channel conversion

sampled response access

flat response layout

array validation and analysis

analysis-point model wrapper

loop extraction

sampled controller blocks

H2/Hinf controller synthesis

fixed-structure tuning

evaluates point-specific loop responses

regulator and estimator assembly

Riccati and Lyapunov validation

System

+Dims()

+Validate()

+Poles()

+IsStable()

+DCGain()

+IsContinuous()

+IsDiscrete()

+Simulate()

+FreqResponse()

+TransferFunction()

TransferFunc

+Dims()

+Eval()

+EvalMulti()

+StateSpace()

+ZPK()

ZPK

+Dims()

+Eval()

+FreqResponse()

+TransferFunction()

+StateSpace()

FRD

+Dims()

+NumFrequencies()

+At()

+Abs()

+SelectFrequencies()

+MapResponse()

+PeakGain()

+FreqResponse()

+Bode()

ModelArray

+Shape()

+Model()

+SelectFlat()

+ModelFlat()

+FreqResponse()

+Step()

GeneralizedModel

+InsertAnalysisPoint()

+AnalysisPoint()

+HasAnalysisPoint()

+CurrentSystem()

GeneralizedClosedLoop

+InsertAnalysisPoint()

+AnalysisPoint()

+OpenLoop()

+ClosedLoop()

+Sensitivity()

+ComplementarySensitivity()

TunableBlock

+CurrentSystem()

+FreeParameters()

+SampleBlock()

descriptorPolicy

+validate()

+poles()

+requireStandard()

+requireRiccatiStandard()

timeDomain

+validateSampleTime()

+frequencyVariable()

+ensureCompatible()

delayTopology

+totalExternal()

+decomposedExternal()

+decomposableExternal()

delayConversionPolicy

+applyDiscreteDelayFields()

+applyContinuousDelayFields()

+replaceDiscreteExternal()

+replaceContinuousExternal()

interconnectionTopology

+seriesDelayPlan()

+parallelDelayPlan()

+seriesRequiresLFT()

+parallelRequiresLFT()

realizationTransformPolicy

+requireStandard()

+requireDelayFree()

+result()

+resultWithOriginalFeedthrough()

+resultWithZeroFeedthrough()

stateSpaceUtilitySeam

+NewDescriptor()

+ToExplicit()

+EliminateStates()

+FixedInputReduction()

+SelectByName()

+SelectByIndex()

+AugmentInternalDelayOutputs()

frequencyEvaluator

+response()

+eval()

timeResponsePlanner

+auto()

+lsim()

simulationDispatcher

+run()

sampledResponseLayout

+offset()

+blockOffset()

generalizedPlantPartition

+validateControllerChannels()

+newController()

+closedLoopPoles()

generalizedTuningSeam

+numericBlockFromAny()

+primaryAnalysisPointName()

+GridTune()

+evaluateTuningGoals()

tuningGoalEvaluator

+evaluateSystem()

+tuningGoalSystem()

+frequencyGainRange()

+maximumComplexSingularValue()

passivitySeam

+SampledPassive()

+Passive()

+FRDPassive()

+SpectralFactor()

physicalAssemblyPlan

+bindPorts()

+bindConnections()

+assemble()

+aggregateMatrices()

modelArraySeam

+validateModelArrayCompatible()

+validateModelArrayHeadersCompatible()

+StackModelArrays()

+ConcatModelArrays()

+flatIndex()

controllerObserverPolicy

+validateNoise()

+regulator()

+estimator()

matrixEquationProblem

+riccatiProblem

+lyapunovProblem

 \ No newline at end of file diff --git a/docs/codebase-public-interface-map.svg b/docs/codebase-public-interface-map.svg index 6c2427e..8425b4b 100644 --- a/docs/codebase-public-interface-map.svg +++ b/docs/codebase-public-interface-map.svg @@ -1 +1 @@ -

Design and synthesis

Transformation and reduction

Analysis interfaces

Representation and domain conversion

Interconnection interfaces

Construction and identification

Model interfaces

External callers

System
fundamental state-space model
A, B, C, D, optional E, delays, names

TransferFunc
polynomial-ratio model

ZPK
zero-pole-gain model

FRD
frequency-response data model

ModelArray
compatible model grid

GeneralizedModel / GeneralizedClosedLoop
analysis-point model interface

TunableBlock
tunable gain, PID, TF, or SS block

FreqResponseMatrix
sampled complex response

TimeResponse
sampled time-domain output

New / NewGain / NewFromSlices
NewWithDelay / Rss / Drss

NewDescriptor / ToExplicit

TransferFunc.StateSpace

NewZPK / NewZPKMIMO

NewFRD

NewModelArray / StackModelArrays

NewGeneralizedModel
NewGeneralizedClosedLoop

ERA
Markov parameters to state-space model

FreqRespEst
sampled input/output to response estimate

Linearize / NewEKF
local nonlinear approximation

NewPhysicalComponent / AssemblePhysical
port-checked component assembly

Series

Parallel

Feedback / SafeFeedback

Append / BlkDiag / Connect
ConnectByName / LFT / SumBlk

FRDSeries / FRDParallel
FRDFeedback / FRDConcat

PadeDelay / ThiranDelay
PullDelaysToLFT / AbsorbDelay
SetDelayModel / GetDelayModel

System.TransferFunction

System.ZPKModel

System.FRD

System.Discretize / DiscretizeWithOpts
DiscretizeZOH / FOH / Impulse / Matched / D2D

System.Undiscretize / System.D2C

StateTransform / EliminateStates
FixedInputReduction / SelectByName / SelectByIndex

Step / Impulse / Initial / Lsim
Simulate / StepInfo / StepInfoForSystem

GenSig
test-signal generation

FreqResponse / Bode / Nyquist / Nichols
Margin / AllMargin / DiskMargin / Bandwidth / Sigma

Poles / Zeros / Damp / IsStable
DCGain / Pzmap

Gram / HSV / Norm
H2Norm / HinfNorm / Covar

Ctrb / Obsv / IsStabilizable / IsDetectable

Loopsens / RootLocus / FRDMargin

Passive / FRDPassive / SpectralFactor

SS2SS / Xperm / Canon

Balreal / Balred / Modred / Sminreal
Reduce / MinimalRealization / ModalTruncate

Stabsep / Modsep / Prescale / Ssbal

Inv / Augstate

Care / Dare / Lyap / DLyap

Lqr / Dlqr / Lqi / Lqrd
Place / Acker

Kalman / Kalmd / Lqe / Lqg
Estim / Reg

H2Syn / HinfSyn

NewPID / NewPIDStd / Pidtune
PID / PID2 / SmithPredictor

Systune / Looptune
tuning goals

 \ No newline at end of file +

Design and synthesis

Transformation and reduction

Analysis interfaces

Representation and domain conversion

Interconnection interfaces

Construction and identification

Model interfaces

External callers

System
fundamental state-space model
A, B, C, D, optional E, delays, names

TransferFunc
polynomial-ratio model

ZPK
zero-pole-gain model

FRD
frequency-response data model

ModelArray
compatible model grid

GeneralizedModel / GeneralizedClosedLoop
analysis-point model interface

TunableBlock
tunable gain, PID, TF, or SS block

FreqResponseMatrix
sampled complex response

TimeResponse
sampled time-domain output

New / NewGain / NewFromSlices
NewWithDelay / Rss / Drss

NewDescriptor / ToExplicit

TransferFunc.StateSpace

NewZPK / NewZPKMIMO

NewFRD

NewModelArray / StackModelArrays
ConcatModelArrays

NewGeneralizedModel
NewGeneralizedClosedLoop

ERA
Markov parameters to state-space model

FreqRespEst
sampled input/output to response estimate

Linearize / NewEKF
local nonlinear approximation

NewPhysicalComponent / AssemblePhysical
constraint-based descriptor assembly

Series

Parallel

Feedback / SafeFeedback

Append / BlkDiag / Connect
ConnectByName / LFT / SumBlk

FRDSeries / FRDParallel
FRDFeedback / FRDConcat

PadeDelay / ThiranDelay
PullDelaysToLFT / AbsorbDelay
SetDelayModel / GetDelayModel

System.TransferFunction

System.ZPKModel

System.FRD

System.Discretize / DiscretizeWithOpts
DiscretizeZOH / FOH / Impulse / Matched / D2D

System.Undiscretize / System.D2C

StateTransform / EliminateStates
FixedInputReduction / SelectByName / SelectByIndex

Step / Impulse / Initial / Lsim
Simulate / StepInfo / StepInfoForSystem

GenSig
test-signal generation

FreqResponse / Bode / Nyquist / Nichols
Margin / AllMargin / DiskMargin / Bandwidth / Sigma

Poles / Zeros / Damp / IsStable
DCGain / Pzmap

Gram / HSV / Norm
H2Norm / HinfNorm / Covar

Ctrb / Obsv / IsStabilizable / IsDetectable

Loopsens / RootLocus / FRDMargin

SampledPassive / FRDPassive
Passive compatibility alias / SpectralFactor

SS2SS / Xperm / Canon

Balreal / Balred / Modred / Sminreal
Reduce / MinimalRealization / ModalTruncate

Stabsep / Modsep / Prescale / Ssbal

Inv / Augstate

Care / Dare / Lyap / DLyap

Lqr / Dlqr / Lqi / Lqrd
Place / Acker

Kalman / Kalmd / Lqe / Lqg
Estim / Reg

H2Syn / HinfSyn

NewPID / NewPIDStd / Pidtune
PID / PID2 / SmithPredictor

GridTune / Systune / Looptune
point-specific tuning goals

 \ No newline at end of file From 9040de6f7b968c494eab3f1895d0f3e18af1af62 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 16:19:37 -0500 Subject: [PATCH 08/10] Fix passivity grid endpoints --- passivity.go | 2 ++ 1 file changed, 2 insertions(+) diff --git a/passivity.go b/passivity.go index 4e11c0b..2ff121a 100644 --- a/passivity.go +++ b/passivity.go @@ -180,6 +180,8 @@ func passivityGrid(opts *PassivityOptions, dt float64) ([]float64, float64, erro lower := math.Min(1e-2, upper*1e-4) omega := make([]float64, 121) copy(omega[1:], logspace(math.Log10(lower), math.Log10(upper), 120)) + omega[1] = lower + omega[len(omega)-1] = upper return omega, tol, nil } From 14029acd0f1c628308bc82c626e8a536079f91a4 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 16:43:36 -0500 Subject: [PATCH 09/10] Address review findings on architectural-correctness PR - untrack .ergo/plans.jsonl and ignore .ergo/ - explain order rejection when it would split a conjugate Schur pair - reject MaxRealPart thresholds that select no modes instead of falling back to the half-order heuristic - give GridTune candidates their own analysis-point map via withSampledController instead of a shallow model copy - reuse weight-response and product buffers across the frequency loop in frequencyGainRange (-30% allocs on the dynamic-weight path) - require strictly increasing passivity grids, matching tuning-goal grid validation, and document the FRDPassive grid contract Co-Authored-By: Claude Fable 5 --- .ergo/plans.jsonl | 48 -------------------------------------- .gitignore | 1 + README.md | 2 +- benchmark_hotspots_test.go | 20 ++++++++++++++++ generalized.go | 19 ++++++++++++++- modal_reduction.go | 21 +++++++++++------ modal_reduction_test.go | 3 +++ passivity.go | 4 ++-- passivity_test.go | 3 +++ tuning.go | 6 ++--- tuning_goal.go | 25 ++++++++++++++++---- 11 files changed, 85 insertions(+), 67 deletions(-) delete mode 100644 .ergo/plans.jsonl diff --git a/.ergo/plans.jsonl b/.ergo/plans.jsonl deleted file mode 100644 index 4db7eae..0000000 --- a/.ergo/plans.jsonl +++ /dev/null @@ -1,48 +0,0 @@ -{"type":"new_task","ts":"2026-07-18T20:09:01.023735Z","data":{"id":"FDNFWB","uuid":"b29e262f-3f32-4561-9dbc-bd7ee8e2b00d","epic_id":"","state":"todo","title":"Architectural and mathematical correctness","body":"","created_at":"2026-07-18T20:09:01.023735Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.02797Z","data":{"id":"2H6QTC","uuid":"02f5b1b3-2fed-42d4-9793-79a169303da0","epic_id":"FDNFWB","state":"todo","title":"Correct ModelArray stacking semantics","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02797Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.027977Z","data":{"id":"XIMOJQ","uuid":"d8ca4cc0-a62e-4010-9046-93c452307a8c","epic_id":"FDNFWB","state":"todo","title":"Make passivity results mathematically honest","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027977Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.02798Z","data":{"id":"KFFJUD","uuid":"9c568d6f-cbb7-425b-9b86-f79988d40d7d","epic_id":"FDNFWB","state":"todo","title":"Implement basis-invariant modal truncation","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02798Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.027982Z","data":{"id":"V7TYI7","uuid":"a43a226c-82ba-4bb1-bbde-18543108fc56","epic_id":"FDNFWB","state":"todo","title":"Assemble physical port connections into dynamics","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027982Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.027985Z","data":{"id":"JA5SHN","uuid":"0417cefc-9a27-4c1d-8b9f-63a952377f8c","epic_id":"FDNFWB","state":"todo","title":"Deepen generalized loops and tuning goals","body":"## Goal\n- Bind analysis points to actual interconnection channels and extract the requested loop.\n- Make MIMO sensitivity and tuning-goal evaluation mathematically correct.\n- Give grid tuning an honest, explicit policy while preserving useful extension seams.\n\n## Acceptance Criteria\n- Different analysis points can identify different loops in one topology.\n- Sensitivity reuses dimension-aware loop-sensitivity behavior.\n- MIMO gain goals use maximum singular value and honor frequency bands/weights.\n- Systune and Looptune either have distinct documented policies or delegate to an explicitly named shared grid tuner.\n- Tests cover multiple points, MIMO loops, custom TunableBlock implementations, and deterministic tuning limits.\n\n## Validation Gates\n- go test -run 'Generalized|AnalysisPoint|Tuning|Systune|Looptune|Loopsens'\n- go test -race -run 'Generalized|Tuning'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027985Z"}} -{"type":"new_task","ts":"2026-07-18T20:09:01.027987Z","data":{"id":"CMLZFF","uuid":"d1dad137-0c41-411b-8c6f-82c3233ac95c","epic_id":"FDNFWB","state":"todo","title":"Integrate architectural correctness changes","body":"## Goal\n- Align package documentation and examples with the corrected mathematical contracts.\n- Run repository-wide static, race, numerical, and benchmark checks before publication.\n\n## Acceptance Criteria\n- README and package docs no longer overstate tracer capabilities.\n- All exported API docs match implemented guarantees and compatibility notes.\n- Full tests, race tests, go fix, go vet, and relevant benchmarks pass.\n- Commits remain grouped by feature and the branch is ready for a draft PR.\n\n## Validation Gates\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n- go test -bench=. -benchmem ./...","created_at":"2026-07-18T20:09:01.027987Z"}} -{"type":"link","ts":"2026-07-18T20:09:11.615311Z","data":{"from_id":"CMLZFF","to_id":"2H6QTC","type":"depends"}} -{"type":"link","ts":"2026-07-18T20:09:11.712916Z","data":{"from_id":"CMLZFF","to_id":"XIMOJQ","type":"depends"}} -{"type":"link","ts":"2026-07-18T20:09:11.862506Z","data":{"from_id":"CMLZFF","to_id":"KFFJUD","type":"depends"}} -{"type":"link","ts":"2026-07-18T20:09:11.955566Z","data":{"from_id":"CMLZFF","to_id":"V7TYI7","type":"depends"}} -{"type":"link","ts":"2026-07-18T20:09:12.049981Z","data":{"from_id":"CMLZFF","to_id":"JA5SHN","type":"depends"}} -{"type":"claim","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","agent_id":"gpt-5@codex","ts":"2026-07-18T20:09:21.963819Z"}} -{"type":"state","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","state":"doing","ts":"2026-07-18T20:09:21.963819Z"}} -{"type":"result","ts":"2026-07-18T20:11:27.026516Z","data":{"task_id":"2H6QTC","summary":"model_array.go","path":"model_array.go","sha256_at_attach":"e305638757980cb032cfa577899f61559b360137c69a46db0c92acf575e33796","mtime_at_attach":"2026-07-18T20:10:13.854733865Z","git_commit_at_attach":"ed73d670c345ec12d9e845cfaceb2128c111cdef","ts":"2026-07-18T20:11:27.026516Z"}} -{"type":"body","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...\n\n## Completion\n- StackModelArrays now creates a leading axis and rejects unequal shapes.\n- ConcatModelArrays preserves the prior flatten-and-append operation under an honest name.\n- Added shape, ordering, void-entry, incompatibility, and concatenation coverage.\n- Validated with gopls diagnostics, focused tests, the full suite, and go vet.\n","ts":"2026-07-18T20:11:27.026516Z"}} -{"type":"state","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","state":"done","ts":"2026-07-18T20:11:27.026516Z"}} -{"type":"claim","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","agent_id":"gpt-5@codex","ts":"2026-07-18T20:11:27.119203Z"}} -{"type":"state","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","state":"doing","ts":"2026-07-18T20:11:27.119203Z"}} -{"type":"result","ts":"2026-07-18T20:15:56.288913Z","data":{"task_id":"XIMOJQ","summary":"passivity.go","path":"passivity.go","sha256_at_attach":"fdbcd7cf267a25e8857ba5145f40483f809088cf224c03d818d4c7ae6fe07d68","mtime_at_attach":"2026-07-18T20:14:46.772346488Z","git_commit_at_attach":"b105aa0cdeecd6c3d5c76b6df28ebe3ae6435f32","ts":"2026-07-18T20:15:56.288913Z"}} -{"type":"body","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...\n\n## Completion\n- Added explicit violated, sampled-pass, and certified result statuses and retained the exact grid and tolerance.\n- Added SampledPassive while preserving Passive as a compatibility alias.\n- Included DC, bounded discrete-time checks at Nyquist, and rejected malformed grids, tolerances, FRD shapes, and non-finite data.\n- Added a persistent restricted-grid counterexample and completed full, race, diagnostic, and vet validation.\n","ts":"2026-07-18T20:15:56.288913Z"}} -{"type":"state","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","state":"done","ts":"2026-07-18T20:15:56.288913Z"}} -{"type":"claim","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","agent_id":"gpt-5@codex","ts":"2026-07-18T20:15:56.38886Z"}} -{"type":"state","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","state":"doing","ts":"2026-07-18T20:15:56.38886Z"}} -{"type":"result","ts":"2026-07-18T20:24:42.638055Z","data":{"task_id":"KFFJUD","summary":"modal_reduction.go","path":"modal_reduction.go","sha256_at_attach":"a8eb1f89dabe851fd369880a046b852679a1ec3299414e46770d521f22f35a5f","mtime_at_attach":"2026-07-18T20:20:16.446837181Z","git_commit_at_attach":"e029f71a069f98e965f0f9778881ac0e184097ad","ts":"2026-07-18T20:24:42.638055Z"}} -{"type":"body","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...\n\n## Completion\n- Replaced source-coordinate deletion with ordered real Schur decomposition and Sylvester block separation.\n- Preserved unstable modes and conjugate pairs and exposed retained poles plus left/right projection bases.\n- Added persistent similarity-invariance, projection, pair-preservation, unstable-mode, and threshold-selection tests.\n- Full tests, focused race tests, diagnostics, and vet pass.\n- Benchstat records 492.8 us/op and 146.8 KiB/op for the corrected N=50/order=10 algorithm versus 5.75 us/op and 25.5 KiB/op for the previous non-modal state slice.\n","ts":"2026-07-18T20:24:42.638055Z"}} -{"type":"state","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","state":"done","ts":"2026-07-18T20:24:42.638055Z"}} -{"type":"new_task","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","uuid":"8ed87628-27b4-4194-800d-8de87d884765","epic_id":"FDNFWB","state":"todo","title":"Update Gonum fork dependency","body":"","created_at":"2026-07-18T20:24:42.731138Z"}} -{"type":"claim","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","agent_id":"gpt-5@codex","ts":"2026-07-18T20:24:42.731138Z"}} -{"type":"state","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","state":"doing","ts":"2026-07-18T20:24:42.731138Z"}} -{"type":"body","ts":"2026-07-18T20:25:03.048658Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the current v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n","ts":"2026-07-18T20:25:03.048658Z"}} -{"type":"link","ts":"2026-07-18T20:25:03.144556Z","data":{"from_id":"CMLZFF","to_id":"4NPDOO","type":"depends"}} -{"type":"result","ts":"2026-07-18T20:27:45.387877Z","data":{"task_id":"4NPDOO","summary":"go.mod","path":"go.mod","sha256_at_attach":"4757c07b5ab37cbaea17c69fd4dc06762d2ed0e6c9216b65e76f4de2dfee3811","mtime_at_attach":"2026-07-18T20:25:07.260545573Z","git_commit_at_attach":"e7552cf47698163a12eb138009db7dd71596b18e","ts":"2026-07-18T20:27:45.387877Z"}} -{"type":"body","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the previous v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n\n## Completion\n- Updated the replace target and checksums to v0.17.7-fork, the latest live fork tag.\n- go mod tidy changed no unrelated modules.\n- Full tests, race tests, go fix, go vet, and diagnostics pass.\n- Vulnerability analysis found only two Go 1.26.4 standard-library advisories and no Gonum/dependency findings.\n- Canon/modal benchmark comparison showed no significant timing or allocation change.\n","ts":"2026-07-18T20:27:45.387877Z"}} -{"type":"state","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","state":"done","ts":"2026-07-18T20:27:45.387877Z"}} -{"type":"claim","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","agent_id":"gpt-5@codex","ts":"2026-07-18T20:27:45.479569Z"}} -{"type":"state","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","state":"doing","ts":"2026-07-18T20:27:45.479569Z"}} -{"type":"result","ts":"2026-07-18T20:40:29.598727Z","data":{"task_id":"V7TYI7","summary":"physical_assembly.go","path":"physical_assembly.go","sha256_at_attach":"0cc166244f3a0d904d0df1cf625ff075c91473983efc13aa639f20a71ad7a876","mtime_at_attach":"2026-07-18T20:38:52.841498285Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:40:29.598727Z"}} -{"type":"body","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...\n\n## Completion\n- Added explicit physical-port input/output channel ownership with compatibility mapping for legacy positional ports.\n- Added an assembly-plan owner for component validation, node unioning, external-channel selection, descriptor constraints, and metadata.\n- Connected nodes enforce across-variable equality and through-variable conservation; grounded nodes preserve reaction variables.\n- Added descriptor frequency response through the regular pencil sE-A and direct scalar-oracle coverage.\n- Added pair, grounded, three-component, descriptor-preservation, and invalid-topology behavior tests.\n- Full and focused race tests, diagnostics, and vet pass.\n- Correct assembly benchmark is 24.9 us/op and 92.2 KiB/op versus 14.9 us/op and 55.1 KiB/op for the prior validation-only append.\n","ts":"2026-07-18T20:40:29.598727Z"}} -{"type":"state","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","state":"done","ts":"2026-07-18T20:40:29.598727Z"}} -{"type":"claim","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","agent_id":"gpt-5@codex","ts":"2026-07-18T20:40:29.696252Z"}} -{"type":"state","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","state":"doing","ts":"2026-07-18T20:40:29.696252Z"}} -{"type":"result","ts":"2026-07-18T20:57:32.859488Z","data":{"task_id":"JA5SHN","summary":"generalized.go","path":"generalized.go","sha256_at_attach":"04348604c2ce16487fcc6c5c920df195c2246f523037c689ead6fe0efd91157d","mtime_at_attach":"2026-07-18T20:56:28.108894409Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:57:32.859488Z"}} -{"type":"state","ts":"2026-07-18T20:57:32.859488Z","data":{"id":"JA5SHN","state":"done","ts":"2026-07-18T20:57:32.859488Z"}} -{"type":"claim","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","agent_id":"codex@gpt-5","ts":"2026-07-18T20:57:52.859874Z"}} -{"type":"state","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","state":"doing","ts":"2026-07-18T20:57:52.859874Z"}} -{"type":"result","ts":"2026-07-18T21:09:30.134851Z","data":{"task_id":"CMLZFF","summary":"docs/codebase-interface-diagram.md","path":"docs/codebase-interface-diagram.md","sha256_at_attach":"926f3544aa5bf37f9e532cb81a2595c257a2d5f0691118000c4cf9348c4fb00b","mtime_at_attach":"2026-07-18T20:59:46.593472404Z","git_commit_at_attach":"9efb609bbe5d07ccc3eab19a74d919ae6caffd4e","ts":"2026-07-18T21:09:30.134851Z"}} -{"type":"state","ts":"2026-07-18T21:09:30.134851Z","data":{"id":"CMLZFF","state":"done","ts":"2026-07-18T21:09:30.134851Z"}} diff --git a/.gitignore b/.gitignore index d8c31bf..5479639 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ .DS_Store pprof_list.txt referance/ +.ergo/ diff --git a/README.md b/README.md index 21b7e65..ea7db08 100644 --- a/README.md +++ b/README.md @@ -180,7 +180,7 @@ func main() { | `Bandwidth` | -3 dB bandwidth | | `RootLocus` | Root locus as a function of loop gain | | `Pzmap` | Poles and transmission zeros for plotting/inspection | -| `SampledPassive` / `FRDPassive` | Passivity evidence on an explicit frequency grid; passing samples are not an analytic certificate | +| `SampledPassive` / `FRDPassive` | Passivity evidence on an explicit frequency grid; passing samples are not an analytic certificate. The grid must be strictly increasing, finite, non-negative, and at or below the Nyquist frequency for sampled models | | `Passive` | Compatibility alias for `SampledPassive` | | `SpectralFactor` | Spectral factor for supported static-gain models | diff --git a/benchmark_hotspots_test.go b/benchmark_hotspots_test.go index 00d0c4b..4e4fb2d 100644 --- a/benchmark_hotspots_test.go +++ b/benchmark_hotspots_test.go @@ -546,6 +546,26 @@ func BenchmarkTuningGoalWeightedGain_MIMO(b *testing.B) { } } +func BenchmarkTuningGoalDynamicWeightedGain_MIMO(b *testing.B) { + sys := benchSysNonSym(8, 3, 3) + goal, err := NewTuningGoal(TuningGoalSpec{ + Name: "weighted", + Type: TuningGoalWeightedGain, + Max: 10, + InputWeight: benchSysNonSym(2, 3, 3), + OutputWeight: benchSysNonSym(2, 3, 3), + }) + if err != nil { + b.Fatal(err) + } + b.ResetTimer() + for i := 0; i < b.N; i++ { + if _, err := goal.Evaluate(sys); err != nil { + b.Fatal(err) + } + } +} + func BenchmarkSystune_SISO(b *testing.B) { plant := benchSysNonSym(2, 1, 1) k, _ := NewTunableReal("K", 0.5, TunableBounds{Lower: 0.1, Upper: 3}) diff --git a/generalized.go b/generalized.go index 3e594d5..7669251 100644 --- a/generalized.go +++ b/generalized.go @@ -1,6 +1,9 @@ package controlsys -import "fmt" +import ( + "fmt" + "maps" +) type NumericBlock interface { CurrentSystem() (*System, error) @@ -171,6 +174,20 @@ func (g *GeneralizedClosedLoop) InsertAnalysisPoint(name string, location Analys return nil } +// withSampledController returns a model bound to sampled without sharing any +// mutable state with the receiver; the plant is shared because it is never +// mutated after construction. +func (g *GeneralizedClosedLoop) withSampledController(sampled TunableBlock) *GeneralizedClosedLoop { + return &GeneralizedClosedLoop{ + name: g.name, + plant: g.plant, + controller: sampled, + tunableController: sampled, + analysisPoints: maps.Clone(g.analysisPoints), + primaryAnalysisPoint: g.primaryAnalysisPoint, + } +} + func (g *GeneralizedClosedLoop) AnalysisPoint(name string) (AnalysisPoint, error) { if g == nil { return AnalysisPoint{}, fmt.Errorf("GeneralizedClosedLoop.AnalysisPoint: nil model: %w", ErrDimensionMismatch) diff --git a/modal_reduction.go b/modal_reduction.go index 95b9e7f..5b3312a 100644 --- a/modal_reduction.go +++ b/modal_reduction.go @@ -54,10 +54,16 @@ func ModalTruncate(sys *System, opts *ModalTruncateOptions) (*ModalReductionResu return nil, err } poles := schurEigenvaluesRaw(t, n) - order := modalReductionOrder(t, poles, n, sys.Dt, opts) - if order < 1 || order > n || splitsSchurBlock(t, n, order) { + order, err := modalReductionOrder(t, poles, n, sys.Dt, opts) + if err != nil { + return nil, err + } + if order < 1 || order > n { return nil, ErrInvalidOrder } + if splitsSchurBlock(t, n, order) { + return nil, fmt.Errorf("ModalTruncate: order %d splits a complex-conjugate mode pair: %w", order, ErrInvalidOrder) + } if order < unstableModalCount(poles, sys.Dt) { return nil, fmt.Errorf("ModalTruncate: order %d would discard an unstable mode: %w", order, ErrInvalidOrder) } @@ -172,18 +178,19 @@ func modalBlockBefore(candidate, current complex128, dt, maxRealPart float64, th return math.Abs(imag(candidate)) < math.Abs(imag(current)) } -func modalReductionOrder(t []float64, poles []complex128, n int, dt float64, opts *ModalTruncateOptions) int { +func modalReductionOrder(t []float64, poles []complex128, n int, dt float64, opts *ModalTruncateOptions) (int, error) { if opts.Order > 0 { - return opts.Order + return opts.Order, nil } if opts.MaxRealPart != 0 { order := 0 for order < n && (real(poles[order]) >= opts.MaxRealPart || modalPoleUnstable(poles[order], dt)) { order += schurBlockSize(t, n, order) } - if order > 0 { - return order + if order == 0 { + return 0, fmt.Errorf("ModalTruncate: no modes satisfy MaxRealPart %g: %w", opts.MaxRealPart, ErrInvalidOrder) } + return order, nil } order := n / 2 if order == 0 { @@ -196,7 +203,7 @@ func modalReductionOrder(t []float64, poles []complex128, n int, dt float64, opt if order < unstable { order = unstable } - return order + return order, nil } func separateSchurBlocks(t []float64, n, order int) ([]float64, error) { diff --git a/modal_reduction_test.go b/modal_reduction_test.go index 3c88162..f04fec9 100644 --- a/modal_reduction_test.go +++ b/modal_reduction_test.go @@ -129,6 +129,9 @@ func TestModalTruncatePreservesUnstableModesAndAutoSelects(t *testing.T) { if _, err := ModalTruncate(sys, &ModalTruncateOptions{Order: 5}); !errors.Is(err, ErrInvalidOrder) { t.Fatalf("invalid order err = %v, want ErrInvalidOrder", err) } + if _, err := ModalTruncate(modalTestSystem(t), &ModalTruncateOptions{MaxRealPart: 5}); !errors.Is(err, ErrInvalidOrder) { + t.Fatalf("empty threshold selection err = %v, want ErrInvalidOrder", err) + } } func modalTestSystem(t *testing.T) *System { diff --git a/passivity.go b/passivity.go index 2ff121a..58f0ddc 100644 --- a/passivity.go +++ b/passivity.go @@ -211,8 +211,8 @@ func validatePassivityGrid(omega []float64, dt float64) error { if frequency > upper*(1+1e-12) { return fmt.Errorf("omega[%d]=%g exceeds Nyquist frequency %g: %w", i, frequency, upper, ErrDimensionMismatch) } - if i > 0 && frequency < omega[i-1] { - return fmt.Errorf("frequency grid is not sorted at index %d: %w", i, ErrDimensionMismatch) + if i > 0 && frequency <= omega[i-1] { + return fmt.Errorf("frequency grid is not strictly increasing at index %d: %w", i, ErrDimensionMismatch) } } return nil diff --git a/passivity_test.go b/passivity_test.go index e98b8c9..79fca1f 100644 --- a/passivity_test.go +++ b/passivity_test.go @@ -83,6 +83,9 @@ func TestSampledPassiveValidatesGridAndUsesDiscreteNyquist(t *testing.T) { if _, err := SampledPassive(discrete, &PassivityOptions{Omega: []float64{1, 0}}); !errors.Is(err, ErrDimensionMismatch) { t.Fatalf("unsorted grid err = %v, want ErrDimensionMismatch", err) } + if _, err := SampledPassive(discrete, &PassivityOptions{Omega: []float64{1, 1}}); !errors.Is(err, ErrDimensionMismatch) { + t.Fatalf("duplicate frequency err = %v, want ErrDimensionMismatch", err) + } if _, err := SampledPassive(discrete, &PassivityOptions{Tol: -1}); !errors.Is(err, ErrDimensionMismatch) { t.Fatalf("negative tolerance err = %v, want ErrDimensionMismatch", err) } diff --git a/tuning.go b/tuning.go index 663d0aa..9b523d9 100644 --- a/tuning.go +++ b/tuning.go @@ -74,14 +74,12 @@ func GridTune(model *GeneralizedClosedLoop, goals []TuningGoal, opts *SystuneOpt if err != nil { return err } - candidate := *model - candidate.controller = sampled - candidate.tunableController = sampled + candidate := model.withSampledController(sampled) closed, err := candidate.ClosedLoop(candidate.primaryAnalysisPointName()) if err != nil { return err } - goalResults, score, pass, err := evaluateTuningGoals(&candidate, closed, goals) + goalResults, score, pass, err := evaluateTuningGoals(candidate, closed, goals) if err != nil { return err } diff --git a/tuning_goal.go b/tuning_goal.go index 2f4e9f2..5938074 100644 --- a/tuning_goal.go +++ b/tuning_goal.go @@ -401,17 +401,28 @@ func frequencyGainRange(sys *System, omega []float64, outputWeight, inputWeight maxGain := 0.0 minGain := math.Inf(1) responseData := make([]complex128, resp.P*resp.M) + var outputData, outputProduct, inputData, inputProduct []complex128 + weightedRows := resp.P + if outputResponse != nil { + outputData = make([]complex128, outputResponse.P*outputResponse.M) + outputProduct = make([]complex128, outputResponse.P*resp.M) + weightedRows = outputResponse.P + } + if inputResponse != nil { + inputData = make([]complex128, inputResponse.P*inputResponse.M) + inputProduct = make([]complex128, weightedRows*inputResponse.M) + } var singularValues complexSingularValueWorkspace for k := range omega { gain := complexResponseAt(resp, k, responseData) if outputResponse != nil { - gain, err = multiplyComplexMatrices(complexResponseAt(outputResponse, k, nil), gain) + gain, err = multiplyComplexMatricesInto(outputProduct, complexResponseAt(outputResponse, k, outputData), gain) if err != nil { return 0, 0, fmt.Errorf("output weight: %w", err) } } if inputResponse != nil { - gain, err = multiplyComplexMatrices(gain, complexResponseAt(inputResponse, k, nil)) + gain, err = multiplyComplexMatricesInto(inputProduct, gain, complexResponseAt(inputResponse, k, inputData)) if err != nil { return 0, 0, fmt.Errorf("input weight: %w", err) } @@ -448,11 +459,17 @@ func complexResponseAt(response *FreqResponseMatrix, frequency int, data []compl return complexMatrix{rows: response.P, cols: response.M, data: data} } -func multiplyComplexMatrices(a, b complexMatrix) (complexMatrix, error) { +func multiplyComplexMatricesInto(dst []complex128, a, b complexMatrix) (complexMatrix, error) { if a.cols != b.rows { return complexMatrix{}, fmt.Errorf("matrix dimensions %dx%d and %dx%d: %w", a.rows, a.cols, b.rows, b.cols, ErrDimensionMismatch) } - result := complexMatrix{rows: a.rows, cols: b.cols, data: make([]complex128, a.rows*b.cols)} + if len(dst) != a.rows*b.cols { + dst = make([]complex128, a.rows*b.cols) + } + result := complexMatrix{rows: a.rows, cols: b.cols, data: dst} + for i := range result.data { + result.data[i] = 0 + } for i := range a.rows { for k := range a.cols { aik := a.data[i*a.cols+k] From a17ba975c72e9607b0ec4120ade10847e51748e5 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 18 Jul 2026 17:10:56 -0500 Subject: [PATCH 10/10] Fix descriptor sweep and implicit port binding --- benchmark_hotspots_test.go | 11 ++++++ frequency.go | 4 ++- physical_assembly.go | 71 ++++++++++++++++++++++++-------------- physical_assembly_test.go | 48 ++++++++++++++++++++++++++ 4 files changed, 108 insertions(+), 26 deletions(-) diff --git a/benchmark_hotspots_test.go b/benchmark_hotspots_test.go index 4e4fb2d..aa01fec 100644 --- a/benchmark_hotspots_test.go +++ b/benchmark_hotspots_test.go @@ -457,6 +457,17 @@ func BenchmarkDescriptorToExplicit_N10(b *testing.B) { } } +func BenchmarkDescriptorFreqResponse_MIMO_N10x200(b *testing.B) { + sys := benchDescriptorSystem(b, 10, 3, 3) + omega := logspace(-2, 3, 200) + b.ResetTimer() + for i := 0; i < b.N; i++ { + if _, err := sys.FreqResponse(omega); err != nil { + b.Fatal(err) + } + } +} + func BenchmarkFixedInputReduction_N10(b *testing.B) { sys := benchSysNonSym(10, 4, 2) fixed := map[int]float64{1: 0.5, 3: -2} diff --git a/frequency.go b/frequency.go index 22d1eab..f1f0e70 100644 --- a/frequency.go +++ b/frequency.go @@ -171,11 +171,13 @@ func (e frequencyEvaluator) response(omega []float64) (*FreqResponseMatrix, erro pm := e.p * e.m data := make([]complex128, nw*pm) if e.sys.IsDescriptor() { + ws := newSSEvalWorkspace(e.n, e.p, e.m) for k, w := range omega { s := e.sAt(w) - if err := e.evalStateSpaceInto(s, data[k*pm:(k+1)*pm]); err != nil { + if err := evalFrSSInto(ws, e.sys, s, e.n, e.p, e.m); err != nil { return nil, err } + copy(data[k*pm:(k+1)*pm], ws.g[:pm]) } applyIODelayPhase(e.sys, omega, data, e.p, e.m, true) return e.matrix(data, omega), nil diff --git a/physical_assembly.go b/physical_assembly.go index fe9fa65..bf7d590 100644 --- a/physical_assembly.go +++ b/physical_assembly.go @@ -135,10 +135,9 @@ func (p *physicalAssemblyPlan) bindPorts() (map[string]int, error) { portIndex := make(map[string]int) for componentIndex, component := range p.components { _, m, outputs := component.System.Dims() - inputCursor := 0 - outputCursor := 0 usedInputs := make([]bool, m) usedOutputs := make([]bool, outputs) + seenPorts := make(map[string]struct{}, len(component.Ports)) for _, port := range component.Ports { if port.Name == "" || port.Dimension <= 0 { return nil, fmt.Errorf("AssemblePhysical: invalid port on %q: %w", component.Name, ErrDimensionMismatch) @@ -147,31 +146,36 @@ func (p *physicalAssemblyPlan) bindPorts() (map[string]int, error) { return nil, fmt.Errorf("AssemblePhysical: invalid port kind on %s.%s: %w", component.Name, port.Name, ErrDimensionMismatch) } key := component.Name + "." + port.Name - if _, exists := portIndex[key]; exists { + if _, exists := seenPorts[key]; exists { return nil, fmt.Errorf("AssemblePhysical: duplicate port %q: %w", key, ErrDimensionMismatch) } - inputs := copyIntSlice(port.Input) - outputsForPort := copyIntSlice(port.Output) - if len(inputs) == 0 && len(outputsForPort) == 0 { - inputs = integerRange(inputCursor, port.Dimension) - outputsForPort = integerRange(outputCursor, port.Dimension) - inputCursor += port.Dimension - outputCursor += port.Dimension + seenPorts[key] = struct{}{} + if len(port.Input) == 0 && len(port.Output) == 0 { + continue } - if len(inputs) != port.Dimension || len(outputsForPort) != port.Dimension { + if len(port.Input) != port.Dimension || len(port.Output) != port.Dimension { return nil, fmt.Errorf("AssemblePhysical: port %q must bind %d input and output channels: %w", key, port.Dimension, ErrDimensionMismatch) } - for _, channel := range inputs { - if channel < 0 || channel >= m || usedInputs[channel] { - return nil, fmt.Errorf("AssemblePhysical: invalid or reused input channel %d on %q: %w", channel, key, ErrDimensionMismatch) - } - usedInputs[channel] = true + if err := reservePhysicalChannels(usedInputs, port.Input, "input", key); err != nil { + return nil, err + } + if err := reservePhysicalChannels(usedOutputs, port.Output, "output", key); err != nil { + return nil, err } - for _, channel := range outputsForPort { - if channel < 0 || channel >= outputs || usedOutputs[channel] { - return nil, fmt.Errorf("AssemblePhysical: invalid or reused output channel %d on %q: %w", channel, key, ErrDimensionMismatch) + } + for _, port := range component.Ports { + key := component.Name + "." + port.Name + inputs := copyIntSlice(port.Input) + outputsForPort := copyIntSlice(port.Output) + if len(inputs) == 0 && len(outputsForPort) == 0 { + inputs = claimUnusedPhysicalChannels(usedInputs, port.Dimension) + if len(inputs) != port.Dimension { + return nil, fmt.Errorf("AssemblePhysical: not enough unused input channels for %q: %w", key, ErrDimensionMismatch) + } + outputsForPort = claimUnusedPhysicalChannels(usedOutputs, port.Dimension) + if len(outputsForPort) != port.Dimension { + return nil, fmt.Errorf("AssemblePhysical: not enough unused output channels for %q: %w", key, ErrDimensionMismatch) } - usedOutputs[channel] = true } binding := physicalPortBinding{key: key, component: componentIndex, port: port, inputs: make([]int, port.Dimension), outputs: make([]int, port.Dimension)} for i := range port.Dimension { @@ -533,12 +537,29 @@ func physicalSignalNames(names []string, count int, prefix, kind string) []strin return out } -func integerRange(start, count int) []int { - values := make([]int, count) - for i := range count { - values[i] = start + i +func reservePhysicalChannels(used []bool, channels []int, signal, key string) error { + for _, channel := range channels { + if channel < 0 || channel >= len(used) || used[channel] { + return fmt.Errorf("AssemblePhysical: invalid or reused %s channel %d on %q: %w", signal, channel, key, ErrDimensionMismatch) + } + used[channel] = true + } + return nil +} + +func claimUnusedPhysicalChannels(used []bool, count int) []int { + channels := make([]int, 0, count) + for channel, reserved := range used { + if reserved { + continue + } + used[channel] = true + channels = append(channels, channel) + if len(channels) == count { + break + } } - return values + return channels } func prefixSystemMetadata(sys *System, prefix string) { diff --git a/physical_assembly_test.go b/physical_assembly_test.go index 47449e5..1548805 100644 --- a/physical_assembly_test.go +++ b/physical_assembly_test.go @@ -107,6 +107,37 @@ func TestPhysicalAssemblyPreservesUnconnectedDescriptorComponent(t *testing.T) { } } +func TestPhysicalAssemblyImplicitPortsSkipExplicitBindings(t *testing.T) { + explicit := PhysicalPort{Name: "external", Kind: PhysicalPortDisplacement, Dimension: 1, Input: []int{0}, Output: []int{0}} + implicit := PhysicalPort{Name: "mount", Kind: PhysicalPortDisplacement, Dimension: 2} + for _, test := range []struct { + name string + ports []PhysicalPort + }{ + {name: "explicit first", ports: []PhysicalPort{explicit, implicit}}, + {name: "explicit last", ports: []PhysicalPort{implicit, explicit}}, + } { + t.Run(test.name, func(t *testing.T) { + component := physicalMixedBindingComponent(t, "mixed", test.ports) + assembled, err := AssemblePhysical("grounded", []PhysicalComponent{component}, []PhysicalConnection{ + {FromComponent: "mixed", FromPort: "mount", Grounded: true}, + }) + if err != nil { + t.Fatalf("AssemblePhysical: %v", err) + } + if n, m, p := assembled.Dims(); n != 3 || m != 1 || p != 1 { + t.Fatalf("dims = (%d,%d,%d), want (3,1,1)", n, m, p) + } + if !sameStrings(assembled.InputName, []string{"mixed.external.force"}) { + t.Fatalf("input names = %v", assembled.InputName) + } + if !sameStrings(assembled.OutputName, []string{"mixed.external.position"}) { + t.Fatalf("output names = %v", assembled.OutputName) + } + }) + } +} + func TestPhysicalAssemblyRejectsInvalidTopologies(t *testing.T) { left := physicalCoupledComponent(t, "left", 1) right := physicalCoupledComponent(t, "right", 2) @@ -167,3 +198,20 @@ func physicalDescriptorComponent(t *testing.T, name string) PhysicalComponent { sys.StateName = []string{"x1", "x2"} return NewPhysicalComponent(name, sys, []PhysicalPort{{Name: "mount", Kind: PhysicalPortDisplacement, Dimension: 1}}) } + +func physicalMixedBindingComponent(t *testing.T, name string, ports []PhysicalPort) PhysicalComponent { + t.Helper() + sys, err := New( + mat.NewDense(1, 1, []float64{-1}), + mat.NewDense(1, 3, []float64{1, 2, 3}), + mat.NewDense(3, 1, []float64{1, 2, 3}), + mat.NewDense(3, 3, nil), + 0, + ) + if err != nil { + t.Fatal(err) + } + sys.InputName = []string{"external.force", "mount.force[0]", "mount.force[1]"} + sys.OutputName = []string{"external.position", "mount.position[0]", "mount.position[1]"} + return NewPhysicalComponent(name, sys, ports) +}