Repair control architecture and numerical contracts

.gitignore

@@ -4,3 +4,4 @@
 .DS_Store
 pprof_list.txt
 referance/
+.ergo/

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
@@ -90,11 +90,12 @@ 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 |
+| `NewGeneralizedClosedLoop` | Build a plant/controller closed-loop model with a plant-output 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 |
@@ -118,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 |
@@ -128,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
 
@@ -177,7 +180,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. 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 |
 
 ### Control Design
@@ -240,7 +244,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 |
@@ -334,7 +338,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 |
@@ -345,8 +350,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
 

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}
@@ -546,6 +557,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})

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}),

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
 
@@ -394,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
 

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<br/>ConcatModelArrays"]
         constructGeneralized["NewGeneralizedModel<br/>NewGeneralizedClosedLoop"]
         identifyERA["ERA<br/>Markov parameters to state-space model"]
         estimateFreqResponse["FreqRespEst<br/>sampled input/output to response estimate"]
         linearizeNonlinear["Linearize / NewEKF<br/>local nonlinear approximation"]
-        assemblePhysical["NewPhysicalComponent / AssemblePhysical<br/>port-checked component assembly"]
+        assemblePhysical["NewPhysicalComponent / AssemblePhysical<br/>constraint-based descriptor assembly"]
     end
 
     subgraph interconnection["Interconnection interfaces"]
@@ -62,7 +62,7 @@ flowchart LR
         energyAnalysis["Gram / HSV / Norm<br/>H2Norm / HinfNorm / Covar"]
         structureAnalysis["Ctrb / Obsv / IsStabilizable / IsDetectable"]
         loopAnalysis["Loopsens / RootLocus / FRDMargin"]
-        passivityAnalysis["Passive / FRDPassive / SpectralFactor"]
+        passivityAnalysis["SampledPassive / FRDPassive<br/>Passive compatibility alias / SpectralFactor"]
     end
 
     subgraph transforms["Transformation and reduction"]
@@ -78,7 +78,7 @@ flowchart LR
         observerDesign["Kalman / Kalmd / Lqe / Lqg<br/>Estim / Reg"]
         robustSynthesis["H2Syn / HinfSyn"]
         pidDesign["NewPID / NewPIDStd / Pidtune<br/>PID / PID2 / SmithPredictor"]
-        fixedStructureTuning["Systune / Looptune<br/>tuning goals"]
+        fixedStructureTuning["GridTune / Systune / Looptune<br/>point-specific tuning goals"]
     end
 
     caller --> constructStateSpace
@@ -229,6 +229,8 @@ classDiagram
     }
 
     class GeneralizedClosedLoop {
+        +InsertAnalysisPoint()
+        +AnalysisPoint()
         +OpenLoop()
         +ClosedLoop()
         +Sensitivity()
@@ -320,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 {
@@ -376,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
@@ -390,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
 ```
@@ -402,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.