From 982472edf88800fb9698d515d3006e73d89939b4 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 4 Apr 2026 08:02:49 -0500 Subject: [PATCH 1/3] fix: Canon modal falls back to Schur for ill-conditioned eigenvectors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit canonModal used eigendecomposition to build the transformation matrix T. When eigenvalues are nearly repeated (e.g. triple pole at ≈-0.01 from 1e6*s³+30000*s²+300*s+1), eigenvectors become nearly parallel, giving κ(T)≈8e13. Inverting T amplifies B entries to ~1e14, destroying the frequency response. Now checks κ(T)>1e12 (matching eigdecomp.go threshold) and falls back to ordered real Schur form (DGEES+DTREXC), which is numerically stable for near-repeated eigenvalues. Reported via python-control/python-control#1077 — the MIMO system there triggers the ill-conditioning. Co-Authored-By: Claude Opus 4.6 (1M context) --- canon.go | 120 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 118 insertions(+), 2 deletions(-) diff --git a/canon.go b/canon.go index f00dea5..10a829a 100644 --- a/canon.go +++ b/canon.go @@ -6,6 +6,9 @@ import ( "math/cmplx" "sort" + "gonum.org/v1/gonum/blas" + "gonum.org/v1/gonum/blas/blas64" + "gonum.org/v1/gonum/lapack" "gonum.org/v1/gonum/mat" ) @@ -49,6 +52,16 @@ func canonModal(sys *System) (*CanonResult, error) { return &CanonResult{Sys: cp, T: &mat.Dense{}}, nil } + res, err := canonModalEig(sys, n, m, p) + if err == nil { + return res, nil + } + return canonModalSchur(sys, n, m, p) +} + +// canonModalEig attempts the eigendecomposition-based modal form. +// Returns error if eigenvector matrix is too ill-conditioned. +func canonModalEig(sys *System, n, m, p int) (*CanonResult, error) { var eig mat.Eigen ok := eig.Factorize(sys.A, mat.EigenRight) if !ok { @@ -124,8 +137,10 @@ func canonModal(sys *System) (*CanonResult, error) { var lu mat.LU lu.Factorize(T) - if luNearSingular(&lu) { - return nil, fmt.Errorf("controlsys: transformation matrix is singular: %w", ErrSingularTransform) + + cond := lu.Cond() + if math.IsNaN(cond) || math.IsInf(cond, 1) || cond > 1e12 { + return nil, fmt.Errorf("controlsys: eigenvector matrix ill-conditioned (κ=%.1e)", cond) } Tinv := mat.NewDense(n, n, nil) @@ -159,6 +174,107 @@ func canonModal(sys *System) (*CanonResult, error) { return &CanonResult{Sys: newSys, T: T}, nil } +// canonModalSchur computes a quasi-modal form via ordered real Schur decomposition. +// Numerically stable even for near-repeated eigenvalues. +// A_modal is quasi-upper-triangular: 1×1 blocks for real eigenvalues, +// 2×2 blocks for complex pairs, ordered by eigenvalue magnitude. +func canonModalSchur(sys *System, n, m, p int) (*CanonResult, error) { + aRaw := sys.A.RawMatrix() + + t := make([]float64, n*n) + 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) + + workQuery := make([]float64, 1) + impl.Dgees(lapack.SchurHess, lapack.SortNone, nil, + n, t, n, wr, wi, z, n, workQuery, -1, bwork) + lwork := int(workQuery[0]) + work := make([]float64, lwork) + + _, ok := impl.Dgees(lapack.SchurHess, lapack.SortNone, nil, + n, t, n, wr, wi, z, n, work, lwork, bwork) + if !ok { + return nil, fmt.Errorf("controlsys: Schur decomposition failed: %w", ErrSchurFailed) + } + + // Sort Schur blocks by eigenvalue magnitude (ascending) + schurSortByMagnitude(t, z, n) + + bRaw := sys.B.RawMatrix() + cRaw := sys.C.RawMatrix() + + zGen := blas64.General{Rows: n, Cols: n, Stride: n, Data: z} + + // B_modal = Z' * B + bNew := make([]float64, n*m) + blas64.Gemm(blas.Trans, blas.NoTrans, 1, zGen, + blas64.General{Rows: n, Cols: m, Stride: bRaw.Stride, Data: bRaw.Data}, + 0, blas64.General{Rows: n, Cols: m, Stride: m, Data: bNew}) + + // C_modal = C * Z + cNew := make([]float64, p*n) + blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, + blas64.General{Rows: p, Cols: n, Stride: cRaw.Stride, Data: cRaw.Data}, + zGen, + 0, blas64.General{Rows: p, Cols: n, Stride: n, Data: cNew}) + + Anew := mat.NewDense(n, n, t) + Bnew := mat.NewDense(n, m, bNew) + Cnew := mat.NewDense(p, n, cNew) + Dnew := denseCopy(sys.D) + + newSys, err := newNoCopy(Anew, Bnew, Cnew, Dnew, sys.Dt) + if err != nil { + return nil, err + } + propagateIONames(newSys, sys) + + T := mat.NewDense(n, n, z) + return &CanonResult{Sys: newSys, T: T}, nil +} + +// schurSortByMagnitude reorders Schur blocks by eigenvalue magnitude (ascending). +func schurSortByMagnitude(t, z []float64, n int) { + trexcWork := make([]float64, n) + nPlaced := 0 + for nPlaced < n { + evals := schurEigenvaluesRaw(t, n) + + bestIdx := nPlaced + bestMag := cmplx.Abs(evals[nPlaced]) + i := nPlaced + for i < n { + blockSize := 1 + if i+1 < n && t[(i+1)*n+i] != 0 { + blockSize = 2 + } + mag := cmplx.Abs(evals[i]) + if mag < bestMag { + bestMag = mag + bestIdx = i + } + i += blockSize + } + + if bestIdx != nPlaced { + _, _, ok := impl.Dtrexc(lapack.UpdateSchur, n, t, n, z, n, bestIdx, nPlaced, trexcWork) + if !ok { + break + } + } + + if nPlaced+1 < n && t[(nPlaced+1)*n+nPlaced] != 0 { + nPlaced += 2 + } else { + nPlaced++ + } + } +} + func canonCompanion(sys *System) (*CanonResult, error) { n, m, p := sys.Dims() if n == 0 { From 44d797105c4833c21c85fef52657e8dc863991e6 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 4 Apr 2026 08:16:46 -0500 Subject: [PATCH 2/3] fix: harden input validation (improper TF, descriptor Riccati, Q PSD) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three defensive checks inspired by python-control edge cases: 1. StateSpace() rejects improper transfer functions (num deg > den deg) before conversion, preventing silent wrong results. (pc#235) 2. Kalman, Lqg, H2Syn, HinfSyn reject descriptor systems (E != nil) since standard CARE/DARE ignore E, producing wrong results. (pc#36) 3. Care/Dare validate Q is positive semi-definite via Cholesky on Q + ε‖Q‖nI. Negative/indefinite Q was silently accepted. (pc#252) R positive-definiteness was already caught by Cholesky. Co-Authored-By: Claude Opus 4.6 (1M context) --- errors.go | 3 ++ h2syn.go | 3 ++ hardening_test.go | 135 ++++++++++++++++++++++++++++++++++++++++++++++ hinfsyn.go | 3 ++ lqg.go | 3 ++ matutil.go | 20 +++++++ observer.go | 3 ++ riccati.go | 6 +++ transfer.go | 4 ++ 9 files changed, 180 insertions(+) create mode 100644 hardening_test.go diff --git a/errors.go b/errors.go index dde128f..5378fde 100644 --- a/errors.go +++ b/errors.go @@ -9,6 +9,7 @@ var ( ErrInvalidSampleTime = errors.New("controlsys: sample time must be positive") ErrSingularDenom = errors.New("controlsys: zero or near-zero leading denominator") ErrOverflow = errors.New("controlsys: coefficient overflow") + ErrImproperTF = errors.New("controlsys: transfer function is improper (numerator degree exceeds denominator degree)") ErrNegativeDelay = errors.New("controlsys: delay must be non-negative") ErrFractionalDelay = errors.New("controlsys: discrete delay must be non-negative integer") @@ -23,6 +24,7 @@ var ( ErrMixedDelayTypes = errors.New("controlsys: InternalDelay and IODelay cannot coexist") ErrNotSymmetric = errors.New("controlsys: matrix must be symmetric") + ErrNotPSD = errors.New("controlsys: Q matrix must be positive semi-definite") ErrSchurFailed = errors.New("controlsys: Schur decomposition failed to converge") ErrSingularEquation = errors.New("controlsys: matrix equation is singular or nearly singular") ErrNoStabilizing = errors.New("controlsys: no stabilizing solution exists") @@ -48,4 +50,5 @@ var ( ErrH2DirectFeedthrough = errors.New("controlsys: H2 synthesis requires D22 = 0") ErrGammaNotAchievable = errors.New("controlsys: no stabilizing controller exists for given gamma") ErrDescriptorSingular = errors.New("controlsys: descriptor matrix E is singular") + ErrDescriptorRiccati = errors.New("controlsys: standard Riccati solvers do not support descriptor systems (E != I)") ) diff --git a/h2syn.go b/h2syn.go index f25f648..f1c9734 100644 --- a/h2syn.go +++ b/h2syn.go @@ -17,6 +17,9 @@ func H2Syn(P *System, nmeas, ncont int) (*H2SynResult, error) { if !P.IsContinuous() { return nil, ErrWrongDomain } + if P.IsDescriptor() { + return nil, ErrDescriptorRiccati + } n, m, p := P.Dims() if n == 0 { diff --git a/hardening_test.go b/hardening_test.go new file mode 100644 index 0000000..f3bde6a --- /dev/null +++ b/hardening_test.go @@ -0,0 +1,135 @@ +package controlsys + +import ( + "errors" + "testing" + + "gonum.org/v1/gonum/mat" +) + +func TestStateSpace_RejectsImproperTF(t *testing.T) { + tf := &TransferFunc{ + Num: [][][]float64{{{1, 0}}}, // s (degree 1) + Den: [][]float64{{1}}, // 1 (degree 0) + } + _, err := tf.StateSpace(nil) + if !errors.Is(err, ErrImproperTF) { + t.Errorf("expected ErrImproperTF, got %v", err) + } +} + +func TestStateSpace_AcceptsProperTF(t *testing.T) { + tf := &TransferFunc{ + Num: [][][]float64{{{1}}}, + Den: [][]float64{{1, 1}}, + } + _, err := tf.StateSpace(nil) + if err != nil { + t.Errorf("proper TF should be accepted, got %v", err) + } +} + +func TestStateSpace_AcceptsBiproperTF(t *testing.T) { + tf := &TransferFunc{ + Num: [][][]float64{{{1, 1}}}, + Den: [][]float64{{1, 2}}, + } + _, err := tf.StateSpace(nil) + if err != nil { + t.Errorf("biproper TF (equal degree) should be accepted, got %v", err) + } +} + +func TestCare_RejectsNegativeDefiniteQ(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0, 1, -2, -3}) + B := mat.NewDense(2, 1, []float64{0, 1}) + Q := mat.NewDense(2, 2, []float64{-1, 0, 0, -1}) + R := mat.NewDense(1, 1, []float64{1}) + + _, err := Care(A, B, Q, R, nil) + if !errors.Is(err, ErrNotPSD) { + t.Errorf("expected ErrNotPSD for negative definite Q, got %v", err) + } +} + +func TestCare_RejectsIndefiniteQ(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0, 1, -2, -3}) + B := mat.NewDense(2, 1, []float64{0, 1}) + Q := mat.NewDense(2, 2, []float64{1, 0, 0, -1}) + R := mat.NewDense(1, 1, []float64{1}) + + _, err := Care(A, B, Q, R, nil) + if !errors.Is(err, ErrNotPSD) { + t.Errorf("expected ErrNotPSD for indefinite Q, got %v", err) + } +} + +func TestCare_AcceptsPSDQ(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0, 1, -2, -3}) + B := mat.NewDense(2, 1, []float64{0, 1}) + Q := mat.NewDense(2, 2, []float64{1, 0, 0, 0}) // PSD but singular + R := mat.NewDense(1, 1, []float64{1}) + + _, err := Care(A, B, Q, R, nil) + if err != nil { + t.Errorf("PSD Q should be accepted, got %v", err) + } +} + +func TestDare_RejectsNegativeDefiniteQ(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0.9, 0.1, 0, 0.8}) + B := mat.NewDense(2, 1, []float64{0, 1}) + Q := mat.NewDense(2, 2, []float64{-1, 0, 0, -1}) + R := mat.NewDense(1, 1, []float64{1}) + + _, err := Dare(A, B, Q, R, nil) + if !errors.Is(err, ErrNotPSD) { + t.Errorf("expected ErrNotPSD for negative definite Q, got %v", err) + } +} + +func TestCare_RejectsNonPDR(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0, 1, -2, -3}) + B := mat.NewDense(2, 1, []float64{0, 1}) + Q := mat.NewDense(2, 2, []float64{1, 0, 0, 1}) + R := mat.NewDense(1, 1, []float64{-1}) + + _, err := Care(A, B, Q, R, nil) + if !errors.Is(err, ErrSingularR) { + t.Errorf("expected ErrSingularR for non-PD R, got %v", err) + } +} + +func TestKalman_RejectsDescriptor(t *testing.T) { + sys, _ := New( + mat.NewDense(2, 2, []float64{-1, 0, 0, -2}), + mat.NewDense(2, 1, []float64{1, 0}), + mat.NewDense(1, 2, []float64{1, 1}), + mat.NewDense(1, 1, []float64{0}), 0) + sys.E = mat.NewDense(2, 2, []float64{1, 0, 0, 2}) + + Qn := mat.NewDense(1, 1, []float64{1}) + Rn := mat.NewDense(1, 1, []float64{1}) + _, err := Kalman(sys, Qn, Rn, nil) + if !errors.Is(err, ErrDescriptorRiccati) { + t.Errorf("expected ErrDescriptorRiccati, got %v", err) + } +} + +func TestLqg_RejectsDescriptor(t *testing.T) { + sys, _ := New( + mat.NewDense(2, 2, []float64{-1, 0, 0, -2}), + mat.NewDense(2, 1, []float64{1, 0}), + mat.NewDense(1, 2, []float64{1, 1}), + mat.NewDense(1, 1, []float64{0}), 0) + sys.E = mat.NewDense(2, 2, []float64{1, 0, 0, 2}) + + Q := mat.NewDense(2, 2, []float64{1, 0, 0, 1}) + R := mat.NewDense(1, 1, []float64{1}) + Qn := mat.NewDense(1, 1, []float64{1}) + Rn := mat.NewDense(1, 1, []float64{1}) + _, err := Lqg(sys, Q, R, Qn, Rn, nil) + if !errors.Is(err, ErrDescriptorRiccati) { + t.Errorf("expected ErrDescriptorRiccati, got %v", err) + } +} diff --git a/hinfsyn.go b/hinfsyn.go index 3d2672a..8251544 100644 --- a/hinfsyn.go +++ b/hinfsyn.go @@ -20,6 +20,9 @@ func HinfSyn(P *System, nmeas, ncont int) (*HinfSynResult, error) { if !P.IsContinuous() { return nil, ErrWrongDomain } + if P.IsDescriptor() { + return nil, ErrDescriptorRiccati + } n, m, p := P.Dims() if n == 0 { diff --git a/lqg.go b/lqg.go index 62e15d8..3e47fa1 100644 --- a/lqg.go +++ b/lqg.go @@ -22,6 +22,9 @@ func Lqg(sys *System, Q, R, Qn, Rn *mat.Dense, opts *RiccatiOpts) (*LqgResult, e if n == 0 { return nil, fmt.Errorf("Lqg: system has no states: %w", ErrDimensionMismatch) } + if sys.IsDescriptor() { + return nil, fmt.Errorf("Lqg: %w", ErrDescriptorRiccati) + } var kRes *RiccatiResult var err error diff --git a/matutil.go b/matutil.go index 3690d79..4b61585 100644 --- a/matutil.go +++ b/matutil.go @@ -3,6 +3,7 @@ package controlsys import ( "math" + "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/mat" ) @@ -89,6 +90,25 @@ func isSymmetric(m *mat.Dense, tol float64) bool { return true } +// isPSD checks if a symmetric matrix is positive semi-definite +// by attempting Cholesky on Q + tol*I. A small shift avoids +// rejecting borderline-zero eigenvalues from roundoff. +func isPSD(m *mat.Dense) bool { + n, _ := m.Dims() + if n == 0 { + return true + } + raw := m.RawMatrix() + tmp := make([]float64, n*n) + copyStrided(tmp, n, raw.Data, raw.Stride, n, n) + symmetrize(tmp, n, n) + tol := eps() * denseNorm(m) * float64(n) + for i := 0; i < n; i++ { + tmp[i*n+i] += tol + } + return impl.Dpotrf(blas.Lower, n, tmp, n) +} + func symmetrize(data []float64, n, stride int) { for i := 0; i < n; i++ { for j := i + 1; j < n; j++ { diff --git a/observer.go b/observer.go index a01375e..de3746d 100644 --- a/observer.go +++ b/observer.go @@ -97,6 +97,9 @@ func Kalman(sys *System, Qn, Rn *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, if n == 0 { return nil, fmt.Errorf("Kalman: system has no states: %w", ErrDimensionMismatch) } + if sys.IsDescriptor() { + return nil, fmt.Errorf("Kalman: %w", ErrDescriptorRiccati) + } qr, qc := Qn.Dims() if qr != m || qc != m { return nil, ErrDimensionMismatch diff --git a/riccati.go b/riccati.go index 1b539c1..1e186c2 100644 --- a/riccati.go +++ b/riccati.go @@ -120,6 +120,9 @@ func Care(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { if !isSymmetric(R, eps()*denseNorm(R)) { return nil, ErrNotSymmetric } + if !isPSD(Q) { + return nil, ErrNotPSD + } var S *mat.Dense if opts != nil && opts.S != nil { @@ -328,6 +331,9 @@ func Dare(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { if !isSymmetric(R, eps()*denseNorm(R)) { return nil, ErrNotSymmetric } + if !isPSD(Q) { + return nil, ErrNotPSD + } var S *mat.Dense if opts != nil && opts.S != nil { diff --git a/transfer.go b/transfer.go index 57beef2..388ea00 100644 --- a/transfer.go +++ b/transfer.go @@ -484,6 +484,10 @@ func (tf *TransferFunc) StateSpace(opts *StateSpaceOpts) (*StateSpaceResult, err return &StateSpaceResult{Sys: sys, MinimalOrder: 0}, nil } + if !tf.Isproper() { + return nil, ErrImproperTF + } + degrees := make([]int, p) totalN := 0 for i := 0; i < p; i++ { From 2392746bfbb2d9d805bc389ce4c88627a7760ddc Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sat, 4 Apr 2026 08:34:49 -0500 Subject: [PATCH 3/3] fix: numerical hardening from python-control PR audit Three fixes from auditing merged python-control PRs: 1. canon.go: Replace explicit Tinv computation with lu.SolveTo() in both canonModalEig and canonCompanion. Avoids forming the inverse matrix, reducing numerical error amplification. (pc PR#101) 2. convert.go: matchedGainFallback used cmplx.Abs ratio which discards sign, producing wrong gain for systems with negative DC gain. Now uses real(contVal/discVal) matching the non-fallback path. (pc PR#951) 3. ctrbobsv.go: IsStabilizable continuous-time check used zero tolerance (real(v) >= 0), so roundoff eigenvalue at Re=1e-16 was falsely classified as unstable. Now uses relative tolerance consistent with the discrete-time check. (pc PR#345) Co-Authored-By: Claude Opus 4.6 (1M context) --- canon.go | 46 ++++++++++++++++++++-------------------------- convert.go | 2 +- ctrbobsv.go | 6 ++++-- 3 files changed, 25 insertions(+), 29 deletions(-) diff --git a/canon.go b/canon.go index 10a829a..498a2ae 100644 --- a/canon.go +++ b/canon.go @@ -143,22 +143,19 @@ func canonModalEig(sys *System, n, m, p int) (*CanonResult, error) { return nil, fmt.Errorf("controlsys: eigenvector matrix ill-conditioned (κ=%.1e)", cond) } - Tinv := mat.NewDense(n, n, nil) - eye := mat.NewDense(n, n, nil) - for i := 0; i < n; i++ { - eye.Set(i, i, 1) - } - if err := lu.SolveTo(Tinv, false, eye); err != nil { - return nil, fmt.Errorf("controlsys: inversion failed: %w", ErrSingularTransform) - } - + // Anew = T⁻¹*A*T: compute AT = A*T, then solve T*Anew = AT + AT := mat.NewDense(n, n, nil) + AT.Mul(sys.A, T) Anew := mat.NewDense(n, n, nil) - tmp := mat.NewDense(n, n, nil) - tmp.Mul(Tinv, sys.A) - Anew.Mul(tmp, T) + if err := lu.SolveTo(Anew, false, AT); err != nil { + return nil, fmt.Errorf("controlsys: solve failed: %w", ErrSingularTransform) + } + // Bnew = T⁻¹*B: solve T*Bnew = B Bnew := mat.NewDense(n, m, nil) - Bnew.Mul(Tinv, sys.B) + if err := lu.SolveTo(Bnew, false, sys.B); err != nil { + return nil, fmt.Errorf("controlsys: solve failed: %w", ErrSingularTransform) + } Cnew := mat.NewDense(p, n, nil) Cnew.Mul(sys.C, T) @@ -328,22 +325,19 @@ func canonCompanion(sys *System) (*CanonResult, error) { return nil, fmt.Errorf("controlsys: transformation matrix singular: %w", ErrSingularTransform) } - Tinv := mat.NewDense(n, n, nil) - eye := mat.NewDense(n, n, nil) - for i := 0; i < n; i++ { - eye.Set(i, i, 1) - } - if err := luT.SolveTo(Tinv, false, eye); err != nil { - return nil, fmt.Errorf("controlsys: inversion failed: %w", ErrSingularTransform) - } - + // Anew = T⁻¹*A*T: compute AT = A*T, then solve T*Anew = AT + AT := mat.NewDense(n, n, nil) + AT.Mul(sys.A, T) Anew := mat.NewDense(n, n, nil) - tmp := mat.NewDense(n, n, nil) - tmp.Mul(Tinv, sys.A) - Anew.Mul(tmp, T) + if err := luT.SolveTo(Anew, false, AT); err != nil { + return nil, fmt.Errorf("controlsys: solve failed: %w", ErrSingularTransform) + } + // Bnew = T⁻¹*B: solve T*Bnew = B Bnew := mat.NewDense(n, 1, nil) - Bnew.Mul(Tinv, sys.B) + if err := luT.SolveTo(Bnew, false, sys.B); err != nil { + return nil, fmt.Errorf("controlsys: solve failed: %w", ErrSingularTransform) + } Cnew := mat.NewDense(1, n, nil) Cnew.Mul(sys.C, T) diff --git a/convert.go b/convert.go index a84fd99..f7698e8 100644 --- a/convert.go +++ b/convert.go @@ -1411,7 +1411,7 @@ func matchedGainFallback(num, den Poly, discZeros, discPoles []complex128, dt fl if cmplx.Abs(discVal) < 1e-14 { return 0, fmt.Errorf("DiscretizeMatched: cannot determine gain: %w", ErrSingularTransform) } - return cmplx.Abs(contVal) / cmplx.Abs(discVal), nil + return real(contVal / discVal), nil } func (sys *System) D2D(newDt float64, opts C2DOptions) (*System, error) { diff --git a/ctrbobsv.go b/ctrbobsv.go index 0f2817b..bd5c55f 100644 --- a/ctrbobsv.go +++ b/ctrbobsv.go @@ -1,6 +1,7 @@ package controlsys import ( + "math" "math/cmplx" "gonum.org/v1/gonum/blas" @@ -193,12 +194,13 @@ func IsStabilizable(A, B *mat.Dense, continuous bool) (bool, error) { vals := eig.Values(nil) for _, v := range vals { + tol := 1e-10 * math.Max(1, cmplx.Abs(v)) if continuous { - if real(v) >= 0 { + if real(v) >= -tol { return false, nil } } else { - if cmplx.Abs(v) >= 1-1e-10 { + if cmplx.Abs(v) >= 1-tol { return false, nil } }