ego_swarm真机代码中的递归最小二乘算法LinearControl::estimateThrustModel

如何计算期望加速度?

在这里插入图片描述
代码求期望加速度如下

 des_acc = des.a + Kv.asDiagonal() * (des.v - odom.v) + Kp.asDiagonal() * (des.p - odom.p);//真实状态的p和v均来自于视觉里程计
 des_acc += Eigen::Vector3d(0,0,param_.gra);
bool 
LinearControl::estimateThrustModel(
    const Eigen::Vector3d &est_a,
    const Parameter_t &param)
{
  ros::Time t_now = ros::Time::now();
  while (timed_thrust_.size() >= 1)
  {
    // Choose data before 35~45ms ago
    std::pair<ros::Time, double> t_t = timed_thrust_.front();
    double time_passed = (t_now - t_t.first).toSec();
    if (time_passed > 0.045) // 45ms
    {
      // printf("continue, time_passed=%f\n", time_passed);
      timed_thrust_.pop();
      continue;
    }
    if (time_passed < 0.035) // 35ms
    {
      // printf("skip, time_passed=%f\n", time_passed);
      return false;
    }

    /***********************************************************/
    /* Recursive least squares algorithm with vanishing memory */
    /***********************************************************/
    double thr = t_t.second;
    timed_thrust_.pop();
    
    /***********************************/
    /* Model: est_a(2) = thr1acc_ * thr */
    /***********************************/
    double gamma = 1 / (rho2_ + thr * P_ * thr);//1/(0.998+10^6*thr^2)
    double K = gamma * P_ * thr;//10^6*thr/(0.998+10^6*thr^2)
    thr2acc_ = thr2acc_ + K * (est_a(2) - thr * thr2acc_);//Ax=b  thr2acc_是x est_a(2)是b thr是A
    P_ = (1 - K * thr) * P_ / rho2_;
    //printf("%6.3f,%6.3f,%6.3f,%6.3f\n", thr2acc_, gamma, K, P_);
    //fflush(stdout);

    // debug_msg_.thr2acc = thr2acc_;
    return true;
  }
  return false;
}

得到了期望的加速度,如何映射为发给PX4的油门的百分比?

Ax=b,其中thr是A,thr2acc_=a/percentage是x, imu估计的est_a(2)是b , t h r = m ( g + a z ) thr= m(g + a_z) thr=m(g+az),但是这里的g和m不准并且存在很多扰动,因此使用参数估计的方法来估计这个映射 t h r 2 a c c thr2acc thr2acc,计算Ax=b, 其中thr2acc_是x, est_a(2)是b, thr是A
在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
x k = ( A k T A k + λ ) − 1 A k T b k x_{k}=\left(A_{k}^{T} A_{k}+\lambda \right)^{-1} A_{k}^{T} b_{k} xk=(AkTAk+λ)1AkTbk

我们不必按照最小二乘法的公式进行计算,首先注意到:

A k T A k = [ A k − 1 T α k T ] [ A k − 1 α k ] = A k − 1 T A k − 1 + α k T α k A_{k}{ }^{T} A_{k}=\left[\begin{array}{ll} A_{k-1}{ }^{T} & \alpha_{k}{ }^{T} \end{array}\right]\left[\begin{array}{c} A_{k-1} \\ \alpha_{k} \end{array}\right]=A_{k-1}{ }^{T} A_{k-1}+\alpha_{k}{ }^{T} \alpha_{k} AkTAk=[Ak1TαkT][Ak1αk]=Ak1TAk1+αkTαk

由上式可得到结论 1 : 求 A k T A k A_{k}{ }^{T} A_{k} AkTAk不需要做 A k T 与 A k A_{k}{ }^{T} 与 A_{k} AkTAk的矩阵乘法,只需要计算 α k T α k ( α k \alpha_{k}{ }^{T} \alpha_{k} ( \alpha_{k} αkTαk(αk的行数比 A k A_{k} Ak小 很多),再加上上一轮的结果 A k − 1 T A k − 1 A_{k-1}{ }^{T} A_{k-1} Ak1TAk1即可,这是对 A k T A k A_{k}{ }^{T} A_{k} AkTAk递归求解
第二,由于:

A k T b k = [ A k − 1 T α k T ] [ b k − 1 β k ] = A k − 1 T b k − 1 + α k T β k ⟶  yields  x k = ( A k T A k ) − 1 ( A k − 1 T b k − 1 + α k T β k ) \begin{array}{c} A_{k}{ }^{T} b_{k}=\left[\begin{array}{ll} A_{k-1}{ }^{T} & \alpha_{k}{ }^{T} \end{array}\right]\left[\begin{array}{c} b_{k-1} \\ \beta_{k} \end{array}\right]=A_{k-1}^{T} b_{k-1}+\alpha_{k}{ }^{T} \beta_{k} \\ \stackrel{\text { yields }}{\longrightarrow} x_{k}=\left(A_{k}{ }^{T} A_{k}\right)^{-1}\left(A_{k-1}{ }^{T} b_{k-1}+\alpha_{k}{ }^{T} \beta_{k}\right) \end{array} AkTbk=[Ak1TαkT][bk1βk]=Ak1Tbk1+αkTβk yields xk=(AkTAk)1(Ak1Tbk1+αkTβk)

由于k-1轮的结果满足: x k − 1 = ( A k − 1 T A k − 1 + λ ) − 1 A k − 1 T b k − 1 x_{k-1}=\left(A_{k-1}^{T} A_{k-1}+\lambda \right)^{-1} A_{k-1}^{T} b_{k-1} xk1=(Ak1TAk1+λ)1Ak1Tbk1 ,可得

( A k − 1 T A k − 1 + λ ) x k − 1 = A k − 1 T b k − 1 x k = ( A k T A k ) − 1 ( ( A k − 1 T A k − 1 + λ ) x k − 1 + α k T β k ) \begin{array}{c} (A_{k-1}{ }^{T} A_{k-1}+\lambda )x_{k-1}=A_{k-1}{ }^{T} b_{k-1} \\ x_{k}=\left(A_{k}{ }^{T} A_{k}\right)^{-1}\left((A_{k-1}^{T} A_{k-1}+\lambda ) x_{k-1}+\alpha_{k}{ }^{T} \beta_{k}\right) \end{array} (Ak1TAk1+λ)xk1=Ak1Tbk1xk=(AkTAk)1((Ak1TAk1+λ)xk1+αkTβk)

再把 A k T A k − α k T α k = A k − 1 T A k − 1 A_{k}^{T} A_{k}-\alpha_{k}{ }^{T} \alpha_{k}=A_{k-1}{ }^{T} A_{k-1} AkTAkαkTαk=Ak1TAk1 代入,得到 x k x_{k} xk 的递归求解:

x k = ( A k T A k + λ ) − 1 ( ( A k T A k − α k T α k + λ ) x k − 1 + α k T β k ) = ( I − ( A k T A k + λ ) − 1 α k T α k ) x k − 1 + ( A k T A k + λ ) − 1 α k T β k = x k − 1 + ( A k T A k + λ ) − 1 α k T ( β k − α k x k − 1 ) \begin{array}{c} x_{k}=\left(A_{k}{ }^{T} A_{k}+\lambda\right)^{-1}\left(\left(A_{k}{ }^{T} A_{k}-\alpha_{k}{ }^{T} \alpha_{k}+\lambda \right) x_{k-1}+\alpha_{k}{ }^{T} \beta_{k}\right) \\ =\left(I-\left(A_{k}{ }^{T} A_{k}+\lambda\right)^{-1} \alpha_{k}{ }^{T} \alpha_{k}\right) x_{k-1}+\left(A_{k}{ }^{T} A_{k}+\lambda\right)^{-1} \alpha_{k}{ }^{T} \beta_{k} \\ =x_{k-1}+\left(A_{k}{ }^{T} A_{k}+\lambda\right)^{-1} \alpha_{k}{ }^{T}\left(\beta_{k}-\alpha_{k} x_{k-1}\right) \end{array} xk=(AkTAk+λ)1((AkTAkαkTαk+λ)xk1+αkTβk)=(I(AkTAk+λ)1αkTαk)xk1+(AkTAk+λ)1αkTβk=xk1+(AkTAk+λ)1αkT(βkαkxk1)

注意, A k T A k A_{k}{ }^{T} A_{k} AkTAk 是递归求解得出的。每一轮都不涉及大规模的矩阵乘法
视频地址

Ax=b thr2acc_是x est_a(2)是b thr是A 给无人机一个推力,实际无人机对应一个什么样的加速度

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值