05月
03

发表于 • 字数统计 3128 • 被 2,408 人看爆

伪随机序列辨识脉冲响应

在经典控制理论中，一个线性时不变系统的输入x(t)与输出y(t)我们可以用如下公式来表示。
image-1651542779909
其中，g(t)为该线性系统的脉冲响应函数。由于系统脉冲响应的拉氏变换即为系统传递函数。因此在理想的情况下，只需要向系统输入理想脉冲，并记录系统输出，即可得到系统传递函数，完成系统建模。如下狄拉克公式所示。
image-1651542845739
然而，在实际应用时，不存在理想的脉冲函数，直接用近似理想的脉冲函数进行激励得到的结果并不准确。
考虑到理想白噪声信号的自相关函数在时移t=0时为无穷，而在时移t≠0的时候为零。因此可以使用白噪声信号替代脉冲信号进行测试。若此时输入白噪声，则由Wiener-Hopf公式，有：
image-1651542871699
因为白噪声的自相关函数为 Rxx(t)=σ2δ(t)R_{xx} (t)=σ^2δ(t) ，因此，上式又可以表示为：
image-1651542908768
因此可以看出，我们可以通过白噪声的响应，求出系统的脉冲响应函数，从而进行系统建模。
不过，白噪声依然是理想信号，并不能用于实际实践中。
m序列与白噪声类似，也是一种随机信号。但是m序列可以由人工生成，适合在实际生产中应用。

MATLAB系统辨识

通过对系统输入M随机序列，得到响应曲线。这里随意选取一个系统举例。
通过自相关函数关系，可以将随机序列生成的输出，转为系统离散脉冲响应。然后可以按照以下流程进行系统辨识。

image-1652065393935

Rxy = zeros(idt_data_num,1);
for k = 1:idt_data_num
  t_temp = 0;
  for i = 1:idt_data_num
      k2 = k+i-1;
      k2 = k2-idt_data_num*(k2>idt_data_num);
      t_temp = t_temp+sys_in(i)*sys_out(k2);
  end
  Rxy(k) = t_temp;
end
Rxy = Rxy/idt_data_num;
markov = Rxy/sys_in(1)^2;
figure;
plot(markov);
title('系统脉冲响应曲线');

得到的结果如下图。
image-1651543308184

接着我们可以开始构造Hankel矩阵。离散脉冲响应序列可以与离散的Markov系数一一对应。因此可以取100个值构造50阶Hankel矩阵，并作平衡实现，进行降阶处理。

matrix_degree = 50;
HANKEL = hankel(markov(2:end));
H_nn = HANKEL(1:matrix_degree,1:matrix_degree);
H_nn1 = HANKEL(1:matrix_degree,2:matrix_degree+1);
[~,S,~]=svd(H_nn,0);
figure;
stem(1:matrix_degree,diag(S.^0.5));
title('Hankel特征值（SVD分解）');
grid on

上述得到的矩阵通过奇异值分解，可以得到奇异值由大到小排列如图。
image-1651543733079
显然可以看出，前3个奇异值显著高于其他奇异值。因此可以考虑设定系统阶次需要高于3阶，经过分别对3-5阶的系统进行辨识，最终选择设置为4阶系统。

由于奇异值分解后可以得到公式：
image-1651543972477
又同样有Hankel矩阵可以分解为能观性矩阵与能控性矩阵乘法，记作T_1。即：
image-1651543994051
将n阶hankel矩阵向后移1位记作T_2=H(n+1,n+1)，有：
image-1651544011998
由此，可以得到我们的能观性矩阵O为：
image-1651544034177
同样，能控性矩阵可以求得为：
image-1651544063100
求得上述系统通过3阶离散系统拟合后的结果，转为连续系统。并可以求得系统传递函数。

hankel_degree=4
[u,s,v]=svd(H_nn,0);
S=s.^0.5;
Obs=u*S
Ctl=S*v'
A=Obs\H_nn1/Ctl;
B=Ctl(:,1);
C=Obs(1,:);
A=A(1:hankel_degree,1:hankel_degree);
B=B(1:hankel_degree,1);
C=C(1,1:hankel_degree);
D=0;
sysd=ss(A,B,C,D,0.01);
sys=d2c(sysd);

绘制出伯德图，幅值和相位均能较好地拟合出实际系统。

Gs=tf(sys)
[mag,phase,tout]=bode(Gs);
figure(2)
% bode(Gs);
subplot(2,1,1)
semilogx(tout,20*log10(abs(mag(:))),'color',[1,0.1,0.1]);
ylabel('幅值','FontSize',11);
grid on
hold on
subplot(2,1,2)
phase1=ones(size(phase,3),1);
for i=1:size(phase,3)
    phase1(i,1)=phase(1,1,i)*phase1(i,1);
end
semilogx(tout,phase1-720,'color',[1,0.1,0.1]);
ylabel('相位','FontSize',11);
grid on
hold on;
[Txy,F]=tfestimate(sys_in,sys_out,[],[],[],100);
angley = 180*angle(Txy)/pi;
gb = size(angley,1);

for i=2:1:gb
    if abs(angley(i)-angley(i-1))>180
        angley(i)= angley(i)-360;
    end
    if abs(angley(i)-angley(i-1))>180
        angley(i)= angley(i)-360;
    end
%     if abs(angley(i)-angley(i-1))>180
%         angley(i)= angley(i)-360;
%     end
%     if abs(angley(i)-angley(i-1))>180
%         angley(i)= angley(i)-360;
%     end
end

figure(2);
subplot(2,1,1);
title('系统伯德图','FontSize',13 );
semilogx(2*pi*F,20*log10(abs(Txy)),'b');
legend('辨识系统','原始系统','FontSize',10)
subplot(2,1,2);
semilogx(2*pi*F,angley,'b');
legend('辨识系统','原始系统','FontSize',10)

image-1651544210362
当然，在频域上的系统拟合不好直接看出效果，因此我们可以在时域系统上进行对比。将一开始的M序列输入拟合的系统中，查看结果。

yy=lsim(c2d(sys,0.01),sys_in(1:300));
plot(yy,'r');
grid on
hold on
plot(sys_out(1:300),'b');
title('系统时域响应图','FontSize',12);
legend('辨识系统','原始系统','FontSize',10)

% xlim([0 182])
% ylim([-37.7 30.0])
figure;
ZPK = zpk(sys)
zplane(ZPK.z{:},ZPK.p{:})
grid on;
save('系统辨识.mat','sys')

margin(ZPK)

image-1651544374086

分享到：