作者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/m.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