系统建模与仿真

Joy6A

于 2023-03-03 21:01:25 发布

阅读量1.3k

点赞数

分类专栏：系统建模与仿真系统建模文章标签： matlab Powered by 金山文档

本文链接：https://blog.csdn.net/ln160864/article/details/128690763

版权

系统建模与仿真同时被 2 个专栏收录

1 篇文章 0 订阅

订阅专栏

系统建模

1 篇文章 0 订阅

订阅专栏

文章介绍了控制系统辨识的基本方法，包括经典辨识的阶跃响应法、脉冲函数法等和现代辨识的最小二乘法、极大似然法。在数据预处理阶段，重点讲述了零均值化和滤波处理。接着，讨论了模型阶次的确定，提出了行列式比法等方法，并展示了如何通过最小二乘法进行模型参数的递推估计。最后，对模型进行了校验，确保了辨识结果的准确性。

摘要由CSDN通过智能技术生成

第一章绪论

经典辨识（非参数模型辨识）：阶跃响应法、脉冲函数法、频率响应法、相关分析法等。

现代辨识方法（参数模型辨识）：最小二乘法、梯度校正法、极大似然法。

最小二乘法：最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

最小二乘格式可用如下方程表示：

$\text{[math]}$

其中 $\text{[math]}$ 是模型的输入变量， $\text{[math]}$ 是模型的输出变量， $\text{[math]}$ 是模型噪声。

第二章数据预处理

输入输出数据通常包含直流成分或低频成分，用任何辨识方法都不好消除对辨识精度的影响，此外，数据中的高频成分对辨识也不利。因此一般需要对输入输出数据进行零均值化和剔除高频成分。

1、零均值化

假设实际观测到的输入输出为 $\text{[math]}$ 和 $\text{[math]}$ ，则零均值化后的数据为：

$\text{[math]}$

其中u0和z0是直流分量，由于实际过程中这俩数是不知道的，因此没有实际意义，通常采用差分法来进行零均值化处理。

代码如下：

Y1=xlsread('F:\桌面\实验输入输出数据.xls',1,'A2:A1001');%读取输入数据
Y2=xlsread('F:\桌面\实验输入输出数据.xls',1,'B2:B1001');%读取输出数据
%%%%%%%%%%%%%%%差分化（零均值化）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y1chafen=diff(Y1);   %对输入输出数据计算向前差分，得到的差分量进行后续辨识
Y2chafen=diff(Y2);
figure(1)
subplot(2,1,1);
plot(Y1chafen);
title('输入数据差分化');
subplot(2,1,2);
plot(Y2chafen);
title('输出数据差分化')；

2、剔除高频

直接利用matlab自带的平滑滤波器进行滤波，去除输入输出数据的“毛刺”。

M=smooth(Y1chafen);           %滑动平均滤波
Z=smooth(Y2chafen);
figure(2)
subplot(2,1,1);
plot(M);
title('输入数据经滤波处理');
subplot(2,1,2);
plot(Z);
title('输出数据经滤波处理');

M序列和白噪声：

M序列是二进制伪随机码序列的一种形式，它的自相关函数接近脉冲函数，且具有输入净扰动小，幅值、周期、时钟节拍容易控制的优点，目前被普遍用于辨识输入信号，下图即为一种周期为15bit的M序列。

但是M序列中含有直流分量，会对辨识系统造成净扰动，而逆M序列可以克服这一缺点。

其与M序列的关系如下：

生产M序列和逆M序列的matlab代码如下：

clc;
clear all;
close all;
L=60;
x1=1;x2=1;x3=1;x4=0;
S=1;

for k=1:L
    M(k)=xor(x3,x4);
    IM=xor(M(k),S);
    if IM==0
        u(k)=-1;
    else
        u(k)=1;
    end

    S=not(S);
    
    x4=x3;x3=x2;x2=x1;x1=M(k);
end
subplot(2,1,1);
stairs(M);grid;
axis([0 L/2 -0.5 1.5]);xlabel('k');ylabel('M序列幅值');title('M序列');
subplot(2,1,2);
stairs(u);grid;
axis([0 L -1.5 1.5]);xlabel('k');ylabel('逆M序列幅值');title('逆M序列');

其输出如下所示：

白噪声（white noise)也称为伪随机序列，它是一种均值为零，谱密度为非零常数的平稳随机过程；

它的方差和平均功率是∞，而且它在任意两个瞬间的取值，都是互不相关的（无论这两个瞬间多近），如下图所示：

注：白噪声在工程上不易实现，因为工业设备不可能按照白噪声的变化规律动作。M序列具有近似白噪声的性质，可保证良好的辨识精度，且在工程上易于实现。

生成均值为0，方差为1 的白噪声可直接使用randn函数，例如：xi=randn(500,1)；即生成一个长度为500的白噪声。如下图所示：

异常值剔除的方法有：肖维勒准则、迪克逊准则、格拉布斯准则、3sigama准则。

本次实验采用肖维勒方法来进行异常值的排查与剔除，肖维勒方法的原理为：在n次测量结果中，如果某误差可能出现的次数小于半次时，就予以剔除。这实质是规定了置信概率为1/2n，根据这一置信概率，可以计算出肖维勒系数，该系数可按下列公式近似计算：

$\text{[math]}$

如果某测量值与平均值之差的绝对值大于标准偏差与肖维勒系数之积，则该测量值被剔除。

$\text{[math]}$

Matlab代码如下：

clc;clear all;close all;
Data=xlsread('F:\桌面\实验输入输出数据.xls',1,'A2:A1001');
Data1=xlsread('F:\桌面\实验输入输出数据.xls',1,'B2:B1001');
Y=[];
w=1+0.4*log(1000);
  for i = 1:1000
   x = Data(i);    
   YiChang = abs(x-mean(x)) > w*std(x);
   Y(i) = YiChang;
  end
u= find(Y() == 1);   % 找出异常值所在的位置
uu = unique(u);    % 剔除重复的行数
Data(uu) = [ ];   %令异常值所在位置为空,即剔除异常值
figure(1)
plot(Data);title('剔除异常值后的输入数据')
figure(2)
plot(Y);title('输入异常值');

第三章模型阶次的确定

常用的方法有：

Hankel矩阵判秩法、行列式比法、F检验法、AIC法和最终预报误差准则法等。

其中Hankel矩阵判秩法需要利用模型脉冲响应组成的Hankel矩阵的秩来判定模型阶次，在

只知道系统的输入输出数据的时候不适用。

而行列式比法是利用输入输出数据所组成的矩阵秩之比来判定模型阶次。在有色噪声的情况下，即模型如下时：

$\text{[math]}$ 比是白噪声时稍微复杂些，具体代码如下（copy的）：

clear all
% num返回的是excel中的数据，txt输出的是文本内容，raw输出的是未处理数据
[num,txt,raw]=xlsread('F:\桌面\实验输入输出数据.xls') ;
u=num(:,1);
z=num(:,2);
N=4;
L=500;
U=u(401:1000);
Z=z(401:1000);
for n=1:N
    for i=1:L        %L行
        for j=1:n    %n列
            Zk(i,j)=Z(n+i-j);
            Uk(i,j)=U(n+i-j);
        end
    end
    Hn{n}=[Zk,Uk];
    H{n}=(    (  (Hn{n})'  )  *     (Hn{n})     ) / L;
end
for n=1:N-1        %求行列式比值
    DR(n)=det(H{n})/det(H{n+1});
end
for i=1:N-1    %比较差值
    if i==1
        E(i)=DR(i);
    else
        E(i)=DR(i)-DR(i-1);
    end
end
[i,n]=max(E);  %确定出行列式比显著增加时的阶次
n
%=======m=1==========
P=500;
a1=[0.01;0.01];b1=[0.01;0.01];c1=0.01;
for k=1:10
    for i=3:500
        estimate_e1(i)=z(i)'+[z(i-1)',z(i-2)']*a1-[u(i-1)',u(i-2)']*b1;
        estimate_v1(i)=estimate_e1(i)+estimate_e1(i-1)*c1;
        Em1(i,1)=estimate_e1(i);
        Vm1(i,1)=estimate_v1(i);   
    end
end
H22_1=[Em1'*Em1,Em1'*Vm1;Vm1'*Em1,Vm1'*Vm1];
%=======m=2==========
a2=[0.01;0.01];b2=[0.01;0.01];c2=[0.01;0.01];
for k=1:10
    for i=3:500
     estimate_e2(i)=z(i)'+[z(i-1)',z(i-2)']*a2-[u(i-1)',u(i-2)']*b2;
     estimate_v2(i)=estimate_e2(i)+[estimate_e2(i-1),estimate_e2(i-2)]*c2;
     Em2(i,1:2)=[estimate_e2(i),estimate_e2(i-1)];
     Vm2(i,1:2)=[estimate_v2(i),estimate_v2(i-1)];
    end
end
H22_2=[Em2'*Em2,Em2'*Vm2;Vm2'*Em2,Vm2'*Vm2];
%=======m=3==========
a3=[0.01;0.01];b3=[0.01;0.01];c3=[0.01;0.01;0.01];
for k=1:10
    for i=4:500
    estimate_e3(i)=z(i)'+[z(i-1)',z(i-2)']*a3-[u(i-1)',u(i-2)']*b3;
    estimate_v3(i)=estimate_e2(i)+[estimate_e2(i-1),estimate_e2(i-2),estimate_e2(i-3)]*c3;
    Em3(i,1:3)=[estimate_e3(i),estimate_e3(i-1),estimate_e3(i-2)];
    Vm3(i,1:3)=[estimate_v3(i),estimate_v3(i-1),estimate_v3(i-2)];
    end
end
H22_3=[Em3'*Em3,Em3'*Vm3;Vm3'*Em3,Vm3'*Vm3];
DR1(1)=det(H22_1)/det(H22_2);
DR1(2)=det(H22_2)/det(H22_3);
for i=1:2    %比较差值
    if i==1
        E1(i)=DR1(i);
    else
        E1(i)=DR1(i)-DR1(i-1);
    end
end
[i,m]=max(E1);  %确定出行列式比显著增加时的阶次
m

注：有色噪声模型的行列式比法在辨别阶次n时和无噪声是一样的。如下图所示：

由于H11(n）只与输入输出数据有关，故辨别n时与无噪声时代码一样。而辨别m时要稍微复杂些，因为H22(m)包含Em和Vm，其中元素是不可测的，需要用估计值代替，Em和Vm元素的估计值可按下式计算：

在进行模型阶次m辨识的时候，需要模型的参数已知，而实际情况不可知时，可对参数进行估计，在知道阶次n的时候，便可确定参数a的个数和b的个数。如若n=2,则模型存在a0,a1,b0,b1四个参数。先对这四个参数进行估计，如我取的是都等于0.01，由于n已经计算出来，故固定n，使m逐次增加，计算DR(m），当DR（m)较DR（m-1）有显著增加时，认为噪声模型的阶次m0=m。

注：目前存在的问题是，在对m进行辨识时，需要模型的其他参数已知（或进行估计），怎么正确的对它进行估计？虽然我尝试过把估计值改成1、100，对最后的结果并没有影响。

看到递推最小二乘法进行参数估计时，需要确定估计的初始参数值 $\text{[math]}$ ，书上写的是 $\text{[math]}$ 为零向量或充分小的正的实向量，所以上面假设a0~b1=0.01是合理的。

模型阶次辨识和参数辨识没有谁先谁后的问题，如果先进行模型参数辨识，在阶次未知时，可以对阶次进行假设，且一般假设 $\text{[math]}$ ，且一般从小往大取，例如可取 $\text{[math]}$ ，若这样辨识得到的参数不满意，可再假设 $\text{[math]}$ 。以此类推。

第四章模型参数辨识

利用最小二乘原理，通过极小化广义误差的平方和函数来确定模型的参数，是最基本、最典型的一种方法。

最小二乘法所用的模型是：

$\text{[math]}$

增广最小二乘法所用的模型为：

$\text{[math]}$

广义最小二乘法所用的模型是：

$\text{[math]}$

由于我大作业的模型是第二种，所以采用（递推）增广最小二乘法来进行辨识。注：利用最小二乘法辨识模型参数的目标是得到无偏一致估计，因此需要结合模型的类型来选择最小二乘法的种类，不然得到的结果就不是无偏一致的，即效果不好。

确定阶次之后的模型： $\text{[math]}$

共有5个参数，a1,a2,b1,b2,d1。其具体递推公式如下：

要注意的是，白噪声是不可观测的，因此也需要估计，在 $\text{[math]}$ 时，v(k)=0;在k>0时，则需要进行估计，公式如下： $\text{[math]}$

具体代码如下：

%递推增广最小二乘参数估计（RELS）
clear all; close all;
load("x1.mat");
load("y1.mat");
na=2; nb=1; nc=1; d=1;%na、nb、nc为A、B、C阶次
L=1000; %仿真长度
% uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)
% yk=zeros(na,1); %输出初值
xik=zeros(nc,1); %噪声初值
xiek=zeros(nc,1); %噪声估计初值
u=randn(L,1); %输入采用白噪声序列
xi=randn(L,1); %白噪声序列
thetae_1=zeros(na+nb+nc+1,1); %na+nb+1+nc为辨识参数个数
P=10^6*eye(na+nb+nc+1);
for k=3:L
    
    phie=[-y1(1:2);x1(d:d+nb);xik]; %组建phie
    
    %递推增广最小二乘法
    K=P*phie/(1+phie'*P*phie);
    thetae(:,k)=thetae_1+K*(y1(k)-phie'*thetae_1);
    P=(eye(na+nb+nc+1)-K*phie')*P;
    xie=y1(k)-phie'*thetae(:,k); %白噪声的估计值
    
    %更新数据
    thetae_1=thetae(:,k);
    
    for i=d+nb:-1:2
        x1(i)=x1(i-1);
    end
    x1(1)=x1(k);
    
    for i=na:-1:2
        y1(i)=y1(i-1);
    end
    y1(1)=y1(k);
    
    for i=nc:-1:2
        xik(i)=xik(i-1);
        xiek(i)=xiek(i-1);
    end
    xik(1)=xi(k);
    xiek(1)=xie;
end
figure(4)
plot([1:L],thetae(1:na,:));
xlabel('k'); ylabel('参数估计a');
legend('a_1','a_2'); axis([0 L -2.5 1.5]);
figure(5)
plot([1:L],thetae(na:na+nb+1,:));
xlabel('k'); ylabel('参数估计b');
legend('b_1','b_2'); axis([0 L -1 1]);
figure(6)
plot([1:L],thetae(na+nb+1:na+nb+nc+1,:));
xlabel('k'); ylabel('参数估计d');
legend('d_1'); axis([0 L -0.2 0.1]);

最后得到的参数为：a1=-1.73 a2=0.7456 b1=0.0129 b2=0.0026 d1=2.8744e-5

第五章模型校验

模型的阶次和参数已经辨识出来，实验只剩最后一步：校验

代码如下：

%%%%%%%%%模型校验%%%%%%%%%%%%%%

L=600;
U=x1(401:1000);
Z=y1(401:1000);
U=U';
Z=Z';
Z1(1)=0;
Z1(2)=0.0026*U(1)+xie(2)-0.000004174*xie(k-1);
for k=3:L
    Z1(k)=0.0129*U(k-1)+0.0026*U(k-2)+1.73*Z1(k-1)-0.7456*Z1(k-2)+xie(k)-0.000004174*xie(k-1);
end
error=Z-Z1;
figure(3) 
hold on
grid on
plot(Z1,'r')
plot(Z,'g')
plot(error,'y')
legend('仿真值','实际值','误差');

输出结果如下：