diff --git a/wpimath/src/main/java/edu/wpi/first/math/StateSpaceUtil.java b/wpimath/src/main/java/edu/wpi/first/math/StateSpaceUtil.java index 35bac8a0e9..8baf401a4b 100644 --- a/wpimath/src/main/java/edu/wpi/first/math/StateSpaceUtil.java +++ b/wpimath/src/main/java/edu/wpi/first/math/StateSpaceUtil.java @@ -87,9 +87,9 @@ public final class StateSpaceUtil { /** * Returns true if (A, B) is a stabilizable pair. * - *

(A,B) is stabilizable if and only if the uncontrollable eigenvalues of A, if any, have + *

(A, B) is stabilizable if and only if the uncontrollable eigenvalues of A, if any, have * absolute values less than one, where an eigenvalue is uncontrollable if rank(λI - A, B) %3C n - * where n is number of states. + * where n is the number of states. * * @param Num representing the size of A. * @param Num representing the columns of B. @@ -103,6 +103,26 @@ public final class StateSpaceUtil { return WPIMathJNI.isStabilizable(A.getNumRows(), B.getNumCols(), A.getData(), B.getData()); } + /** + * Returns true if (A, C) is a detectable pair. + * + *

(A, C) is detectable if and only if the unobservable eigenvalues of A, if any, have absolute + * values less than one, where an eigenvalue is unobservable if rank(λI - A; C) %3C n where n is + * the number of states. + * + * @param Num representing the size of A. + * @param Num representing the rows of C. + * @param A System matrix. + * @param C Output matrix. + * @return If the system is detectable. + */ + @SuppressWarnings("MethodTypeParameterName") + public static boolean isDetectable( + Matrix A, Matrix C) { + return WPIMathJNI.isStabilizable( + A.getNumRows(), C.getNumRows(), A.transpose().getData(), C.transpose().getData()); + } + /** * Convert a {@link Pose2d} to a vector of [x, y, theta], where theta is in radians. * diff --git a/wpimath/src/main/java/edu/wpi/first/math/WPIMathJNI.java b/wpimath/src/main/java/edu/wpi/first/math/WPIMathJNI.java index 1e2ccaf4e6..54445d3600 100644 --- a/wpimath/src/main/java/edu/wpi/first/math/WPIMathJNI.java +++ b/wpimath/src/main/java/edu/wpi/first/math/WPIMathJNI.java @@ -79,9 +79,9 @@ public final class WPIMathJNI { /** * Returns true if (A, B) is a stabilizable pair. * - *

(A,B) is stabilizable if and only if the uncontrollable eigenvalues of A, if any, have + *

(A, B) is stabilizable if and only if the uncontrollable eigenvalues of A, if any, have * absolute values less than one, where an eigenvalue is uncontrollable if rank(lambda * I - A, B) - * < n where n is number of states. + * < n where n is the number of states. * * @param states the number of states of the system. * @param inputs the number of inputs to the system. diff --git a/wpimath/src/main/java/edu/wpi/first/math/controller/LinearQuadraticRegulator.java b/wpimath/src/main/java/edu/wpi/first/math/controller/LinearQuadraticRegulator.java index 1c957b7069..cf4b26dea8 100644 --- a/wpimath/src/main/java/edu/wpi/first/math/controller/LinearQuadraticRegulator.java +++ b/wpimath/src/main/java/edu/wpi/first/math/controller/LinearQuadraticRegulator.java @@ -5,6 +5,7 @@ package edu.wpi.first.math.controller; import edu.wpi.first.math.Drake; +import edu.wpi.first.math.MathSharedStore; import edu.wpi.first.math.Matrix; import edu.wpi.first.math.Nat; import edu.wpi.first.math.Num; @@ -90,7 +91,7 @@ public class LinearQuadraticRegulator A, Matrix B, @@ -101,6 +102,18 @@ public class LinearQuadraticRegulator& A, const Eigen::Ref& B) { - // This function checks if (A,B) is a stabilizable pair. - // (A,B) is stabilizable if and only if the uncontrollable eigenvalues of - // A, if any, have absolute values less than one, where an eigenvalue is - // uncontrollable if Rank[lambda * I - A, B] < n. - int n = B.rows(), m = B.cols(); - Eigen::EigenSolver es(A); - for (int i = 0; i < n; i++) { + int states = B.rows(); + int inputs = B.cols(); + Eigen::EigenSolver es{A}; + for (int i = 0; i < states; ++i) { if (es.eigenvalues()[i].real() * es.eigenvalues()[i].real() + es.eigenvalues()[i].imag() * es.eigenvalues()[i].imag() < 1) { continue; } - Eigen::MatrixXcd E(n, n + m); - E << es.eigenvalues()[i] * Eigen::MatrixXcd::Identity(n, n) - A, B; - Eigen::ColPivHouseholderQR qr(E); - if (qr.rank() != n) { + Eigen::MatrixXcd E{states, states + inputs}; + E << es.eigenvalues()[i] * Eigen::MatrixXcd::Identity(states, states) - A, + B; + Eigen::ColPivHouseholderQR qr{E}; + if (qr.rank() < states) { return false; } } @@ -196,7 +204,7 @@ Java_edu_wpi_first_math_WPIMathJNI_isStabilizable Eigen::Matrix> B{nativeB, states, inputs}; - bool isStabilizable = check_stabilizable(A, B); // NOLINT + bool isStabilizable = check_stabilizable(A, B); env->ReleaseDoubleArrayElements(aSrc, nativeA, 0); env->ReleaseDoubleArrayElements(bSrc, nativeB, 0); diff --git a/wpimath/src/main/native/include/frc/StateSpaceUtil.h b/wpimath/src/main/native/include/frc/StateSpaceUtil.h index 2d41915169..3300a48ea0 100644 --- a/wpimath/src/main/native/include/frc/StateSpaceUtil.h +++ b/wpimath/src/main/native/include/frc/StateSpaceUtil.h @@ -78,7 +78,7 @@ bool IsStabilizableImpl(const Eigen::Matrix& A, Eigen::ColPivHouseholderQR< Eigen::Matrix, States, States + Inputs>> - qr(E); + qr{E}; if (qr.rank() < States) { return false; } @@ -258,9 +258,9 @@ Eigen::Vector PoseTo4dVector(const Pose2d& pose); /** * Returns true if (A, B) is a stabilizable pair. * - * (A,B) is stabilizable if and only if the uncontrollable eigenvalues of A, if + * (A, B) is stabilizable if and only if the uncontrollable eigenvalues of A, if * any, have absolute values less than one, where an eigenvalue is - * uncontrollable if rank(λI - A, B) < n where n is number of states. + * uncontrollable if rank(λI - A, B) < n where n is the number of states. * * @param A System matrix. * @param B Input matrix. @@ -271,6 +271,23 @@ bool IsStabilizable(const Eigen::Matrix& A, return detail::IsStabilizableImpl(A, B); } +/** + * Returns true if (A, C) is a detectable pair. + * + * (A, C) is detectable if and only if the unobservable eigenvalues of A, if + * any, have absolute values less than one, where an eigenvalue is unobservable + * if rank(λI - A; C) < n where n is the number of states. + * + * @param A System matrix. + * @param C Output matrix. + */ +template +bool IsDetectable(const Eigen::Matrix& A, + const Eigen::Matrix& C) { + return detail::IsStabilizableImpl(A.transpose(), + C.transpose()); +} + // Template specializations are used here to make common state-input pairs // compile faster. template <> diff --git a/wpimath/src/main/native/include/frc/controller/LinearQuadraticRegulator.h b/wpimath/src/main/native/include/frc/controller/LinearQuadraticRegulator.h index c324702b81..a10f1ee0ac 100644 --- a/wpimath/src/main/native/include/frc/controller/LinearQuadraticRegulator.h +++ b/wpimath/src/main/native/include/frc/controller/LinearQuadraticRegulator.h @@ -4,6 +4,10 @@ #pragma once +#include + +#include + #include #include @@ -16,6 +20,7 @@ #include "frc/system/LinearSystem.h" #include "units/time.h" #include "unsupported/Eigen/MatrixFunctions" +#include "wpimath/MathShared.h" namespace frc { namespace detail { @@ -82,6 +87,16 @@ class LinearQuadraticRegulatorImpl { Eigen::Matrix discB; DiscretizeAB(A, B, dt, &discA, &discB); + if (!IsStabilizable(discA, discB)) { + std::string msg = fmt::format( + "The system passed to the LQR is uncontrollable!\n\nA =\n{}\nB " + "=\n{}\n", + discA, discB); + + wpi::math::MathSharedStore::ReportError(msg); + throw std::invalid_argument(msg); + } + Eigen::Matrix S = drake::math::DiscreteAlgebraicRiccatiEquation(discA, discB, Q, R); diff --git a/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h b/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h index 985f2c14d9..8191f82b25 100644 --- a/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h +++ b/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h @@ -71,10 +71,7 @@ class ExtendedKalmanFilter { Eigen::Matrix discR = DiscretizeR(m_contR, dt); - // IsStabilizable(Aᵀ, Cᵀ) will tell us if the system is observable. - bool isObservable = - IsStabilizable(discA.transpose(), C.transpose()); - if (isObservable && Outputs <= States) { + if (IsDetectable(discA, C) && Outputs <= States) { m_initP = drake::math::DiscreteAlgebraicRiccatiEquation( discA.transpose(), C.transpose(), discQ, discR); } else { @@ -137,10 +134,7 @@ class ExtendedKalmanFilter { Eigen::Matrix discR = DiscretizeR(m_contR, dt); - // IsStabilizable(Aᵀ, Cᵀ) will tell us if the system is observable. - bool isObservable = - IsStabilizable(discA.transpose(), C.transpose()); - if (isObservable && Outputs <= States) { + if (IsDetectable(discA, C) && Outputs <= States) { m_initP = drake::math::DiscreteAlgebraicRiccatiEquation( discA.transpose(), C.transpose(), discQ, discR); } else { diff --git a/wpimath/src/main/native/include/frc/estimator/KalmanFilter.h b/wpimath/src/main/native/include/frc/estimator/KalmanFilter.h index 363753e5a3..437b76eb16 100644 --- a/wpimath/src/main/native/include/frc/estimator/KalmanFilter.h +++ b/wpimath/src/main/native/include/frc/estimator/KalmanFilter.h @@ -4,7 +4,10 @@ #pragma once +#include + #include +#include #include #include @@ -65,14 +68,14 @@ class KalmanFilterImpl { const auto& C = plant.C(); - // IsStabilizable(Aᵀ, Cᵀ) will tell us if the system is observable. - bool isObservable = - IsStabilizable(discA.transpose(), C.transpose()); - if (!isObservable) { - wpi::math::MathSharedStore::ReportError( - "The system passed to the Kalman filter is not observable!"); - throw std::invalid_argument( - "The system passed to the Kalman filter is not observable!"); + if (!IsDetectable(discA, C)) { + std::string msg = fmt::format( + "The system passed to the Kalman filter is " + "unobservable!\n\nA =\n{}\nC =\n{}\n", + discA, C); + + wpi::math::MathSharedStore::ReportError(msg); + throw std::invalid_argument(msg); } Eigen::Matrix P = diff --git a/wpimath/src/test/java/edu/wpi/first/math/StateSpaceUtilTest.java b/wpimath/src/test/java/edu/wpi/first/math/StateSpaceUtilTest.java index 795e76aa9d..456506c70c 100644 --- a/wpimath/src/test/java/edu/wpi/first/math/StateSpaceUtilTest.java +++ b/wpimath/src/test/java/edu/wpi/first/math/StateSpaceUtilTest.java @@ -76,6 +76,33 @@ public class StateSpaceUtilTest { assertTrue(StateSpaceUtil.isStabilizable(A, B)); } + @Test + @SuppressWarnings("LocalVariableName") + public void testIsDetectable() { + Matrix A; + Matrix C = Matrix.mat(Nat.N1(), Nat.N2()).fill(0, 1); + + // First eigenvalue is unobservable and unstable. + // Second eigenvalue is observable and stable. + A = Matrix.mat(Nat.N2(), Nat.N2()).fill(1.2, 0, 0, 0.5); + assertFalse(StateSpaceUtil.isDetectable(A, C)); + + // First eigenvalue is unobservable and marginally stable. + // Second eigenvalue is observable and stable. + A = Matrix.mat(Nat.N2(), Nat.N2()).fill(1, 0, 0, 0.5); + assertFalse(StateSpaceUtil.isDetectable(A, C)); + + // First eigenvalue is unobservable and stable. + // Second eigenvalue is observable and stable. + A = Matrix.mat(Nat.N2(), Nat.N2()).fill(0.2, 0, 0, 0.5); + assertTrue(StateSpaceUtil.isDetectable(A, C)); + + // First eigenvalue is unobservable and stable. + // Second eigenvalue is observable and unstable. + A = Matrix.mat(Nat.N2(), Nat.N2()).fill(0.2, 0, 0, 1.2); + assertTrue(StateSpaceUtil.isDetectable(A, C)); + } + @Test public void testMakeWhiteNoiseVector() { var firstData = new ArrayList(); diff --git a/wpimath/src/test/native/cpp/StateSpaceUtilTest.cpp b/wpimath/src/test/native/cpp/StateSpaceUtilTest.cpp index 3710cffd6c..57b93bb42b 100644 --- a/wpimath/src/test/native/cpp/StateSpaceUtilTest.cpp +++ b/wpimath/src/test/native/cpp/StateSpaceUtilTest.cpp @@ -134,3 +134,27 @@ TEST(StateSpaceUtilTest, IsStabilizable) { EXPECT_TRUE((frc::IsStabilizable<2, 1>( Eigen::Matrix{{0.2, 0}, {0, 1.2}}, B))); } + +TEST(StateSpaceUtilTest, IsDetectable) { + Eigen::Matrix C{0, 1}; + + // First eigenvalue is unobservable and unstable. + // Second eigenvalue is observable and stable. + EXPECT_FALSE((frc::IsDetectable<2, 1>( + Eigen::Matrix{{1.2, 0}, {0, 0.5}}, C))); + + // First eigenvalue is unobservable and marginally stable. + // Second eigenvalue is observable and stable. + EXPECT_FALSE((frc::IsDetectable<2, 1>( + Eigen::Matrix{{1, 0}, {0, 0.5}}, C))); + + // First eigenvalue is unobservable and stable. + // Second eigenvalue is observable and stable. + EXPECT_TRUE((frc::IsDetectable<2, 1>( + Eigen::Matrix{{0.2, 0}, {0, 0.5}}, C))); + + // First eigenvalue is unobservable and stable. + // Second eigenvalue is observable and unstable. + EXPECT_TRUE((frc::IsDetectable<2, 1>( + Eigen::Matrix{{0.2, 0}, {0, 1.2}}, C))); +}