作者dennisshyu12 (路人)
看板MATLAB
标题[问题] 照着论文做状态分析但就是无法得出一样
时间Mon Feb 22 00:27:52 2021
想请教版上的大神，是否是小弟程式有误导至结果不同，或者是公式的部分打错了...
请大神帮帮忙
[状况说明]
四种状态数据分数阶导数後，做可拓的物质元素分析 ，皆照着步骤与公式做，
但结果就是与论文不同，达不到百分之百，
论文名:Fractional-Order Chaotic Self-Synchronization Based Tracking Faults Diagnosis
of Ball Bearing System
以下是我的程式，附上步骤与说明
%-------------------------------------------------------------------------------
% % 步骤
% training
% -------------------------------
% 1.从测试集中，随机取一段距离(一千笔数据)，将4种状态的数据代到分数阶混沌方程
，并把结果作欧式(norm)存到一个矩阵，重复作直到回圈结束(测试也一样，不过只要代
入一种数据即可)。
%
% 2.寻找各种状态得最大最小值，储存为经典域，并各往外扩长10%
% testing
% ----------------------------
% 3.将每笔测试数据，计算出与可拓域与经典域之距离
%
% 4.计算每笔测试数据与4种状态得关联度，并进行正规化
%
% 5.比较结果，最大值即判断为该种状态。
%4种DATA合并
S(:,1) = X097_DE_time(1:240000); %normal
S(:,2) = X122_DE_time(1:240000); %ball
S(:,3) = X109_DE_time(1:240000); %inner
S(:,4) = X135_DE_time(1:240000); %outer
%参数固定
a = 5;
b = -10;
c = -3.8;
alpha = 0.3;
% % trainning
%
-----------------------------------------------------------------------------------------------------------------
R=randi([48000,134000],1,50); % 随机几段 %几次
for s=1:4
if s==1 % s 是状态
for i = 1:50
for T = R(i) : (R(i)+1000)
%以下是 分数阶混沌方程 公式
E(1,s) = S(T,s);
E(2,s) = S(T+1,s);
E(3,s) = S(T+2,s);
g1= [ a , -E(3,s) , 0 ; E(3,s) , b , 0 ; 1/3*E(2,s) , 0 , c ];
h1= ( gamma(2) / gamma( 2 - alpha )) ;
h2= [ E(1,s)^(1- alpha ) ; E(2,s)^(1-alpha) ;
E(3,s)^(1-alpha) ];
e1(:,i)=g1*h2.*h1;
end
O1(1,i) = norm(e1(1,:));
O1(2,i) = norm(e1(2,:));
end
end
if s==2
for i = 1:50
for T = R(i) : (R(i)+1000)
E(1,s) = S(T,s);
E(2,s) = S(T+1,s);
E(3,s) = S(T+2,s);
g1= [ a , -E(3,s) , 0 ; E(3,s) , b , 0 ; 1/3*E(2,s) , 0 , c ];
h1= ( gamma(2) / gamma( 2 - alpha )) ;
h2= [ E(1,s)^(1- alpha ) ; E(2,s)^(1-alpha) ;
E(3,s)^(1-alpha) ];
e2(:,i)=g1*h2.*h1;
end
O2(1,i) = norm(e2(1,:));
O2(2,i) = norm(e2(2,:));
end
end
if s==3
for i = 1:50
for T = R(i) : (R(i)+1000)
E(1,s) = S(T,s);
E(2,s) = S(T+1,s);
E(3,s) = S(T+2,s);
g1= [ a , -E(3,s) , 0 ; E(3,s) , b , 0 ; 1/3*E(2,s) , 0 , c ];
h1= ( gamma(2) / gamma( 2 - alpha )) ;
h2= [ E(1,s)^(1- alpha ) ; E(2,s)^(1-alpha) ;
E(3,s)^(1-alpha) ];
e3(:,i)=g1*h2.*h1;
end
O3(1,i) = norm(e3(1,:));
O3(2,i) = norm(e3(2,:));
end
end
if s==4
for i = 1:50
for T = R(i) : (R(i)+1000)
E(1,s) = S(T,s);
E(2,s) = S(T+1,s);
E(3,s) = S(T+2,s);
g1= [ a , -E(3,s) , 0 ; E(3,s) , b , 0 ; 1/3*E(2,s) , 0 , c ];
h1= ( gamma(2) / gamma( 2 - alpha )) ;
h2= [ E(1,s)^(1- alpha ) ; E(2,s)^(1-alpha) ;
E(3,s)^(1-alpha) ];
e4(:,i)=g1*h2.*h1;
end
O4(1,i) = norm(e4(1,:));
O4(2,i) = norm(e4(2,:));
end
end
end
% % 找各种状态大小 存为经典域 并扩大10%当做可拓
%
-----------------------------------------------------------------------------------------------
% % 经典
N1(1,:) = [min(O1(1,:)), max(O1(1,:))];
N1(2,:) = [min(O1(2,:)) max(O1(2,:))];
N2(1,:) = [min(O2(1,:)), max(O2(1,:))];
N2(2,:) = [min(O2(2,:)) max(O2(2,:))];
N3(1,:) = [min(O3(1,:)), max(O3(1,:))];
N3(2,:) = [min(O3(2,:)) max(O3(2,:))];
N4(1,:) = [min(O4(1,:)), max(O4(1,:))];
N4(2,:) = [min(O4(2,:)) max(O4(2,:))];
% % 可拓
EX1(:,1) = N1(:,1).*0.9;
EX1(:,2) = N1(:,2).*1.1;
EX2(:,1) = N2(:,1).*0.9;
EX2(:,2) = N2(:,2).*1.1;
EX3(:,1) = N3(:,1).*0.9;
EX3(:,2) = N3(:,2).*1.1;
EX4(:,1) = N4(:,1).*0.9;
EX4(:,2) = N4(:,2).*1.1;
%
----------------------------------------------------------------------------------
%每笔测试数据
m=10;
R=randi([48000,134000],1, m);
s=3%第几种状态
if s==3
for i = 1: m
for T = R(i) : (R(i)+1000)
E(1,s) = S(T,s);
E(2,s) = S(T+1,s);
E(3,s) = S(T+2,s);
g1= [ a , -E(3,s) , 0 ; E(3,s) , b , 0 ; 1/3*E(2,s) , 0 , c ];
h1= ( gamma(2) / gamma( 2 - alpha )) ;
h2= [ E(1,s)^(1- alpha ) ; E(2,s)^(1-alpha) ;
E(3,s)^(1-alpha) ];
e1(:,i)=g1*h2.*h1;
end
O5(1,i) = norm(e1(1,:));
O5(2,i) = norm(e1(2,:));
end
end
% % --------------------------------------------------
%计算出待测物(O5)与可拓经典域的距离
for i = 1 : m
%经典
D1(1,i) = abs( O5(1,i)-(N1(1,1)+N1(1,2))/2 ) - ((N1(1,2)-N1(1,1))/2) ;
D1(2,i) = abs( O5(2,i)-(N1(2,1)+N1(2,2))/2 ) - ((N1(2,2)-N1(2,1))/2) ;
D2(1,i) = abs( O5(1,i)-(N2(1,1)+N2(1,2))/2 ) - ((N2(1,2)-N2(1,1))/2) ;
D2(2,i) = abs( O5(2,i)-(N2(2,1)+N2(2,2))/2 ) - ((N2(2,2)-N2(2,1))/2) ;
D3(1,i) = abs( O5(1,i)-(N3(1,1)+N3(1,2))/2 ) - ((N3(1,2)-N3(1,1))/2) ;
D3(2,i) = abs( O5(2,i)-(N3(2,1)+N3(2,2))/2 ) - ((N3(2,2)-N3(2,1))/2) ;
D4(1,i) = abs( O5(1,i)-(N4(1,1)+N4(1,2))/2 ) - ((N4(1,2)-N4(1,1))/2) ;
D4(2,i) = abs( O5(2,i)-(N4(2,1)+N4(2,2))/2 ) - ((N4(2,2)-N4(2,1))/2) ;
% 可拓
E1(1,i) = abs( O5(1,i) - (EX1(1,1)+EX1(1,2))/2 ) -
((EX1(1,2)-EX1(1,1))/2) ;
E1(2,i) = abs( O5(2,i) - (EX1(2,1)+EX1(2,2))/2 ) -
((EX1(2,2)-EX1(2,1))/2) ;
E2(1,i) = abs( O5(1,i) - (EX2(1,1)+EX2(1,2))/2 ) -
((EX2(1,2)-EX2(1,1))/2) ;
E2(2,i) = abs( O5(2,i) - (EX2(2,1)+EX2(2,2))/2 ) -
((EX2(2,2)-EX2(2,1))/2) ;
E3(1,i) = abs( O5(1,i) - (EX3(1,1)+EX3(1,2))/2 ) -
((EX3(1,2)-EX3(1,1))/2) ;
E3(2,i) = abs( O5(2,i) - (EX3(2,1)+EX3(2,2))/2 ) -
((EX3(2,2)-EX3(2,1))/2) ;
E4(1,i) = abs( O5(1,i) - (EX4(1,1)+EX4(1,2))/2 ) -
((EX4(1,2)-EX4(1,1))/2) ;
E4(2,i) = abs( O5(2,i) - (EX4(2,1)+EX4(2,2))/2 ) -
((EX4(2,2)-EX4(2,1))/2) ;
end
% %
%
---------------------------------------------------------------------------------------
%关联函式 判断待测目标(O5)与四种状态的关系
% 状态一
for i = 1 : m
%状态1 e1
if O5(1,i) >=N1(1,1) && O5(1,i) <=N1(1,2)
K1(1,i) = -D1(1,i) / abs( N1(1,1)-N1(1,2) );
else
K1(1,i) = D1(1,i) / ( E1(1,i) - D1(1,i));
end
%状态1 e1
if O5(2,i) >=N1(2,1) && O5(2,i) <=N1(2,2)
K1(2,i) = -D1(2,i) / abs(N1(2,1)-N1(2,2));
else
K1(2,i) = D1(2,i) / ( E1(2,i) - D1(2,i));
end
Kk1(1,i) = K1(1,i)*0.5+ K1(2,i)*0.5;
end
Kk1;
% 状态二
for i = 1 : m
%状态2 e1
if O5(1,i) >=N2(1,1) && O5(1,i) <=N2(1,2)
K2(1,i) = -D2(1,i) / abs( N2(1,1)-N2(1,2) );
else
K2(1,i) = D2(1,i) / ( E2(1,i) - D2(1,i));
end
%状态2 e1
if O5(2,i) >=N2(2,1) && O5(2,i) <=N2(2,2)
K2(2,i) = -D2(2,i) / abs(N2(2,1)-N2(2,2));
else
K2(2,i) = D2(2,i) / ( E2(2,i) - D2(2,i));
end
Kk2(1,i) = K2(1,i)*0.5+ K2(2,i)*0.5;
end
Kk2;
% 状态三
for i = 1 : m
%状态3 e1
if O5(1,i) >=N3(1,1) && O5(1,i) <=N3(1,2)
K3(1,i) = -D3(1,i) / abs( N3(1,1)-N3(1,2) );
else
K3(1,i) = D3(1,i) / ( E3(1,i) - D3(1,i));
end
%状态3 e2
if O5(2,i) >=N3(2,1) && O5(2,i) <=N3(2,2)
K3(2,i) = -D3(2,i) / abs(N3(2,1)-N3(2,2));
else
K3(2,i) = D3(2,i) / ( E3(2,i) - D3(2,i));
end
Kk3(1,i) = K3(1,i)*0.5+ K3(2,i)*0.5;
end
Kk3;
for i = 1 : m
%状态4 e1
if O5(1,i) >=N4(1,1) && O5(1,i) <=N4(1,2)
K4(1,i) = -D4(1,i) / abs( N4(1,1)-N4(1,2) );
else
K4(1,i) = D4(1,i) / ( E4(1,i) - D4(1,i));
end
%状态4 e2
if O5(2,i) >=N4(2,1) && O5(2,i) <=N4(2,2)
K4(2,i) = -D4(2,i) / abs(N4(2,1)-N4(2,2));
else
K4(2,i) = D4(2,i) / ( E4(2,i) - D4(2,i));
end
Kk4(1,i) = K4(1,i)*0.5+ K4(2,i)*0.5;
end
Kk4;
F=[Kk1;Kk2;Kk3;Kk4]';
%正规化
for j= 1:m
for i=1:4
Q(j,i)=(2*F(j,i)-max(F(j,1:4))-min(F(j,1:4)))/(max(F(j,1:4))-min(F(j,1:4)));
end
end
Q
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 111.250.212.32 (台湾)
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/MATLAB/M.1613924874.A.2E1.html
1F：推 book5601: 对於这种请别人debug 的文章，好像都乏人问津。劝你各位 02/27 11:08
2F：→ book5601: 做论文还是要多主动思考把大问题分割成数个小问题，确认 02/27 11:08
3F：→ book5601: 每个环节都在自己的掌握之後再连结各个小问题，问题应该 02/27 11:08
4F：→ book5601: 就解决了。但不能排除有许多条件是论文没有提及的，最快 02/27 11:08
5F：→ book5601: 且最不出错的方法是问本人。 02/27 11:08
6F：嘘 jatj: ...... 02/27 23:01
7F：→ guteres: 推一楼，我想作者都有留连络方式，我们以前都是直接问 04/08 21:38
8F：推 el3qu4: 偷偷跟你说 论文都会暗藏小地方不可能全部都给你看....问 04/09 12:47
9F：→ el3qu4: 原作者可能好一点 04/09 12:47