系统辨识和最小二乘参数估计Matlab仿真
一、系统辨识
在控制系统的分析中，首先要建立系统的数学模型，控制系统的数学模型是定量描述系统或过程内部物理量（或变量）之间关系的数学表达式。
一般来说，建立控制系统数学模型有两种基本方法：
（1）机理建模（白箱模型）：即根据系统内在运行机制、物料和能量守恒等物理学、化学规律建立系统的数学模型，一般步骤如下：
Step1：根据系统工作原理及其在控制系统中的作用，确定输入和输出；
Step2：根据物料和能量守恒等关系列写基本方程式；
Step3：消去中间量；
Step4：获得系统模型；
（2）实验法建模（黑箱模型）：即对于机理尚不清楚或机理过于复杂的系统，可以人为的对其施加某种测试信号，并记录其输出响应，或者记录正常运行时的输入输出数据，然后利用这些输入输出数据确定系统模型结构和参数。
多年来，系统辨识已经发展为一门独立学科分支，通过系统辨识建立一个对象的数学模型，通常包括两方面的工作：一是模型结构的确定（模型的类型、阶次），二是模型参数估计。
根据时间是否连续，参数模型又可以分为连续时间系统模型和离散时间系统参数模型，这两类模型均可采用输入输出模型和状态空间模型描述，离散系统采用差分方程描述，以单输入单输出（SISO）离散系统参数模型为例。
1.确定性模型
SISO系统确定性模型可表示为：
u(k)和y(k)分别为输入和输出，d为纯延时。
2.随机性模型
如果受到随机扰动，则式子可写为：
为系统随机扰动，其结构如图：
系统辨识的一般步骤如图：
从图中可以看出，利用辨识的方法建立系统数学模型，从实验设计到模型获得，需要这些步骤。
二、最小二乘参数估计
1.批处理最小二乘
考虑以下CAR模型：
为白噪声，结构参数na、nb和d已知，参数估计的任务就是根据可测量
的输入输出，确定如下参数：
仿真实例：
式中，为方差为1的白噪声，选用幅值为1的逆M序列作为输入，LS 算法进行参数估计，仿真结果如图：
仿真程序(Matlab):
%批处理最小二乘参数估计（LS）
clear all;
a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数
na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次
L=500; %数据长度
uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)
yk=zeros(na,1); %输出初值
x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值
xi=randn(L,1); %白噪声序列
theta=[a(2:na+1);b]; %对象参数真值
for k=1:L
phi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi矩阵
y(k)=phi(k,:)*theta+xi(k); %采集输出数据
IM=xor(S,x4); %产生逆M序列
if IM==0
u(k)=-1;
else
u(k)=1;
end
S=not(S); M=xor(x3,x4); %产生M序列
%更新数据
x4=x3; x3=x2; x2=x1; x1=M;
for i=d+nb:-1:2
uk(i)=uk(i-1);
end
uk(1)=u(k);
for i=na:-1:2
yk(i)=yk(i-1);
end
yk(1)=y(k);
end
thetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae
2.递推最小二乘
在批处理最小二乘法时，由于每次处理的数据量较大，而且常常要求对象参数能够在线实时估计，解决的方法就是将其化成递推算法，其基本思想为：
算法介绍：
仿真实例:
式中，为方差为0.1的白噪声，取初值，选择方差为1的白噪声作为输入信号u(k)，采用RLS算法进行参数估计，仿真结果如图：
仿真程序（Matlab）：
%递推最小二乘参数估计（RLS）
clear all; close all;
a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数
na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=400; %仿真长度
uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)
yk=zeros(na,1); %输出初值
u=randn(L,1); %输入采用白噪声序列
xi=sqrt(0.1)*randn(L,1); %白噪声序列
theta=[a(2:na+1);b]; %对象参数真值
thetae_1=zeros(na+nb+1,1); %thetae初值
P=10^6*eye(na+nb+1);
for k=1:L
phi=[-yk;uk(d:d+nb)]; %此处phi为列向量
y(k)=phi'*theta+xi(k); %采集输出数据
%递推最小二乘法
K=P*phi/(1+phi'*P*phi);
thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);
P=(eye(na+nb+1)-K*phi')*P;
%更新数据
thetae_1=thetae(:,k);
for i=d+nb:-1:2
uk(i)=uk(i-1);
end
uk(1)=u(k);
for i=na:-1:2
yk(i)=yk(i-1);
end
yk(1)=y(k);
end
plot([1:L],thetae); %line([1,L],[theta,theta]); xlabel('k'); ylabel('参数估计a、b');
legend('a_1','a_2','b_0','b_1'); axis([0 L -2 2]);
系统辨识和最小二乘参数估计Matlab仿真 一、系统辨识 在控制系统的分析中，首先要建立系统的数学模型，控制系统的数学模型是定量描述系统或过程内部物理量（或变量）之间关系的数学表达式。 一般来说，建立控制系统数学模型有两种基本方法： （1）机理建模（白箱模型）：即根据系统内在运行机制、物料和能量守恒等物理学、化学规律建立系统的数学模型，一般步骤如下： Step1：根据系统工作原理及其在控制系统中的作用，确定输入和输出； Step2：根据物料和能量守恒等关系列写基本方程式； Step3：消去中间量； Step4：获得系统模型； （2）实验法建模（黑箱模型）：即对于机理尚不清楚或机理过于复杂的系统，可以人为的对其施加某种测试信号，并记录其输出响应，或者记录正常运行时的输入输出数据，然后利用这些输入输出数据确定系统模型结构和参数。 多年来，系统辨识已经发展为一门独立学科分支，通过系统辨识建立一个对象的数学模型，通常包括两方面的工作：一是模型结构的确定（模型的类型、阶次），二是模型参数估计。 根据时间是否连续，参数模型又可以分为连续时间系统模型和离散时间系统参数模型，这两类模型均可采用输入输出模型和状态空间模型描述，离散系统采用差分方程描述，以单输入单输出（SISO）离散系统参数模型为例。 1.确定性模型 SISO系统确定性模型可表示为： u(k)和y(k)分别为输入和输出，d为纯延时。 2.随机性模型 如果受到随机扰动，则式子可写为： 为系统随机扰动，其结构如图：
系统辨识的一般步骤如图： 从图中可以看出，利用辨识的方法建立系统数学模型，从实验设计到模型获得，需要这些步骤。 二、最小二乘参数估计 1.批处理最小二乘 考虑以下CAR模型：
用matlab实现最小二乘递推算法辨识系统参 数 自动化系统仿真实验室指导教师： 学生姓名班级计082-2 班学号撰写时间： 全文结束》》-3-1 成绩评定： 一．设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作； 3、了解并掌握Matlab中一些函数的作用与使用；二．设计要求最小二乘递推算法辨识系统参数，利用matlab编程实现，设初始参数为零。z(k)-1、5*z(k-1)+0、7*z(k-2)=1*u(k-1)+0、5*u(k-2)+v(k); 选择如下形式的辨识模型：z(k)+a1*z(k- 1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k);三．实验程序 m=3;N=100;uk=rand(1,N);for i=1:Nuk(i)=uk(i)*(-1)^(i-1);endyk=zeros(1,N); for k=3:N yk(k)=1、5*yk(k-1)-0、 7*yk(k-2)+uk(k-1)+0、5*uk(k-2); end%j=100;kn=0;%y=yk(m:j);%psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j- 2)];%pn=inv(psi*psi);%theta=(inv(psi*psi)*psi*y);theta=[0 ;0;0;0];pn=10^6*eye(4);for t=3:Nps=([yk(t-1);yk(t-
2);uk(t-1);uk(t-2)]);pn=pn- pn*ps*ps*pn*(inv(1+ps*pn*ps));theta=theta+pn*ps*(yk(t)-ps*theta);thet=theta;a1=thet(1);a2=thet(2);b1=thet(3);b2= thet(4); a1t(t)=a1;a2t(t)=a2;b1t(t)=b1;b2t(t)=b2;endt=1:N;plot(t,a 1t(t),t,a2t(t),t,b1t(t),t,b2t(t));text(20,1、 47,a1);text(20,-0、67,a2);text(20,0、97,b1);text(20,0、47,b2);四．设计实验结果及分析实验结果图：仿真结果表明，大约递推到第步时，参数辨识的结果基本到稳态状态，即a1=1、5999，b1=1，c1=0、5，d1=-0、7。五、设计感受这周的课程设计告一段落了，时间短暂，意义重大。通过这次次练习的机会，重新把matlab课本看了一遍，另外学习了系统辨识的有关内容，收获颇丰。对matlab的使用更加纯熟，也锻炼了自己在课本中搜索信息和知识的能力。在设计过程中虽然遇到了一些问题，但经过一次又一次的思考，一遍又一遍的检查终于找出了原因所在，也暴露出了前期我在这方面的知识欠缺和经验不足。同时我也进一步认识了matlab软件强大的功能。在以后的学习和工作中必定有很大的用处。
最小二乘法在系统辨识中的应用 王文进 控制科学与控制工程学院 控制理论与控制工程专业 2009010211 摘要：在实际的工程中，经常要对一个系统建立数学模型。很多时候，要面对一个未知的系统，对于这些未知系统，我们所知道的仅仅是它们的一些输入输出数据，我们要根据这些测量的输入输出数据，建立系统的数学模型。由此诞生了系统辨识这门科学，系统辨识就是研究怎样利用对未知系统的输入输出数据建立描述系统的数学模型的科学。系统辨识在工程中的应用非常广泛，系统辨识的方法有很多种，最小二乘法是一种应用及其广泛的系统辨识方法。本文主要讲述了最小二乘估计在系统辨识中的应用。 首先，为了便于介绍，用一个最基本的单输入单输出模型来引入系统辨识中的最小二乘估计。 例如：y = ax + （1） 其中：y、x 可测，为不可测的干扰项，a未知参数。通过N 次实验，得到测量数据y k和x k ，其中k=1、2、3、…，我们所需要做的就是通过这N次实验得到的数据，来确定未知参数a 。在忽略不可测干扰项的前提下，基本的思想就是要使观测点y k和由式（1）确定的估计点y的差的平方和达到最小。用公式表达出来就是要使J最小： 确定未知参数a的具体方法就是令： J a = 0 , 导出 a 通过上面最基本的单输入单输出模型，我们对系统辨识中的最小二乘法有了初步的了解，但在实际的工程中，系统一般为多输入系统，下面就用一个实际的例子来分析。在接下来的表述中，为了便于区分，向量均用带下划线的字母表示。 水泥在凝固过程中，由于发生了一系列的化学反应，会释放出一定的热量。若水泥成分及其组成比例不同，释放的热量也会不同。 水泥凝固放热量与水泥成分的关系模型如下： y = a0+ a1x1+…+ a n x n + 其中，y为水泥凝固时的放热量（卡/克）；x1~x2为水泥的几种成分。
最小二乘参数估计 摘要： 最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T l ΦΦΦ-∧=1θ。 最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。 关键词： 最小二乘（Least-squares ），系统辨识（System Identification ） 目录： 1.目的 (1) 2.设备 (1) 3引言 (1) 3.1 课题背景 (1) 4数学模型的结构辨识 (2) 5 程序 (3) 5.1 M 序列子函数 ................................................................................. 错误！未定义书签。 5.2主程序............................................................................................... 错误！未定义书签。 6实验结果： ................................................................................................................................... 3 7参考文献： ................................................................................................. 错误！未定义书签。 1.目的 1.1掌握系统辨识的理论、方法及应用 1.2熟练Matlab 下最小二乘法编程 1.3掌握M 序列产生方法 2.设备 PC 机1台（含Matlab 软件） 3引言 3.1 课题背景 最小二乘理论是有高斯（K.F.Gauss ）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。 最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最
最小二乘法拟合 在科学实验的统计方法研究中,往往要从一组实验数据中寻找出自变量x 和因变量y之间的函数关系y=f(x) 。由于观测数据往往不够准确，因此并 不要求y=f(x)经过所有的点,而只要求在给定点上误差 按照某种标准达到最小，通常采用欧氏范数作为误差量度的标准。这就是所谓的最小二乘法。在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行。 函数polyfit是指用一个多项式函数来对已知数据进行拟合，我们以下列数据为例介绍这个函数的用法： >> x=0:0.1:1； >> y=[ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2 ] 为了使用polyfit，首先必须指定我们希望以多少阶多项式对以上数据进行拟合，如果我们指定一阶多项式，结果为线性近似，通常称为线性回归。我们选择二阶多项式进行拟合。 >> P= polyfit (x, y, 2) P= -9.8108 20.1293 -0.0317 函数返回的是一个多项式系数的行向量，写成多项式形式为： 为了比较拟合结果，我们绘制两者的图形： >> xi=linspace (0, 1, 100); %绘图的X-轴数据。 >> Z=polyval (p, xi); %得到多项式在数据点处的值。 当然，我们也可以选择更高幂次的多项式进行拟合，如10阶： >> p=polyfit (x, y, 10); >> xi=linspace (0, 1,100);
>> z=ployval (p, xi); 读者可以上机绘图进行比较，曲线在数据点附近更加接近数据点的测量值了，但从整体上来说，曲线波动比较大，并不一定适合实际使用的需要，所以在进行高阶曲线拟合时，“越高越好”的观点不一定对的。
《系统辨识》基于MATLAB的最小二乘法（一阶）的仿真 clc clear % ①白噪声的生成过程如下：e=randn(1,500); e=e/std(e); e=e-mean(e); A=0; %白噪声的均值为0 B=sqrt(0.1); %白噪声的方差为0.1 e=A+B*e; %绘制白噪声图 k=1:500; subplot(4,1,1) %画四行一列图形窗口中的第一个图形 plot(k,e,'r'); xlabel('k'), ylabel('e');title('(0,1)均匀分布的随机序列') % ②生成M序列的过程如下：X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初始状态（0101）， Yi寄存器的各级输出 m=500; %M序列的总长度 for i=1:m Y4=X4; Y3=X3; Y2=X2; Y1=X1; X4=Y3; X3=Y2; X2=Y1; X1=xor(Y3,Y4); %异或运算 if Y4==0 U(i)=-1; else U(i)=Y4; end end M=U; u=U; %绘制M序列图? i1=i k=1:1:i1; subplot(4,1,2) %画四行一列图形窗口中的第二个图形 plot(k,U,k,U,'rx') stem(M) xlabel('k') ylabel('M序列') title('移位寄存器产生的M序列') % ③参数估计的过程如下： %绘制参数估计的相关图形 z=zeros(1,500); %定义输出观测值的长度 for k=2:500 z(k)=0.9*z(k-1)+u(k-1)+e(k);%用理想输出值作为观测值 end subplot(4,1,3) %画四行一列图形窗口中的第三个图形 i=1:1:500; %横坐标的范围从1到500，步长为1 plot(i,z) %图形的横坐标是采样时刻i，纵坐标是输出观测值Z， 图形格式为连续曲线
Matlab系统辨识尝试之详细过程1 前面介绍了Matlab系统辨识工具箱的一些用法，这里拿一个直观的例子来尝试工具箱的具体用法。比较长，给个简单目录吧： 1.辨识的准备 2.辨识数据结构的构造 3.GUI辨识 4.辨识效果 5.对固有频率的辨识 6.结构化辨识 7.灰箱辨识 8.加入kalman滤波的灰箱辨识 1.辨识的准备 在辨识前，首先要根据自己辨识的情况，确定要辨识的状态空间模型的一些特点，如连续还是离散的；有无直通 分量（即从输入直通到输出的分量）；输入延迟；初始状态等。了解了这些情况就可以更快速的配置辨识时的一些设 置选项。 2.辨识数据结构的构造 使用原始数据构造iddata结构： data=iddata(y,u,Ts); 这里以一个弹簧质量系统的仿真为例 代码如下，其中用到了函数MDOFSolve，这在之前的博文介绍过（https://www.doczj.com/doc/0310213856.html,/?p=183），拿来用即可。如果发现运行有错误，可以将MDOFSolve函数开头的一句 omega2=real(eval(omega2)); 注释掉。 %弹簧质量系统建模 clc clear close all m=200; k=980*1000;
c=1.5*1000; m1=1*m; m2=1.5*m; k1=1*k; k2=2*k; k3=k1; %%由振动力学知识求固有频率 M=[m10;0m2]; K=[k1+k2-k2;-k2k3+k2]; [omega,phi,phin]=MDOFSolve(M,K); fprintf('固有频率:%fHz\n',subs(omega/2/pi)); %%转化到状态空间 innum=2; outnum=2; statenum=4; A=[0100; -(k1+k2)/m10k2/m10; 0001; k2/m20-(k3+k2)/m20]; B=[00; 1/m10; 00; 01/m2]; C=[1000; 0010]; D=zeros(outnum,innum); K=zeros(statenum,innum); mcon=idss(A,B,C,D,K,'Ts',0);%连续时间模型 figure impulse(mcon) %%信号仿真，构造数据供辨识 n=511;%输入信号长度 Ts=0.001; t=0:Ts:(n-1)*Ts; u1=idinput(n,'prbs');%输入1为伪随机信号 u2=zeros(n,1);%输入2为空 u=[u1u2]; simdat=iddata([],u,Ts);%形成输入数据对象 e=randn(n,2)*1e-7; simopt=simOptions('AddNoise',true,'NoiseData',e);%添加噪声yn=sim(mcon,simdat,simopt);%加噪声仿真 y=sim(mcon,simdat);%无噪声仿真
例1、考虑仿真对象 )()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- 其中，)(k v 是服从正态分布的白噪声N )1,0(。输入信号采用4阶M 序列，幅度为1。选择如下形式的辨识模型 )()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ 设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LS θ?为LS θ?=(T T -ΦΦΦ1)z 。其中，被辨识参数LS θ?、观测矩阵Φ的表达式为: ????? ???????=2121?b b a a LS θ (3)(4)(16)z z z ??????=???? ??z (2)(1)(2)((3)(2)(3)(2)(15)(14)(15)(14)z z u u z z u u z z u u --????--??Φ=????--?? 程序框图如图1所示。Matlab 仿真程序如下： %二阶系统的最小二乘一次完成算法辨识程序，文件名：LS.m
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画输入信号u的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线 subplot(3,1,3) %画三行一列图形窗口中的第三个图形 stem(z),grid on %画出输出观测值z的径线图形，并显示坐标网格u,z %显示输入信号和输出观测信号 %L=14 %数据长度 HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵 赋值
【1】随机序列产生程序 【2】白噪声产生程序 【3】M序列产生程序 【4】二阶系统一次性完成最小二乘辨识程序 【5】实际压力系统的最小二乘辨识程序 【6】递推的最小二乘辨识程序 【7】增广的最小二乘辨识程序 【8】梯度校正的最小二乘辨识程序 【9】递推的极大似然辨识程序 【10】Bayes辨识程序 【11】改进的神经网络MBP算法对噪声系统辨识程序【12】多维非线性函数辨识程序的Matlab程序【13】模糊神经网络解耦Matlab程序 【14】F-检验法部分程序 【1】随机序列产生程序 A=6; x0=1;M=255; for k=1:100 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=v1; x0=x1; v0=v1; end v2=v k1=k; %grapher k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');title('(0,1)均匀分布的随机序列') 【2】白噪声产生程序 A=6; x0=1; M=255; f=2; N=100; for k=1:N x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*f; x0=x1;
v0=v1; end v2=v k1=k; %grapher k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');title('(-1,+1)均匀分布的白噪声') 【3】M序列产生程序 X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101），Yi为移位寄存器各级输出m=60; %置M序列总长度 for i=1:m %1# Y4=X4; Y3=X3; Y2=X2; Y1=X1; X4=Y3; X3=Y2; X2=Y1; X1=xor(Y3,Y4); %异或运算 if Y4==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M序列') title('移位寄存器产生的M序列') 【4】二阶系统一次性完成最小二乘辨识程序 %FLch3LSeg1 u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号u的经线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线
设一非线性系统如下所示: ()()() ()()111y k y k e u k k βαε--=-+-+ 0.75α= 0.35β=0.25γ=,()k ε是零均值，方差为0.01的噪声序列（均匀白噪声）。 （1）试设计一种激励信号能持续激励系统的各工作点（平衡点） （2）用适当的方法辨识出系统的等价模型（用另一组数据来检验模型的泛化性） 说明：下面讨论的都是离散系统，所以时间坐标均采用离散时间节点k 。 解： (1) 线性化处理寻找系统的合理输入信号 可以求得系统的平衡点为： ()0.75y k α== (1.1) 按题意要求最后系统必须收敛于平衡点附近，即： ()lim 0.75k y k →∞ = (1.2) 为了找出系统的合理输入信号，使得系统最终工作在平衡点附近，这里可以将系统线性化处理，将上述非线性系统进行泰勒展开得： ()()()()()()()2323111111112!3!! n n y k y k y k y k y k u k k n αβαβαβαβγε--+---++-=-+ 因为 ，后面()1y k -的高阶项都可以扔掉（只作为寻找输入信号使用）， 所以系统可以化为下式： 此时不妨设系统输出()y k 的最后的极限为A , 从式1.2得0.75A =。 那么应该满足 ()()lim lim 1k k y k y k A →∞ →∞ =-= (1.5) 从而有 ()()()11A u k k αβγε-=-+ (1.6) 同时为了抵消系统的部分噪声，这里采用MA TLAB 软件编程产生另一服从同一分布的均匀噪声()1k ε，将式1.6变形得： ()()()0111 1u k u k k αβ εγ γ --= - (1.7) 式1.7中()0u k 是一个最后收敛于系统平衡点0.75的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号()0u k ，同时考虑系统的平衡点，即设计为： ()/00.75k T u k e -=- (1.8) T 是一阶线性系统对应的时间常数， 反应到输入基本信号()0u k 上就是过零点作()0u k 对应()()()()11y k y k u k k αβγε--=-+01αβ<（1.3） (1.4)
最小二乘最优问题（转） 默认分类2009-05-21 14:56:33 阅读62 评论1 字号：大中小 1.约束线性最小二乘 有约束线性最小二乘的标准形式为 sub.to 其中：C、A、Aeq 为矩阵；d、b、beq、lb、ub、x 是向量。 在MA TLAB5.x 中，约束线性最小二乘用函数conls 求解。 函数lsqlin 格式x = lsqlin(C,d,A,b) %求在约束条件下，方程Cx = d 的最小二乘解x。 x = lsqlin(C,d,A,b,Aeq,beq) %Aeq、beq 满足等式约束，若没有不等式约束，则设A=[ ]，b=[ ]。 x = lsqlin(C,d,A,b,Aeq,beq,lb,ub) %lb、ub 满足，若没有等式约束，则Aeq=[ ]，beq=[ ]。 x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0) % x0 为初始解向量，若x 没有界，则lb=[ ]，ub=[ ]。 x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options) % options 为指定优化参 数 [x,resnorm] = lsqlin(...) % resnorm=norm(C*x-d)^2，即2-范数。 [x,resnorm,residual] = lsqlin(...) %residual=C*x-d，即残差。 [x,resnorm,residual,exitflag] = lsqlin(...) %exitflag 为终止迭代的条 件 [x,resnorm,residual,exitflag,output] = lsqlin(...) % output 表示输出
西北工业大学系统辩识大作业 题目：最小二乘法系统辨识
一、 问题重述： 用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数 离散化有 z^4 - 3.935 z^3 + 5.806 z^2 - 3.807 z + 0.9362 ---------------------------------------------- = z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187 噪声的成形滤波器 离散化有 4.004e-010 z^3 + 4.232e-009 z^2 + 4.066e-009 z + 3.551e-010 ----------------------------------------------------------------------------- = z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187 采样时间0.01s 要求：1.用Matlab 写出程序代码； 2.画出实际模型和辨识得到模型的误差曲线； 3.画出递推算法迭代时各辨识参数的变化曲线； 最小二乘法: 在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态 ,静态 , 线性 ,非线性系统。在使用最小二乘法进行参数估计时 ,为了实现实时控制 ,必须优化成参数递推算法 ,即最小二乘递推算法。这种辨识方法主要用于在线辨识。MATLAB 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。对 4324326.51411.5320120232320 Y s s s s G U s s s s ++++== ++++432 120120232320 E N W s s s s == ++++
最小二乘递推算法的MATLAB仿真 针对辨识模型，有z(k)-+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k)模型结构，对其进行最小二乘递推算法的MATLAB仿真，对比真值与估计值。更改a1、a2、b1、b2参数，观察结果。 仿真对象：z(k)-1.5*z(k-1)+0.7*z(k-2)=u(k-1)+0.5*u(k-2)+v(k) 程序如下： L=15; y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的初始值 for i=1:L; %移位循环 x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; %取出作为输出信号，即M序列 if y(i)>0.5,u(i)=-0.03; %输入信号 else u(i)=0.03; end y1=x1;y2=x2;y3=x3;y4=x4; end figure(1); stem(u),grid on z(2)=0;z(1)=0; for k=3:15; z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号 end c0=[0.001 0.001 0.001 0.001]'; %直接给出被识别参数的初始值 p0=10^6*eye(4,4); %直接给出初始状态P0 E=0.000000005; c=[c0,zeros(4,14)]; e=zeros(4,15); for k=3:15; %开始求k h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %开始求k的值 d1=z(k)-h1'*c0;c1=c0+k1*d1; e1=c1-c0; e2=e1./c0; %求参数的相对变化 e(:,k)=e2; c0=c1; c(:,k)=c1; p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出P(k)的值 p0=p1;
方法一、最小二乘一次性算法： 首先对最小二乘法的一次性辨识算法做简要介绍如下： 过程的黑箱模型如图所示： 其中u(k)和z(k)分别是过程的输入输出，)(1-z G 描述输入输出关系的模型，成为过程模型。 过程的输入输出关系可以描述成以下最小二乘格式： )()()(k n k h k z T +=θ （1） 其中z(k)为系统输出，θ是待辨识的参数，h(k)是观测数据向量，n(k) 是均值为0的随机噪声。 利用数据序列{z （k ）}和{h （k ）}极小化下列准则函数： ∑=-=L k T k h k z J 12])()([)(θθ （2） 使J 最小的θ的估计值^ θ，成为最小二乘估计值。 具体的对于时不变SISO 动态过程的数学模型为 )()()()()(11k n k u z B k z z A +=-- （3） 应该利用过程的输入、输出数据确定)(1-z A 和 )(1-Z B 的系数。 对于求解θ的估计值^θ，一般对模型的阶次 a n ， b n 已定，且b a n n >；其次将（3）模 型写成最小二乘格式 )()()(k n k h k z T +=θ （4） 式中 ?????=------=T n n T b a b a b b b a a a n k u k u n k z k z k h ],,,,,,,[)](,),1(),(,),1([)(2121 θ （5）
L k ,,2,1 = 因此结合式（4）（5）可以得到一个线性方程组 L L L n H Z +=θ （6） 其中 ???==T L T L L n n n n L z z z z )](),2(),1([)](),2(),1([ （7） 对此可以分析得出，L H 矩阵的行数为),max(b a n n L -，列数b a n n +。 在过程的输入为2n 阶次，噪声为方差为1，均值为0的随机序列，数据长度)(b a n n L +>的情况下，取加权矩阵L Λ为正定的单位矩阵I ，可以得出： L T L L T L z H H H 1^ )(-=θ （8） 其次，利用在Matlab 中编写M 文件，实现上述算法。 此次算法的实现，采用6阶M 序作为过程黑箱的输入；噪声采用方差为1，均值为0的随机数序列；黑箱模型假设为：y(k)-1.5y(k-1)+0.7y(k-2)=2u(k-1)+0.5u(k-2)，则系统输出为Z(k)-1.5Z(k-1)+0.7Z(k-2)=2U(k-1)+0.5U(k-2)+n （k ）；模型的阶次2,2==b a n n ；数据长度取L=200。 程序清单如下见附录：最小二乘一次性算法Matlab 程序 运行结果如下： 图1 最小二乘一次性算法参数真值与估计值 其中re 为真值，ans 为估计值^ θ 结果发现辨识出的参数与真值之间存在细微误差，这是由于系统噪声以及数据长度L 的限制引起的，最小二乘辨识法是一种无偏估计方法。 方法二、最小二乘递推算法： 最小二乘一次性算法计算量大，并且浪费存储空间，不利于在线应用，由此引出最小
R1=6378137; C=3e8; DeltaT=1e-3; SatellitePosition=[177461757273651;12127-9774210911;13324-181********;14000 -130********;19376-15756-73651;zeros(19,4);6400200-12000]; SatellitePosNew=ones(1,3); VisSatNum=0; %VisSatNum:可见的卫星个数; CalculateOK=1; for k=1:24 if(SatellitePosition(k,4)==1) VisSatNum=VisSatNum+1; SatellitePosNew=[SatellitePosNew;SatellitePosition(k,1:3)]; end end SatellitePosNew(1,:)=[];%把SatellitePosNew的第一行去掉 if(VisSatNum<4) CalculateOK=0; CalUserPosition=[000]; return end Prange=ones(1,VisSatNum); UserPos=SatellitePosition(25,1:3); for n=1:VisSatNum Prange(1,n)= sqrt((SatellitePosNew(n,:)-UserPos)*(SatellitePosNew(n,:)-UserPos)')+C*DeltaT; end CalculateRecord=[111]; XYZ0=[6400100-1200]; DeltaT0=0; Wxyz=SatellitePosNew; Error=1000; ComputeTime=0; while((Error>1)&(ComputeTime<1000)) ComputeTime=ComputeTime+1; R=ones(1,VisSatNum); for n=1:VisSatNum R(1,n)=sqrt((Wxyz(n,:)-XYZ0)*(Wxyz(n,:)-XYZ0)')+DeltaT*C;%DeltaT0???? end
系统辨识大作业 最小二乘法及其相关估值方法应用 学院：自动化学院 学号： 姓名：日期：
基于最小二乘法的多种系统辨识方法研究 一、实验原理 1.最小二乘法 在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。 设单输入-单输出线性定长系统的差分方程为 (5.1.1) 式中:为随机干扰；为理论上的输出值。只有通过观测才能得到，在观测过程中往往附加有随机干扰。的观测值可表示为 (5.1.2) 式中:为随机干扰。由式(5.1.2)得 (5.1.3) 将式(5.1.3)带入式(5.1.1)得 (5.1.4) 我们可能不知道的统计特性，在这种情况下，往往把看做均值为0的白噪声。 设 (5.1.5) 则式(5.1.4)可写成 (5.1.6) 在观测时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定不仅包含了的测量误差，而且包含了的测量误差和系统内部噪声。假定是不相关随机序列(实际上是相关随机序列)。 现分别测出个随机输入值，则可写成个方程，即 上述个方程可写成向量-矩阵形式 (5.1.7) 设 则式(5.1.7)可写为
(5.1.8) 式中：为维输出向量；为维噪声向量；为维参数向量；为测量矩阵。因此式(5.1.8)是一个含有个未知参数，由个方程组成的联立方程组。如果，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。如果，方程组正好与未知数数目相等，当噪声时，就能准确地解出 (5.1.9) 如果噪声，则 (5.1.10) 从上式可以看出噪声对参数估计是有影响的，为了尽量较小噪声对估值的影响。在给定输出向量和测量矩阵的条件下求系统参数的估值，这就是系统辨识问题。可用最小二乘法来求的估值，以下讨论最小二乘法估计。 2.最小二乘法估计算法 设表示的最优估值，表示的最优估值，则有 (5.1.11) 写出式(5.1.11)的某一行，则有 (5.1.12) 设表示与之差，即 - (5.1.13)式中 成为残差。把分别代入式(5.1.13)可得残差。设 则有 (5.1.14) 最小二乘估计要求残差的平方和为最小，即按照指数函数 (5.1.15) 为最小来确定估值。求对的偏导数并令其等于0可得 (5.1.16) (5.1.17)
M序列发生程序 function [out]=Mseries(order) %M-series function % by Chen Hao,2011/12/16 for i=1:order x(i)=1; end for n=1:252 y(1:2)=x(order-1:order); for i=order:-1:2 x(i)=x(i-1); end x(1)=xor(y(1),y(2)); out(n)=1-2*x(order); end stem(out); end 计算程序 function [out]=addfunction( A,B,order,n ) %addfunction % by Chen Hao,2011/12/9 order=n*order; for i=1:order out(i)=0; end for k=0:order-1 for i=0:order-1 out(k+1)=out(k+1)+A(i-k+order+1)*B(i+1); end end for k=1:order out(k)=(out(k)-out(order))/(order+1); end end 主程序 function main() %Main function % by Chen Hao,2011/12/9 out=Mseries(6); time=1:252; simin=[time' out']; assignin('base','simin',simin); sim('Mseries_distinguish'); z=simout.signals.values(2:253); out=[out out]; g=addfunction(out,z,63,4); g=g(1:63); plot(g); num=[1 2]; den=[3 1 1]; s=tf(num,den); t=0:62; h=impulse(s,t);
最小二乘法参数辨识 1 引言 系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。现代控制理论中的一个分支。通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类μ={M}（即给定一类已知结构的模型），一类输入信号u和等价准则J=L(y，yM)(一般情况下，J是误差函数，是过程输出y和模型输出yM的一个泛函)；然后选择使误差函数J达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。 2 系统辨识的目的 在提出和解决一个辨识问题时，明确最终使用模型的目的是至关重要的。它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。通过辨识建立数学模型通常有四个目的。 ①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。 ②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。用于系统分析的仿真模型要求能真实反映系统的特性。用于系统设计的仿真，则强调设计参数能正确地符合它本身的物理意义。 ③预测这是辨识的一个重要应用方面，其目的是用迄今为止系统的可测量的输入和输出去预测系统输出的未来的演变。例如最常见的气象预报，洪水预报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。预测模型辨识的等价准则主要是使预测误差平方和最小。只要预测误差小就是好的预测模型，对模型的结构及参数则很少再有其他要求。这时辨识的准则和模型应用的目的是一致的，因此可以得到较好的预测模型。 ④控制为了设计控制系统就需要知道描述系统动态特性的数学模型，建立这些模型的目的在于设计控制器。建立什么样的模型合适，取决于设计的方法和准备采用的控制策略。 3 系统辨识的方法 经典方法： 经典的系统辨识方法的发展已经比较成熟和完善，他包括阶跃响应法、脉冲
基于最小二乘法的系统参数辨识 吴令红，熊晓燕，张涛 太原理工大学机械电子研究所，太原 (030024) E-mail lhwu0818@https://www.doczj.com/doc/0310213856.html, 摘要：系统辨识是自动控制学科的一个重要分支，由于其特殊作用，已经广泛应用于各种领域，尤其是复杂系统或参数不容易确定的系统的建模。过去，系统辨识主要用于线性系统的建模，经过多年的研究，已经形成成熟的理论。但随着社会、科学的发展，非线性系统越来越受到人们的关注，其控制与模型之间的矛盾越来越明显，因而非线性系统的辨识问题也越来越受到重视，其辨识理论不断发展和完善本。文重点介绍了系统参数辨识中最小二乘法的基本原理，并通过悬臂梁模型的辨识实例，具体说明了基于最小二乘法参数辨识在Matlab 中的实现方法。结果表明基于最小二乘法具有算法简单、精度较高等优点。 关键词：系统辨识；参数辨识；滑动平均模型(ARX)；最小二乘法；Matlab 中图分类号：TH-9 1. 引言 所谓辨识就是通过测取研究对象在人为输入作用下的输出响应，或正常运行时的输入输出数据记录，加以必要的数据处理和数学计算，估计出对象的数学模型。这是因为对象的动态特性被认为必然表现在它的变化着的输入输出数据之中，辨识只不过是利用数学的方法从数据序列中提炼出对象的数学模型而已[1]。 最小二乘法是系统参数辨识中最基本最常用的方法。最小二乘法因其算法简单、理论成熟和通用性强而广泛应用于系统参数辨识中。本文基于悬臂梁的实测数据，介绍了最小二乘法的参数辨识在Matlab中的实现。 2. 系统辨识 一般而言，建立系统的数学模型有两种方法：激励分析法和系统辨识法。前者是按照系统所遵循的物化（或社会、经济等）规律分析推导出模型。后者则是从实际系统运行和实验数据处理获得模型。如图1所示，系统辨识就是从系统的输入输出数据测算系统数学模型的理论和方法。更进一步的定义是L.A.Zadeh曾经与1962年给出的，即“系统辨识是在输入和输出的基础上，从系统的一类系统范围内，确立一个与所实验系统等价的系统”。另外，系统辨识还应该具有3个基本要素，即模型类、数据和准则[5]。被辨识系统模型根据模型形式可分为参数模型和非参数模型两大类。所谓参数模型是指微分方程、差分方程、状态方程等形式的数学模型；而非参数模型是指频率响应、脉冲响应、传递函数等隐含参数的数学模型。在辨识工程中，模型的确定主要根据经验对实际对象的特性进行一定程度上的假设，如对象的模型是线性的还是非线性的、是参数模型还是非参数模型等。在模型确定之后，就可以根据对象的输入输出数据，按照一定的辨识算法确定模型的参数[4]。 y 图1 被研究的动态系统