非线性模型参数估计的理论知识补充

三轴运动平台

---

1.背景知识

什么是参数估计？
参数估计是指对有一定统计分布的观测子样建立起某种函数模型,在某一准则函数约束下采用一定的数值计算方法来估计参数的过程。
测量领域中参数估计的过程大概分为三步
1、在观测子样和估计参数间建立起符合统计规律的函数模型
2、根据评估目标选择合适的估计准则
3、采用适当的数值解算方法,确保估计参数解的精确性和稳定性。

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

2.非线性最小二乘估计的定义和存在性定理

非线性模型响应的误差方程为：

$$\begin{align}V =f(\hat X)-L\end{align}$$ 残差平方和为：
$$\begin{align}V'V&=||V||^2=||f(\hat X)-L||^2\\&=(f(\hat X)-L)'(f(\hat X)-L)\end{align}$$ $\qquad$根据等价观测理论，所有观测值都可以变换为独立观测值，可以采用同精度观测讨论。

---定义
非线性模型中参数$X$的一个估计量$\hat X$,如果满足关系：

$$\begin{align}V'V=min\end{align}$$
$\qquad$则称$\hat X$是$X$的一个非线性最小二乘估计，用$\hat X_{NLS}$表示。这个定义与线性模型最小二乘估计的定义是完全一致的，在几何意义上就是观测空间至解空间的距离最短，或者说$f(\hat X)$是解轨迹$\pi$上离观测值$L$最近的点，L到$\pi$的距离就是$||V||$。

$\qquad$上式定义了非线性最小二乘的估计量$\hat X$，而在参数空间是否存在这样的量？

定理 假设参数空间$\Theta$为$R^t$上的紧子集，$f(X)$关于$X$在参数空间$\Theta$上连续，则必存在$R^n$上的可测函数$\hat X=\hat X(L)$,使得：

$$\begin{align}||f(\hat X(L))-L||^2=\min_{X\in\Theta}||f(X)-L||^2,L\in R^n\end{align}$$ $\qquad$由于$||f(X)-L||^2$在$\Theta\times R^n$上关于$X$和$L$都是连续函数，又因为$\Theta$为紧子集，所以对于任意的$L$,$||f(X)-L||^2$在参数空间$\Theta$上必存在极小值$\hat X(L)$

3.非线性最小二乘估计的近似解法[1][2]

3.1.线性近似

$\qquad$当非线性模型的非线性强度较弱时，可以将非线性模型在$X^{(0)}$处线性化，并用线性模型的求解理论来解算非线性模型式，得误差方程：

$$\begin{align} \because L&=f(X)+\Delta \\ & =BX+{{B}_{0}}+\Delta \\ \therefore V&=B(X^{(0)})dX-(L-f(X^{(0)})) \end{align}$$
根据最小二乘原理可以解得：
$$\begin{align}dX=(B'(X^{(0)})B(X^{(0)}))^{-1}B'(X^{(0)})(L-f(X^{(0)}))\end{align}$$
于是参数$X$的最小二乘估计量为：
$$\begin{align}\hat X=X^{(0)}+dX\end{align}$$
$\qquad$当然还需要考虑线性近似所引起的模型误差对最小二乘估计量$\hat X$的影响，只有当其小于等于观测误差对参数估值$\hat X$的影响时可以忽略不计。

3.2.数值迭代解法

$\qquad$对于非线性强度很强的非线性模型，由于线性近似将产生大于观测误差的模型误差，所以一般采用迭代的方法求解。

$$\begin{align}min\quad V'V(\hat X)&=(f(\hat X)-L)'(f(\hat X)-L)\\ &=f'(\hat X)f(\hat X)-2f'(\hat X)L+L'L\end{align}$$
$\qquad$由于$L'L$是一常量，所以上式等价于目标函数为：
$$\begin{align}min \quad R(\hat X)=f'(\hat X)f(\hat X)-2f'(\hat X)L\end{align}$$
此问题转化为非线性无约束最优化问题。
$\qquad$因为$f(\hat X)$是$\hat X$的非线性函数，所以对上式求一阶偏导数，并令其为零，并不能得到$\hat X$的显式表达式，所以求不出$\hat X$的解析解，因此，只能设法寻找某一近似解$X^*$使得：
$$\begin{align}R(X^*)\le R(\hat X)\end{align}$$ 成立，一般这个过程只有采用迭代的方法。

---

---常见的迭代方法有：

3.2.1 - 牛顿法

$\qquad$可以假设$R(\hat X)$的极小值$X^*$的一个近似值为$X^{(k)}$,在$X^{(k)}$附近将$R(X^*)$泰勒展开，取至二次项得：

$$\begin{align} R({{X}^{*}})&=R({{X}^{(k)}}+d{{X}^{(k)}}) \\ & =R({{X}^{(k)}})+{{g}^{(k)}}d{{X}^{(k)}}+\frac{1}{2}d{{X}^{'(k)}}{{G}_{k}}d{{X}^{(k)}} \end{align}$$
式中：

$$\begin{align} {{g}^{(k)}} &=\left( \begin{matrix} g_{1}^{(k)} & g_{2}^{(k)} & \cdots & g_{t}^{(k)} \\ \end{matrix} \right) \\ & ={{\left. \left( \begin{matrix} \frac{\partial R}{\partial {{x}_{1}}} & \frac{\partial R}{\partial {{x}_{2}}} & \cdots & \frac{\partial R}{\partial {{x}_{t}}} \\ \end{matrix} \right) \right|}_{X={{X}^{(k)}}}} \end{align}$$

$${{G}_{k}}= \left. \begin{pmatrix} \frac{{{\partial }^{2}}R}{\partial x_{1}^{2}} & \frac{{{\partial }^{2}}R}{\partial {{x}_{1}}\partial {{x}_{2}}} & \cdots &\frac{{{\partial }^{2}}R}{\partial {{x}_{1}}\partial {{x}_{t}}}\\ \frac{{{\partial }^{2}}R}{\partial {{x}_{2}}\partial {{x}_{1}}} & \frac{{{\partial }^{2}}R}{\partial x_{2}^{2}} & \cdots & \frac{{{\partial }^{2}}R}{\partial {{x}_{1}}\partial {{x}_{t}}}\\ \vdots &\vdots&\ddots &\vdots\\ \frac{{{\partial }^{2}}R}{\partial {{x}_{t}}\partial {{x}_{1}}} & \frac{{{\partial }^{2}}R}{\partial {{x}_{t}}\partial {{x}_{2}}} & \cdots & \frac{{{\partial }^{2}}R}{\partial x_{t}^{2}} \\ \end{pmatrix} \right|_{X={{X}^{(k)}}}={{G}_{k}}'$$
$G_k$称为$X^{(k)}$处的$Hessian$矩阵
$$\begin{align}d{{X}^{(k)}}={{X}^{*}}-{{X}^{(k)}}\end{align}$$
$g^{(k)}$是$R(X)$在$X^{(k)}$处的梯度方向。
$\qquad$由于$X^{(k)}$是$X^*$的一个已知的近似值，所以式（16）只是$dX^{(k)}$的函数，为了求得使得式（16）成立的$dX^{(k)}$，将该式对$dX^{(k)}$求偏导，并令其为0，得：
$$\begin{align} {{g}^{(k)}}+d{{X}^{'(k)}}{{G}_{k}}&=0 \\ {{G}_{k}}d{{X}^{(k)}}& =-g{{'}^{(k)}} \end{align}$$
$\qquad$当$G_k$非奇异时，由上式可以解得使式（16）成立的$dX^{(k)}$，即:
$$\begin{align} d{{X}^{(k)}}=-G_{k}^{-1}g{{'}^{(k)}} \end{align}$$
$\qquad$由假设知当$dX^{(k)}$充分小时，$dX^{(k)}$能够使得式（16）成立，但由于$X^*$未知，故$dX^{(k)}$不能充分小，需不断地迭代，直至$dX^{(k)}$充分小，其迭代公式为：

$$\begin{align} {{X}^{(k+1)}}&={{X}^{(k)}}+d{{X}^{(k)}} \\ & ={{X}^{(k)}}-G_{k}^{-1}g{{'}^{(k)}} \end{align}$$
$\qquad$上式就是牛顿迭代的基本公式，其终止迭代条件就是
$$\begin{align} R({{X}^{(k+1)}})&=R({{X}^{(k)}}) \\ {{g}^{(k)}}&=0 \end{align}$$
$\qquad$由于$R(X^{(k)})$是一个绝对值较大的数，而$g^{(k)}$的各元素的绝对值都很小，因此，由于计算机有效数字的限制，以第一终止迭代条件作为迭代收敛条件的收敛速度要比第二终止迭代条件作为迭代收敛条件快。
牛顿法的迭代步骤：

1. 选取初值$X^{(0)}$,并且令$k=0$
2. 按式（18）计算梯度方向$g^{(k)},若g^{(k)}=0则终止迭代$
3. 计算矩阵$G_k$
4. 求解线性方程组式（22），得到$dX^{(k)}$
5. 按式（24）计算新的近似值$X^{(k+1)}$
6. 计算目标函数值$R(X^{(k+1)})$，若$R({{X}^{(k+1)}})\neq R({{X}^{(k)}})$则转至2继续迭代。
7. 终止迭代，输出$X^{(k+1)}$和$R(X^{(k+1)})$，结束

总结：牛顿法虽然有很快的收敛速度，但是它总是局部收敛，而且对初值十分敏感。

---

3.2.2 - 信赖域法

$\qquad$综上,可以发现牛顿法的基本思想是用二次函数：

$$\begin{align} Q({{X}^{(k)}})=R({{X}^{(k)}})+{{g}^{(k)}}dX+\frac{1}{2}dX{{'}^{(k)}}{{G}_{k}}d{{X}^{(k)}} \end{align}$$
去逼近$R(X^*)$。只有当$dX^{(k)}$充分小时，$Q(X^{(k)})$才能很好的逼近$R（X^{(k)}）$。
$\qquad$既然只有当$dX^{(k)}$充分小时，$Q(X^{(k)})$才能很好的逼近$R（X^{(k)}）$，那么可以对$dX$加限制，然后再限制条件下寻求$R(X^*)$的极小值，将无约束最优化问题转换为约束最优化问题:
目标函数：
$$\min Q({{X}^{(k+1)}})=R({{X}^{(k)}})+{{g}^{(k)}}d{{X}^{(k)}}+\frac{1}{2}dX{{'}^{(k)}}{{G}_{k}}d{{X}^{(k)}}$$
约束条件：
$$\left\| d{{X}^{(k)}} \right\|\le {{h}_{k}}$$
式中，$h_k$是一个随迭代而变化的正数。

$\qquad$约束条件$\left\| d{{X}^{(k)}} \right\|\le {{h}_{k}}$限制了$d{X}^{(k)}$,使得$\left\| d{{X}^{(k)}} \right\|$不大于${{h}_{k}}$,这样$\left\| d{{X}^{(k)}} \right\|$总在一个给定的小区域内活动，这个区域是可信赖的，称该方法为信赖域法。
$\qquad$常数$h_k$取决于$Q(X^{(k+1)})$对$R(X^{(k+1)})$的逼近程度：

$$\begin{align} {{r}_{k}}=\frac{R({{x}^{(k+1)}})}{Q({{x}^{(k+1)}})} \end{align}$$
$\qquad{r}_{k}$越接近1，$Q(X^{(k+1)})$对$R(X^{(k+1)})$的逼近程度越好：
$$\begin{align} {{h}_{k+1}}=\left\{ \begin{matrix} \frac{\left\| d{{X}^{(k)}} \right\|}{4}, & {{r}_{k}}<0.25 \\[2ex] 2\left\| d{{X}^{(k)}} \right\|,& {{r}_{k}}>0.75 \\[2ex] {{h}_{k}}, &0.25<{{r}_{k}}<0.75 \\[2ex] \end{matrix} \right. \end{align}$$

信赖域法的迭代步骤：

1. 选取初值$X^{(0)}$和$h_0,k=0$。

2. 按照式（17）和式（18）计算梯度方向$g^{(k)}$和矩阵$G_k$,若$g^{(k)}=0$,则停止迭代

3. 按式（20）计算$dX^{(k)}$,并检查$\left\|dX^{(k)}\right\|$是否满足约束条件，若不满足，则采取适当的方法对$\left\|dX^{(k)}\right\|$予以压缩，然后再区域$h_k$内求使得$Q(X^{(k+1)})=min$的$dX^{(k)}$。

4. 计算$X^*$的新的近似值 $X^{(k+1)}=X^{(k)}+dX^{(k)}$
5. 按式（28）计算$r_k$，并按式（29）确定$h_{k+1}$
6. 检查$R(X^{k+1})=R(X^{k})$是否成立，若不成立，则继续迭代

7. 终止迭代，输出$X^{k+1}$和$R(X^{k})$,结束

   总结：信赖域法也与初值有关，但仍然是局部收敛，并不像想象的那样全局收敛。

---

3.2.3 - 拟牛顿法

$\qquad$拟牛顿法是在牛顿法的基础之上，用一个仅包含一阶偏导数信息的对称矩阵$Q_k$去逼近二阶偏导数矩阵$G_k$（Hessian）矩阵,然后再按牛顿法予以迭代。其中有一种数值方法可以确定$Q_k$阵：

$$\begin{align}{{G}_{k}}=\left( \begin{matrix} \frac{{{\partial }^{2}}R}{\partial x_{1}^{2}} & \frac{{{\partial }^{2}}R}{\partial {{x}_{1}}\partial {{x}_{2}}} & \cdots & \frac{{{\partial }^{2}}R}{\partial {{x}_{1}}\partial {{x}_{t}}} \\ \frac{{{\partial }^{2}}R}{\partial {{x}_{2}}\partial {{x}_{1}}} & \frac{{{\partial }^{2}}R}{\partial x_{2}^{2}} & \cdots & \frac{{{\partial }^{2}}R}{\partial {{x}_{2}}\partial {{x}_{t}}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{{{\partial }^{2}}R}{\partial {{x}_{t}}\partial {{x}_{1}}} & \frac{{{\partial }^{2}}R}{\partial {{x}_{t}}\partial {{x}_{2}}} & \cdots & \frac{{{\partial }^{2}}R}{\partial x_{t}^{2}} \\ \end{matrix} \right)=\left( \begin{matrix} \frac{\partial {{g}_{1}}}{\partial {{x}_{1}}} & \frac{\partial {{g}_{1}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{g}_{1}}}{\partial {{x}_{t}}} \\ \frac{\partial {{g}_{2}}}{\partial {{x}_{1}}} & \frac{\partial {{g}_{2}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{g}_{2}}}{\partial {{x}_{t}}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial {{g}_{t}}}{\partial {{x}_{1}}} & \frac{\partial {{g}_{t}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{g}_{t}}}{\partial {{x}_{t}}} \\ \end{matrix} \right)\end{align}$$
根据多元函数偏导数的定义：
$$\begin{align} \frac{\partial f({{x}_{1}},{{x}_{2}}\cdots {{x}_{t}})}{\partial {{x}_{i}}}=\underset{\Delta {{x}_{i}}\to 0}{\mathop{\lim }}\,\frac{f({{x}_{1}},\cdots ,{{x}_{i}}+\Delta {{x}_{i}},\cdots ,{{x}_{t}})-f({{x}_{1}},\cdots ,{{x}_{t}})}{\Delta {{x}_{i}}} \end{align}$$

$$\begin{align} {{G}_{k}}=\left( \begin{matrix} \underset{d{{x}_{1}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{1}}({{x}_{1}}+d{{x}_{1}},{{x}_{2}},\cdots ,{{x}_{t}})-{{g}_{1}}(x)}{d{{x}_{1}}} & \cdots & \underset{d{{x}_{t}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{1}}({{x}_{1}},\cdots ,{{x}_{t}}+d{{x}_{t}})-{{g}_{1}}(x)}{d{{x}_{t}}} \\ \underset{d{{x}_{1}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{2}}({{x}_{1}}+d{{x}_{1}},{{x}_{2}},\cdots ,{{x}_{t}})-{{g}_{2}}(x)}{d{{x}_{1}}} & \cdots & \underset{d{{x}_{t}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{2}}({{x}_{1}},\cdots ,{{x}_{t}}+d{{x}_{t}})-{{g}_{2}}(x)}{d{{x}_{t}}} \\ \vdots & \ddots & \vdots \\ \underset{d{{x}_{1}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{t}}({{x}_{1}}+d{{x}_{1}},{{x}_{2}},\cdots ,{{x}_{t}})-{{g}_{t}}(x)}{d{{x}_{1}}} & \cdots & \underset{d{{x}_{t}}\to 0}{\mathop{\lim }}\,\frac{{{g}_{t}}({{x}_{1}},\cdots ,{{x}_{t}}+d{{x}_{t}})-{{g}_{t}}(x)}{d{{x}_{t}}} \\ \end{matrix} \right) \end{align}$$

令：

$$\begin{align} {{\overline{g}}^{(k)}}=\left( \begin{matrix} \overline{g}_{1}^{(k)} \\ \overline{g}_{2}^{(k)} \\ \vdots \\ \overline{g}_{t}^{(k)} \\ \end{matrix} \right)=\left( \begin{matrix} \overline{g}_{1}^{(k)}(x_{1}^{(k+1)},x_{2}^{(k)},\cdots ,x_{t}^{(k)}) \\ \overline{g}_{2}^{(k)}(x_{1}^{(k)},x_{2}^{(k+1)},\cdots ,x_{t}^{(k)}) \\ \vdots \\ \overline{g}_{t}^{(k)}(x_{1}^{(k)},x_{2}^{(k)},\cdots ,x_{t}^{(k+1)}) \\ \end{matrix} \right) \end{align}$$
去掉极限后，就可以得到$G_k$的近似矩阵：
$$\begin{align} {{G}_{k}}\approx {{Q}_{k}}=\left( \begin{matrix} \frac{\overline{g}_{1}^{(k)}-g_{1}^{(k)}}{dx_{1}^{(k)}} & \frac{1}{2}(\frac{\overline{g}_{1}^{(k)}-g_{1}^{(k)}}{dx_{2}^{(k)}}+\frac{\overline{g}_{2}^{(k)}-g_{2}^{(k)}}{dx_{1}^{(k)}}) & \cdots & \frac{1}{2}(\frac{\overline{g}_{1}^{(k)}-g_{1}^{(k)}}{dx_{t}^{(k)}}+\frac{\overline{g}_{t}^{(k)}-g_{t}^{(k)}}{dx_{t}^{(k)}}) \\ {} & \frac{\overline{g}_{2}^{(k)}-g_{2}^{(k)}}{dx_{2}^{(k)}} & \cdots & \frac{1}{2}(\frac{\overline{g}_{2}^{(k)}-g_{2}^{(k)}}{dx_{t}^{(k)}}+\frac{\overline{g}_{t}^{(k)}-g_{t}^{(k)}}{dx_{2}^{(k)}}) \\ * & {} & \ddots & \vdots \\ {} & {} & {} & \frac{\overline{g}_{t}^{(k)}-g_{t}^{(k)}}{dx_{t}^{(k)}} \\ \end{matrix} \right) \end{align}$$
$\qquad$用上式定义的对称矩阵$Q_k$能够较为准确地逼近$G_k$，又只包含$R$的一阶偏导数信息，不需要求二阶Hessian矩阵，有了$Q_k$之后，一切迭代均按照牛顿法进行，由于拟牛顿法一开始就要计算$Q_k$阵，所以计算前除了给定$X$的初值$X^{(0)}$之外，还必须给定$dX$的初值$dX^{(0)}$。这里可以这样确定：当给定$X^{(0)}$之后，将$X^{(0)}$减去一个很接近的$X^{(0)}$的向量$\overline X^{(0)}$,则差值就是$dX^{(0)}$，即：
$$\begin{align} d{{X}^{(0)}}={{X}^{(0)}}-{{\overline{X}}^{(0)}} \end{align}$$

总结：拟牛顿法也与初值有关，仍然是局部收敛，只是相对于牛顿法而言降低的计算量但计算复杂度依旧很高。

---

3.2.4 - 最速下降法

$\qquad$综上三种方法，都需要计算$R(X)$的二阶偏导数矩阵，而最速下降法只假定$R(X)$在解$\hat X$附近具有二阶连续偏导数，且此二阶偏导数矩阵的行列式大于、等于0，并不涉及到具体的Hessian矩阵计算。
$\qquad$最速下降法的基本思想是基于这样的事实：目标函数$R(\hat X)$在沿着$X^{(k)}$处的梯度方向$g^{(k)}$上数值增加最快。既然我们现在是要求目标函数$R(\hat X)$的最小值，因此，若是在寻找$R(\hat X)$的最小值点$X^*$的过程中，采用的是沿着$R(\hat X)$在$X^{(k)}$处的负梯度方向上寻找，势必使$R(\hat X)$的数值下降最快，所以称按负梯度方向搜寻$X^*$的方法为最速下降法，即：

$$\begin{align} {{X}^{(k+1)}}={{X}^{(k)}}+d{{X}^{(k)}}={{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}} \end{align}$$
式中：$\lambda$为实数，称步长；梯度方向$g^{k}$仍按式（17）计算

目标函数：

$$\begin{align} & R({{X}^{(k+1)}})=\min (R({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}})) \\ & {{V}^{'}}V({{X}^{(k+1)}})=\min ({{V}^{'}}V({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}})) \end{align}$$
将${{V}^{'}}V({{X}^{(k+1)}})$在$X^{(k)}$处泰勒展开，取至二次项，但略去$V$对$X$的二阶偏导数
$$\begin{align} {{V}^{'}}V({{X}^{(k+1)}})&=({{V}^{'}}V({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}}))=v_{1}^{2}({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}}) \\ & +v_{2}^{2}({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}})+\cdots +v_{n}^{2}({{X}^{(k)}}-{{\lambda }^{(k)}}{{g}^{'(k)}}) \\ & \approx v_{1}^{2}({{X}^{(k)}})-2{{\lambda }^{(k)}}{{v}_{1}}({{X}^{(k)}})\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\tfrac{\partial {{f}_{1}}}{\partial {{x}_{j}}}+{{({{\lambda }^{(k)}})}^{2}}\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\sum\limits_{l=1}^{t}{g_{l}^{(k)}}\frac{\partial {{f}_{1}}}{\partial {{x}_{j}}}\frac{\partial {{f}_{1}}}{\partial {{x}_{l}}} \\ & +v_{2}^{2}({{X}^{(k)}})-2{{\lambda }^{(k)}}{{v}_{2}}({{X}^{(k)}})\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\tfrac{\partial {{f}_{2}}}{\partial {{x}_{j}}}+{{({{\lambda }^{(k)}})}^{2}}\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\sum\limits_{l=1}^{t}{g_{l}^{(k)}}\frac{\partial {{f}_{2}}}{\partial {{x}_{j}}}\frac{\partial {{f}_{2}}}{\partial {{x}_{l}}} \\ & +\cdots +v_{n}^{2}({{X}^{(k)}})-2{{\lambda }^{(k)}}{{v}_{n}}({{X}^{(k)}})\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\tfrac{\partial {{f}_{n}}}{\partial {{x}_{j}}}+{{({{\lambda }^{(k)}})}^{2}}\sum\limits_{j=1}^{t}{g_{j}^{(k)}}\sum\limits_{l=1}^{t}{g_{l}^{(k)}}\frac{\partial {{f}_{n}}}{\partial {{x}_{j}}}\frac{\partial {{f}_{n}}}{\partial {{x}_{l}}} \\ & ={{V}^{'}}V({{X}^{(k)}})-2{{\lambda }^{(k)}}{{V}^{'}}({{X}^{(k)}})B({{X}^{(k)}}){{g}^{'(k)}}+{{({{\lambda }^{(k)}})}^{2}}\sum\limits_{i=1}^{n}{g_{i}^{(k)}}W_{i}^{(k)}g_{i}^{'(k)} \end{align}$$
式中：
$$\begin{align} & V({{X}^{(k)}})={{(\begin{matrix} {{f}_{1}}(X{}^{(k)})-{{L}_{1}} & {{f}_{2}}(X{}^{(k)})-{{L}_{2}} & \cdots & {{f}_{n}}(X{}^{(k)})-{{L}_{n}} \\ \end{matrix})}^{'}} \\ & B({{X}^{(k)}})={{\left( \begin{matrix} \frac{\partial {{f}_{1}}}{\partial {{x}_{1}}} & \frac{\partial {{f}_{1}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{f}_{1}}}{\partial {{x}_{t}}} \\ \frac{\partial {{f}_{2}}}{\partial {{x}_{1}}} & \frac{\partial {{f}_{2}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{f}_{2}}}{\partial {{x}_{t}}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial {{f}_{n}}}{\partial {{x}_{1}}} & \frac{\partial {{f}_{n}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{f}_{n}}}{\partial {{x}_{t}}} \\ \end{matrix} \right)}_{X={{X}^{(k)}}}} \\ & W_{i}^{(k)}={{\left( \begin{matrix} {{(\frac{\partial {{f}_{1}}}{\partial {{x}_{1}}})}^{2}} & \frac{\partial {{f}_{i}}}{\partial {{x}_{1}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{2}}} & \cdots & \frac{\partial {{f}_{i}}}{\partial {{x}_{1}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{t}}} \\ \frac{\partial {{f}_{i}}}{\partial {{x}_{2}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{1}}} & {{(\frac{\partial {{f}_{2}}}{\partial {{x}_{2}}})}^{2}} & \cdots & \frac{\partial {{f}_{i}}}{\partial {{x}_{2}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{t}}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial {{f}_{i}}}{\partial {{x}_{t}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{1}}} & \frac{\partial {{f}_{i}}}{\partial {{x}_{t}}}\frac{\partial {{f}_{i}}}{\partial {{x}_{2}}} & \cdots & {{(\frac{\partial {{f}_{i}}}{\partial {{x}_{t}}})}^{2}} \\ \end{matrix} \right)}_{X={{X}^{(k)}}}} \\ \end{align}$$

$\qquad$为了确定满足条件的$\lambda$，将式（44）对$\lambda$求导数，令其为0，得：

$$\begin{align} & {{\lambda }^{(k)}}=\frac{{{V}^{'}}({{X}^{(k)}})B({{X}^{(k)}}){{g}^{'(k)}}}{\sum\limits_{i=1}^{n}{g_{i}^{(k)}W_{i}^{(k)}g_{i}^{'(k)}}} \\ & =\frac{{{(f({{X}^{(k)}})-L)}^{'}}B({{X}^{(k)}}){{B}^{'}}({{X}^{(k)}})(f({{X}^{(k)}})-L)}{2{{(f({{X}^{(k)}})-L)}^{'}}B({{X}^{(k)}}){{B}^{'}}({{X}^{(k)}})B({{X}^{(k)}}){{B}^{'}}({{X}^{(k)}})(f({{X}^{(k)}})-L)} \end{align}$$

最速下降法的迭代步骤如下：

1. 选取初值$X^{(0)}$。
2. 按照式（17）计算梯度方向$g^{(k)}$，如果$g^{(k)}=0$，则终止迭代
3. 按照式（46）和式（47）计算矩阵$B(X^{(k)})$和$W_i^{(k)}$
4. 按照式（48）计算$\lambda^{(k)}$
5. 按照式（36）计算新的近似值$X^{(k+1)}$,并且计算$R(X^{(k+1)})$
6. 若$R(X^{(k+1)}) \neq R(X^{(k)})$
7. 终止迭代，输出$X^{(k+1)}、R(X^{(k+1)})$，结束

总结：一般来说，最速下降法对任意初值都能收敛，但收敛速度并不像该方法的名称那样是最速的，相反它是最慢的，这是因为最速下降法在接近最小值点时会产生拉锯现象。

---

3.2.5 - 高斯-牛顿法

$\qquad$前面一些方法都是求目标函数$R(X^*)$最小化的非线性最优化算法，与我们测量平差中已掌握的方法相差很多。因此，测量平差中的一些软件都不能直接应用。而高斯-牛顿法则不同，几乎可以完全不改变原测量平差程序。

测量平差：
$\qquad$由于测量仪器的精度不完善和人为因素及外界条件的影响，测量误差总是不可避免的。为了提高成果的质量，处理好这些测量中存在的误差问题，观测值的个数往往要多于确定未知量所必须观测的个数，也就是要进行多余观测。有了多余观测，势必在观测结果之间产生矛盾，测量平差的目的就在于消除这些矛盾而求得观测量的最可靠结果并评定测量成果的精度。测量平差采用的原理就是“最小二乘法”。

-
$\qquad$高斯-牛顿法的基本出发点就是在初值$X^{(0)}$处对非线性模型进行线性近似，并按照传统的平差方法求出一次近似值$X^{(1)}$,然后反复迭代，直到前后两次$V^{'}V$值相等，迭代步骤如下：
$\qquad$假设非线性模型存在一阶连续偏导数，且参数$X$之间相互独立，则在近似值$X^{(0)}$处线性化，得误差方程：

$$\begin{align} V=B(X^{(0)})dX-(L-f(X^{(0)})) \end{align}$$
$\qquad$根据最小二乘原理，有：
$$\begin{align} X^{(1)}=X^{(0)}+(B^{'}(X^{(0)})B(X^{(0)}))^{-1}B^{'}(X^{(0)})(L-f（X^{(0)})) \end{align}$$
$\qquad$求得$X^{(1)}$之后，再以$X^{(1)}$为近似值迭代，其迭代公式为:
$$\begin{align} X^{(k+1)}=X^{(k)}+(B^{'}(X^{(k)})B(X^{(k)}))^{-1}B^{'}(X^{(k)})(L-f（X^{(k)})) \end{align}$$
终止迭代条件为：
$$\begin{align} R(X^{(k+1)})=R(X^{(k)}) \end{align}$$
$\qquad$高斯-牛顿法具有一定的合理性，因为若$L=f(X)+\Delta$是线性模型，则有$B(X^{(0)})=B$,$f(X^{(0)})=BX^{(0)}$,于是:
$$\begin{align} X^{(1)} &=X^{(0)}+(B^{'}B)^{-1}B^{'}(L-BX^{0})\\ & = (B^{'}B)^{-1}B^{'}L\\ \end{align}$$
$\qquad$上式表明:若$L=f(X)+\Delta$是线性模型，则由高斯-牛顿法从任意初值出发，经一次迭代就可以得到最小二乘估计的精确解。当非线性模型$L=f(X)+\Delta$的非线性强度较弱时，高斯-牛顿法就是较好的方法

总结：虽然高斯-牛顿法有一定的合理性，但在具体执行时，可能会产生一些问题。首先是对初值的依赖性较大，当初值较差时，会出现迭代发散现象，使迭代无法进行下去，好在我们在实际计算时，总是用观测值算出$X^{(0)}$,$X^{(0)}$与$X$很接近，故一般可迭代收敛。

---

3.2.6 - 改进的高斯-牛顿法

$\qquad$高斯-牛顿法对初值的依赖性较强，当初值较差时，会出现迭代发散现象。为了克服这个缺点，需要对高斯-牛顿法进行下列改进：

定理 设$X^{(k)}$是$X^*$的近似值，则$V^{'}V(X^{(k)})$一定不能达到最小，于是有：

$$\begin{align} \left.\frac{\partial V^{'}V(X^{(k)})}{\partial X}\right|_{X=X^{(k)}}\neq 0 \end{align}$$
$\qquad$设$dX^{(k)}=(B^{'}(X^{(k)})B(X^{(k)}))^{-1}B^{'}(X^{(k)})(L-f（X^{(k)}))$,那么必存在$\lambda^*>0$，使$\lambda \in [0,\lambda^*]$时，有
$$\begin{align} V^{'}V(X^{(k)}+\lambda dX^{(k)})<V^{'}V(X^{(k)}) \end{align}$$

-
$\qquad$此定理说明，当我们用高斯-牛顿法求出$dX^{(k)}$后，若适当选取$\lambda^{(k)}$,使:

$$\begin{align} X^{(k+1)}=X^{(k)}+\lambda dX^{(k)} \end{align}$$
$\qquad$则一定有：$V^{'}V(X^{(k+1)})<V^{'}V(X^{(k)})$,这样就能保证$V^{'}V(X^{(k)})$逐步向$V^{'}V$的极小值靠近，避免迭代过程的波动性，从而保证得到收敛的非线性最小二乘估计。根据这个思想，可构成如下得迭代算法：

• 根据初值$X^{(0)}$,对非线性模型式线性化，并按线性的最小二乘估计求解$dX^{(0)}$。
• 适当选取$\lambda^{(0)}$,并按式（58）计算$X^{(1)}$,同时计算$R(X^{(1)})$
• 若$R(X^{(1)})=R(X^{(0)})$,则停止迭代，否则令$X^{(0)}=X^{(1)}$,

• 输出$X^{(1)}$和$R^{(1)}$，结束

  $\qquad$该算法的关键是要计算$\lambda^{(k)}$,韦博成提出了对$R(X)$采用三点抛物线近似方法。即分别求出$\lambda=0,\lambda=0.5$和$\lambda=1$时$R(X)$的值，这三个值为:$R(X^{(k)}),R(X^{(k)}+\frac 12 dX^{(k)}),R(X^{(k)}+dX^{(k)})$,则：
  

  $$\begin{align} \lambda^{(k)}=\frac{\frac 12 + \frac 14(R(X^{(k)})-R(X^{(k)}+dX^{(k)}))}{R(X^{(k)})-2R(X^{(k)}+\frac 12 dX^{(k)}+R(X^{(k)}+dX^{(k)}))} \end{align}$$

  总结：改进的高斯-牛顿法除具备高斯-牛顿法的全部优点之外，还在很大程度上克服了高斯-牛顿法强烈依赖初值的缺点，是一种很有实用价值的方法。

---

3.2.7 - 阻尼最小二乘法

$\qquad$高斯-牛顿法和改进的高斯-牛顿法有解的必要条件是$B(X^{(k)})$矩阵列满秩，但在非线性秩亏自由网平差中，由于缺少基准，$B(X^{(k)})$矩阵总是列降秩的。在这种情况下，高斯-牛顿法和改进的高斯-牛顿法都不能使用。另外，当非线性模型中存在复共线关系时，尽管$B(X^{(k)})$矩阵列满秩，但$B^{'}(X^{(k)})B(X^{(k)})$的条件数很大，使得$B^{'}(X^{(k)})B(X^{(k)})$呈病态。这时候高斯-牛顿法和改进的高斯-牛顿法也无法使用。为了克服这个缺点，可以采用增大$B^{'}(X^{(k)})B(X^{(k)})$的主对角线元素的办法:

$$\begin{align} X^{(k+1)}=X^{(k)}+(B^{'}(X^{(k)})B(X^{(k)})+\alpha^{(k)}I)^{-1}B^{'}(X^{(k)})(L-f(X^{(k)})) \end{align}$$
式中：$\alpha^{(k)}$为大于等于0的任意常数，称为阻尼因子。因此称上式迭代的算法为阻尼最小二乘法，又称Levenberg - Marquardt（LM）算法，后来又有学者Fletcher对之前的LM算法进行改进称LMF算法。显然引进阻尼因子后，矩阵$B^{'}(X^{(k)})B(X^{(k)})+\alpha^{(k)}I$对任何正数$\alpha$，总具有对称正定的性质，因为此时矩阵$B^{'}(X^{(k)})B(X^{(k)})+\alpha^{(k)}I$的所有特征值均为正。
$\qquad$如果在上式的基础上再引入步长因子$\lambda^{(k)}$,则:
$$\begin{align} X^{(k+1)}=X^{(k)}+\lambda^{(k)}(B^{'}(X^{(k)})B(X^{(k)})+\alpha^{(k)}I)^{-1}B^{'}(X^{(k)})(L-f(X^{(k)})) \end{align}$$
$\qquad$至于$\alpha^{(k)}$的确定，对于非线性秩亏自由网平差问题，一般取：
$$\begin{align} \alpha^{(k)}=tr(B^{'}(X^{(k)})B(X^{(k)}))\sqrt k \end{align}$$
$\qquad$ 然后根据$V^{(k)}$和$V^{(k+1)}$来调整$\alpha^{(k+1)}$：
$$\begin{align} \alpha^{(k+1)}=\alpha^{(k)}\frac{V^{'}V^{(k+1)}}{V^{'}V^{(k)}} \end{align}$$

适用于非线性秩亏自由网平差的阻尼最小二乘法的具体算法如下：
1.选取初值$X^{(k)}$(此时$k=0$)，计算$V^{'}V^{(k)}$。
2.将非线性模型在$X^{(k)}$处线性化。
3.按一般间接平差组成法方程。若$k=0$,则计算$\alpha^{(k)}=tr(B^{'}(X^{(k)})B(X^{(k)}))$。
4.按照式（59）计算$\lambda^{(k)}$。
5.按照式（60）计算$X^{(k+1)}$,并计算$V^{'}V^{(k+1)}$，若$V^{'}V^{(k+1)}=V^{'}V^{(k)}$，则停止迭代。
6.按照式（62）计算$\lambda^{(k+1)}$。继续从2开始迭代
7.输出$X^{(k+1)}$,$V^{'}V^{(k+1)}$，结束

4、平台应用

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

$$\begin{align}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}\end{align}$$
$\qquad$因此，需要的辨识的模型可以为：
$$\begin{align}G(s)=\frac{a_0}{b_3s^3+b_2s^2+b_1s+1}\end{align}$$
$\qquad$得到的仍为高阶系统模型，在系统速度处于1/20 ~ 1/5额定转速时机电一体化伺服系统可以近似用二阶模型近似：

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

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

$$\begin{align} G(s)=\frac{Y(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) \end{align}$$
$\qquad$这里的$T$为伺服器反馈信息时间差。
$$\begin{align} \min_{K,T_p}Q(K,T_p)&=\sum_{k=1}^m[y(k)-f(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\\ \end{align}$$

5.应用MATLAB工具箱辨识平台电机模型参数

什么是过程模型？
$\qquad$过程模型的结构是简单连续时间传递函数描述了动态线性系统中的一项或者多项元素：

• 静态增益$K_p$
• 一个或多个时间常数$T_{pk}$,对于复杂极点模型，时间常数被称为$T_\omega$(等于所述固有频率的倒数)以及阻尼系数$\zeta$
• 过程零点$T_z$
• 系统可能存在的响应输出之前的时间延迟Td（死区时间）。
• 纯积分
  $\qquad$过程模型普遍用于描述动态系统在许多行业，适用于不同的生产环境。这些模型的优点是它们非常简单，支持传输延迟估计，以及模型系数有一个简单的解释为极点和零点。
  $\qquad$我们可以创建不同的模型结构通过改变极点的数量，增加纯积分环节或者添加系统的时延和零点个数。并且可以指定一阶、二阶或者三阶模型，当然极点可能都是实根或者是复根（欠阻尼模型）。
  $\qquad$例如，下面的模型结构是一阶连续时间过程模型，其中$K$是静态增益，$T_{P1}$是时间常数，和$T_d$是输入 - 输出延迟：
  $$G(s)=\frac{K_p}{1+sT_{p1}}e^{-sT_d}$$
  $\qquad$这样，$Y(s)=G(s)U(s)+E(s)$,这里$Y(s),U(s),E(s)$分别表示经过拉普拉斯变换的输出，输入和输出误差。输出误差$e(t)$是方差为$\lambda$的高斯白噪声。也可以通过在输出上添加一个干扰模型$H(s)$输出有色噪声,像是$Y(s)=G(s)U(s)+H(s)E(s)$

5.1使用system identification app辨识过程模型[3]

• 将数据导入系统识别的应用程序
• 对需要的数据进行一系列数据预处理操作，如果需要非零偏，比如模型包含纯积分则不消除趋势。在其它情况下，为提高模型的精确性，应该消除数据趋势。

1、在系统识别应用中，选择估值>过程模型打开过程模型对话框。

2、如果模型包含多个输入或多个输出，你可以指定是否估计相同的传递函数为所有输入 - 输出对，或各自具有不同的传递函数。选择输入和输出领域的输入和输出通道，字段只有当有多个输入或输出出现。
3、在模型传递函数区，可以使用下列选项指定模型的结构：

• 在Poles下面的选择框可以选择极点的个数，边上可以选择极点都是实根或者都是复根
• 选择Zero复选框则传递函数模型包含一个零点
• 选择Delay复选框则传递函数模型包含一个滞后环节
• 选择Integrator复选框则传递函数模型包含一个纯积分环节

4、在Initial Guess区域，选择Auto-selected自动抉择来计算估计的初始参数值。在参数表中的Initial Guess则显示为Auto。如果没有较好的初始猜测值，自动工作要优于特设值。

5、如果我们大概知道参数的值，可以在initial Guess列中输入此值，估计算法就会使用该值作为我们的起始点。类似的，如果我们完全知道该参数值，就在initial Guess列中输入该值并勾选Known复选框来修正该值。如果我们知道该参数的可能值所在的范围，那么我们可以在相应的Bounds字段中输入此范围，使得算法的搜索范围缩小，加速算法估计。
6、在Disturbance Model列表，我们可以选择一个可用选项

7、在Focus列表，选择如何衡量不同频率下拟合的相对重要性
8、在初始状态列表中，指定算法如何处理初始状态。
9、在协方差列表中，如果希望算法计算参数不确定度，请选择估计。 这些不确定性的影响在图上显示为模型置信区域。如果要省略估计不确定度，请选择无。 跳过不确定性计算可能会减少复杂模型和大型数据集的计算时间。
10、在“Model Name”字段中，编辑模型的名称或保留默认值。模型的名称在模型板中应该是唯一的。
11、要查看估计进度，请选中显示进度“Display progress”复选框。这将打开一个进度查看器窗口，其中报告估计进度。
12、点击“Regularization”正则化获得模型参数的正则化估计。在“正则化选项”对话框中指定正则化常量。
13、单击“Estimate”估计将此模型添加到系统识别应用程序中的模型板。
14、要在当前迭代完成后停止搜索并保存结果，请单击停止迭代“Stop Iterations”。要从当前模型继续迭代，请单击“Continue”继续按钮将当前参数值分配为下一次搜索的初始猜测。
得到模型参数估计后：
1、通过在“系统识别”应用程序的“Model Views模型视图”区域中选择相应的复选框来验证模型。
2、通过单击Value —> Initial Guess (值->初始猜测)按钮来优化模型，将当前参数值分配为下一次搜索的初始猜测，编辑名称字段，然后单击估计。
3、将模型导出到MATLAB工作区，以进一步分析，将其拖动到系统标识应用程序中的工作区矩形。

5.2使用system identification command line辨识过程模型

clear all;close all;
N=1000; %设置试验数据长度
A=[1,-1.5, 0.7]; B=[0,1,0.5]; C=[1,0.5];D=[1,0.5]; %模型参数
Model=idpoly(A,B); %理想系统模型
% figure(1);step(Model); axis([0 100 0 10]);grid; %绘制阶跃响应曲线
%产生输入输出数据
Model_Wnoise=idpoly(A, B ,1); %模型中加入白噪声
U=iddata([],idinput(N,'prbs')); %伪随机序列
E=iddata([],idinput(N,'rgs')); %白噪声序列
Y=sim(Model_Wnoise, [U,E]); %产生输出数据
Model_Noise1=idpoly(A,B,C); %有色噪声模型 1
Y1=sim(Model_Noise1, [U,E]);
Model_Noise2=idpoly(A,B,1,D); %有色噪声模型 2
Y2=sim(Model_Noise2, [U,E]);
data=iddata(Y,U); %输入输出数据组
opt = procestOptions;
opt.Display='on';
opt.DisturbanceModel='ARMA2';
opt.SearchMethod='lm';
%搜索方式有'lm','au','gr','gn'
sys=procest(data,'p1d',opt);
compare(data,sys);

该工具箱过程模型辨识模块中可以采用3种搜索方式：
‘lm’---Levenberg-Marquardt line search
‘gr’---Gradient-descent line search
‘gn’---Gauss-Netwon line search
‘au’---Nonlinear least squares with automatically chosen line search method

---

---