精华内容
下载资源
问答
  • 姿态篇:.非线性最小二乘与飞控传感器校准

    千次阅读 多人点赞 2018-08-12 01:12:52
    [深入浅出多旋翼飞控开发][高级篇][一][非线性最小二乘与飞控传感器校准] 作者:BlueSky QQ : 352707983 Github 前言 搞好了传感器,那意味着飞控已经完成了一半。 不用猜了,这句话正是鄙人说的:)...

    [深入浅出多旋翼飞控开发][姿态篇][四][非线性最小二乘与飞控传感器校准]

    作者:王伟韵
    QQ : 352707983
    Github

    前言

    搞好了传感器,那意味着飞控已经完成了一半。

    不用猜了,这句话正是鄙人说的:)。飞控的软硬件相关工作,绝大部分是围绕传感器展开的,它们的重要性,可能比你想象中的更大。

    硬件上,为了寻找性能优秀又符合产品实际需求的传感器如陀螺仪、加速度计、磁力计等,我们经常需要对来自不同厂商的样品进行实际性能测试与评估,即使实际上常用的型号就那么几个。为了提升飞控的性能,许多传感器问题我们还需要从硬件上进行解决,如降低电源噪声、隔离机身震动甚至还要为传感器提供恒定温度的工作环境。

    软件上,传感器校准、滤波与融合更是飞控程序中的重点,作为一名渣渣飞控工程师,浑浑噩噩工作数年,发现到头来大部分时间都是花在了传感器数据的处理上面,而且未来需要解决的问题,依然和传感器有关。

    本篇文章将结合实际飞控代码,以循序渐进的方式向读者讲述,如何使用最优化算法,解决飞控中最基本的传感器校准问题。

    一.传感器校准简述

    飞控上会使用各种各样的传感器,但对于多旋翼而言,最核心的是陀螺仪无疑。离开了陀螺仪,你甚至无法让一架多旋翼飞行器正常离地。其次则是加速度计,提供姿态观测,让飞行器可以得到一个较为客观的自身姿态观测信息,使得自稳飞行成为可能(能够自动保持自身姿态水平或一定角度)。除此之外,还有许许多多的的其它传感器,共同为飞控提供稳定飞行所需要的数据。

    其中,陀螺仪、加速度计和磁力计在飞行前(出厂前),均需要进行校准,才能达到正常工作所需要的性能。这便引出了一个问题:什么是传感器校准(标定)。

    答案很简单,所有事物都不是非黑即白,而是时刻处于一个不可预知的混沌状态。工厂按照特定的技术及工艺标准制造生产出一批传感器,这些传感器中所有个体的实际状态处于一个随机分布且均值为期望误差的序列上。简单来说,我们买来10个传感器,这些传感器均存在着特定误差且这10个传感器的误差各不一致。这就要求我们在实际使用这些传感器时,需要对每一个进行单独的校准,得到特定的校准参数,从而正常发挥传感器的性能。

    二.传感器误差类型

    接下来我们会面临下面两个问题:

    • 传感器有哪些误差
    • 如何测量传感器的误差

    如果从我们目前所使用的MEMS传感器的自身原理及制造工艺出发,分析传感器的误差来源,那可能需要长篇大论一番了,但那是超出作者能力的事情,也不是本篇文章所讨论的重点。

    根据实际经验,我们可以把传感器误差分为两大类:

    • 零偏误差(offset)
    • 刻度误差(scale)

    这两个是我们最容易理解,也是飞控中三大传感器均需要解决的误差类型,其中

    零偏误差: 指传感器在零激励下的输出误差。如陀螺仪在静止状态下,其角速度测量值理论为0,而实际上我们将飞控静止放置,陀螺仪依然会输出诸如1°/s之类的数据,这便是该陀螺仪传感器的零偏误差。

    刻度误差: 指传感器在一定激励输入下,其输出值与输入值的比值。比如我们将飞控放在一个以500°/s恒定角速度运动的转台上,而陀螺仪实际测量得到的角速度输出为495°/s,500 / 495 ≈ 1.01 便是该陀螺仪的刻度误差。

    从上述说明中我们可以看出,要正确测量传感器的误差,必须向传感器施加特定的激励信号,通过比较传感器的实际输出与激励输入量,计算得到传感器的误差。

    不同类型的传感器,所需要的激励信号源也不一样。如陀螺仪,要提供精确的激励信号,必须依靠专业的转台设备,这是个人及小公司均不具备的条件。好在对于大多数应用而言,陀螺仪的刻度误差在一个可以忍受的范围内,因此可被我们忽略。但对于精益求精的产品而言,要提升飞控性能到极致,刻度误差校准还是必须的,如某行业龙头,其飞控产品及整机产品,出厂前均需使用转台进行校准。缺少这些条件的个人爱好者,就只能暂时忽略这一项,只校准陀螺仪的零偏误差了,而大多数情况下,陀螺仪的零偏校准方式很简单,使用均值校准即可。因此,陀螺仪的校准不在本篇文章的讨论范围内。

    而本文中讨论的重点:加速度计磁力计,其信号激励源均来自于自然界,前者为地球的重力场,后者为地球的地磁场。

    三.传感器校准的简单方式

    1.加速度计

    假设有下面的情形,我们将飞控水平放置于水平面,测得加速度计z轴输出值g1=0.97gg_1= 0.97g,然后将飞控翻转并同样水平置于水平面,再测得加速度计z轴输出值g2=0.99gg_2 = -0.99g

    设z轴零偏误差为scalezscale_z,刻度误差为offsetzoffset_z,已知实际重力加速度为g,可以得到下列两条等式:

    • (g1offsetz)scalez=g(g_1 - offset_z) * scale_z = g
    • (g2offsetz)scalez=g(g_2 - offset_z) * scale_z = -g

    从而可以计算得到:

    • scalez1.02scale_z ≈ 1.02
    • offsetz=0.01goffset_z = -0.01g

    将这个方法推广到其余两轴,便能得到其校准参数scalex,scaley,offsetx,offsetyscale_x, scale_y,offset_x,offset_y

    缺陷:

    这个校准方法较为简单,但是其可靠性建立于一个理想条件下:加速度计在采集每一个方向的数据时,必须严格贴合水平面。否则只能得到一个不太准确的校准数值。实际应用中,我们很难保证校准时让飞控的每一个面保持水平,因为飞行器本身就是一个不规则物体,更多时候我们采集到的加速度数据为[0.05g,0.08g,0,95g][0.05g,0.08g,0,95g],而非[0,0,0.98g][0,0,0.98g]。或者说,这种方式的实现成本过高,难以实际应用。

    2.磁力计

    磁力计的简单校准方法原理等同于加速度计,即采集得到磁力计的三个轴的测量极值(最大与最小值),然后单独计算出每个轴的零偏误差与刻度误差。实现区别在于由于传感器激励信号源不一样,我们需要将磁力计在三维空间中进行各个方向的360°旋转,才能尽可能准确采集到极值。

    缺陷:

    实际应用中,磁力计校准对于飞行器使用者来说可能颇为频繁,而飞行器大小又各不一致,在很多情况下,我们是很难将飞行器翻来覆去地在各个方向上进行旋转,其可操作性非常低,且不方便普通用户的使用。如果旋转不全面,便无法采集到磁力计的真实极值,校准效果会大打折扣,甚至校准值严重偏离真实值。

    四.建立传感器误差模型

    对于程序员而言,解决实际问题的标准姿势如下:

    提出实际问题—>建立数学模型—>转换为程序—>求解

    而上一节中,我们描述了一个较为简单的用于加速度计与磁力计的校准方法,但是其前提条件较为苛刻,且校准效果不佳。于是,在这里我们提出了一个待解决的问题:找到一种使用便捷并且效果更好的校准算法。使用便捷,意味着约束条件越少越好,即在校准算法对传感器的采样点的数量和位置要求不能过于严格的情况下依然可以取得良好的校准效果。

    为了解决这个问题,得到传感器的误差参数,首先我们需要建立传感器的误差模型,并将其转换为数学方程,从而进行求解。

    以加速度计为例建立传感器误差模型

    设三轴加速度计的零偏误差与刻度误差分别为
    ox,oy,oz,sx,sy,szo_x ,o_y,o_z,s_x,s_y,s_z
    某一时刻该三轴加速度计的测量值分别为
    x,y,z x,y,z
    设当前时刻真实加速度向量为
    a=(ax,ay,az)a=(a_x,a_y,a_z)
    有以下关系:
    {ax=(xox)sxay=(yoy)syaz=(zoz)sz\begin{cases} a_x=(x-o_x)s_x\\ a_y=(y-o_y)s_y\\ a_z=(z-o_z)s_z\\ \end{cases}

    静止状态下,不管加速度计方向如何,其测量得到的加速度为重力加速度在各轴上的投影。而在地球上,重力加速度向量模值总是等于1g1g。假设ax,ay,aza_x,a_y,a_z单位为gg,便有以下等式成立:
    ax2+ay2+az2=1 a_x^2+a_y^2+a_z^2=1

    (xox)2sx2+(yoy)2sy2+(zoz)2sz2=1 (x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2=1

    至此,我们将加速度计的误差求解转化成为了一个数学问题: 已知 x,y,zx,y,z ,由上一条等式,求解出 ox,oy,oz,sx,sy,szo_x ,o_y,o_z,s_x,s_y,s_z

    五.求解误差方程

    上一节中,我们得到了传感器的误差方程:(xox)2sx2+(yoy)2sy2+(zoz)2sz2=1(x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2=1
    求解方程,对于经历过九年义务教育的我们来说是一件非常熟悉的事情,但这里存在着几个棘手的问题:

    • 该方程为六元二次方程,存在着6个未知数,如果只存在x,y,zx,y,z一组观测值,那将会有无穷多解。因此我们需要获得多组观测值xi,yi,zix_i,y_i,z_i,以构建方程组,从而求出唯一解,满足方程组里的每一个方程

    • 对于六元二次方程而言,使用六组独立的观测值,构建出包含6个方程的方程组,理论上可以求出方程组的唯一解。可是事实上,我们从传感器上获取到的观测值,都是带有一定随机噪声的,即观测值并非完全准确,从而导致求得的方程解未必精确

    • 使用计算机,我们可以轻松求出线性方程组的唯一精确解(如高斯消元法),而非线性方程问题,无论是从理论上还是计算方法上,都要复杂得多,大部分时候,对于无特定形式的非线性方程f(x)=0f(x)=0,是没有办法直接求得精确解的。

    为了解决上述问题,这里引进了一种数学方法:最小二乘法

    最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。

    对于现实中的绝大部分问题,我们都是无法求得精确解的,于是退而求其次,使用最小二乘法,我们可以得到其近似解。

    以传感器误差方程为例,设rr为某特定解下方程的误差,有
    r=1(xox)2sx2+(yoy)2sy2+(zoz)2sz2 r=1 - (x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2
    对于方程组而言,其误差的平方和函数为
    S=i=0nri2 S=\sum_{i=0}^{n} r_i^2
    所谓最小二乘法,便是求得一组特定解,使得SS最小(求函数极值),该条件下的解,便是方程组的最优近似解。

    对于线性方程组,可将误差函数S=0S=0转换为矩阵方程并进行求解即可,可惜现实问题大部分均为非线性问题,如本文中的传感器误差方程便是一个非线性方程,是无法变换为矩阵方程形式的,因此我们必须换一个思路来解决问题。

    假定该函数存在着局部最小值,那给定一个初始解,通过不断地迭代计算,从而找到误差平方和函数SS的局部最优解,这便是非线性最小二乘的基本思想。

    接下来,我们将讨论求解非线性最小二乘的具体算法。

    六.高斯牛顿法

    非线性最小二乘的求解算法有许多种,比较常见的有高斯牛顿法Levenberg-Marquardt法(LM),DogLeg法等。其中高斯牛顿法是最为经典的一种算法,APM飞控中便是使用了该算法,而LM法可以看成是高斯牛顿法的改良版,鲁棒性更好,在实际工程中应用更广泛,天穹飞控中则使用了LM法。而LM法的基础部分和高斯牛顿法完全一致,因此我们从高斯牛顿法开始讲起。

    由于算法原理细节部分较为复杂,需要用到大量公式,为了方便理解,将结合天穹飞控的实际代码进行讲解,高斯牛顿法和LM法的具体实现大部分是一样的,所以可以结合着来看。对于初学者来说彻底弄懂这一部分的原理对自己能力的提升帮助较大,因为此类最优化算法在飞控之外的非常多工程领域都是有着大量应用的,如计算机视觉和机器学习等。

    天穹飞控中的LM法完整代码

    要搞明白高斯牛顿法,我们要先来了解一下牛顿法

    1.牛顿法

    为求得函数最优解(极值点),牛顿法将问题转换为了求解f(x)=0f^{'} (x)=0

    首先给定一个初始解x0x_0,然后在x0x_0处进行一阶泰勒展开,得到
    f(x)f(x0)+(xx0)f(x0) f^{'} (x) \approx f^{'} (x_0) + (x-x_0)f^{''} (x_0)

    f(x0)+(xx0)f(x0)=0 f^{'} (x_0) + (x-x_0)f^{''} (x_0) = 0
    求解得到xx,相比于x0x_0,有f(x)<f(x0)f^{'} (x) < f^{'} (x_0),即新解xx比初始解更接近最优解。当然该解和最优解之间还有一定未知的距离,于是我们将新解作为初始解,重复这一步骤,经过不断地迭代计算,最终会收敛于f(x)=0f^{'} (x)=0

    牛顿法寻找极值点的动画示意图:
    牛顿法寻找极值点的动画示意图
    牛顿法的优点:

    • 收敛速度快,且可用于求解任意连续函数的最优化问题

    牛顿法的缺点:

    • 可以看出,牛顿法实际上等同于寻找函数的局部极值点,因此必须选取一个相对接近的初始点,否则会导致无法正确收敛

    • 需要求解目标函数的Hessian矩阵(二阶导)的逆矩阵,计算比较复杂

    2.高斯牛顿法

    高斯牛顿法是一种非线性最小二乘最优化方法。 其利用了目标函数的泰勒展开式把非线性函数的最小二乘化问题化为每次迭代的线性函数的最小二乘化问题

    高斯牛顿法实际上是牛顿法的简化形式,或者说改进版。多维情况下,当维数较大时,牛顿法中的Hessian矩阵求逆是往往行不通的,而高斯牛顿法对牛顿法中的Hessian矩阵进行了化简,使得问题可解。

    下面讲解高斯牛顿法的具体实现

    首先定义目标函数
    S(β)=i=05ri2(β) S(\beta)= \sum_{i=0}^{5}r _i^2( \beta)
    其中
    β=[β0,β1,β2,β3,β4,β5] \beta=[\beta_0,\beta_1,\beta_2,\beta_3,\beta_4,\beta_5]
    对应了加速度计误差方程中的六个变量,即
    [ox,oy,oz,sx,sy,sz] [o_x ,o_y,o_z,s_x,s_y,s_z]
    ri2(β)r _i^2( \beta)为传感器误差方程的误差平方和函数,即
    ri2(β)=(1(xiox)2sx2+(yioy)2sy2+(zioz)2sz2)2 r _i^2( \beta)=(1 - (x_i-o_x)^2s_x^2+(y_i-o_y)^2s_y^2+(z_i-o_z)^2s_z^2)^2

    设定一个初解βα\beta_\alpha,将目标函数S(β)S(\beta)βα\beta_\alpha附近泰勒展开,并舍弃高阶无穷项,近似有:
    S(β)=S(βα)+[S(βα)]T(ββα)+12(ββα)THS(βα)(ββα) S(\beta)= S(\beta_\alpha)+[\bigtriangledown S(\beta_\alpha)]^T(\beta-\beta_\alpha)+\frac{1}{2}(\beta-\beta_\alpha)^TH_S(\beta_\alpha)(\beta-\beta_\alpha)
    我们的目标是求出特定解,使S(β)S(\beta)有最小值。当函数有最小值(属于极小值)时,其一阶导为0,即
    S(β)=S(βα)+HS(βα)(ββα)=0 \bigtriangledown S(\beta)=\bigtriangledown S(\beta_\alpha)+H_S(\beta_\alpha)(\beta-\beta_\alpha)=0
    从而得到了高斯牛顿法的迭代公式
    β=βα[HS(βα)]1S(βα) \beta=\beta_\alpha - [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha)

    =[HS(βα)]1S(βα)\bigtriangleup = [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha)
    则有
    β=βα\beta=\beta_\alpha - \bigtriangleup
    这里便是迭代的关键了,给定初解后,每一次计算出\bigtriangleup,便得到了一个新解,然后使用新解不断进行迭代计算,使其逐渐逼近最优解。

    可能看到这里,你依然是一脸懵逼,但没有关系,你只需要知道,为了求解出函数的最优解,也就是传感器的最优校准参数,我们只需要不断去计算\bigtriangleup就行了。

    接下来便讲解如何计算\bigtriangleup (这里是整个算法实现的重点)

    \bigtriangleup中的HSH_SSS函数的Hessian矩阵(实际上是S函数的二阶导)
    HS=[2Sβ022Sβ0β52Sβ5β02Sβ52] H_S=\left[ \begin{matrix} \frac{\partial^2 S}{\partial \beta_0^2}&\cdots&\frac{\partial^2 S}{\partial \beta_0 \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial^2 S}{\partial \beta_5 \beta_0}&\cdots&\frac{\partial^2 S}{\partial \beta_5^2}&\\ \end{matrix} \right]
    前文提到,高斯牛顿法对牛顿法所做的改进,便是对Hessian矩阵进行了简化,其中有
    2Sβiβj=2k(rkβirkβj+rk2rkβiβj) \frac{\partial^2 S}{\partial \beta_i \beta_j}=2\sum_{k}(\frac{\partial r_k}{\partial \beta_i} \frac{\partial r_k}{\partial \beta_j} + r_k\frac{\partial^2 r_k}{\partial \beta_i \beta_j} )
    忽略二阶微分项,近似有
    2Sβiβj=2krkβirkβj \frac{\partial^2 S}{\partial \beta_i \beta_j}=2\sum_{k}\frac{\partial r_k}{\partial \beta_i} \frac{\partial r_k}{\partial \beta_j}
    于是有
    HS=2[r0β0r0β5r5β0r5β5]T[r0β0r0β5r5β0r5β5]=2JrTJr H_S=2 \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right]^T \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right] =2J_r^TJ_r
    这里的JrJ_r实际上就是函数rr(传感器误差方程的误差函数)的雅可比矩阵

    \bigtriangleup中的另外一个成员S\bigtriangledown SSS函数的一阶导的列向量,即
    S=[Sβ0Sβ5] \bigtriangledown S = \left[ \begin{matrix} \frac{\partial S}{\partial \beta_0}\\ \vdots\\ \frac{\partial S}{\partial \beta_5}\\ \end{matrix} \right]
    又有
    Sβi=k(rk2)βi=2krkβirk \frac{\partial S}{\partial \beta_i}=\sum_{k}\frac{\partial (r_k^2)}{\partial \beta_i} =2 \sum_{k}\frac{\partial r_k}{\partial \beta_i} r_k
    故$\bigtriangledown S $可改写为
    S=[Sβ0Sβ5]=2[r0β0r0β5r5β0r5β5]T[r0r5]=2JrTr \bigtriangledown S =\left[ \begin{matrix} \frac{\partial S}{\partial \beta_0}\\ \vdots\\ \frac{\partial S}{\partial \beta_5}\\ \end{matrix} \right] =2 \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right]^T \left[ \begin{matrix} r_0\\ \vdots\\ r_5\\ \end{matrix} \right] =2J_r^Tr

    于是我们终于得到了\bigtriangleup的最终形式
    =[Jr(βα)TJr(βα)]1Jr(βα)Tr(βα) \bigtriangleup = [J_r(\beta_\alpha)^T J_r(\beta_\alpha)]^{-1} J_r(\beta_\alpha)^T r(\beta_\alpha)
    或者换一个形式
    Jr(βα)TJr(βα)=Jr(βα)Tr(βα) J_r(\beta_\alpha)^T J_r(\beta_\alpha) \bigtriangleup =J_r(\beta_\alpha)^T r(\beta_\alpha)
    这就成了一个线性方程组,未知量是\bigtriangleup,使用高斯消元法,便能求解出\bigtriangleup的唯一解。

    至此,高斯牛顿法的所有计算步骤均结束,我们只需要不断进行迭代计算:求解出\bigtriangleup,从而得到新解β\beta,然后进行下一次迭代,不断逼近函数最优解。

    随着迭代次数的增加,最终会得到一个无限接近最优的解,但是考虑到时间成本和实际使用需求,在校准程序中我们设置一个阈值epseps,当方程解达到一定精度时(定义为迭代步长的平方2<eps\bigtriangleup ^2 < eps),停止迭代,结束此次校准。

    看完了,还是有点懵,怎么办?

    下面贴上高斯牛顿法的完整实现代码,结合代码进行理解和分析,事半功倍,如果还没看懂,那就多看几遍。

    主程序

    /**********************************************************************************************************
    *函 数 名: GaussNewton
    *功能说明: 高斯牛顿法求解传感器误差方程,得到校准参数
    *形    参: 传感器采样数据(6组) 零偏误差指针 比例误差指针 数据向量长度
    *返 回 值: 无
    **********************************************************************************************************/
    void GaussNewton(Vector3f_t inputData[6], Vector3f_t* offset, Vector3f_t* scale, float length)
    {
        uint32_t cnt    = 0;
        double   eps    = 0.000000001;
        double   change = 100.0;
        float    data[3];  
        float    beta[6];      //方程解
        float    delta[6];     //迭代步长
        float    JtR[6];       //梯度矩阵
        float    JtJ[6][6];    //Hessian矩阵
       
        //设定方程解初值
        beta[0] = beta[1] = beta[2] = 0;
        beta[3] = beta[4] = beta[5] = 1 / length;
    
        //开始迭代,当迭代步长小于eps时结束计算,得到方程近似最优解
        while(change > eps) 
        {
            //矩阵初始化
            ResetMatrices(JtR, JtJ);
    
            //计算误差方程函数的梯度JtR和Hessian矩阵JtJ
            for(uint8_t i=0; i<6; i++) 
            {
                data[0] = inputData[i].x;
                data[1] = inputData[i].y;
                data[2] = inputData[i].z;
                UpdateMatrices(JtR, JtJ, beta, data);
            }
    
            //高斯消元法求解方程:JtJ * delta = JtR,得到delta
            GaussEliminateSolveDelta(JtR, JtJ, delta);
    
            //计算迭代步长
            change = delta[0]*delta[0] +
                     delta[0]*delta[0] +
                     delta[1]*delta[1] +
                     delta[2]*delta[2] +
                     delta[3]*delta[3] / (beta[3]*beta[3]) +
                     delta[4]*delta[4] / (beta[4]*beta[4]) +
                     delta[5]*delta[5] / (beta[5]*beta[5]);
    
            //更新方程解
            for(uint8_t i=0; i<6; i++) 
            {
                beta[i] -= delta[i];
            }
                
            //限制迭代次数
            if(cnt++ > 100)
                break;
        }
    
        //更新校准参数
        scale->x  = beta[3] * length;
        scale->y  = beta[4] * length;
        scale->z  = beta[5] * length;
        offset->x = beta[0];
        offset->y = beta[1];
        offset->z = beta[2];
    }
    

    主程序中调用了三个子函数,其中:

    • ResetMatrices,用于初始化矩阵变量
    static void ResetMatrices(float JtR[6], float JtJ[6][6])
    {
        int16_t j,k;
        for(j=0; j<6; j++) 
        {
            JtR[j] = 0.0f;
            for(k=0; k<6; k++) 
            {
                JtJ[j][k] = 0.0f;
            }
        }
    }
    
    • UpdateMatrices,用于计算求解\bigtriangleup所用到的Jr(βα)TJr(βα)J_r(\beta_\alpha)^T J_r(\beta_\alpha)(Hessian矩阵)和Jr(βα)Tr(βα)J_r(\beta_\alpha)^T r(\beta_\alpha)这两个矩阵
    static void UpdateMatrices(float JtR[6], float JtJ[6][6], float beta[6], float data[3])
    {
        int16_t j, k;
        float dx, b;
        float residual = 1.0;
        float jacobian[6];
        
        for(j=0; j<3; j++) 
        { 
            b = beta[3+j];
            dx = (float)data[j] - beta[j];
            //计算残差 (传感器误差方程的误差)
            residual -= b*b*dx*dx;
            
            //计算雅可比矩阵
            jacobian[j] = 2.0f*b*b*dx;
            jacobian[3+j] = -2.0f*b*dx*dx;
        }
        
        for(j=0; j<6; j++) 
        {
            //计算函数梯度
            JtR[j] += jacobian[j]*residual;
            
            for(k=0; k<6; k++) 
            {
                //计算Hessian矩阵(简化形式,省略二阶偏导),即雅可比矩阵与其转置的乘积
                JtJ[j][k] += jacobian[j]*jacobian[k];
            }
        }
    }
    
    • GaussEliminateSolveDelta,使用高斯消元法求解\bigtriangleup
    static void GaussEliminateSolveDelta(float JtR[6], float JtJ[6][6], float delta[6])
    {
        int16_t i,j,k;
        float mu;
       
        //逐次消元,将线性方程组转换为上三角方程组
        for(i=0; i<6; i++)
        {
            //若JtJ[i][i]不为0,将该列在JtJ[i][i]以下的元素消为0
            for(j=i+1; j<6; j++)
            {
                mu = JtJ[i][j] / JtJ[i][i];
                if(mu != 0.0f) 
                {
                    JtR[j] -= mu * JtR[i];
                    for(k=j; k<6; k++)
                    {
                        JtJ[k][j] -= mu * JtJ[k][i];
                    }
                }
            }
        }
    
        //回代得到方程组的解
        for(i=5; i>=0; i--)
        {
            JtR[i] /= JtJ[i][i];
            JtJ[i][i] = 1.0f;
            
            for(j=0; j<i; j++)
            {
                mu = JtJ[i][j];
                JtR[j] -= mu * JtR[i];
                JtJ[i][j] = 0.0f;
            }
        }
    
        for(i=0; i<6; i++)
        {
            delta[i] = JtR[i];
        }
    }
    

    事实上,在实际应用中使用高斯牛顿法,会出现许许多多的问题,因此下面介绍一种基于高斯牛顿改良优化而来的算法。

    七.Levenberg-Marquardt法

    实际应用中,高斯牛顿法有着如下缺点:

    • 初始点选取较为严苛,若初始点距离极小值点过远,迭代步长过大会导致迭代下一代的函数值不一定小于上一代的函数值
    • 高斯牛顿法中的Hessian矩阵在某些情况下可能是非正定的(不可逆)
    • 残差rr过大时可能导致高斯牛顿法无法收敛

    而Levenberg-Marquardt法(LM法,又称阻尼最小二乘法)改善了高斯牛顿法的不足之处,并结合了高斯牛顿法和梯度下降法的优点,而实际上所做的改动非常小。

    在高斯牛顿法中,迭代步长\bigtriangleup被定义为
    =[HS(βα)]1S(βα) \bigtriangleup = [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha)
    而在LM法中,加入了阻尼因子μ\mu,有
    =[HS(βα)+μI]1S(βα) \bigtriangleup = [H_S(\beta_\alpha) + \mu I]^{-1}\bigtriangledown S(\beta_\alpha)
    从而消除了Hessian矩阵非正定的情况,更重要的是迭代过程中阻尼因子μ\mu是可变的,用于调整下降方向,当:

    • μ\mu趋近于0时,算法退化为高斯牛顿法
    • μ\mu很大时,效果相当于梯度下降法

    使用LM法时, 需要加入一些μ\mu的调整策略,使得距离解较远的时候,使用梯度下降法,距离极小值较近的时候,使μ\mu减小,这样就近似高斯牛顿法,能得到比梯度下降法更快的收敛速度。

    算法的具体实现这里就不贴出来了,可以到这里查看:天穹飞控中的LM法

    八.传感器数据采集

    前几节介绍了如何使用非线性最小二乘法校准传感器,下面简单说明校准时采集所需传感器数据的具体步骤。

    1.加速度计

    在天穹飞控中,加速度校准使用标准的六面采集法,需要将飞控按照上、下、前、后、左、右六个方向静止放置一段时间(2s左右),顺序随意且飞控方向大致准确即可,飞控自动检测当前加速度计方向,并采集多组数据取平均值,最后得到6组方向各异的数据,导入LM校准函数,便能求出加速度计的6个校准参数。

    加速度校准采集代码

    天穹飞控兼容mavlink,可以使用QGroundControl地面站完成加速度计的校准步骤,如下图:
    使用QGC校准加速度计

    点击左侧的Accelerometer按钮,确认OK,飞控便会进入加速度校准流程,同时QGC会显示出当前校准进度以及图片提示,当对应方向采集数据完毕时,对应图片外框会变绿,红色的则是还没有采集数据的方向。

    2.磁力计

    磁力计的数据采集方式采用比较通用的两圈旋转方式,即:将飞控水平旋转一圈,然后再竖直旋转一圈。在这个过程中,记录下各轴的最大和最小值数据。

    磁力计校准采集代码

    同样可以使用QGC进行磁力计校准,点击Compass并确认OK。开始按照上述步骤旋转飞控,完成校准。

    然而相对加速度计校准而言,磁力计校准对算法的要求更高。由于加速度计校准环境通常比较单一,且静止状态下加速度计数据基本不受外界干扰,所以采集到的数据都是比较准确的,使得算法很容易收敛,实测大概只需4-5次的迭代运算,就可以达到所需精度。

    而磁力计数据采集在实际应用过程中,可能存在以下问题:

    • 校准环境比较复杂,而且校准频率非常频繁。用户可能是在室内校准,也可能是在城市的建筑间校准,或者是在充满了未知磁场干扰的环境下校准。这便导致了在校准期间,可能由于各种未知的外部干扰,使得磁力计数据采集出现“野值”

    • 为了简化校准步骤,我们采用了较为简单的两圈旋转方式,但由于地磁场的分布方向原因,实际上这样我们采集到的磁力计数据并不是在整个球面上分布均匀的,如下图所示(图片来自正在开发中的天穹地面站):
      磁力计数据分布示意图
      为了解决问题1,我们对采集到的数据做了一定的滤波平滑处理,同时实时计算数据的一个平均模值,当新值的大小超过平均模值一定比例时,抛弃这个值,避免在构建误差方程组时,引入包含较大干扰量的观测值,导致算法无法收敛或最终计算出来的传感器校准参数存在误差。

    问题2考验了算法的性能,即在样本分布不均匀的情况下是否依然可以收敛,经过测试和验证,天穹飞控中使用的LM法是能够经受得住这个考验的。

    九.总结

    本篇文章介绍了飞控中最基本的传感器误差问题,并提出了有效的解决方法。其中所使用到的非线性最小二乘算法,在现代计算机工程领域有着非常广泛的应用,所以即使对飞控技术兴趣不大的童鞋,本文也是值得一看的。

    而对于我们的飞控而言,完成了传感器基本误差的校准,只是相当于迈出了第一步,后面还有无数难题待解决,再接下来的教程中,会陆续针对不同的问题寻求最优的解决方案,当然由于本人水平有限,只能起到抛砖引玉的作用,希望有更多的有识之士,加入这个开源飞控项目,一起完善,共同进步。

    最后再留一下传送门:

    天穹飞控
    飞控交流论坛
    飞控交流Q群:472648354

    展开全文
  • 二、误差函数—最小二乘法 三、引出案例函数曲线 、目标函数 五、优化目标函数 六、优化目标函数—梯度下降法 七、优化目标函数—求解线性方程组 八、python编程实现拟合曲线函数 九、结果分析 一、函数...

    原文地址:https://www.jianshu.com/p/af0a4f71c05a

    目录

    一、函数的近似表示—高次多项式

    二、误差函数—最小二乘法

    三、引出案例函数曲线

    四、目标函数

    五、优化目标函数

    六、优化目标函数—梯度下降法

    七、优化目标函数—求解线性方程组

    八、python编程实现拟合曲线函数

    九、结果分析


     

    一、函数的近似表示—高次多项式

    为了研究一些复杂的函数,我们希望对函数自变量进行有限的加、减、乘法三种算数运算,便可以求出原函数值,因此我们通常使用高次多项式来近似表示函数

                                                   高次多项式近似表示f(x)

    可以看到,我们可以用n阶多项式的系数线性组合来近似表示f(x)

    特别说明,由泰勒展开式可知,当pn(x)的各阶导数和f(x)的各阶导数都相等,则f(x)和pn(x)的误差只是(x-x0)^n的高阶无穷小

                                                              泰勒公式

    二、误差函数—最小二乘法

    我们用f(x)的真实函数值减去多项式函数的结果的平方,来表示f(x)和多项式函数的误差关系,即

                                                      最小二乘法表示误差

    我们令x0=0,则最小二乘法表示的误差为

                                                   最小二乘法表示误差

    三、引出案例函数曲线

    有了上述数学知识,下面我们用一个案例来介绍最小二乘法拟合非线性函数的算法步骤

    案例:求一个多项式来拟合下列函数的函数值

                                                                    案例函数

    其中我们加入了噪点数据noise,使得函数在定义域内随机的震荡

    我们画出f(x)在[-1, 1]之间的图像,即

                                                                                                  案例函数图像

    案例问题即为:已知上述N个数据点,求其函数f(x)的表达式?

    显然,f(x)也就是机器学习要学习的目标

    四、目标函数

    下面我们开始推导机器学习的过程

    机器学习的目标为:

    (1) 学习一个f(x)多项式,可以拟合真实数据的变化趋势

    (2)f(x)的目标:使每一个真实数值到f(x)的拟合数值的距离之和最小

    翻译为数学语言,即

                             学习到的f(x)

                                                   目标函数

    五、优化目标函数

    为了使目标函数最优,我们对每个系数求其偏导数为0,即

                                                     优化目标函数

    至此,我们需要根据上述k的等式,求出所有的系数ak,有两种方法可以求解

    (1)梯度下降法

    (2)求解线性方程

    六、优化目标函数—梯度下降法

    采用梯度下降法,可以避免我们用数学方法直接一步来求解上述k个方程的最优解,而是采用迭代逼近最优解的思路,其具体步骤为:

    (1)初始化k个系数值,开始迭代

    (2)每次迭代,分别求出各个系数ak对应的梯度值

    (3)用梯度值和学习率来更新每一个系数ak

    (4)保证每次更新完所有系数,对应的损失函数值都在减少(说明梯度一直在下降)

    翻译为数学语言为:

    计算ak对应的梯度值

                                                     计算ak的梯度值

    更新ak

                                  更新ak

    七、优化目标函数—求解线性方程组

    求解线性方程组让我们直接一步求出偏导数最优解,其具体步骤为:

    (1)将最优化偏导数的方程组写为矩阵乘法形式:XA=Y

    (2)利用数学消元算法来求解XA=Y

    我们对k个偏导数=0的等式进行化简处理,即

                                                        k个偏导数等式

    下面我们将其写为矩阵相乘的形式,将所有只含有x的系数写为第一个矩阵X,即

                                                 含有x的系数矩阵

    将只含有待求解的多项式系数ak写为第二个矩阵A,即

          含有a的系数矩阵

    最后将含有y的系数写为第三个矩阵Y,即

                  含有y的系数矩阵

    此时,我们就可以用矩阵的乘法来表示最初的k个偏导数为0的等式,即

             k个偏导数等式的矩阵形式

    注意:其中X是k*k的矩阵,A是k*1的矩阵,Y是k*1的矩阵,符合矩阵的乘法规则

    从矩阵表达式可以看出,X和Y矩阵是已知的,A矩阵只包含待求解的f(x)的所有多项式系数,则问题即转化为线性方程组:XA=Y

    求解线性方程组,大概有如下算法

    (1)高斯消元法(全主元、列主消元)

    (2)LU三角分解法

    (3)雅克比迭代法

    (4)逐次超松弛(SOR)迭代法

    (5)高斯-赛德尔迭代法

    这里我们采用高斯消元法来求解A,具体算法步骤较为复杂,我们在单独的章节介绍高斯消元算法

    八、python编程实现拟合曲线函数

    生成带有噪点的待拟合的数据集合

                                                                           生成带有噪点的待拟合的数据集合

    计算最小二乘法当前的误差

                                                                                计算最小二乘法当前的误差

    迭代解法:最小二乘法+梯度下降法

                                                                           迭代解法:最小二乘法+梯度下降法

    数学解法:最小二乘法+求解线性方程组

                                                                           数学解法:最小二乘法+求解线性方程组

    九、结果分析

    我们来拟合N=10次幂的多项式,分别用梯度下降法和求解线性方程组的算法来比较拟合结果

    (1)使用梯度下降法来优化目标函数,其拟合结果为:

                                                                                         梯度下降法拟合结果

    其误差变化为

                                                                                               梯度下降法优化的误差

    (2)使用求解线性方程组来优化目标函数,其拟合结果为:

                                                                                        求解线性方程组拟合结果

    其误差变化为

                               求解线性方程组优化的误差

    比较上述两种解法,以及多项式拟合的思想,我们可以总结出:

    • 利用函数的近似多项式表示,和最小二乘法来表示目标函数,我们可以高效的拟合非线性函数
    • 梯度下降法可以一步步迭代逼近目标函数的极值,而不用直接求出偏导数最优解
    • 求解方程组,使用数学消元的算法,直接一步优化完成了目标函数的偏导数最优解
    • 从结果上比较,求解方程组的效率和拟合结果要比梯度下降法好

    案例代码见:最小二乘法—多项式拟合非线性函数

    示例论文预测:基于非线性最小二乘曲线拟合法的电子商务预测模型

    展开全文
  • 例题:下列给定程序中,函数fun的功能是:统计一无符合整数中各位数字值为0的个数,通过形参传回主函数,并把该整数中各位上最大的数字值作为函数值返回。 例如。若输入无符号整数10080,则数字值为0的个数为3,...

    例题:下列给定程序中,函数fun的功能是:统计一个无符合整数中各位数字值为0的个数,通过形参传回主函数,并把该整数中各位上最大的数字值作为函数值返回。

    例如。若输入无符号整数10080,则数字值为0的个数为3,各位上数字值最大的是8。
    注意:不要改动main函数,不能增行或删行,也不能更改程序的结构。

    代码如下:

    #include<stdio.h>
    int fun(unsigned m,int*z)
    {
    	int n=0,max=0,t;
    	do
    	{
    		t=m%10;
    		if(t==0)
    			n++;
    		if(max<t)
    			max=t;
    		m=m/10;
    	}while(m);
    	*z=n;
    	return max;
    }
    main()
    {
    	unsigned m;
    	int z,max;
    	printf("\nInput m(unsigned):");
    	scanf("%d",&m);
    	max=fun(m,&z);
    	printf("\nThe result:max=%d z=%d\n",max,z);
    }
    

    输出运行窗口如下:
    在这里插入图片描述
    本周其他练习

    C语言程序设计专栏

    C语言编程>第二十三周 ① 下列给定程序中,函数fun的功能是:求n!(n<20),所求阶乘的值作为函数值返回。例如,若n=5,则应输出120。

    C语言编程>第二十三周 ② 请补充fun函数,该函数的功能是:交换数组a中最大和最小两个元素的位置,结果重新保存在原数组中,其它元素位置不变。注意数组a中没有相同元素。

    C语言编程>第二十三周 ③ 下列给定程序中,函数fun的功能是:利用插入排序法对字符串中的字符按从小到大的顺序进行排序。插入法的基本算法是:先对字符串中的头两个元素进行排序;然后把第三字符插入到前两个字符中,插入后前三个字符依然有序;再把第四个字符插入到前三个字符中……待排序的字符串已在主函数中赋予。

    C语言编程>第二十三周 ④ 请补充fun 函数,该函数的功能是:删除字符数组中比指定字符小的字符,指定字符从键盘输入,结果仍保存在原数组中。

    C语言编程>第二十三周 ⑤ 请补充main函数,该函数的功能是:求1~100(不包括100)以内所有素数的平均值。

    C语言编程>第二十三周 ⑥ 下列给定程序中函数fun的功能是:删除字符串s中的所有空白字符(包括Tab字符、回车符及换行符)。输入字符串时用 “#”结束输入。

    C语言编程>第二十三周 ⑦ 请补充main函数,该函数的功能是:求n!。

    C语言编程>第二十三周 ⑧ 下列给定程序中,函数fun的功能是:统计一个无符合整数中各位数字值为0的个数,通过形参传回主函数,并把该整数中各位上最大的数字值作为函数值返回。

    越努力越幸运!
    加油,奥力给!!!

    展开全文
  • python最小二乘法拟合实例

    千次阅读 2019-11-05 14:50:04
    最小二乘法主要解决的是插值中出现的三问题: ① 高多项式会龙格现象; ② 定函数关系不合理(数据误差); ⑶评估逼近效果: 假设“逼近”规律的近似函数为y=f(x), 即有yi﹡= f(xi)(i=1,2,…,m).它与观测值yi之差...

    1)曲线拟合的最小二乘法
    ⑴两种逼近概念:
    ①插值:在节点处函数值相同;
    ②拟合:在数值上误差最小;
    ⑵最小二乘法主要解决的是插值中出现的三个问题:
    ① 高多项式会龙格现象;
    ② 定函数关系不合理(数据误差);
    ⑶评估逼近效果:
    假设“逼近”规律的近似函数为y=f(x), 即有yi﹡= f(xi)(i=1,2,…,m).它与观测值yi之差 δi = yi﹡ - yi = f(xi) – yi,i = 1, 2, …,m。这称为残差。残差大大小可以作为衡量近似函数好坏的标准。按照使残差的平方和∑δi²最小的规则求得近似函数y=f(x)的方法称为最佳平方逼近,也称为曲线拟合的最小二乘法。
    举例:
    1)假定某天的气温变化记录如下,使用最小二乘法找出这一天的气温变化规律.
    t/h 0 1 2 3 4 5 6 7 8 9 10 11 12
    T/℃ 15 14 14 14 14 15 16 18 20 22 23 25 28
    t/h 13 14 15 16 17 18 19 20 21 22 23 24
    T/℃ 31 32 31 29 27 25 24 22 20 18 17 16
    要求:编程实现考虑下列类型的拟合函数,计算误差平方和,并作图比较效果.
    (1)二次、三次、四次多项式拟合函数;
    (2)形如 在这里插入图片描述
    计算各个多项式的二范数的平方、残差的方差和残差绝对值的均值如表1-1所示.
    在这里插入图片描述
    残差的方差的意义在于评估各个偏差的离散程度,避免只是部分拟合的效果比较好,而部分偏差又非常大。残差绝对值的均值的意义在于评估总体水平上残差的大小。通过比较我们可以发现总体上是四次多项式的拟合效果比较好。

    接下我们分别作出各个多项式在样本数据上的拟合效果,结果如图1-1所示。
    图1-1
    通过图1-1我们也可以看出四次项的拟合效果比较好。
    但通过增加多项式的最高次数之后,测试发现随着次数增加拟合的效果越好,出现了过拟合的现象,所以目前的评价体系还有问题。参考机器学习中的常见防止过拟合的方法是加入L1或者L2正则化项。
    《目前还在学习阶段有些问题还请谅解》

    代码如下:

    import matplotlib.pyplot as plt
    from matplotlib.font_manager import FontProperties
    import numpy as np
    from scipy.optimize import leastsq
    
    x = np.arange(0,25,1,dtype='float')
    y = np.array([15,14,14,14,14,15,16,18,20,22,23,25,
         28,31,32,31,29,27,25,24,22,20,18,17,16], dtype='float')
    
    # param:最小化的目标| *_fun:拟合函数
    # 二次
    param0 = [0, 0, 0]
    def quadratic_fun(a, x):
        k1, k2, b = a
        return k1 * x **2 + k2 * x + b
    # 三次
    param1 = [0, 0, 0, 0]
    def cubic_fun(a, x):
        k1, k2, k3, b = a
        return k1 * x ** 3 + k2 * x**2 + k3 * x + b
    # 四次
    param2 = [0, 0, 0, 0, 0]
    def fpower_fun(a, x):
        k1, k2, k3, k4, b = a
        return k1 * x ** 4 + k2 * x**3 + k3 * x**2 + k4 * x + b
    # 自定义函数
    param3 = [0, 0, 0]
    def myfuns(s, x):
        a, b, c = s
        return a * np.exp(-b*(x-c))
    param4 = [0, 0, 0, 0, 0, 0, 0, 0]
    def testfun(a, x):
        k1, k2, k3, k4, k5, k6, k7, b = a
        return  k7 * x ** 7 + k6 * x ** 6 + k1 * x ** 5 + k2 * x**4 + k3 * x**3 + k4 * x ** 2 + k5 * x + b
    # 偏差
    def dist(a, fun, x, y):
        return fun(a, x) - y
    # 作图
    myfont = FontProperties(fname="C:\Windows\Fonts\msyh.ttc")
    plt.figure()
    plt.title(u'某天的气温变化', fontproperties=myfont)
    plt.xlabel(u't/h')
    plt.ylabel(u'T/摄氏度', fontproperties=myfont)
    # 坐标轴的范围xmin, xmax, ymin, ymax
    plt.axis([0, 24, 10, 35])
    plt.grid(True)
    plt.plot(x, y, 'k.')
    
    
    funs = [quadratic_fun, cubic_fun, fpower_fun, myfuns, testfun]
    params = [param0, param1, param2, param3, param4]
    colors = ['blue', 'red', 'black', 'green', 'yellow']
    fun_name = ['quadratic_fun', 'cubic_fun', '4power_fun', 'a*exp(-b*(t-c))', 'testfun']
    
    for i, (func, param, color, name) in enumerate(zip(funs, params, colors, fun_name)):
        var = leastsq(dist, param, args=(func, x, y))
        plt.plot(x, func(var[0], x), color)
        print('[%s] 二范数: %.4f, abs(bias): %.4f, bias-std: %.4f' % (name,
              ((y-func(var[0], x))**2).sum(),
              (abs(y-func(var[0], x))).mean(), 
              (y-func(var[0], x)).std())
        )
    plt.legend(['sample data', 'quadratic_fun', 'cubic_fun', '4power_fun', 'a*exp(-b*(t-c))', 'fivpower_fun'], loc='upper left') 
    plt.show()
    

    参考书目:
    [M] 郑继明,朱伟,刘勇等 数值分析

    展开全文
  • 例题:给定程序中,函数fun 的功能是:将形参student所指结构体数组中年龄最小者的数据作为函数值返回,并在main函数中输出。 注意:请勿改动主函数main与其它函数中的任何内容,仅在fun函数的横线上填写所需的若干...
  • 例题:下列给定程序中,函数fun的功能是:求n!(n<20),所求阶乘的值作为函数值返回。例如,若n=5,则应输出120。 请修改程序中的错误,得出正确的结果。 注意:不要改动main函数,不能增行或删行,也不能更改程序...
  • 寻找最小的k

    千次阅读 2011-11-13 20:27:18
    因为后续,我们还专门写了程序进行测试,即针对1000w的数据寻找其中最小的k数的问题,采取两种实现,一是采取常规的建立k元素最大堆后通过比较寻找最小的k数的方案,一是采取建立n元素的最小堆,然后取其...
  • 下面,我试图用最清晰易懂,最易令人理解的思维或方式阐述有关寻找最小的k数这问题(这几天一直在想,除了计数排序外,这题到底还有没有其它的O(n)的算法? )。希望,有任何问题,欢迎不吝指正。谢谢。 寻找...
  • 《Python》SciPy——数值计算库

    千次阅读 2018-07-31 18:47:55
    例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。作为入门介绍,让我们看看如何用SciPy进行插值处理、信号滤波以及用C语言加速计算。 一、最小二乘拟合     假设有一组实验数据(x[i], ...
  • 数值分析实习作业 院系 学号 姓名 数据科学与计算机学院 16339053 张志宏 1.对函数f(x)f(x)f(x)进行插值 f(x)=11+x2f(x)=11+x2f(x)=\frac{1}{1+x^{2}} (1)令插值节点为等距节点−5,−4,−3,...
  • 最小费用流及其求法

    万次阅读 多人点赞 2019-05-04 08:33:48
    【3】树:基本概念与最小生成树 【4】匹配问题: 匈牙利算法 、最优指派、相等子图 【5】Euler 图和 Hamilton 图 【6】计划评审方法和关键路线法【统筹方法】:广泛地用于系统分析和项 目管理 【7】最小费用流...
  • 例题:下列给定程序中,函数fun的功能是:实现两整数的交换。 例如,给x和y分别输入60和65,输出为:x=65 y=60。 注意:不要改动main函数,不能增行或删行,也不能更改程序的结构。 代码如下: #include<stdio....
  • 最小二乘支持向量机LSSVM

    千次阅读 2020-03-13 21:54:39
    Suykens等人提出的最小二乘支持向量机(Least Squares Support Vector Machines,LS-SVM)从机器学习损失函数着手,在其优化问题的目标函数中使用二范数,并利用等式约束条件代替SVM标准算法中的不等式约束条件,使得...
  • 例题:请补充main函数,该函数的功能是:计算四个学生各科的平均分。 例如,当score[N][M]={{83,65,63},{89,93,95},{90,63,80},{56,75,77}}时,则平均分为:79.5,74.0,78.8 注意:仅在横线上填写所需的若干表达式...
  • 注:N表示数字型,C表示字符型,D表示日期型,[]表示内中参数可被忽略,fmt表示格式...数值型函数(Number Functions) 数值型函数输入数字型参数并返回数值型的值。多数该类函数的返回值支持38位小数点,诸如:COS, COS
  • 偏微分方程的数值解系列博文: 偏微分方程的数值解(一):定解问题 &...偏微分方程的数值解(): 化工应用————扩散系统之浓度分布 偏微分方程的数值解(五): 二维状态空间的偏微分方程的 MATLA...
  • 在数学中,辗转相除法,又称欧几里得算法,是求最大公约数的算法...辗转相除法基于如下原理:两整数的最大公约数等于其中较小的数和两数的相除余数的最大公约数。例如,252和105的最大公约数是21(252 = 21 × 12;10
  • 文档介绍的还不错,我估计任何一本数值分析教材上讲的都非常清楚。 推导就不再写了,我主要参考下面两页PPT,公式和例子讲的比较清楚。 公式: 例子: matlab代码如下: clear all; close...
  • 深度学习之数学基础(数值计算)

    万次阅读 2017-10-26 23:27:53
    信息论是应用数学的一分支,主要研究的是对一信号能够提供信息的多少进行量化。如果说概率使我们能够做出不确定性的陈述以及在不确定性存在的情况下进行推理,那信息论就是使我们能够量化概率分布中不确定性的...
  • 令,为这些随机变量的均值,对于任意有: 在机器学习中,引理2称为Chernoff边界(Chernoff bound),它说明:假设我们用随机变量的均值去估计参数,估计的参数和实际参数的差超过一特定数值的概率有一确定的上界,...
  • 第2章 基本数据类型与数值表达式2.1 知识要点计算机的基本功能是进行数据处理。在C++语言中,数据处理的基本对象是常量和变量。运算是对各种形式的数据进行处理。数据在内存中存放的情况由数据类型所决定。数据的...
  • CSI:信息的表示与处理-数值陷阱(二) 前言 本篇继续上一篇,进行浮点数的介绍,浮点数的表示并不像整型那样简单,其在计算机中的运算也会使用更多地时钟周期。我们都知道计算机并不能绝对正确的表示浮点数,都是在...
  • 线性回归实例——投资额(python、OLS最小二乘) 一、问题描述:     建立投资额模型,研究某地区实际投资额与国民生产值(GNP)及物价指数(PI)的关系,根据对未来GNP及PI的估计,预测未来投资额。以下是该...
  • 本篇文章中,我们一起探讨了OpenCV中霍夫变换相关的知识点...此博文一共有四个配套的简短的示例程序,其详细注释过的代码都在文中贴出,且文章最后提供了综合示例程序的下载。 先尝鲜一下其中一个示例程序的运行截图:
  • 再把第四个字符插入到前三个字符中……待排序的字符串已在主函数中赋予。 请修改程序中的错误,得出正确的结果。 注意:不要改动main函数,不能增行或删行,也不能更改程序的结构。 代码如下: #include<stdio.h&...
  • 如果数据集是线性可分的,那么感知机获得的模型可能有很多,而支持向量机选择的是间隔最大的那一。 支持向量机还支持核技巧,从而使它成为实质上的非线性分类器。 支持向量机支持处理线性可分...
  • SVM解释:、线性不可分的情况

    万次阅读 多人点赞 2018-07-23 08:41:42
    之前的博客介绍了在数据为线性可分的情况下,如何用SVM对数据集训练,从而得到一线性分类器,也就是超平面WX+b=0WX+b=0WX + b = 0. 但是我已经强调过多次,线性可分的情况有相当的局限,所以SVM的终极目标还是要...
  • 例题:请补充fun函数,该函数的功能是:交换数组a中最大和最小元素的位置,结果重新保存在原数组中,其它元素位置不变。注意数组a中没有相同元素。 例如,输入 “20,40,78,60,52”,则输出结果为 “78,40,20,...
  • 本文代码实现基本按照《数据结构》课本目录顺序,外加大量的复杂算法实现,一篇文章足够。能换你一收藏了吧?
  •  这里的一元/二元指的就是输出标签的情况,这个具体的例子我还没看到,example文件夹中四个例子,也都是只用了Unigram,没有用Bigarm,因此感觉一般Unigram feature就够了。 5.3.3 模板例子  这是CoNLL 2000...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 27,040
精华内容 10,816
关键字:

下列四个数值最小