非线性最小二乘法的理论

三轴运动平台

---

1、 背景

$\qquad$非线性最小二乘算法的应用大约于上世纪80年代初，非线性最小二乘法除可直接用于估计静态非线性模型的参数外，在时间序列建模、连续动态模型的参数估计中，也往往遇到求解非线性最小二乘问题。其算法主要有两类：一类是搜索算法，另一类是迭代算法
$\qquad$搜索算法的思路是：按一定的规则选择若干组参数值，分别计算它们的目标函数值并比较大小；选出使目标函数值最小的参数值，同时舍弃其他的参数值；然后按规则补充新的参数值，再与原来留下的参数值进行比较，选出使目标函数达到最小的参数值。如此继续进行，直到选不出更好的参数值为止。以不同的规则选择参数值，即可构成不同的搜索算法。常用的方法有单纯形搜索法、复合形搜索法、随机搜索法等。
$\qquad$迭代算法是从参数的某一初始猜测值$\theta$出发，然后产生一系列的参数点，如果这个参数序列收敛到使目标函数极小的参数点娈，那么对充分大的$N$就可用$N$作为娈。迭代算法的一般步骤是：
1、给出初始猜测值$\theta$，并置迭代步数$i=1$
2、确定一个向量$v$作为第$i$步的迭代方向
3、用寻优的方法决定一个标量步长$\rho$
4、检查停机规则是否满足，如果不满足，则将$i$加1再从2开始重复；如果满足，则取$\theta$为值

典型的迭代算法有：

• 牛顿－拉夫森法
• 高斯-牛顿迭代法
• 修正高斯-牛顿迭代算法
• 阻尼最小二乘法
• Levenberg-Marquardt（莱文伯-麦夸脱）算法
• 变尺度法
  ...
  在这些算法中初值得选取至关重要。

2、 算法核心

一般非线性最小二乘法及其典型迭代算法：
提出问题：---量测关系式

$$y_k=f(x_k,B)+\Delta_k\qquad (k=1,\cdots,m;B=(b_1,b_2,\cdots,b_m);m\geq n)$$
其中，$y_k$可以是量测量，物理量测量，数学量测量
$\qquad B$是待估时常参数，可以是单变量，也可以是多变量
$\qquad x_k$可以是时间变量，也可以是其他已知函数，$x_k$可以是单变量，也可以是多变量
$\qquad \Delta_k$是对量测$y_k$的测量误差，可以是偏倚误差，也可以是随机误差，可以是平稳随机过程，也可以是非平稳随机过程；
$\qquad f(x_k,B)$是$B$的非线性函数。
$\qquad$一般非线性最小二乘问题是：给定量测$y_k$与关系式，要求确定参数$B$，使：
$$\min_{B} Q(B)=\sum_{k=1}^m \Delta_k^2=\sum_{k=1}^{m}[y_k-f(x_k,B)]^2$$
$\qquad$该式是一个无条件极值问题，根据微积分学教程的方法，可以求：
$$\frac{\partial Q(B)}{\partial B}=0$$
的n个方程的解，找出平衡值（为极值存在的必要条件），再来验证是极大值还是极小值。
$\qquad$由上式可以看出：一般非线性最小二乘问题的解，是转化为非线性方程组的求解问题。而非线性方程组并无解析解的一般方法，为此发展了一系列的数值解法[1]。
1、高斯-牛顿迭代法
$\qquad$选给定$b_i$的初始近似值$b_i^{(0)}$,记录真值与初始$b_i^{(0)}$之差（未知）为$\delta_i$
$$b_i=b_i^{(0)}+\delta_i\qquad(i=1,2,\cdots n)$$
$\qquad \delta_i$由下面联立方程组给出：
$$\begin{pmatrix} a_{11} & a_{12} & \cdots &a_{1n}\\ a_{21} & a_{22} & \cdots &a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots &a_{nn}\\ \end{pmatrix} \begin{pmatrix} \delta_1\\ \delta_2\\ \vdots\\ \delta_n \end{pmatrix} \begin{pmatrix} a_{1y}\\ a_{2y}\\ \vdots\\ a_{ny} \end{pmatrix}$$
其中：
$$\begin{cases} a_{ij}=\sum_{k=1}^{m}\frac{\partial f_{ko}}{\partial b_i}\frac{\partial f_{ko}}{\partial b_j}, \qquad\qquad\qquad\qquad\qquad\quad (i,j=1,2\cdots,n)\\[2ex] a_{iy}=\sum_{k=1}^{m}\frac{\partial f_{ko}}{\partial b_i}(y_k-f_{ko}),\qquad\qquad\qquad\qquad\quad (i=1,2,\cdots,n)\\[2ex] f_{ko}=f(x_k,b_1^{(0)},\cdots,b_n^{(0)})\\[2ex] \left. \frac{\partial f_{ko}}{\partial b_i}=\frac{\partial f_{ko}=f(x_k,b_1,\cdots,b_n)}{\partial b_i}\right|_{x=x_k,b_i=b_i^{(0)}}\qquad (i=1,2,\cdots,n)\\[2ex] \end{cases}$$
$\quad$当$|\delta_i|$值较大时，可令当前的$b_i$值代替原来的近似值$b_i^{(0)}$，重复计算$a_{ij}、a_{iy}$并求解联立方程得到新的$\delta_i$，直至$|\delta_i|$值小到可以忽略不计为止。

2、麦夸脱（Levenberg-Marquardt）迭代法
$\qquad$该法与高斯-牛顿法不同，在于确定$\delta_i$方程组的不同，此处取联立方程组为：

$$\begin{pmatrix} a_{11}+d & a_{12} & \cdots &a_{1n}\\ a_{21} & a_{22}+d & \cdots &a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots &a_{nn}+d\\ \end{pmatrix} \begin{pmatrix} \delta_1\\ \delta_2\\ \vdots\\ \delta_n \end{pmatrix} \begin{pmatrix} a_{1y}\\ a_{2y}\\ \vdots\\ a_{ny} \end{pmatrix}$$
$\qquad$上式中，$d=0$则退化为高斯-牛顿迭代法，$d$的选取方法为：
1）算出初始值$b_i^{(0)}$所对应的残差$Q^{(0)}$，指定一个常数$C$（$C>1$,例如$C=10$），并给出$d$的一个初始值$d^{(0)}$（例如：当$a_{11}=1$时，令$d^{(0)}=0.01$）。

2）进行下一步迭代时，令

$$d=c^\alpha d^{(0)}\qquad(\alpha = -1,0,1,2,\cdots)$$
其中$\alpha$值尽可能取得小些，但需要保证由上面方程组解出的$\delta_i$（进而是$b_i$）相应的残差平方和$Q$小于$Q_{(0)}$,即：
$$Q<Q^{(0)}$$
先令$d=c^{(-1)} d^{(0)}$,若上式成立，则此次迭代完成；否则令d=c^{(0)} d^{(0)}$, 若上式成立，则这一次迭代完成，否则再令d=c^{(1)} d^{(0)}\cdots$直至上式成立。

3）以$d,b_i$和$Q$的当前值代替$d^{(0)}、b_i^{(0)}$和$Q^{(0)}$，重复第2）步，做下次迭代，直至$|\delta_i|$可以忽略。

3、平台应用

$\qquad$建立速度输入（$U(s)$）速度输出（$V(s)$）的数学模型：

$$G(s)=\frac{V(s)}{U(s)}=\frac{K_v^g\cdot K_a\cdot K_L}{s^3[(1+K_a)J_MJ_LN_G^2]+s^2[N_G^2(J_LK_v^gK_a+J_MD_L+J_MK_aD_L)]+\\s(J_LK_L+N_G^2K_v^gK_aD_L+K_aK_LJ_L)+K_aK_LD_L}$$
$\qquad$因此，需要的辨识的模型可以为：
$$G(s)=\frac{a_0}{b_3s^3+b_2s^2+b_1s+1}$$
$\qquad$得到的仍为高阶系统模型，在系统速度处于1/20 ~ 1/5额定转速时机电一体化伺服系统可以近似用二阶模型近似：

$\qquad$因此可以建立输入输出降阶模型为：
$$G(s)=\frac{Y(s)}{U(s)}=\frac{K}{s(T_ps+1)}$$

$\qquad$记录输入量（给定目标速度）、输出量（实际瞬间转速，实际电机位置），可以通过辨识实际输出转速（物理量）与输入量（数字量）之间的输入输出模型进而获得实际输出电机位置（物理量）与输入量（数字量）之间的输入输出模型，即：

$$G(s)=\frac{V(s)}{U(s)}=\frac{K}{T_ps+1}\\ \iff T_p\cdot \frac{y(t)-y(t-1)}{T}+y(t)=K\cdot u(t)\\ \iff y(t)=\frac{T_p}{T+T_p}\cdot y(t-1)+\frac{T}{T+T_p}\cdot K\cdot u(t)$$
$\qquad$这里的$T$为伺服器反馈信息时间差。
$$\min_{K,T_p}Q(K,T_p)=\sum_{k=1}^m[y(k)-f(x_k,y_{k-1}K,T_p)]^2\\ \min_{K,T_p}Q(K,T_p)=\sum_{k=1}^m[y(k)-(\frac{T_p}{T+T_p}\cdot y(k-1)+\frac{T}{T+T_p}\cdot K\cdot u(k))]^2\\ $$

N 参考文献

---

1 Lawson C. L.，Hanson R. J. Solving Least Squares
Problems[M]. Prentice-Hall Inc,1974.
2 刘钦圣. 最小二乘问题计算方法[M]. 北京：北京工业
大学出版社，1989
[3]王新洲、非线性模型参数估计理论与应用.武汉大学出版社,2002