diff --git a/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/ExtendedKalmanFilter.java b/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/ExtendedKalmanFilter.java index 68a84e89aa..7474b02f74 100644 --- a/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/ExtendedKalmanFilter.java +++ b/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/ExtendedKalmanFilter.java @@ -39,13 +39,13 @@ public class ExtendedKalmanFilter, Matrix, Matrix> m_h; private final Matrix m_contQ; - private Matrix m_discR; private final Matrix m_initP; private final Matrix m_contR; @SuppressWarnings("MemberName") private Matrix m_xHat; @SuppressWarnings("MemberName") private Matrix m_P; + private double m_dtSeconds; /** * Constructs an extended Kalman filter. @@ -78,10 +78,11 @@ public class ExtendedKalmanFilter(inputs, Nat.N1())); @@ -92,13 +93,13 @@ public class ExtendedKalmanFilter(states, states); } @@ -223,7 +224,7 @@ public class ExtendedKalmanFilter u, Matrix y) { - correct(m_outputs, u, y, m_h, m_discR); + correct(m_outputs, u, y, m_h, m_contR); } /** @@ -261,8 +262,9 @@ public class ExtendedKalmanFilter R ) { final var C = NumericalJacobian.numericalJacobianX(rows, m_states, h, m_xHat, u); + final var discR = Discretization.discretizeR(R, m_dtSeconds); - final var S = C.times(m_P).times(C.transpose()).plus(R); + final var S = C.times(m_P).times(C.transpose()).plus(discR); // We want to put K = PC^T S^-1 into Ax = b form so we can solve it more // efficiently. diff --git a/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/UnscentedKalmanFilter.java b/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/UnscentedKalmanFilter.java index bac4511427..ca99153297 100644 --- a/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/UnscentedKalmanFilter.java +++ b/wpimath/src/main/java/edu/wpi/first/wpilibj/estimator/UnscentedKalmanFilter.java @@ -45,8 +45,8 @@ public class UnscentedKalmanFilter m_P; private final Matrix m_contQ; private final Matrix m_contR; - private Matrix m_discR; private Matrix m_sigmasF; + private double m_dtSeconds; private final MerweScaledSigmaPoints m_pts; @@ -61,7 +61,7 @@ public class UnscentedKalmanFilter states, Nat outputs, @@ -71,7 +71,7 @@ public class UnscentedKalmanFilter> h, Matrix stateStdDevs, Matrix measurementStdDevs, - double nominalDtSeconds) { + double dtSeconds) { this.m_states = states; this.m_outputs = outputs; @@ -81,7 +81,7 @@ public class UnscentedKalmanFilter(states); @@ -238,7 +238,7 @@ public class UnscentedKalmanFilter u, Matrix y) { - correct(m_outputs, u, y, m_h, m_discR); + correct(m_outputs, u, y, m_h, m_contR); } /** @@ -272,6 +272,8 @@ public class UnscentedKalmanFilter y, BiFunction, Matrix, Matrix> h, Matrix R) { + final var discR = Discretization.discretizeR(R, m_dtSeconds); + // Transform sigma points into measurement space Matrix sigmasH = new Matrix<>(new SimpleMatrix( rows.getNum(), 2 * m_states.getNum() + 1)); @@ -287,7 +289,7 @@ public class UnscentedKalmanFilter Pxy = new Matrix<>(m_states, rows); diff --git a/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h b/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h index 8d60ba8271..6f0fc85482 100644 --- a/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h +++ b/wpimath/src/main/native/include/frc/estimator/ExtendedKalmanFilter.h @@ -49,6 +49,7 @@ class ExtendedKalmanFilter { : m_f(f), m_h(h) { m_contQ = MakeCovMatrix(stateStdDevs); m_contR = MakeCovMatrix(measurementStdDevs); + m_dt = dt; Reset(); @@ -63,14 +64,15 @@ class ExtendedKalmanFilter { Eigen::Matrix discQ; DiscretizeAQTaylor(contA, m_contQ, dt, &discA, &discQ); - m_discR = DiscretizeR(m_contR, dt); + Eigen::Matrix discR = + DiscretizeR(m_contR, dt); // IsStabilizable(A^T, C^T) will tell us if the system is observable. bool isObservable = IsStabilizable(discA.transpose(), C.transpose()); if (isObservable && Outputs <= States) { m_initP = drake::math::DiscreteAlgebraicRiccatiEquation( - discA.transpose(), C.transpose(), discQ, m_discR); + discA.transpose(), C.transpose(), discQ, discR); } else { m_initP = Eigen::Matrix::Zero(); } @@ -139,6 +141,8 @@ class ExtendedKalmanFilter { * @param dt Timestep for prediction. */ void Predict(const Eigen::Matrix& u, units::second_t dt) { + m_dt = dt; + // Find continuous A Eigen::Matrix contA = NumericalJacobianX(m_f, m_xHat, u); @@ -150,7 +154,6 @@ class ExtendedKalmanFilter { m_xHat = RungeKutta(m_f, m_xHat, u, dt); m_P = discA * m_P * discA.transpose() + discQ; - m_discR = DiscretizeR(m_contR, dt); } /** @@ -161,7 +164,7 @@ class ExtendedKalmanFilter { */ void Correct(const Eigen::Matrix& u, const Eigen::Matrix& y) { - Correct(u, y, m_h, m_discR); + Correct(u, y, m_h, m_contR); } /** @@ -187,8 +190,9 @@ class ExtendedKalmanFilter { const Eigen::Matrix& R) { const Eigen::Matrix C = NumericalJacobianX(h, m_xHat, u); + const Eigen::Matrix discR = DiscretizeR(R, m_dt); - Eigen::Matrix S = C * m_P * C.transpose() + R; + Eigen::Matrix S = C * m_P * C.transpose() + discR; // We want to put K = PC^T S^-1 into Ax = b form so we can solve it more // efficiently. @@ -222,7 +226,7 @@ class ExtendedKalmanFilter { Eigen::Matrix m_P; Eigen::Matrix m_contQ; Eigen::Matrix m_contR; - Eigen::Matrix m_discR; + units::second_t m_dt; Eigen::Matrix m_initP; }; diff --git a/wpimath/src/main/native/include/frc/estimator/UnscentedKalmanFilter.h b/wpimath/src/main/native/include/frc/estimator/UnscentedKalmanFilter.h index f07d480115..8c2c31f0e6 100644 --- a/wpimath/src/main/native/include/frc/estimator/UnscentedKalmanFilter.h +++ b/wpimath/src/main/native/include/frc/estimator/UnscentedKalmanFilter.h @@ -50,8 +50,7 @@ class UnscentedKalmanFilter { : m_f(f), m_h(h) { m_contQ = MakeCovMatrix(stateStdDevs); m_contR = MakeCovMatrix(measurementStdDevs); - - m_discR = DiscretizeR(m_contR, dt); + m_dt = dt; Reset(); } @@ -119,6 +118,8 @@ class UnscentedKalmanFilter { * @param dt Timestep for prediction. */ void Predict(const Eigen::Matrix& u, units::second_t dt) { + m_dt = dt; + // Discretize Q before projecting mean and covariance forward Eigen::Matrix contA = NumericalJacobianX(m_f, m_xHat, u); @@ -141,7 +142,6 @@ class UnscentedKalmanFilter { m_P = std::get<1>(ret); m_P += discQ; - m_discR = DiscretizeR(m_contR, dt); } /** @@ -152,7 +152,7 @@ class UnscentedKalmanFilter { */ void Correct(const Eigen::Matrix& u, const Eigen::Matrix& y) { - Correct(u, y, m_h, m_discR); + Correct(u, y, m_h, m_contR); } /** @@ -176,6 +176,8 @@ class UnscentedKalmanFilter { const Eigen::Matrix&)> h, const Eigen::Matrix& R) { + const Eigen::Matrix discR = DiscretizeR(R, m_dt); + // Transform sigma points into measurement space Eigen::Matrix sigmasH; Eigen::Matrix sigmas = @@ -188,7 +190,7 @@ class UnscentedKalmanFilter { // Mean and covariance of prediction passed through UT auto [yHat, Py] = UnscentedTransform(sigmasH, m_pts.Wm(), m_pts.Wc()); - Py += R; + Py += discR; // Compute cross covariance of the state and the measurements Eigen::Matrix Pxy; @@ -224,8 +226,8 @@ class UnscentedKalmanFilter { Eigen::Matrix m_P; Eigen::Matrix m_contQ; Eigen::Matrix m_contR; - Eigen::Matrix m_discR; Eigen::Matrix m_sigmasF; + units::second_t m_dt; MerweScaledSigmaPoints m_pts; };