基于最小二乘法的阶梯式双容水箱系统辨识

前言

这个是本科时候做的一个小小课设（轻喷），不是很复杂，主要是利用最小二乘法、带遗忘因子的最小二乘法和增广最小二乘法来拟合阶梯式双容水箱水位与进水量之间的关系，在Simulink中建模，绘制了一个App Designer界面~

最小二乘法

最小二乘法是个有趣的方法，我记得初中的时候就提到了最小二乘法，那个时候针对一个拟合直线$y=kx+b$给出的公式是：
$$k=\frac{\sum xy-n\bar{x}\bar{y}}{\sum x^2-n{\bar{x}}^2}$$
我记得也没有很深入的解释，只是从平均之类的角度讲了讲，不过到了系统辨识之中，又对最小二乘法有了更深一步的学习。

基础最小二乘法

假设变量y与n维变量$X=[x_1,x_2,...,x_n{]}^T$有线性关系：$y={\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n=X^T\Theta$
我们现在拥有N个时刻的观测值 Y N = { y ( k ) } Y_N=\{y(k)\} YN​={y(k)}，则有如下的关系：
$$\begin{gathered} Y_N=[y(1),y(2),...,y(N){]}^T,X_N=[X^T(1),X^T(2),...,X^T(N){]}^T={\left[\begin{array}{ccc}{x}_1(1)& \dots & x_n(1)\\ ⋮& & ⋮\\ x_1(N)& \dots & x_n(N)\end{array}\right]}_{N\times n} \\ \Theta ={\left[\begin{array}{cccc}{\theta }_1& {\theta }_2& \dots & {\theta }_n\end{array}\right]}^T,Y_N=X_N\Theta \end{gathered}$$
推导的话网上有很多~我就不献丑了。
这时存在几种情况：
①$N=n$且$X_N^{-1}$存在时，$\Theta =X_N^{-1}{Y}_N$存在且唯一；
②$N>n$且$(X_N^T{X}_N{)}^{-1}$存在时，$\Theta =(X_N^T{X}_N{)}^{-1}{X}_N^T{Y}_N$；
上述两种情况都是精准的等式，而实际过程中会存在很多噪声以及建模的误差。
来一个实际的例子计算一下：

X=[0.3,1.1,2,2.9],Y=[0.1,1.2,4,8.4]，假设存在$y=a{x}^2+bx+c$，求出这三个参数：
设$X=[x^2,x,1{]}^T,\theta =[a,b,c{]}^T,y=X^T\theta$
有：$\left[\begin{array}{c}0.1\\ 1.2\\ 4\\ 8.4\end{array}\right]=\left[\begin{array}{ccc}0.09& 0.3& 1\\ 1.21& 1.1& 1\\ 4& 2& 1\\ 8.41& 2.9& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{c}0.09a+0.3b+c\\ 1.21a+1.1b+c\\ 4a+2b+c\\ 8.41a+2.9b+c\end{array}\right]$
$\theta ={\left(\left[\begin{array}{cccc}0.09& 1.21& 4& 8.41\\ 0.3& 1.1& 2& 2.9\\ 1& 1& 1& 1\end{array}\right]\left[\begin{array}{ccc}0.09& 0.3& 1\\ 1.21& 1.1& 1\\ 4& 2& 1\\ 8.41& 2.9& 1\end{array}\right]\right)}^{-1}\left[\begin{array}{cccc}0.09& 1.21& 4& 8.41\\ 0.3& 1.1& 2& 2.9\\ 1& 1& 1& 1\end{array}\right]\left[\begin{array}{c}0.1\\ 1.2\\ 4\\ 8.4\end{array}\right]=\left[\begin{array}{c}1.0032\\ -0.016\\ 0.0116\end{array}\right]$

带遗忘因子的最小二乘法

加权最小二乘法

在说带遗忘因子的最小二乘法之前，我们说一下加权最小二乘法，其实加权最小二乘法就是在我们求解中增加了一个权重矩阵，这个权重矩阵告诉了我们每一个观测值的重要性，即：
$${\hat{\Theta }}_{WLS}={\left(X_N^{T}W{X}_N\right)}^{-1}{X}_N^T{Y}_N,W={\left[\begin{array}{cccc}{\alpha }_1& 0& \dots & 0\\ 0& {\alpha }_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_N\end{array}\right]}_{N\times N}$$
显然当权值矩阵为单位阵时，退化为基础的最小二乘。
那什么是带遗忘因子的最小二乘法呢？我们还要说另一个最小二乘法的成员——递推最小二乘法。

递推最小二乘法

在工业现场中，我们往往需要对数据进行在线的分析计算，而且实际的数据是非常庞大的，如果都离线运算，一次性计算大量的矩阵逆，这非常的消耗计算资源，因此我们需要一种在线的算法，也就是递推最小二乘法。
递推最小二乘法利用上一时刻的估计值加上一个修正项来得到当前时刻新的估计值，这个推导过程也比较繁琐，我这里也直接给出结果：
对于：$$A(y^{-1})y(k)=B(y^{-1})u(k)\Rightarrow y(k)=h^T(k)\Theta$$
其中：$$\begin{gathered} h(k)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m){]}^T \\ \Theta =[a_1,...,a_n,b_1,...,b_m{]}^T \end{gathered}$$
则递推最小二乘的公式为：$$\begin{gathered} \hat{\theta }(k)=\hat{\theta }(k-1)+K(k)\left[y(k)-h^T(k)\hat{\theta }(k-1)\right] \\ K(k)=\frac{P(k-1)h(k)}{{\alpha }^{-1}(k)+h^T(k)P(k-1)h(k)} \\ P(k)=\left(I-K(k)h^T(k)\right)P(k-1) \end{gathered}$$
其中$P(k)$初值赋为较大的值（单位阵乘以较大的数），$\theta$赋初值为较小的值。

带遗忘因子的最小二乘法

终于我们讲到了带遗忘因子的最小二乘法，其实也很简单，就是在实时计算中，我们附加了每一个观测值一个权值，而这个权值与我们运行的时长是相关的，这是为了让我们不要受到太久远之前数值的影响，而主要分析比较新的数值的作用，即将我们的权重矩阵改为：
$$\begin{gathered} W=\left[\begin{array}{cccc}{\beta }^{2(N-1})& & & \\ & {\beta }^{2(N-2)}& & \\ & & \ddots & \\ & & & 1\end{array}\right]=\left[\begin{array}{cccc}{\alpha }^{N-1}& & & \\ & {\alpha }^{N-2}& & \\ & & \ddots & \\ & & & 1\end{array}\right] \\ \alpha (k)=\frac{1}{\alpha }\alpha (k-1)\in (0,1) \end{gathered}$$
结合之前的公式，越早的数据乘以更高次幂的遗忘因子，也就会变得更小。这样让我们可以逐渐减少历史数据的影响。

增广最小二乘法

正如之前所说，我们上述的最小二乘法都是基于精确的等式的，而忽略了噪声的影响，即实际的公式应为：
$$A(y^{-1})y(k)=B(y^{-1})u(k)+C(y^{-1})e(k)$$
其中$C(y^{-1})e(k)$为有色噪声，$e(k)$为白噪声。
经过一系列调整后，我们得到增广最小二乘的公式：
$$\begin{gathered} \hat{\theta }(k)=\hat{\theta }(k-1)+K(k)\left[y(k)-{\hat{h}}^T(k)\hat{\theta }(k-1)\right] \\ K(k)=\frac{P(k-1)h(k)}{{\alpha }^{-1}(k)+h^T(k)P(k-1)h(k)} \\ P(k)=\left(I-K(k)h^T(k)\right)P(k-1) \end{gathered}$$
其中$\hat{h}(k)=h(k)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),\hat{e}(k-1),...,\hat{e}(k-s){]}^T$，即我们将噪声模型的参数也考虑了进去一并进行估计，因此称之为增广最小二乘法。

至此，概念与最小二乘的内容就告一段落了。

阶梯式双容水箱建模

数学建模

如图为双容水箱的结构简图：

这个结构图比较简单，其中的一些参数：

$H$表示扬程，单位$m$，当然这个是简化了的，即相对于基准线的高度差；
$h_1$和$h_2$分别表示对应水箱的液位，单位$m$；
$R_1$,$R_2$,$R_3$分别表示液阻，单位$s/m^2$；
$A_1$和$A_2$分别表示两个水箱的底面积，单位$m^2$，这里把两个水箱都视为立方体。

根据双容水箱特性进行机理建模：
第一个水箱：$$动态流量=输入流量-输出流量\Rightarrow \frac{d{h}_1}{dt}{A}_1=\frac{H-h_1}{R_1}-\frac{h_1}{R_2}$$
第二个水箱：$$动态流量=输入流量-输出流量\Rightarrow \frac{d{h}_2}{dt}{A}_2=\frac{h1}{R_2}-\frac{h_2}{R_3}$$
看这两个式子是不是很好理解？停留在中间的等于进入的减出去的即可。
（顺带一提，这也对于小时候的一个问题：一个水池又进水又放水。在控制中，这个问题很有意义，因为在实际工况中我们需要保证液位稳定在一个范围之内才能保证正常的运行。）

回归正题，为了获得我们需要的下位水箱高度和输入之间的关系，我们需要联立两个式子得到如下的式子：
$${\ddot{h}}_2=\frac{1}{R_1{R}_2{R}_3}\left[R_{3}H-(R_1+R_2)h_2-(R_2{R}_3{A}_2+R_1{R}_3{A}_2+R_1{R}_2{A}_1){\dot{h}}_2\right]$$
但在计算机处理时，我们实际面对的是离散的信号，因此需要把上式离散化，得到我们真正应用的模型：
$$h_2(k)=-a_1{h}_2(k-1)-a_2{h}_2(k-2)+bH(k-2)$$
其中：
$$\begin{gathered} a_1=\frac{1}{K}(R_2{R}_3{A}_2+R_1{R}_3{A}_2+R_1{R}_2{A}_1-2K)h_2(k-1) \\ {a}_2=\frac{1}{K}(R_1+R_2-R_2{R}_3{A}_2-R_1{R}_3{A}_2-R_1{R}_{2}A1+K)h_2(k-1) \\ b=\frac{R_3}{K} \\ K=R_1{R}_2{R}_3{A}_1{A}_2 \end{gathered}$$
推导花了好久，总是在小地方出错…感兴趣的朋友可以一起推一下。

至此，我们完成了对阶梯式双容水箱的数学建模。

Simulink建模

如图为双容水箱的Simulink建模，这里输入DeltaH即为输入水量/扬程，中间分为两个水箱，并且要注意，水箱的高度是有上限的，这要求我们增加一个限幅器，保证水箱内液位不高于最高位置，也与后续我们的识别有关。

辨识目标

本设计的目标是通过最小二乘法求解出双容水箱的三个液阻值，通过数学机理部分的公式可以求解，但是三个方程中存在高次项，存在复数解，利用matlab的solve函数会因为极小的误差而产生难以想象的错解，因此为了避免该情况，我们可以通过先辨识上位水箱再辨识整体二阶水箱的方法，利用上位水箱的参数再来求解二阶水箱的参数：
$$\frac{d{h}_1}{dt}{A}_1=\frac{H-h_1}{R_1}-\frac{h_1}{R_2}$$
转化为差分方程：
$$h(k)=-(\frac{1}{R_1}+\frac{1}{R_2})h(k-1)+\frac{1}{R_1}H$$
此时差分方程的参数：
$$a=-(\frac{1}{R_1}+\frac{1}{R_2}),b=\frac{1}{R_1}$$
这样求解出a和b两个数值且是唯一解，再应用到双容水箱模型中即可。

仿真实验

操作界面

如图为我画的一个操作界面，里面主要设置了一些仿真的参数，假设我们的水箱底面积是已知的，我们的识别对象是参数$a_1,a_2,b,c_1,c_2$，也就是水箱的三个液阻以及有色噪声的两个系数。

实验一 无噪声

	扬程	液阻1	液阻2	液阻3
数值	5	35	25	28

如上图所示，蓝色曲线为上位水箱液位，红色曲线为下位水箱液位。可以观察到，由于液阻1较大，液阻2较小导致稳定后水箱1中的液位低于水箱2的液位。
此时辨识的上位水箱的参数为：

其中R1=34.97，较实际设定值35有0.085%的误差，可以接受；R2=25，与设定值相同。此时得到两个液阻阻值，可以进行下一步辨识。

这样我们再根据辨识的参数来计算液阻3：

如图所示，计算得到的液阻R3=28.02，与实际设定值28有0.07%的误差。

实验二 有噪声+溢出情况

	扬程	液阻1	液阻2	液阻3	噪声1	噪声2
数值	5	20	20	20	1	2

如图所示，水箱由于高度只有2m，液位超过2m会直接进入饱和，即转折后的水平线。此时权值均为1的增广最小二乘法拟合的结果从形状上和误差上看都比较接近系统输出，但显然图像中丢失了由于噪声而引起的波动。
观察一下上位水箱：

发生了明显的强制转折，就是因为达到液位上限了。
看一下辨识的参数：

跟我们预设的参数不能说一模一样，只能说毫不相关！其实我们只是拟合出了发生饱和的水箱模型而非真正的模型。

我们调整遗忘因子的大小，请注意！$这里的遗忘因子我们选择超过1的数值$，这是因为后续的饱和值与液阻大小没有关系，是因为水位超过水箱高度了，因此与之前所讲的遗忘因子不同，我们反而需要加大过去数据对结果的影响，即增大历史数据的权重，而减小新数据的权重，这样我们得到如下结果：

注意我们拟合的结果，液位高度高于水箱高度的2m，并且存在肉眼可见的噪声波动，再来看看液阻参数的辨识结果：

这样，我们得到的液阻计算结果与我们预先设定的比较接近啦。

总结

通过带遗忘因子的增广最小二乘法，我们实现了对双容水箱模型参数的辨识，并且根据最终结果可以看出，其识别效果较好，数据比较贴合实际情况，但毕竟本次设计是基于Simulink仿真，少了很多实际模型中的干扰、硬件特性等情况，因此如果要针对现实情况中一个双容水箱的辨识，还需要进一步分析机理、分析干扰才能更好地辨识系统。
希望和大家多多学习！
ps：分享一下本科做课设的经验，高等学府我不敢说，都是大佬，但是一般来说课设不会非常难，它的目的还是为了能让我们把课上学的知识能够实际应用一下，原来我也对这些很苦恼很抓狂，但是真的，有一些想法，一些背景还是很好做的。