移动地平线估算器（MHE）和 cEKF 应用于 4 个圆柱形坦克（Matlab代码实现）

最新推荐文章于 2023-12-31 10:28:26 发布

荔枝科研社

最新推荐文章于 2023-12-31 10:28:26 发布

阅读量396

点赞数

文章标签：数学建模

本文链接： https://blog.csdn.net/weixin_46039719/article/details/131197596

版权

👨‍🎓个人主页：研学社的博客

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

移动地平线估算器（MHE）和 cEKF 应用于 4 个圆柱形坦克（Matlab代码实现）

disp(' ')
disp('***************************************************************************')
disp('Plot comparison between "MHE" and "cEKF" estimators!!!')
disp('***************************************************************************')
disp(' ')

if exist('MHE_plot', 'var') & exist('MHE_plot', 'var')
disp('... saving data to ''MHE_cEKF_data.mat''!!!')
save('MHE_cEKF_data', 'MHE_plot', 'cEKF_plot', 'ymea_plot')
disp(' ')
disp('... ploting the just simulated data!!!')

else

disp('... ploting the previous saved data!!!')
disp(' ')
end

clear all

load('MHE_cEKF_data')
load('prbs_signal')

t = MHE_plot.time;

ym1 = ymea_plot.signals.values(:,1);
ym2 = ymea_plot.signals.values(:,2);

h1p               =   MHE_plot.signals.values(:,5);
h2p               =   MHE_plot.signals.values(:,6);
h3p               =   MHE_plot.signals.values(:,7);
h4p               =   MHE_plot.signals.values(:,8);

h1m               =   MHE_plot.signals.values(:,1);
h2m               =   MHE_plot.signals.values(:,2);
h3m               =   MHE_plot.signals.values(:,3);
h4m               =   MHE_plot.signals.values(:,4);

h1c               =   cEKF_plot.signals.values(:,1);
h2c               =   cEKF_plot.signals.values(:,2);
h3c               =   cEKF_plot.signals.values(:,3);
h4c               =   cEKF_plot.signals.values(:,4);

f1 = F1_prbs ;
f2 = F2_prbs ;

op1 = menu('Plot Mode:', ...
'Plot: N=2, All',...
'Plot: N=2, h1, h2',...
'Subplot.english: MHE, cEKF, PRBS',...
'Sub.portuguese: MHE, cEKF, PRBS',...
'Subplot: N=2',...
'Subplot: F1, F2');

switch op1
case 1
caso   =   1;
case 2
caso   =   2;
case 3
caso   =   3;
case 4
caso   =   4;
case 5
caso   =   5;
case 6
caso   =   6;
end

figure(1)

if caso== 4,
% Portuguese

subplot(3,1,1)
plot(t, ym1, t, ym2, t, h1m,'.', t, h2m,'.');
axis([0 100 6 15])
ylabel('N韛el [cm]')
title('Compara玢o de Estimadores (MHE x cEKF): modelo 4 tanques cil韓dricos')
legend('ym1','ym2','h1_M_H_E','h2_M_H_E' )
text(40,12, 'P_0 = diag[1 1 1 1]*10^3 cm^2')
text(40,11, 'Q = diag[1 1 1 1] cm^2')
text(40,10, 'R = diag[1 1]*10 cm^2')
text(40, 9, 'NA = [1.0 1.0] cm')

text(20, 9,'Janela MHE (N=2)')

subplot(3,1,2)
plot(t, ym1, t, ym2, t, h1c,'.', t, h2c,'.');
axis([0 100 6 15])
ylabel('N韛el [cm]')
legend('h1','h2','h1_c_E_K_F','h2_c_E_K_F' )

subplot(3,1,3)
plot( t, f1, t, f2);
axis([0 100 3 11])
xlabel('Tempo [min]')
ylabel('Vaz鉶 [cm^3/min]')
legend('F1_i_n_p_u_t','F2_i_n_p_u_t' )
text(10,10,'Sinal PRBS')
end

figure(2)

switch caso
case 1
plot(t, h1p, t, h2p, t, h3p, t, h4p,...
t, h1m, t, h2m, t, h3m, t, h4m, t, h1c, t, h2c, t, h3c, t, h4c);
axis([0 100 0 13])
xlabel('Time [min]')
ylabel('Tank level [cm]')
title('Estimators comparison (MHE x cEKF) applied to 4 cylindrical tanks')
legend('h1','h2','h3','h4','h1_M_H_E','h2_M_H_E','h3_M_H_E',...
'h4_M_H_E', 'h1_c_E_K_F','h2_c_E_K_F','h3_c_E_K_F','h4_c_E_K_F' )

case 2
plot(t, h1p, t, h2p, t, h1m,'.', t, h2m,'.', t, h1c,':', t, h2c,':');
axis([0 100 6 13])
xlabel('Time [min]')
ylabel('Level [cm]')
title('Estimators comparison (MHE x cEKF) applied to 4 cylindrical tanks')
legend('h1','h2','h1_M_H_E','h2_M_H_E','h1_c_E_K_F','h2_c_E_K_F' )
text(60,12,'MHE Horizon (N=2)')

case 3
% English
subplot(3,1,1)
plot(t, h1p, t, h2p, t, h1m,'.', t, h2m,'.');
axis([0 100 6 15])
ylabel('N韛el [cm]')
title('Estimators comparison (MHE x cEKF) applied to 4 cylindrical tanks')
legend('h1','h2','h1_M_H_E','h2_M_H_E' )
text(40,12, 'P_0 = diag[1 1 1 1]*10^3 cm^2')
text(40,11, 'Q = diag[1 1 1 1] cm^2')
text(40,10, 'R = diag[1 1]*10 cm^2')
text(40, 9, 'NA = [1.0 1.0] cm')

text(20, 9,'MHE Horizon (N=2)')

subplot(3,1,2)
plot(t, h1p, t, h2p, t, h1c,'.', t, h2c,'.');
axis([0 100 6 15])
ylabel('Level [cm]')
legend('h1','h2','h1_c_E_K_F','h2_c_E_K_F' )

subplot(3,1,3)
plot( t, f1, t, f2);
axis([0 100 3 11])
xlabel('Time [min]')
ylabel('Flow [cm^3/min]')
legend('F1_i_n_p_u_t','F2_i_n_p_u_t' )
text(10,10,'PRBS signal')

case 4
% Portuguese

subplot(3,1,1)
plot(t, h1p, t, h2p, t, h1m,'.', t, h2m,'.');
axis([0 100 6 15])
ylabel('N韛el [cm]')
title('Compara玢o de Estimadores (MHE x cEKF): modelo 4 tanques cil韓dricos')
legend('h1','h2','h1_M_H_E','h2_M_H_E' )
text(40,12, 'P_0 = diag[1 1 1 1]*10^3 cm^2')
text(40,11, 'Q = diag[1 1 1 1] cm^2')
text(40,10, 'R = diag[1 1]*10 cm^2')
text(40, 9, 'NA = [1.0 1.0] cm')

text(20, 9,'Janela MHE (N=2)')

subplot(3,1,2)
plot(t, h1p, t, h2p, t, h1c,'.', t, h2c,'.');
axis([0 100 6 15])
ylabel('N韛el [cm]')
legend('h1','h2','h1_c_E_K_F','h2_c_E_K_F' )

subplot(3,1,3)
plot( t, f1, t, f2);
axis([0 100 3 11])
xlabel('Tempo [min]')
ylabel('Vaz鉶 [cm^3/min]')
legend('F1_i_n_p_u_t','F2_i_n_p_u_t' )
text(10,10,'Sinal PRBS')

case 5
subplot(2,1,1)
plot(t, h1p, t, h2p, t, h1m,'.', t, h2m,'.');
axis([0 100 6 13])
ylabel('Level [cm]')
title('Estimators comparison (MHE x cEKF) applied to 4 cylindrical tanks')
legend('h1','h2','h1_M_H_E','h2_M_H_E' )
text(60,12,'MHE Horizon (N=2)')

subplot(2,1,2)
plot(t, h1p, t, h2p, t, h1c,'.', t, h2c,'.');
axis([0 100 6 13])
xlabel('Time [min]')
ylabel('Level [cm]')
legend('h1','h2','h1_c_E_K_F','h2_c_E_K_F' )

case 6
subplot(3,1,1)
plot(t, h1p, t, h2p, t, h1m,'.', t, h2m,'.');
axis([0 100 6 13])
ylabel('Level [cm]')
title('Estimators comparison (MHE x cEKF) applied to 4 cylindrical tanks')
legend('h1','h2','h1_M_H_E','h2_M_H_E' )
text(60,12,'MHE Horizon (N=2)')

subplot(3,1,2)
plot(t, h1p, t, h2p, t, h1c,'.', t, h2c,'.');
axis([0 100 6 13])
ylabel('Level [cm]')
legend('h1','h2','h1_c_E_K_F','h2_c_E_K_F' )

subplot(3,1,3)
plot( t, f1, t, f2);
axis([0 100 6 8.5])
xlabel('Time [min]')
ylabel('Flow [cm^3/min]')
legend('F1_i_n_p_u_t','F2_i_n_p_u_t' )
text(40,7.5,'PRBS signal')

end