[問題] 照著論文做狀態分析但就是無法得出一樣

看板MATLAB作者dennisshyu12 (路人)時間3年前 (2021/02/22 00:27)推噓1(2推 1噓 6→)

留言9則, 4人參與討論串1/1

想請教版上的大神，是否是小弟程式有誤導至結果不同，或者是公式的部分打錯了... 請大神幫幫忙 [狀況說明] 四種狀態數據分數階導數後，做可拓的物質元素分析，皆照著步驟與公式做，但結果就是與論文不同，達不到百分之百，論文名: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://www.ptt.cc/bbs/MATLAB/M.1613924874.A.2E1.html