Akademia Górniczo-Hutnicza
Katedra Robotyki i Mechatroniki
Identyfikacja i analiza sygnałów
Projekt
Damian GÅ‚Ä…b
Gr. II
Mechatronika
Ad.1 Schemat przyjętego układu
Ad.2 Równania dynamiki założonego układu w postaci równań czasowych
Ad.3 Równania w postaci macierzowej
gdzie:
· [M]- jest macierzÄ… mas
· [C]  jest macierzÄ… tÅ‚umieÅ„
· [K]  jest macierzÄ… sztywnoÅ›ci
· {f(t)}  jest wektorem siÅ‚
· {x(t)}  jest wektorem odpowiedzi
Ad.4
Przyjęte parametry:
Kod programu:
clear all
clc
% Wartości współczynników
%mas
m1=8;
m2=3;
m3=6;
m4=4;
m5=5;
m6=6;
%sztywności
k1=1000;
k2=20500;
k3=1700;
k4=25000;
k5=9000;
k6=5000;
k7=28000;
k8=7200;
k9=10000;
%tłumienia
c1=5;
c2=3;
c3=7;
c4=8;
c5=4;
c6=2;
c7=3;
c8=4;
c9=3;
M=[m1 0 0 0 0 0;
0 m2 0 0 0 0;
0 0 m3 0 0 0;
0 0 0 m4 0 0;
0 0 0 0 m5 0;
0 0 0 0 0 m6];
C=[(c1+c2) -c1 -c2 0 0 0;
-c1 (c1+c3+c4) 0 -c3 -c4 0;
-c2 0 (c2+c5+c6) 0 -c5 -c6;
0 -c3 0 (c3+c7) 0 0;
0 -c4 -c5 0 (c4+c5+c8) 0;
0 0 -c6 0 0 (c6+c9)];
K=[k1+k2 -k1 -k2 0 0 0;
-k1 k1+k3+k4 0 -k3 -k4 0;
-k2 0 k2+k5+k6 0 -k5 -k6;
0 -k3 0 k3+k7 0 0;
0 -k4 -k5 0 k4+k5+k8 0;
0 0 -k6 0 0 k6+k9];
% Budowanie macierzy do równań stanu
ZER = zeros(size(M));
A = [ZER,M;M,C];
B = [-M,ZER;ZER,K];
% Rozwiązanie uogólnionego zagadnienia wlasnego
[PHI,LAMBDA]=eig(-B,A);
% Czestotliwosci drgan wlasnych tlumionych [Hz]
WD=imag(diag(LAMBDA))/2/pi;
% Czestotliwosci drgan wlasnych [Hz]
WW=sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)/2/pi
% Tlumienie modalne
KSI=-real(diag(LAMBDA))./sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)
Otrzymane wyniki:
Częstotliwości drgań własnych układu symulacyjnego wynoszą więc odpowiednio:
19.7331 Hz, 13.9738 Hz, 13.6894 Hz, 3.2862 Hz, 7.7903 Hz, 7.9998 Hz
Współczynniki tłumienia modalnego wynoszą:
0.0294, 0.0126, 0.0126, 0.0085, 0.0252, 0.0120
Ad.5 Symulacja układu przy wymuszeniu białym szumem na masie M1
Kod programu:
clear all
clc
% Wartości współczynników
%mas
m1=8;
m2=3;
m3=6;
m4=4;
m5=5;
m6=6;
%sztywności
k1=1000;
k2=20500;
k3=1700;
k4=25000;
k5=9000;
k6=5000;
k7=28000;
k8=7200;
k9=10000;
%tłumienia
c1=5;
c2=3;
c3=7;
c4=8;
c5=4;
c6=2;
c7=3;
c8=4;
c9=3;
A=[0 1 0 0 0 0 0 0 0 0 0 0;
-((k1+k2)/m1) -((c1+c2)/m1) k1/m1 c1/m1 k2/m1 c2/m1 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0;
k1/m2 c1/m2 -((k1+k3+k4)/m2) -((c1+c3+c4)/m2) 0 0 k3/m2 c3/m2 k4/m2 c4/m2 0 0;
0 0 0 0 0 1 0 0 0 0 0 0;
k2/m3 c2/m3 0 0 -((k2+k5+k6)/m3) -((c2+c5+c6)/m3) 0 0 k5/m3 c5/m3 k6/m3 c6/m3;
0 0 0 0 0 0 0 1 0 0 0 0;
0 0 k3/m4 c3/m4 0 0 -((k3+k7)/m4) -((c3+c7)/m4) 0 0 0 0;
0 0 0 0 0 0 0 0 0 1 0 0;
0 0 k4/m5 c4/m5 k5/m5 c5/m5 0 0 -((k4+k5+k8)/m5) -((c4+c5+c8)/m5) 0 0;
0 0 0 0 0 0 0 0 0 0 0 1;
0 0 0 0 k6/m6 c6/m6 0 0 0 0 -((k6+k9)/m6) -((c6+c9)/m6)];
B=[0; 1/m1; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
C=[1 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 0 0 1 0];
D=[0; 0; 0; 0; 0; 0];
sys=ss(A,B,C,D); % wyznaczenie parametrów modalnych modelu
t=0:0.01:50; % czas
F=2*rand(1,length(t)); % wymuszenie jako bialy szum
Y_out=lsim(sys,F,t); % odpowiedz na biały szum
x1=Y_out(:,1);
x2=Y_out(:,2);
x3=Y_out(:,3);
x4=Y_out(:,4);
x5=Y_out(:,5);
x6=Y_out(:,6);
figure(1)
plot(t,x1)
title('Przemieszzcenie masy M1 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
figure(2)
plot(t,x2)
title('Przemieszzcenie masy M2 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
figure(3)
plot(t,x3)
title('Przemieszzcenie masy M3 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
figure(4)
plot(t,x4)
title('Przemieszzcenie masy M4 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
figure(5)
plot(t,x5)
title('Przemieszzcenie masy M5 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
figure(6)
plot(t,x6)
title('Przemieszzcenie masy M6 w czasie')
xlabel('Czas [sec]')
ylabel('przemieszczenie [mm]')
Odpowiedzi w postaci wykresów przemieszczeń w czasie mas:
Ad.6 Identyfikacja nieparametryczna układu.
Kod:
%nieparametryczna
window=hanning(512);
nfft=512;
noverlap=128;
dflag='none';
figure(1)
subplot(3,1,1)
psd(x1,nfft, Fs, window, noverlap, dflag);
subplot(3,1,2)
psd(x2,nfft, Fs, window, noverlap, dflag);
subplot(3,1,3)
psd(x3,nfft, Fs, window, noverlap, dflag);
figure(2)
subplot(3,1,1)
psd(x4,nfft, Fs, window, noverlap, dflag);
subplot(3,1,2)
psd(x5,nfft, Fs, window, noverlap, dflag);
subplot(3,1,3)
psd(x6,nfft, Fs, window, noverlap, dflag);
Ad.7 Identyfikacja parametryczna układu.
Wymuszenie na masie 1 odpowiedz
na masie 5.
Kod programu:
[Wn,Z] = damp(sys);
dat = IDDATA(x5,F',1/Fs);
model = armax([dat],[20 20 2 1]);
model_z = zpk(model)
[Wn_es, Z_es P] = damp(model_z);
[Wn/2/pi Z]
[Wn_es/2/pi Z_es]
Odpowiedzi:
ans =
3.2862 0.0085
3.2862 0.0085
7.7903 0.0252
7.7903 0.0252
7.9998 0.0120
7.9998 0.0120
13.6894 0.0126
13.6894 0.0126
13.9738 0.0126
13.9738 0.0126
19.7331 0.0294
19.7331 0.0294
ans =
3.2862 0.0085
3.2862 0.0085
7.7903 0.0252
7.7903 0.0252
7.9998 0.0120
7.9998 0.0120
13.6894 0.0126
13.6894 0.0126
13.9738 0.0126
13.9738 0.0126
19.7331 0.0294
19.7331 0.0294
35.1097 0.0246
35.1097 0.0246
42.6818 0.0230
42.6818 0.0230
49.7358 0.0211
49.7358 0.0211
56.6022 0.0191
56.6022 0.0191
Ad.8 Skalogram układu wymuszonego białym szumem (dla jednego punktu odpowiedzi)
wykorzystując ciągłą transformatę falkową.
Kod programu:
% dobór wektora skali
a=scal2frq(20,'cmor 6-5',1/Fs)
b=scal2frq(400,'cmor 6-5',1/Fs)
ws=20:1:400;%wektor skali
f = cwt(x5,ws, 'cmor 6-5');
c=size(f);
time = repmat(t,c(1,1),1);
fr = scal2frq(ws,'cmor 6-5',1/Fs);
freq = repmat(fr,length(t),1);
surf(time,freq',abs(f))
shading interp
view([65 45])
ylabel('Czestotliwosc','FontSize',16)
xlabel('Czas','FontSize',16)
Skalogram:
Ad.9 Wnioski:
Dzięki zastosowaniu 3ch metod identyfikacji otrzymujemy informacje na temat układu. Identyfikacja
nieparametryczna podaje nam informacje o widmie mocy poszczególnych częstotliwości w analizowanym
sygnale. Układ zawierał 2 pary częstotliwości własnych różniących się o mniej niż 0.5 Hz. Te częstotliwości
zostały zidentyfikowane jako jedna ( zlały się ). Dokładność metody nieparametrycznej zależy od długości i
częstotliwości próbkowania sygnału. Metoda parametryczna (zastosowanie modelu ARMAX) dała dokładne
częstotliwości własne układu a także współczynniki tłumienia układu. Ciągła transformata falkowa podaje
bardzo interesujacą informacje o rozkładzie częstotliwosci w czasie trwania sygnału. Częstotliwości
minimalnie różniące się są praktycznie nie do rozróżnienia.