系统辨识最小二乘Matlab仿真

系统辨识和最小二乘参数估计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]);。