00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "robot_localization/ekf.h"
00034 #include "robot_localization/filter_common.h"
00035
00036 #include <XmlRpcException.h>
00037
00038 #include <iomanip>
00039 #include <limits>
00040 #include <sstream>
00041 #include <vector>
00042
00043 namespace RobotLocalization
00044 {
00045 Ekf::Ekf(std::vector<double>) :
00046 FilterBase()
00047 {
00048 }
00049
00050 Ekf::~Ekf()
00051 {
00052 }
00053
00054 void Ekf::correct(const Measurement &measurement)
00055 {
00056 FB_DEBUG("---------------------- Ekf::correct ----------------------\n" <<
00057 "State is:\n" << state_ << "\n"
00058 "Topic is:\n" << measurement.topicName_ << "\n"
00059 "Measurement is:\n" << measurement.measurement_ << "\n"
00060 "Measurement topic name is:\n" << measurement.topicName_ << "\n\n"
00061 "Measurement covariance is:\n" << measurement.covariance_ << "\n");
00062
00063
00064
00065
00066
00067
00068 std::vector<size_t> updateIndices;
00069 for (size_t i = 0; i < measurement.updateVector_.size(); ++i)
00070 {
00071 if (measurement.updateVector_[i])
00072 {
00073
00074 if (std::isnan(measurement.measurement_(i)))
00075 {
00076 FB_DEBUG("Value at index " << i << " was nan. Excluding from update.\n");
00077 }
00078 else if (std::isinf(measurement.measurement_(i)))
00079 {
00080 FB_DEBUG("Value at index " << i << " was inf. Excluding from update.\n");
00081 }
00082 else
00083 {
00084 updateIndices.push_back(i);
00085 }
00086 }
00087 }
00088
00089 FB_DEBUG("Update indices are:\n" << updateIndices << "\n");
00090
00091 size_t updateSize = updateIndices.size();
00092
00093
00094 Eigen::VectorXd stateSubset(updateSize);
00095 Eigen::VectorXd measurementSubset(updateSize);
00096 Eigen::MatrixXd measurementCovarianceSubset(updateSize, updateSize);
00097 Eigen::MatrixXd stateToMeasurementSubset(updateSize, state_.rows());
00098 Eigen::MatrixXd kalmanGainSubset(state_.rows(), updateSize);
00099 Eigen::VectorXd innovationSubset(updateSize);
00100
00101 stateSubset.setZero();
00102 measurementSubset.setZero();
00103 measurementCovarianceSubset.setZero();
00104 stateToMeasurementSubset.setZero();
00105 kalmanGainSubset.setZero();
00106 innovationSubset.setZero();
00107
00108
00109 for (size_t i = 0; i < updateSize; ++i)
00110 {
00111 measurementSubset(i) = measurement.measurement_(updateIndices[i]);
00112 stateSubset(i) = state_(updateIndices[i]);
00113
00114 for (size_t j = 0; j < updateSize; ++j)
00115 {
00116 measurementCovarianceSubset(i, j) = measurement.covariance_(updateIndices[i], updateIndices[j]);
00117 }
00118
00119
00120
00121
00122 if (measurementCovarianceSubset(i, i) < 0)
00123 {
00124 FB_DEBUG("WARNING: Negative covariance for index " << i <<
00125 " of measurement (value is" << measurementCovarianceSubset(i, i) <<
00126 "). Using absolute value...\n");
00127
00128 measurementCovarianceSubset(i, i) = ::fabs(measurementCovarianceSubset(i, i));
00129 }
00130
00131
00132
00133
00134
00135
00136
00137 if (measurementCovarianceSubset(i, i) < 1e-9)
00138 {
00139 FB_DEBUG("WARNING: measurement had very small error covariance for index " << updateIndices[i] <<
00140 ". Adding some noise to maintain filter stability.\n");
00141
00142 measurementCovarianceSubset(i, i) = 1e-9;
00143 }
00144 }
00145
00146
00147
00148 for (size_t i = 0; i < updateSize; ++i)
00149 {
00150 stateToMeasurementSubset(i, updateIndices[i]) = 1;
00151 }
00152
00153 FB_DEBUG("Current state subset is:\n" << stateSubset <<
00154 "\nMeasurement subset is:\n" << measurementSubset <<
00155 "\nMeasurement covariance subset is:\n" << measurementCovarianceSubset <<
00156 "\nState-to-measurement subset is:\n" << stateToMeasurementSubset << "\n");
00157
00158
00159 Eigen::MatrixXd pht = estimateErrorCovariance_ * stateToMeasurementSubset.transpose();
00160 Eigen::MatrixXd hphrInv = (stateToMeasurementSubset * pht + measurementCovarianceSubset).inverse();
00161 kalmanGainSubset.noalias() = pht * hphrInv;
00162
00163 innovationSubset = (measurementSubset - stateSubset);
00164
00165
00166 for (size_t i = 0; i < updateSize; ++i)
00167 {
00168 if (updateIndices[i] == StateMemberRoll ||
00169 updateIndices[i] == StateMemberPitch ||
00170 updateIndices[i] == StateMemberYaw)
00171 {
00172 while (innovationSubset(i) < -PI)
00173 {
00174 innovationSubset(i) += TAU;
00175 }
00176
00177 while (innovationSubset(i) > PI)
00178 {
00179 innovationSubset(i) -= TAU;
00180 }
00181 }
00182 }
00183
00184
00185 if (checkMahalanobisThreshold(innovationSubset, hphrInv, measurement.mahalanobisThresh_))
00186 {
00187
00188 state_.noalias() += kalmanGainSubset * innovationSubset;
00189
00190
00191 Eigen::MatrixXd gainResidual = identity_;
00192 gainResidual.noalias() -= kalmanGainSubset * stateToMeasurementSubset;
00193 estimateErrorCovariance_ = gainResidual * estimateErrorCovariance_ * gainResidual.transpose();
00194 estimateErrorCovariance_.noalias() += kalmanGainSubset *
00195 measurementCovarianceSubset *
00196 kalmanGainSubset.transpose();
00197
00198
00199 wrapStateAngles();
00200
00201 FB_DEBUG("Kalman gain subset is:\n" << kalmanGainSubset <<
00202 "\nInnovation is:\n" << innovationSubset <<
00203 "\nCorrected full state is:\n" << state_ <<
00204 "\nCorrected full estimate error covariance is:\n" << estimateErrorCovariance_ <<
00205 "\n\n---------------------- /Ekf::correct ----------------------\n");
00206 }
00207 }
00208
00209 void Ekf::predict(const double referenceTime, const double delta)
00210 {
00211 FB_DEBUG("---------------------- Ekf::predict ----------------------\n" <<
00212 "delta is " << delta << "\n" <<
00213 "state is " << state_ << "\n");
00214
00215 double roll = state_(StateMemberRoll);
00216 double pitch = state_(StateMemberPitch);
00217 double yaw = state_(StateMemberYaw);
00218 double xVel = state_(StateMemberVx);
00219 double yVel = state_(StateMemberVy);
00220 double zVel = state_(StateMemberVz);
00221 double rollVel = state_(StateMemberVroll);
00222 double pitchVel = state_(StateMemberVpitch);
00223 double yawVel = state_(StateMemberVyaw);
00224 double xAcc = state_(StateMemberAx);
00225 double yAcc = state_(StateMemberAy);
00226 double zAcc = state_(StateMemberAz);
00227
00228
00229 double sp = ::sin(pitch);
00230 double cp = ::cos(pitch);
00231
00232 double sr = ::sin(roll);
00233 double cr = ::cos(roll);
00234
00235 double sy = ::sin(yaw);
00236 double cy = ::cos(yaw);
00237
00238 prepareControl(referenceTime, delta);
00239
00240
00241 transferFunction_(StateMemberX, StateMemberVx) = cy * cp * delta;
00242 transferFunction_(StateMemberX, StateMemberVy) = (cy * sp * sr - sy * cr) * delta;
00243 transferFunction_(StateMemberX, StateMemberVz) = (cy * sp * cr + sy * sr) * delta;
00244 transferFunction_(StateMemberX, StateMemberAx) = 0.5 * transferFunction_(StateMemberX, StateMemberVx) * delta;
00245 transferFunction_(StateMemberX, StateMemberAy) = 0.5 * transferFunction_(StateMemberX, StateMemberVy) * delta;
00246 transferFunction_(StateMemberX, StateMemberAz) = 0.5 * transferFunction_(StateMemberX, StateMemberVz) * delta;
00247 transferFunction_(StateMemberY, StateMemberVx) = sy * cp * delta;
00248 transferFunction_(StateMemberY, StateMemberVy) = (sy * sp * sr + cy * cr) * delta;
00249 transferFunction_(StateMemberY, StateMemberVz) = (sy * sp * cr - cy * sr) * delta;
00250 transferFunction_(StateMemberY, StateMemberAx) = 0.5 * transferFunction_(StateMemberY, StateMemberVx) * delta;
00251 transferFunction_(StateMemberY, StateMemberAy) = 0.5 * transferFunction_(StateMemberY, StateMemberVy) * delta;
00252 transferFunction_(StateMemberY, StateMemberAz) = 0.5 * transferFunction_(StateMemberY, StateMemberVz) * delta;
00253 transferFunction_(StateMemberZ, StateMemberVx) = -sp * delta;
00254 transferFunction_(StateMemberZ, StateMemberVy) = cp * sr * delta;
00255 transferFunction_(StateMemberZ, StateMemberVz) = cp * cr * delta;
00256 transferFunction_(StateMemberZ, StateMemberAx) = 0.5 * transferFunction_(StateMemberZ, StateMemberVx) * delta;
00257 transferFunction_(StateMemberZ, StateMemberAy) = 0.5 * transferFunction_(StateMemberZ, StateMemberVy) * delta;
00258 transferFunction_(StateMemberZ, StateMemberAz) = 0.5 * transferFunction_(StateMemberZ, StateMemberVz) * delta;
00259 transferFunction_(StateMemberRoll, StateMemberVroll) = transferFunction_(StateMemberX, StateMemberVx);
00260 transferFunction_(StateMemberRoll, StateMemberVpitch) = transferFunction_(StateMemberX, StateMemberVy);
00261 transferFunction_(StateMemberRoll, StateMemberVyaw) = transferFunction_(StateMemberX, StateMemberVz);
00262 transferFunction_(StateMemberPitch, StateMemberVroll) = transferFunction_(StateMemberY, StateMemberVx);
00263 transferFunction_(StateMemberPitch, StateMemberVpitch) = transferFunction_(StateMemberY, StateMemberVy);
00264 transferFunction_(StateMemberPitch, StateMemberVyaw) = transferFunction_(StateMemberY, StateMemberVz);
00265 transferFunction_(StateMemberYaw, StateMemberVroll) = transferFunction_(StateMemberZ, StateMemberVx);
00266 transferFunction_(StateMemberYaw, StateMemberVpitch) = transferFunction_(StateMemberZ, StateMemberVy);
00267 transferFunction_(StateMemberYaw, StateMemberVyaw) = transferFunction_(StateMemberZ, StateMemberVz);
00268 transferFunction_(StateMemberVx, StateMemberAx) = delta;
00269 transferFunction_(StateMemberVy, StateMemberAy) = delta;
00270 transferFunction_(StateMemberVz, StateMemberAz) = delta;
00271
00272
00273
00274 double xCoeff = 0.0;
00275 double yCoeff = 0.0;
00276 double zCoeff = 0.0;
00277 double oneHalfATSquared = 0.5 * delta * delta;
00278
00279 yCoeff = cy * sp * cr + sy * sr;
00280 zCoeff = -cy * sp * sr + sy * cr;
00281 double dFx_dR = (yCoeff * yVel + zCoeff * zVel) * delta +
00282 (yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00283 double dFR_dR = 1 + (yCoeff * pitchVel + zCoeff * yawVel) * delta;
00284
00285 xCoeff = -cy * sp;
00286 yCoeff = cy * cp * sr;
00287 zCoeff = cy * cp * cr;
00288 double dFx_dP = (xCoeff * xVel + yCoeff * yVel + zCoeff * zVel) * delta +
00289 (xCoeff * xAcc + yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00290 double dFR_dP = (xCoeff * rollVel + yCoeff * pitchVel + zCoeff * yawVel) * delta;
00291
00292 xCoeff = -sy * cp;
00293 yCoeff = -sy * sp * sr - cy * cr;
00294 zCoeff = -sy * sp * cr + cy * sr;
00295 double dFx_dY = (xCoeff * xVel + yCoeff * yVel + zCoeff * zVel) * delta +
00296 (xCoeff * xAcc + yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00297 double dFR_dY = (xCoeff * rollVel + yCoeff * pitchVel + zCoeff * yawVel) * delta;
00298
00299 yCoeff = sy * sp * cr - cy * sr;
00300 zCoeff = -sy * sp * sr - cy * cr;
00301 double dFy_dR = (yCoeff * yVel + zCoeff * zVel) * delta +
00302 (yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00303 double dFP_dR = (yCoeff * pitchVel + zCoeff * yawVel) * delta;
00304
00305 xCoeff = -sy * sp;
00306 yCoeff = sy * cp * sr;
00307 zCoeff = sy * cp * cr;
00308 double dFy_dP = (xCoeff * xVel + yCoeff * yVel + zCoeff * zVel) * delta +
00309 (xCoeff * xAcc + yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00310 double dFP_dP = 1 + (xCoeff * rollVel + yCoeff * pitchVel + zCoeff * yawVel) * delta;
00311
00312 xCoeff = cy * cp;
00313 yCoeff = cy * sp * sr - sy * cr;
00314 zCoeff = cy * sp * cr + sy * sr;
00315 double dFy_dY = (xCoeff * xVel + yCoeff * yVel + zCoeff * zVel) * delta +
00316 (xCoeff * xAcc + yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00317 double dFP_dY = (xCoeff * rollVel + yCoeff * pitchVel + zCoeff * yawVel) * delta;
00318
00319 yCoeff = cp * cr;
00320 zCoeff = -cp * sr;
00321 double dFz_dR = (yCoeff * yVel + zCoeff * zVel) * delta +
00322 (yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00323 double dFY_dR = (yCoeff * pitchVel + zCoeff * yawVel) * delta;
00324
00325 xCoeff = -cp;
00326 yCoeff = -sp * sr;
00327 zCoeff = -sp * cr;
00328 double dFz_dP = (xCoeff * xVel + yCoeff * yVel + zCoeff * zVel) * delta +
00329 (xCoeff * xAcc + yCoeff * yAcc + zCoeff * zAcc) * oneHalfATSquared;
00330 double dFY_dP = (xCoeff * rollVel + yCoeff * pitchVel + zCoeff * yawVel) * delta;
00331
00332
00333 transferFunctionJacobian_ = transferFunction_;
00334 transferFunctionJacobian_(StateMemberX, StateMemberRoll) = dFx_dR;
00335 transferFunctionJacobian_(StateMemberX, StateMemberPitch) = dFx_dP;
00336 transferFunctionJacobian_(StateMemberX, StateMemberYaw) = dFx_dY;
00337 transferFunctionJacobian_(StateMemberY, StateMemberRoll) = dFy_dR;
00338 transferFunctionJacobian_(StateMemberY, StateMemberPitch) = dFy_dP;
00339 transferFunctionJacobian_(StateMemberY, StateMemberYaw) = dFy_dY;
00340 transferFunctionJacobian_(StateMemberZ, StateMemberRoll) = dFz_dR;
00341 transferFunctionJacobian_(StateMemberZ, StateMemberPitch) = dFz_dP;
00342 transferFunctionJacobian_(StateMemberRoll, StateMemberRoll) = dFR_dR;
00343 transferFunctionJacobian_(StateMemberRoll, StateMemberPitch) = dFR_dP;
00344 transferFunctionJacobian_(StateMemberRoll, StateMemberYaw) = dFR_dY;
00345 transferFunctionJacobian_(StateMemberPitch, StateMemberRoll) = dFP_dR;
00346 transferFunctionJacobian_(StateMemberPitch, StateMemberPitch) = dFP_dP;
00347 transferFunctionJacobian_(StateMemberPitch, StateMemberYaw) = dFP_dY;
00348 transferFunctionJacobian_(StateMemberYaw, StateMemberRoll) = dFY_dR;
00349 transferFunctionJacobian_(StateMemberYaw, StateMemberPitch) = dFY_dP;
00350
00351 FB_DEBUG("Transfer function is:\n" << transferFunction_ <<
00352 "\nTransfer function Jacobian is:\n" << transferFunctionJacobian_ <<
00353 "\nProcess noise covariance is:\n" << processNoiseCovariance_ <<
00354 "\nCurrent state is:\n" << state_ << "\n");
00355
00356 Eigen::MatrixXd *processNoiseCovariance = &processNoiseCovariance_;
00357
00358 if (useDynamicProcessNoiseCovariance_)
00359 {
00360 computeDynamicProcessNoiseCovariance(state_, delta);
00361 processNoiseCovariance = &dynamicProcessNoiseCovariance_;
00362 }
00363
00364
00365 state_(StateMemberVroll) += controlAcceleration_(ControlMemberVroll) * delta;
00366 state_(StateMemberVpitch) += controlAcceleration_(ControlMemberVpitch) * delta;
00367 state_(StateMemberVyaw) += controlAcceleration_(ControlMemberVyaw) * delta;
00368
00369 state_(StateMemberAx) = (controlUpdateVector_[ControlMemberVx] ?
00370 controlAcceleration_(ControlMemberVx) : state_(StateMemberAx));
00371 state_(StateMemberAy) = (controlUpdateVector_[ControlMemberVy] ?
00372 controlAcceleration_(ControlMemberVy) : state_(StateMemberAy));
00373 state_(StateMemberAz) = (controlUpdateVector_[ControlMemberVz] ?
00374 controlAcceleration_(ControlMemberVz) : state_(StateMemberAz));
00375
00376
00377 state_ = transferFunction_ * state_;
00378
00379
00380 wrapStateAngles();
00381
00382 FB_DEBUG("Predicted state is:\n" << state_ <<
00383 "\nCurrent estimate error covariance is:\n" << estimateErrorCovariance_ << "\n");
00384
00385
00386 estimateErrorCovariance_ = (transferFunctionJacobian_ *
00387 estimateErrorCovariance_ *
00388 transferFunctionJacobian_.transpose());
00389 estimateErrorCovariance_.noalias() += delta * (*processNoiseCovariance);
00390
00391 FB_DEBUG("Predicted estimate error covariance is:\n" << estimateErrorCovariance_ <<
00392 "\n\n--------------------- /Ekf::predict ----------------------\n");
00393 }
00394
00395 }