观测信号(包括异常值)的状态估计方法(Matlab代码实现)

简介: 观测信号(包括异常值)的状态估计方法(Matlab代码实现)

💥1 概述

传统的系统状态估计方法只用到连续信号,而离散测量信号所包含的信息没有得到利用.提出一种基于混合信号(包括连续和离散)的系统状态估计方法,既利用了连续信号,也用到离散信号的信息。该方法将离散信号的变化视作系统的离散事件,提取其准确的信息并参与系统状态估计,构成具有混合系统特性的新型状态估计器。还讨论了该估计器的稳定性条件和设计方法。仿真实验证明这种所提出的状态估计方法可以有效地改善系统的状态估计性能。本文模拟 MCV 观察器,使用多个候选信号的中位数的状态估计方法来进行观测信号(包括异常值)。


📚2 运行结果

部分代码:

set(0,'DefaultAxesLinewidth',2,'DefaultLineLineWidth',2);
set(0,'defaultAxesFontSize',14);
set(0,'defaultAxesFontName','arial');
set(0,'defaultTextFontName','arial');
close all
clear
A = [0.7 0.5 -0.1;0 0.7 0.1;-0.3 0 0.9];
Bd = [0.1;0.1;0.2];
Bu = [-1.2;-0.8;1.4];
C = [1 2 -1;0 -5 -0.2];
D = [0.01;0.02];
E = [1 1 1];
A0 = A;
A1 = A^2;
A2 = A^3;
%A3 = A^4;
%A4 = A^5;
B0 = sqrt(3)*[zeros(3,1) Bd];
D0 = sqrt(3)*[D zeros(2,1)];
B1 = sqrt(4)*[zeros(3,1) Bd A*Bd];
D1 = sqrt(4)*[D zeros(2,1) zeros(2,1)];
B2 = sqrt(5)*[zeros(3,1) Bd A*Bd A^2*Bd];
D2 = sqrt(5)*[D zeros(2,1) zeros(2,1) zeros(2,1)];
gammaoptimal = 100000;
for i = 1:30
    for j = 1:30
        for k = 1:30
alpha0 = 0.03*i;
alpha1 = 0.03*j;
alpha2 = 0.03*k;
setlmis([])
[gamma,n,sgamma] = lmivar(1,[1 1]);
[P,n,sP] = lmivar(1,[3 1]);
[Y0,n,sY0] = lmivar(2,[3 2]);
[Y1,n,sY1] = lmivar(2,[3 2]);
[Y2,n,sY2] = lmivar(2,[3 2]);
S0 = newlmi;
lmiterm([-S0 1 1 P],(1-alpha0),1)
lmiterm([-S0 3 1 P],1,A0)
lmiterm([-S0 3 1 Y0],eye(3),C)
lmiterm([-S0 3 2 P],1,B0)
lmiterm([-S0 2 2 0],alpha0)
lmiterm([-S0 3 2 Y0],eye(3),D0)
lmiterm([-S0 3 3 P],1,1)
S1 = newlmi;
lmiterm([-S1 1 1 P],(1-alpha1),1)
lmiterm([-S1 3 1 P],1,A1)
lmiterm([-S1 3 1 Y1],1,C)
lmiterm([-S1 3 2 P],1,B1)
lmiterm([-S1 2 2 0],alpha1)
lmiterm([-S1 3 2 Y1],1,D1)
lmiterm([-S1 3 3 P],1,1)
S2 = newlmi;
lmiterm([-S2 1 1 P],(1-alpha2),1)
lmiterm([-S2 3 1 P],1,A2)
lmiterm([-S2 3 1 Y2],1,C)
lmiterm([-S2 3 2 P],1,B2)
lmiterm([-S2 2 2 0],alpha2)
lmiterm([-S2 3 2 Y2],1,D2)
lmiterm([-S2 3 3 P],1,1)
Sn = newlmi;
lmiterm([-Sn 1 1 P],1,1)
lmiterm([-Sn 2 1 0],E)
lmiterm([-Sn 2 2 gamma],1,1)
LMIs = getlmis;
c = [1 zeros(1,24)];
[copt,xopt] = mincx(LMIs,c);
if isempty(xopt)
else
Popt = dec2mat(LMIs,xopt,P);
Y0opt = dec2mat(LMIs,xopt,Y0);
Y1opt = dec2mat(LMIs,xopt,Y1);
Y2opt = dec2mat(LMIs,xopt,Y2);
gammaopt2 = dec2mat(LMIs,xopt,gamma);
L0p = -Popt^(-1)*Y0opt;
L1p = -Popt^(-1)*Y1opt;
L2p = -Popt^(-1)*Y2opt;
al0=abs(max(eig(A0-L0p*C)))
al1=abs(max(eig(A1-L1p*C)))
al2=abs(max(eig(A2-L2p*C)))
if (gammaopt2 < gammaoptimal)&&(alpha0<1-al0^2)&&(alpha1<1-al1^2)&&(alpha2<1-al2^2)
    gammaoptimal = gammaopt2;
    L0 = L0p;
    L1 = L1p;
    L2 = L2p;
    iopti = i;
    jopti = j;
    kopti = k;
end
end
dellmi(LMIs,S0);
dellmi(LMIs,S1);
dellmi(LMIs,S2);
dellmi(LMIs,Sn);
        end
    end
end
x(1,1) = 0;
x(2,1) = 0;
x(3,1) = 0;
xa(1,1) = 0;
xa(2,1) = 0;
xa(3,1) = 0;
x0(1,1) = 0;
x0(2,1) = 0;
x0(3,1) = 0;
x1(1,1) = 0;
x1(2,1) = 0;
x1(3,1) = 0;
x2(1,1) = 0;
x2(2,1) = 0;
x2(3,1) = 0;
xtotal(1,1) = 0;
xtotal(2,1) = 0;
xtotal(3,1) = 0;
Y(1,1) = x1(1,1);
Y(2,1) = x1(2,1);
Ts = 0.01;
t = 0:1:600;
a = 0:0.01:2*pi;
u = sin(a);
F = [];
e_proposed = 0;
for k = 1:600;
    x(:,k+1) = A * x(:,k) + Bu * u(:,k);
    y(:,k) = C * x(:,k);
    xa(:,k+1) = A * xa(:,k) + Bu * u(:,k) + Bd*(rand-0.5)*2;
    ya(:,k) = C * xa(:,k) + D*(rand-0.5)*2;
    if k < 4
    x0(:,k+1) = (A-L0*C)* xtotal(:,k) + Bu * u(:,k) + L0 * ya(:,k);
    x1(:,k+1) = (A-L0*C)* xtotal(:,k) + Bu * u(:,k) + L0 * ya(:,k);
    x2(:,k+1) = (A-L0*C)* xtotal(:,k) + Bu * u(:,k) + L0 * ya(:,k);
            F(1,1) = norm(x0(:,k+1));
            F(2,1) = norm(x1(:,k+1));
            F(3,1) = norm(x2(:,k+1));
            [I1,I2] = sort(F);
            if I2(2) == 1
                xtotal(:,k+1) = x0(:,k+1);
            elseif I2(2) == 2
                xtotal(:,k+1) = x1(:,k+1);
            elseif I2(2) == 3
                xtotal(:,k+1) = x2(:,k+1);
            end
    else
            if rem(k,400) == 0
                ya(1,k) = 100;
                figure(2)
                plot(k,100,'gs',...
    'MarkerSize',10,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor',[0.5,0.5,0.5]);
hold on;
            elseif rem(k,200) == 0
                ya(1,k) = -100;
                figure(2)
                plot(k,-100,'gs',...
    'MarkerSize',10,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor',[0.5,0.5,0.5]);
hold on;
            elseif rem(k,100) == 0
                ya(1,k) = 0;
                figure(2)
                plot(k,0,'gs',...
    'MarkerSize',10,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor',[0.5,0.5,0.5]);
hold on;
            end
            if rem(k,60) == 5
                ya(2,k) = 80;
                                figure(2)
                plot(k,80,'gs',...
    'MarkerSize',10,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor',[0.5,0.5,0.5]);
hold on;
            elseif rem(k,30) == 5
                ya(2,k) = -80;
                                figure(2)
                plot(k,-80,'gs',...
    'MarkerSize',10,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor',[0.5,0.5,0.5]);
hold on;
            end
            x0(:,k+1) = (A-L0*C)* xtotal(:,k) + Bu * u(:,k) + L0 * ya(:,k);
            x1(:,k+1) = (A^2-L1*C)*xtotal(:,k-1) + L1*ya(:,k-1) + A*Bu*u(:,k-1) + Bu*u(:,k);
            x2(:,k+1) = (A^3-L2*C) * xtotal(:,k-2) + L2*ya(:,k-2) + A^2*Bu*u(:,k-2) + A*Bu*u(:,k-1) + Bu*u(:,k);
            F(1,1) =  [1 1] * ya(:,k) - [1 1] * C * xtotal(:,k-1) + [1 1] * C * Bu * u(:,k-1);
            F(2,1) = [1 1] * ya(:,k-1) - [1 1] * C * xtotal(:,k-2) + [1 1] * C * Bu * u(:,k-2);
            F(3,1) = [1 1] * ya(:,k-2) - [1 1] * C * xtotal(:,k-3) + [1 1] * C * Bu * u(:,k-3);
            [I1,I2] = sort(F);
            if I2(2) == 1
                xtotal(:,k+1) = x0(:,k+1);
            elseif I2(2) == 2
                xtotal(:,k+1) = x1(:,k+1);
            elseif I2(2) == 3
                xtotal(:,k+1) = x2(:,k+1);
            end
  % xtotal(:,k+1) = (A-L0*C)* xtotal(:,k) + B * u(:,k) + L0 * ya(:,k);
    end
         e_proposed = e_proposed + norm(xtotal(:,k)-xa(:,k));
         F2(k) = I2(2);
end
figure(1),stairs(t,xa(1,:),'k','LineStyle','--','LineWidth',2),hold on,stairs(t,xa(2,:),'k','LineStyle','--','LineWidth',2),hold on,stairs(t,xa(3,:),'k','LineStyle','--','LineWidth',2)
figure(1),stairs(t,xtotal(1,:),'r'),hold on,stairs(t,xtotal(2,:),'r'),hold on,stairs(t,xtotal(3,:),'r'),axis([0 600 -30 30]);
%figure(2),stairs(t,x2(1,:),'b'), hold on
%figure(3),stairs(t,x3(1,:),'k'), hold on
%figure(4),stairs(t,xtotal(1,:),'k'), hold on
ylabel('x_p,x')
xlabel('k')
ya(1,601) = 0;
ya(2,601) = 0;
figure(2)
stairs(t,ya(1,:),'k')
hold on,axis([0 600 -105 105]);
stairs(t,ya(2,:),'k'),axis([0 600 -105 105]);
hold on
ylabel('y')
xlabel('k')
figure(3) %z,zQ1
stairs(t,E(1)*(xa(1,:)-xtotal(1,:))+E(2)*(xa(2,:)-xtotal(2,:))+E(3)*(xa(3,:)-xtotal(3,:)),'k','LineWidth',0.5);
axis([0 600 -2.5 2.5])
%set(gca,'fontsize',14);
xlabel('k')
ylabel('z_e')
figure(4) %z,zQ1
subplot(3,1,1)
stairs(t,(xa(1,:)-xtotal(1,:)),'k','LineWidth',0.5);
axis([0 600 -.5 .5])
ylabel('e_{1}')
xlabel('k')
%set(gca,'fontsize',14);
subplot(3,1,2)
stairs(t,xa(2,:)-xtotal(2,:),'k','LineWidth',0.5);
axis([0 600 -.5 .5])
%set(gca,'fontsize',14);
ylabel('e_{2}')
xlabel('k')
subplot(3,1,3)
stairs(t,xa(3,:)-xtotal(3,:),'k','LineWidth',0.5);
axis([0 600 -.5 .5])
%set(gca,'fontsize',14);
ylabel('e_{3}')
xlabel('k')
gammaopt = sqrt(gammaoptimal)
    iopti
    jopti
    kopti
    L0
    L1
    L2

🎉3 参考文献

部分理论来源于网络,如有侵权请联系删除。

[1]Hiroshi Okajima (2023) State estimation method using median of multiple candidates for observation signals including outliers [Source Code].


[2]莫以为,萧德云.基于混合信号的状态估计方法[J].控制理论与应用,2004(03):327-334.


🌈4 Matlab代码实现

相关文章
|
2月前
|
算法 数据安全/隐私保护 计算机视觉
基于二维CS-SCHT变换和LABS方法的水印嵌入和提取算法matlab仿真
该内容包括一个算法的运行展示和详细步骤,使用了MATLAB2022a。算法涉及水印嵌入和提取,利用LAB色彩空间可能用于隐藏水印。水印通过二维CS-SCHT变换、低频系数处理和特定解码策略来提取。代码段展示了水印置乱、图像处理(如噪声、旋转、剪切等攻击)以及水印的逆置乱和提取过程。最后,计算并保存了比特率,用于评估水印的稳健性。
|
1天前
|
算法 安全 数据库
基于结点电压法的配电网状态估计算法matlab仿真
**摘要** 该程序实现了基于结点电压法的配电网状态估计算法,旨在提升数据的准确性和可靠性。在MATLAB2022a中运行,显示了状态估计过程中的电压和相位估计值,以及误差随迭代变化的图表。算法通过迭代计算雅可比矩阵,结合基尔霍夫定律解决线性方程组,估算网络节点电压。状态估计过程中应用了高斯-牛顿或莱文贝格-马夸尔特法,处理量测数据并考虑约束条件,以提高估计精度。程序结果以图形形式展示电压幅值和角度估计的比较,以及估计误差的演变,体现了算法在处理配电网状态估计问题的有效性。
|
5天前
|
机器学习/深度学习 自然语言处理 算法
m基于深度学习的OFDM+QPSK链路信道估计和均衡算法误码率matlab仿真,对比LS,MMSE及LMMSE传统算法
**摘要:** 升级版MATLAB仿真对比了深度学习与LS、MMSE、LMMSE的OFDM信道估计算法,新增自动样本生成、复杂度分析及抗频偏性能评估。深度学习在无线通信中,尤其在OFDM的信道估计问题上展现潜力,解决了传统方法的局限。程序涉及信道估计器设计,深度学习模型通过学习导频信息估计信道响应,适应频域变化。核心代码展示了信号处理流程,包括编码、调制、信道模拟、降噪、信道估计和解调。
27 8
|
9天前
|
传感器 算法
ANC主动降噪理论及Matlab代码实现
ANC主动降噪理论及Matlab代码实现
|
2月前
|
算法 数据安全/隐私保护 C++
基于二维CS-SCHT变换和扩频方法的彩色图像水印嵌入和提取算法matlab仿真
该内容是关于一个图像水印算法的描述。在MATLAB2022a中运行,算法包括水印的嵌入和提取。首先,RGB图像转换为YUV格式,然后水印通过特定规则嵌入到Y分量中,并经过Arnold置乱增强安全性。水印提取时,经过逆过程恢复,使用了二维CS-SCHT变换和噪声对比度(NC)计算来评估水印的鲁棒性。代码中展示了从RGB到YUV的转换、水印嵌入、JPEG压缩攻击模拟以及水印提取的步骤。
|
1天前
|
算法
m基于GA遗传优化的高斯白噪声信道SNR估计算法matlab仿真
**MATLAB2022a模拟展示了遗传算法在AWGN信道中估计SNR的效能。该算法利用生物进化原理全局寻优,解决通信系统中复杂环境下的SNR估计问题。核心代码执行多代选择、重组和突变操作,逐步优化SNR估计。结果以图形形式对比了真实SNR与估计值,并显示了均方根误差(RMSE),体现了算法的准确性。**
8 0
基于高通滤波器的ECG信号滤波及心率统计matlab仿真
**摘要:** 使用MATLAB2022a,实施高通滤波对ECG信号预处理,消除基线漂移,随后分析心率。系统仿真展示效果,核心代码涉及IIR HPF设计,如二阶滤波器的差分方程。通过滤波后的信号,检测R波计算RR间期,从而得到心率。滤波与R波检测是心电生理研究的关键步骤,平衡滤波性能与计算资源是设计挑战。
|
10天前
|
机器学习/深度学习 算法 语音技术
基于语音信号MFCC特征提取和GRNN神经网络的人员身份检测算法matlab仿真
**语音识别算法概览** MATLAB2022a中实现,结合MFCC与GRNN技术进行说话人身份检测。MFCC利用人耳感知特性提取语音频谱特征,GRNN作为非线性映射工具,擅长序列学习,确保高效识别。预加重、分帧、加窗、FFT、滤波器组、IDCT构成MFCC步骤,GRNN以其快速学习与鲁棒性处理不稳定数据。适用于多种领域。
|
13天前
|
资源调度 SoC
基于UKF无迹卡尔曼滤波的电池Soc估计matlab仿真
**摘要:** 使用MATLAB2022a,基于UKF的电池SOC估计仿真比较真实值,展示非线性滤波在电动车电池管理中的效用。电池电气模型描述电压、电流与SoC的非线性关系,UKF利用无迹变换处理非线性,通过预测和更新步骤实时估计SoC,优化状态估计。尽管UKF有效,但依赖准确模型参数。
|
2月前
|
数据安全/隐私保护
地震波功率谱密度函数、功率谱密度曲线,反应谱转功率谱,matlab代码
地震波格式转换、时程转换、峰值调整、规范反应谱、计算反应谱、计算持时、生成人工波、时频域转换、数据滤波、基线校正、Arias截波、傅里叶变换、耐震时程曲线、脉冲波合成与提取、三联反应谱、地震动参数、延性反应谱、地震波缩尺、功率谱密度