参数辨识示例 报告

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

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

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

虽然目前参数辨识的领域己经发展了多种算法，但是用于气动参数估计的算法主要有：极大似然法(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 --和是多项式，定义为：12111212z ...z ...n n mm 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),...,()Tn m Ta a ab 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),...,()TL TL 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 ----⎡⎤⎡⎤⎢⎥⎢⎥----⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥------⎣⎦⎣⎦LL L LM M MM 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.0211.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)]Tk k K k y k k k θθϕθ=-+--()[()()](1)T P k I K k k P k ϕ=--(1)()()1()(1)()TP 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 仿真图像输出 endtheta1 %输出参数估计值plot([1:20],theta2) %输出参数逐步递推估计的轨迹图像 xlabel('k'); %设置横坐标为步长k ylabel('参数估计a ，b ，c'); %纵坐标为估计参数a ，b ，c legend('a','b','c'); %标示相应曲线对应的参数得到结果：参数估计变化曲线：2468101214161820-6-4-2246k参数估计a ，b ，ca b c可以计算得到各辨识参数误差分别为：2.47%，0.438%，3.74%。