参数辨识示例 报告

参数辨识
参数辨识的步骤
飞行器气动参数辨识是一个系统工程，包括四部分：①试验设计，使试验能为辨识提供含有足够信息量且信息分布均匀的试验数据；②气动模型结果确定，即从候选模型集中，根据一定的准则和经验，选出最优的气动模型构式；③气动参数辨识，根据辨识准则和数据求取模型中待定参数，这是气动辨识定量研究的核心阶段；④模型检验，确认所得气动模型是否确实反映了飞行器动力学系统中气动力的本质属性。

这四个部分环环相扣，缺一不可，要反复进行，直到对所得气动模型满意为止。

参数辨识的方法
参数辨识方法主要有最小二乘算法、极大似然法、集员辨识法、贝叶斯法、岭估计法、超椭球法和鲁棒辨识法等多种辨识方法。

虽然目前参数辨识的领域己经发展了多种算法，但是用于气动参数估计的算法主要有：极大似然法(ML)，广义Kalman滤波(EKF)法，模型估计法(EBM )、分割及多分割算法(PIA及MPIA)、最小二乘法，微分动态规划法等。

因为最小二乘法和极大似然法是两种经典的算法，目前己经发展得相当成熟。

最小二乘法适于线性模型的参数辨识，可以用于飞行器系统辨识中很多的线性模型，如惯性仪表误差系数的辨识，线性时变离散系统初始状态的辨识及多项式曲线拟合等。

目前最小二乘法已经广泛应用于工程实际中。

而极大似然算法因其具有渐进一致性、估计的无偏性、良好的收敛特性等特点而被广泛应用于飞行器参数辨识领域。

最小二乘法大约是1975年高斯在其著名的星体运动轨道预报研究工作中提出来的。

后来，最小二乘法就成了估计理论的奠基石。

由于最小二乘法原理简单，编程容易，所以它颇受人们重视，应用相当广泛。

极大似然估计算法在实践中不断地被加以改进，这种改进主要表现在三个方
面，即算法的优化、是否考虑过程噪声以及灵敏度的计算等。

极大似然法具有集中处理数据的特点，即对一段时间历程上的数据集中分析，直接得到所需的气动参数，此方法也具有较高的精度，其估计具有一致性和无偏性。

但是，它也有缺点，主要是计算所需时间较长，易于进行集中处理而不适合在线实时处理。

递推最小二乘法
设时不变动态过程的数学模型为
11()()()()()A z z k B z u k v k --=+
式中()u k 和()z k 为系统的输入量和输出量；()v k 是噪声；11()()A z B z --和是多项式，定义为：
1211
1
2
12z ...z ...n n m
m A a z a z B b z b z b z
-------+++++（）=1+（）=
现在问题是如何利用系统的{}(),();1,...,u k z k k L =序列，估计多项式
11()()A z B z --和中的未知系数。

定义：
[]
[]
1212,,...,,,,...,()(1),...,(),(1),...,()T
n m T
a a a
b b b h k z k z k n u k u k m θ==------
则可将数学模型化为一个标准的最小二乘格式
()()()T z k h k v k θ=+ 1,2,...,k L =
上式可写成一个线性方程组，形式如下：
L l l Z H V =+
其中
[][]
(1),(2),...,()(1),(2),...,()T
L T
L Z z z z L V v v v L ==
12(0)(1)(0)(1)(1)(1)(2)
(1)
(2)(2)(1)()(1)()()L T z z n u u m h z z n u u m h H z L z L n u L u L m h L ----⎡⎤⎡⎤⎢⎥⎢
⎥----⎢⎥⎢⎥==
⎢⎥⎢⎥⎢⎥⎢⎥
------⎣⎦
⎣⎦L
L L L
M M M
M L L
最小二乘的准则函数()J θ可取为：
()()()T L l L l J Z H Z H θθθ=--
极小化()J θ，可求出系统系数θ的估计值^
θ。

下面给出最小二乘估计一般估值计算公式：
^
11()T l l l L H H H z θ-=
对应的方法叫做最小二乘法。

随着待估计的参数和观测数据的增多，最小二乘一次完成算法的计算量增加计算困难不利于在线辨识。

现在把它化成递推算法，其格式为：
新的估计值^
()k θ=老的估计值^
(1)k θ-+修正项
这样的优点是：
(1)每一步的计算量比较小，因而能够用比较少的计算量完成较大的任务。

(2)具有跟踪时变参数的能力，能辨识含有时变参数的系数模型。

(3)在自适应控制系统和基本诊断目的的故障检测中，需要递推算法。

(4)在参数估计达到给定的精确度时算法可以给出收敛中止的判据。

示例：
一模型形式为
1122y a u a u e =++
取得3组实验数据为
12():1,1,3;():1,2,1;():0,1,2y k x k x k -
求出参数12,a a 的最小二乘估计。

该模型中，1210(,),(1,1,3),2112T T a a y θ⎡⎤
⎢⎥==Φ=-⎢⎥⎢⎥⎣⎦
，所以 606,055T T y ⎛⎫⎛⎫
ΦΦ=Φ= ⎪ ⎪⎝⎭⎝⎭
解正则方程组
T T y θΦΦ=Φ
即
12606055a a ⎛⎫⎛⎫⎛⎫= ⎪ ⎪ ⎪⎝⎭⎝⎭
⎝⎭ 可得参数的最小二乘估计为：121,1a a ==。

MATLAB 算例
通过上例我们可初步了解应用最小二乘法进行了解，下面通过MATLAB 编程实现递推最小二乘法求解系统参数辨识，这里选用单输入单输出系统
2y au bu c =++作为待辨识模型，要求对,,a b c 进行参数辨识，首先需要取得若干
组输入输出数据，故我们事先选定实际系统为20.2345 4.5678 5.1234y u u =-+，通过对系统进行测量我们得到以下初始数据。

[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]u =
[ 0.7901,-3.0742,-6.4695,-9.3958,-11.8531,-13.8414,-15.3607,-16.4110,-16.9923,-17.1046,
-16.7479,-15.9222,-14.6275,-12.8638,-10.6311,-7.9294,-4.7587,-1.1190, 2.9897,7.5674]
y =加入随机误差，我们得到将要使用的初始数据。

u=[1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10.02
11.02 12.02 13.02 14.02 15.02 16.03 17.02 18.01 19.02 20.04];
y=[0.791 -3.072 -6.465 -9.398 -11.83 -13.414 -15.607 -16.10 -16.23 -17.046
-16.479 -15.222 -14.275 -12.638 -10.311 -7.294 -4.87 -1.90 2.97 7.74]
对该模型参数辨识采用递推最小二乘法，即RLS （ recurisive least square ），它是一种能够对模型参数进行在线实时估计的辨识方法。

其基本思想可以概括为：新的估计值^
()k θ=旧的估计值^
(1)k θ-+修正项
下面将批处理最小二乘法改写为递推形式即递推最小二乘参数估计的计算方法。

批处理最小二乘估计^
θ为^
1()T T Y θ-=ΦΦΦ，设k 时刻的批处理最小二乘估计为：
^
1()T T k k k k Y θ-=ΦΦΦ
令111()()[(1)()(1)]T T k k P k P k k k φφ---=ΦΦ=-+- K 时刻的最小二乘估计可以表示为
11ˆˆˆ()()()[()()](1)()[()()(1)]T T T k k k k k P k Y P k Y k y k k K k y k k k θφθφθ--=Φ=Φ+=-+--；式中()()()K k P k k φ=，因为要推导出()P k 和()K k 的递推方程，因此这里使用矩阵求逆引理：设1,(),A A BC I CA B -++均为非奇异方阵，则
111111()()A BC A A B I CA B CA ------+=-+通过运用矩阵求逆引理，把复杂的矩阵求
逆转化为标量求倒数，大大减小了计算量。

1N P +与N P 间的递推关系。

最终得到递推最小二乘参数递推估计公式如下：
^^
^
()(1)()[()()(1)]T
k k K k y k k k θθϕθ=-+--
()[()()](1)T P k I K k k P k ϕ=--
(1)()
()1()(1)()
T
P k k K k k P k k ϕϕϕ-=
+- 程序流程图
Step 1 ：设置初值)0(ˆθ
和P(0)，输入初始数据； Step2 ：采样当前输出y(k)、和输入u(k)
Step3 ：利用递推公式计算 )(k K 、)(ˆk θ和)(k P ； Step4 ：k=k+1，返回step2，继续循环。

程序代码：
X=[1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10.02 11.02 12.02 13.02 14.02 15.02 16.03 17.02 18.01 19.02 20.04];
Y=[0.791 -3.072 -6.465 -9.398 -11.83 -13.414 -15.607 -16.10 -16.23 -17.046 -16.479 -15.222 -14.275 -12.638 -10.311 -7.294 -4.87 -1.90 2.97 7.74];
%实验输入数据、实验输出数据 syms a b c % 定义待辨识参数
theta=[a;b;c]; %theta 包含待辨识参数a ，b ，c theta1=zeros(3,1); %对象参数初始化 P=10^6*eye(3); %构造初始P 阵 for k=1:20 %仿真步长范围1到20 phi=[X(k)*X(k);X(k);1];
%y=aX*X+bX+c=phi'*theta
%theta=[a ；b ；c]；phi=[X(k)*X(k);X(k);1] K=P*phi/(1+phi'*P*phi); %递推最小二乘法K 阵的递推公式 theta=theta1+K*(Y(k)-phi'*theta1); %theta 的递推公式 P=(eye(3)-K*phi')*P; %递推最小二乘法P 阵的递推公式 theta1=theta; %theta 的最终估计向量 theta2(:,k)=theta %theta 估计向量矩阵化，目的是为了 %下面的plot 仿真图像输出 end
theta1 %输出参数估计值
plot([1:20],theta2) %输出参数逐步递推估计的轨迹图像 xlabel('k'); %设置横坐标为步长k ylabel('参数估计a ，b ，c'); %纵坐标为估计参数a ，b ，c legend('a','b','c'); %标示相应曲线对应的参数
得到结果：
参数估计变化曲线：
2468
10
1214161820
-6-4
-2
2
4
6
k
参数估计a ，b ，c
a b c
可以计算得到各辨识参数误差分别为：2.47%，0.438%，3.74%。

考虑在初始参数中引入了随机误差，故模型参数的最小二乘估计在误差允许的范围内基本符合要求。

2468
101214161820
-20-15
-10
-5
5
10
x
y
上图为辨识出的系统与真实系统曲线比对，可观察到曲线拟合程度较好，说明此方法对系统的辨识结果是可以接受的。

算例扩展
通过上例我们对应用递推最小二乘法进行系统参数辨识有了一定了解，下面对此算例进行扩展，选用多输入单输出系统2112y au bu cu d =+++作为待辨识模型，要求对,,,a b c d 进行参数辨识，首先需要取得若干组输入输出数据，故我们事先选定实际系统为2112234y u u u =-+-，通过对系统进行测量并加入随机误差我们得到以下初始数据。

1u =[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]; 2u =[1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
y=[-2.01 -4.06 2.07 3.91 14.02 19.91 34.01 44.04 62.08 76.04 98.92 115.943 142.98 164 194.1 220.91 253.88 284.1 322.91 356]
对该模型参数辨识采用递推最小二乘法。

得到结果
theta1 =
0.9984 -1.9391 3.2436 -4.3414
参数估计变化曲线：
2468
101214161820
-5-4-3-2-101234k
参数估计a ，b ，c ，d
可以计算得到各辨识参数误差分别为：0.16%，3.045%，8.12%，8.535%。

考虑在初始参数中引入了随机误差，故模型参数的最小二乘估计在误差允许的范围内基本符合要求。

由此可见，递推最小二乘法在单输入单输出和多输入单输出系统的参数辨识中，均可达到较好的辨识效果。

本周总结
本周对系统辨识方法进行了更深入的学习，使用了递推最小二乘法进行了系统参数辨识的学习应用，并编写了MATLAB程序对实例系统进行了辨识，得到了较为满意的结果。

不过也有些许不足：参数变化曲线不够平滑，分析原因可能是因为我引入的是随机误差，而不是高斯白噪声，对结果有了一定影响，也有可能是采样数据个数的限制，只有20个输入输出数据对于递推最小二乘法来说有些不足，如果要是能扩展到100个输入输出数据时，相信辨识精度会有进一步提高。

在本周学习中，也发现了自己在基础知识方面的欠缺，因为本科时候学习的专业是自动化，对飞行器基础知识的了解有限，所以这次系统辨识的模型没有直接选用飞行器模型。

并且发现自己对于MATLAB，C语言等工具的使用不够熟练，还需要多多做些练习来进一步熟悉。

而在下一步的学习中，我将把重点放在飞行器基础知识的学习，并编写对飞行器气动模型参数辨识的程序，学习有关闭环系统辨识的知识。