动态系统数据驱动验证案例研究（Matlab代码实现）

最新推荐文章于 2024-09-12 15:02:06 发布

荔枝科研社

最新推荐文章于 2024-09-12 15:02:06 发布

阅读量193

点赞数

分类专栏：物理应用/数据驱动文章标签：数学建模

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

版权

物理应用/数据驱动专栏收录该内容

116 篇文章 7 订阅

订阅专栏

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

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

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

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

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

本代码为动态系统数据驱动验证案例研究，共四个案例。

第一个案例是从以前的非线性分析和验证工作中获取的范德波尔振荡器，并为比较提供了基线。该问题的目标是估计系统能够成功返回原点（0，0）的一组初始条件。

第二个案例是并发学习模型参考自适应控制器（CL-MRAC）实现。在这个示例问题中，我们希望估计完整的参数集（W1，W2），自适应系统仍然能够从中跟踪±1.0范围内的参考轨迹。

第三个案例是飞机回归安全控制器，该控制器旨在从单独的控制器接管并将系统返回到标称状态。在这个问题中，我们感兴趣的是估计回归控制器可以成功接管的扰动状态集。

第四个案例是基于强化学习的双连杆机器人控制器，即机器人问题。控制器使用Sarsa的先前强化学习步骤进行训练，旨在将两个链接相互重叠，类似于体操运动员在高杆上挥杆。这里感兴趣的问题是验证训练的策略是否可以在 1000 个模拟步骤（\Delta t = 0.5）内从各种初始条件开始稳定 acrobot。

📚2 运行结果

主函数部分代码：

clear all;
close all;
clc;


%% Load training and  testing data
load('vdp_data.mat'); %both training and testing data


%% Train the SVM Classifier
tic

% Setup training data
data1 = dataX( dataY == -1 , : );
data2 = dataX( dataY == 1, : );
data3 = [data1;data2];
theclass = ones( length(data3) ,1);
theclass(1:length(data1)) = -1;

dpoly = 4; % polynomial kernel order

% Train the SVM Classifier
cl = fitcsvm(data3,theclass,'KernelFunction','polynomial',...
    'PolynomialOrder',dpoly,'BoxConstraint',10,'ClassNames',[-1,1],...
    'Solver','SMO','Standardize',true,'ScoreTransform','sign',...
    'Cost',[0 1; 1 0]);

disp('------------------------------------------------------------------');

%Test on the training data
disp('Testing on training data:');
[labels1,scores] = predict(cl,data3);
errfrac = sum(labels1 ~= theclass) / size(theclass,1);
disp(['error = ' num2str(errfrac) ' %']);
istabInc = find(theclass == 1 & theclass ~= labels1); % false positive
disp(['FP = ' num2str(length(istabInc)) ' of ' num2str(length(find(theclass == 1)))]);
iunstInc = find(theclass == -1 & theclass ~= labels1); % false negative
disp(['FN = ' num2str(length(iunstInc)) ' of ' num2str(length(find(theclass == -1)))]);
toc

% Generalization error
tic
for ii = 1:100
    M = crossval(cl, 'Kfold', 10);
    L(ii) = kfoldLoss(M);
end
disp('------------------------------------------------------------------');
disp(['Generalization error = ' num2str(mean(L))]);
toc


%% Predict scores over the testing grid
tic
disp('------------------------------------------------------------------');
disp('Empirical testing:');
[labels1,scores] = predict(cl,xGrid);
errfrac = sum(labels1 ~= xGridLabels) / size(xGridLabels,1);
disp(['error = ' num2str(errfrac) ' %']);
istabInc = find(xGridLabels == 1 & xGridLabels ~= labels1); % false positive
disp(['FP = ' num2str(length(istabInc)) ' of ' num2str(length(find(xGridLabels == 1)))]);
iunstInc = find(xGridLabels == -1 & xGridLabels ~= labels1); % false negative
disp(['FN = ' num2str(length(iunstInc)) ' of ' num2str(length(find(xGridLabels == -1)))]);
toc


%% Plot test results:
figure; hold on;
plot(xGrid(xGridLabels==1 & labels1~=xGridLabels,1),xGrid(xGridLabels==1 & labels1~=xGridLabels,2),'r.')
plot(xGrid(xGridLabels==-1 & labels1~=xGridLabels,1),xGrid(xGridLabels==-1 & labels1~=xGridLabels,2),'b.')
legend('False Positive', 'False Negative');
hold off;
figure; hold on;
plot(xGrid(xGridLabels==1 & labels1==xGridLabels,1),xGrid(xGridLabels==1 & labels1==xGridLabels,2),'b.')
plot(xGrid(xGridLabels==-1 & labels1==xGridLabels,1),xGrid(xGridLabels==-1 & labels1==xGridLabels,2),'r.')
legend('Correct Unstable', 'Correct Stable');
hold off;

% Plot the data and the decision boundary:
figure; hold on;
h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.',16);
h(3) = plot(data3(cl.IsSupportVector,1),data3(cl.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'Stab','Unstab','Support Vectors'});
xlim([-3 3]);
ylim([-3 3]);
xlabel('x_1');
ylabel('x_2');
hold off;

% Plot the decision boundary along with the true limit cycle:
figure; hold on;
set(gcf,'renderer','zbuffer');
[CC HH] = contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'--k');
set(HH,'LineWidth',3)
vdp1 = @(t,z) [z(2), (1 - z(1)^2)*z(2) - z(1)]';
[T,Z] = ode45(vdp1,[0 20],[2 0]);
plot(Z(:,1),Z(:,2),'
r.
')
title(['
SVM Classification: Polynomial kernel degree 
' num2str(dpoly)]);
xlabel('
x_1
');
ylabel('
x_2
');
hold off;