|
| 1 | +// -*- mode:c++; fill-column: 100; -*- |
| 2 | + |
| 3 | +#include "ekf_fusion.h" |
| 4 | +#include <cmath> |
| 5 | +#include <iostream> |
| 6 | + |
| 7 | +namespace vesc_driver |
| 8 | +{ |
| 9 | + |
| 10 | +EKFFusion::EKFFusion(float wheelbase) |
| 11 | + : wheelbase_(wheelbase) |
| 12 | +{ |
| 13 | + // Initialize state to zero |
| 14 | + state_.setZero(); |
| 15 | + |
| 16 | + // Initialize covariance matrix with moderate uncertainty |
| 17 | + P_.setIdentity(); |
| 18 | + P_(0, 0) = 0.01; // x position variance |
| 19 | + P_(1, 1) = 0.01; // y position variance |
| 20 | + P_(2, 2) = 0.01; // theta variance |
| 21 | + P_(3, 3) = 0.1; // velocity variance |
| 22 | + P_(4, 4) = 0.1; // angular velocity variance |
| 23 | + |
| 24 | + // Process noise covariance - tuned for vehicle dynamics |
| 25 | + Q_.setIdentity(); |
| 26 | + Q_(0, 0) = 0.001; // x position process noise |
| 27 | + Q_(1, 1) = 0.001; // y position process noise |
| 28 | + Q_(2, 2) = 0.005; // theta process noise |
| 29 | + Q_(3, 3) = 0.1; // velocity process noise |
| 30 | + Q_(4, 4) = 0.1; // angular velocity process noise |
| 31 | + |
| 32 | + // IMU measurement noise - gyroscope noise |
| 33 | + R_imu_ = 0.01; // rad/s variance |
| 34 | +} |
| 35 | + |
| 36 | +void EKFFusion::predict(double dt, float rpm, float steering_angle, |
| 37 | + float speed_to_erpm_gain, float speed_to_erpm_offset) |
| 38 | +{ |
| 39 | + if (dt <= 0.0 || dt > 1.0) { |
| 40 | + // Skip invalid time steps |
| 41 | + return; |
| 42 | + } |
| 43 | + |
| 44 | + // Convert RPM to linear velocity |
| 45 | + float lin_vel = (rpm - speed_to_erpm_offset) / speed_to_erpm_gain; |
| 46 | + |
| 47 | + // Clamp velocity to zero for minuscule values |
| 48 | + if (std::fabs(lin_vel) < 0.01) { |
| 49 | + lin_vel = 0.0; |
| 50 | + } |
| 51 | + |
| 52 | + // Calculate angular velocity from steering angle |
| 53 | + float rot_vel = 0.0; |
| 54 | + if (std::fabs(steering_angle) > 1e-6) { |
| 55 | + float turn_radius = wheelbase_ / std::tan(steering_angle); |
| 56 | + if (std::fabs(turn_radius) > 1e-6) { |
| 57 | + rot_vel = lin_vel / turn_radius; |
| 58 | + } |
| 59 | + } |
| 60 | + |
| 61 | + // Extract current state |
| 62 | + float x = state_(0); |
| 63 | + float y = state_(1); |
| 64 | + float theta = state_(2); |
| 65 | + |
| 66 | + // Predict new state using odometry motion model |
| 67 | + float cos_theta = std::cos(theta); |
| 68 | + float sin_theta = std::sin(theta); |
| 69 | + |
| 70 | + state_(0) = x + lin_vel * dt * cos_theta; // x |
| 71 | + state_(1) = y + lin_vel * dt * sin_theta; // y |
| 72 | + state_(2) = theta + rot_vel * dt; // theta |
| 73 | + state_(3) = lin_vel; // v |
| 74 | + state_(4) = rot_vel; // omega |
| 75 | + |
| 76 | + // Normalize theta to [-pi, pi] |
| 77 | + while (state_(2) > M_PI) state_(2) -= 2.0 * M_PI; |
| 78 | + while (state_(2) < -M_PI) state_(2) += 2.0 * M_PI; |
| 79 | + |
| 80 | + // Compute Jacobian of motion model F |
| 81 | + Eigen::Matrix<float, 5, 5> F = Eigen::Matrix<float, 5, 5>::Identity(); |
| 82 | + F(0, 2) = -lin_vel * dt * sin_theta; // dx/dtheta |
| 83 | + F(0, 3) = dt * cos_theta; // dx/dv |
| 84 | + F(1, 2) = lin_vel * dt * cos_theta; // dy/dtheta |
| 85 | + F(1, 3) = dt * sin_theta; // dy/dv |
| 86 | + F(2, 4) = dt; // dtheta/domega |
| 87 | + |
| 88 | + // Predict covariance: P = F * P * F^T + Q |
| 89 | + P_ = F * P_ * F.transpose() + Q_; |
| 90 | +} |
| 91 | + |
| 92 | +void EKFFusion::updateIMU(float angular_velocity_z, bool available) |
| 93 | +{ |
| 94 | + if (!available) { |
| 95 | + return; // No IMU data available |
| 96 | + } |
| 97 | + |
| 98 | + // Measurement model: z = H * x + v |
| 99 | + // We're measuring omega (angular velocity) from IMU |
| 100 | + // H = [0, 0, 0, 0, 1] - selects the omega state |
| 101 | + |
| 102 | + Eigen::Matrix<float, 1, 5> H; |
| 103 | + H << 0, 0, 0, 0, 1; |
| 104 | + |
| 105 | + // Innovation (measurement residual) |
| 106 | + float y = angular_velocity_z - state_(4); |
| 107 | + |
| 108 | + // Innovation covariance |
| 109 | + float S = H * P_ * H.transpose() + R_imu_; |
| 110 | + |
| 111 | + // Kalman gain |
| 112 | + Eigen::Matrix<float, 5, 1> K = P_ * H.transpose() / S; |
| 113 | + |
| 114 | + // Update state estimate |
| 115 | + state_ = state_ + K * y; |
| 116 | + |
| 117 | + // Normalize theta |
| 118 | + while (state_(2) > M_PI) state_(2) -= 2.0 * M_PI; |
| 119 | + while (state_(2) < -M_PI) state_(2) += 2.0 * M_PI; |
| 120 | + |
| 121 | + // Update covariance estimate |
| 122 | + Eigen::Matrix<float, 5, 5> I = Eigen::Matrix<float, 5, 5>::Identity(); |
| 123 | + P_ = (I - K * H) * P_; |
| 124 | +} |
| 125 | + |
| 126 | +void EKFFusion::getState(float& x, float& y, float& theta, float& v, float& omega) const |
| 127 | +{ |
| 128 | + x = state_(0); |
| 129 | + y = state_(1); |
| 130 | + theta = state_(2); |
| 131 | + v = state_(3); |
| 132 | + omega = state_(4); |
| 133 | +} |
| 134 | + |
| 135 | +void EKFFusion::reset() |
| 136 | +{ |
| 137 | + state_.setZero(); |
| 138 | + P_.setIdentity(); |
| 139 | + P_(0, 0) = 0.01; |
| 140 | + P_(1, 1) = 0.01; |
| 141 | + P_(2, 2) = 0.01; |
| 142 | + P_(3, 3) = 0.1; |
| 143 | + P_(4, 4) = 0.1; |
| 144 | +} |
| 145 | + |
| 146 | +} // namespace vesc_driver |
0 commit comments