Ich ermittle die Impulsantwort eines PT1-Gliedes mit Hilfe einer numerischen Simulation. Aus dem Vektor der Impulsantwort bestimme ich über die diskrete Fourier Transformation die Übertragungsfunktion. Soweit so gut.
Da ich die Übertragungsfunktion nur im Frequenzbereich F_min bis F_max bestimmen muss, und ich den Speicher- und Rechenaufwand möglichst minimieren möchte, habe ich mir für die erforderliche zeitliche Auflösung "dt" und die Aufnahmedauer "T" der Impulsantwort Folgendes überlegt:
dt = 1/(2*F_max) gemäß Nyquist-Theorem
T = 1/F_min
Lasse ich mir die resultierende Übertragungsfunktion plotten erhalte ich ein merkwürdiges Ergebnis:
Code:
% Konstanten
A = 1; % Amplitudenverhältnis des PT1-Gliedes
B = 15; % Zeitkonstante des PT1-Gliedes
F_min = 1E-4; % Frequenzauflösung bzw. kleinste Frequenz der Ü-Funktion
F_max = 1E+0; % Maximale Frequenz der Übertragungsfunktion
% Zeit- und Frequenzvektoren
dt = 1/(2*F_max); % Zeitabstand zwischen 2 Werten der Impulsantwort
T = 1/F_min; % Gesamt-Zeitdauer der Impulsantwort
t = 0:dt:T; % Zeitvektor der Impulsantwort
dF = 1/T; % Frequenzauflösung der Übertragungsfunktion
F = 0:dF:1/dt; % Frequenzvektor der Übertragungsfunktion
% Rekonstruktion der Übertragungsfunktion aus Impulsantwort
h = (A/B)*exp(-t./B); % Impulsantwort des PT1-Gliedes
H = fft(h)*dt; % Übertragungsfunktion des PT1-Gliedes
Der Phasengang müsste ohne Umwege asymptotisch zu -pi/2 gehen und das Amplitudenverhältnis sollte bei niedriger Frequenz genau 1 sein. Erhöht man die zeitliche Auflösung der Impulsantwort um das 100-fache (dt = 1/(200*F_max)) ergibt sich ein korrektes, aber ressourcenhungriges Ergebnis:
Code:
% Konstanten
A = 1; % Amplitudenverhältnis des PT1-Gliedes
B = 15; % Zeitkonstante des PT1-Gliedes
F_min = 1E-4; % Frequenzauflösung bzw. kleinste Frequenz der Ü-Funktion
F_max = 1E+0; % Maximale Frequenz der Übertragungsfunktion
% Zeit- und Frequenzvektoren
dt = 1/(200*F_max); % Zeitabstand zwischen 2 Werten der Impulsantwort
T = 1/F_min; % Gesamt-Zeitdauer der Impulsantwort
t = 0:dt:T; % Zeitvektor der Impulsantwort
dF = 1/T; % Frequenzauflösung der Übertragungsfunktion
F = 0:dF:1/dt; % Frequenzvektor der Übertragungsfunktion
% Rekonstruktion der Übertragungsfunktion aus Impulsantwort
h = (A/B)*exp(-t./B); % Impulsantwort des PT1-Gliedes
H = fft(h)*dt; % Übertragungsfunktion des PT1-Gliedes
Du hast zwar das Abtasttheorem beachtet, aber in der Praxis ist der Faktor 2 eben viel zu wenig. Besonders da du eine sehr hohe Auflösung im niederfreq. Bereich haben möchtest. Außerdem ist deine Frequenzachse sowie die Skalierung der Amplitude falsch. Schau mal in der Skriptecke unter "umfassendes FFT Bsp." wie es richtig sein müsste. Du kannst mit der FFT auch immer nur bis zur halben Abtastfreq. darstellen, was du hier ebenfalls mißachtest.
Hast du die Control System und Signal Processing Toolbox für Matlab? Da würde es nämlich mit tf, bode etc. oder dem sptool deutlich einfacher gehen.
Ebenfalls sollte dir bewusst sein, dass einem kontinuierlichen PT1 entsprichst, du hier aber mit einem diskreten System arbeiten möchtest/musst. Da gibt es Unterschiede
Du hast zwar das Abtasttheorem beachtet, aber in der Praxis ist der Faktor 2 eben viel zu wenig.
Soweit sogut, aber nach welchem Kriterium lege ich die Abtastrate fest? Ist da Ausprobieren angesagt? (wäre ja iterativ möglich, bis Konvergenz erreicht)
Zitat:
Außerdem ist deine Frequenzachse sowie die Skalierung der Amplitude falsch.
Da sehe ich den Fehler nicht, aber darum gehts ja jetzt auch nicht primär...
Zitat:
Hast du die Control System und Signal Processing Toolbox für Matlab? Da würde es nämlich mit tf, bode etc. oder dem sptool deutlich einfacher gehen.
Ja, allerdings will ich das Programm ohnehin später in Fortran umsetzen, sodass ich das Problem irgendwie lösen muss.
Hallo Eric,
habe gerade im Buch "Regelungstechnik 2" von Jan Lunze(S. 420-S.421) reingesehen.
Er schreibt, dass sich in der Praxis Abtastfrequenz= 6...20 mal Grenzfrequenz als sinnvoll erwiesen hat.
Die Grenzfrequenz ist die größte im Regelkreis auftretende Frequenz.
Da du nur eine Regelstrecke hast ist in deinem Fall die Zeitkonstante deines Pt1-Systems die einzige Zeitkonstante die hierbei berücksichtigt werden muss.
Es hat natürlich auch einen Grund, warum das so lange dauert. Ist die Anzahl deiner Messwerte keine 2er Potenz (wie auch hier in dem Bsp.) so wird der deutlich langsamere DFT Algorith. bei dem fft() Befehl verwendet.
Ich habe mal verschiedene Sachen probiert, aber der Grund für den abweichenden Phasengang kann ich nicht finden. Wenn ich allerdings die Impulsantwort hd der diskreten Übertragungsfkt. an die fft() übergebe, erhalte ich das gleiche Ergebnis, wie auch beim Bodeplot. Hier stimmen dann Amplituden- und Phasengang überein. Verwendet man aber h oder hs, stimmt der Amplitudengang immer mit der diskreten Impulsantwort hd überein. Der Phasengang über fft(h) oder fft(hs) weicht bei beiden aber gegenüber dem Bodebefehl ab. Ebenfalls verstehe ich das hier mit der Abtastfreq. noch nicht ganz. Der Default beim Bodeplot ist Fs = 1, wobei dann alle 3 Amplitudengänge bis weit über die Eckfreq. hinaus übereinstimmen. Wenn man Fs für das diskrete System verringert, nähert sich der Phasengang immer mehr dem kontinuierlichen PT1 an, was ich schon nicht nachvollziehen kann. Es müsste doch umgekehrt sein.
Code:
clear;
% Konstanten
A = 1; % Verstärkung des PT1-Gliedes
B = 15; % Zeitkonstante des PT1-Gliedes
F_min = 1E-4; % Frequenzauflösung bzw. kleinste Frequenz der Ü-Funktion
F_max = 1E+0; % Maximale Frequenz der Übertragungsfunktion
N = 2^14; % Anzahl Abtastwerte, immer 2er Potenz nehmen da sonst der langsame DFT Alg. genutzt wird
Fs = 1; % Abtastfreq.
Fn = 0.5*Fs; % Nyquistfreq.
df = Fs/N; % Frequenzauflösung des Spektrums
% Zeit- und Frequenzvektoren
Ts = 1/Fs; % Abtastrate
t = 0:Ts:(N-1)*Ts; % Zeitvektor der Impulsantwort
% Rekonstruktion der Übertragungsfunktion aus Impulsantwort
h = (A/B)*(exp(-t./B)); % Impulsantwort des PT1-Gliedes
h = h/sum(h); % auf 1 normieren % Test: Impulsantwort ist die Ableitung der Sprungantwort
%sprung = A*(1 - exp(-t./B));
% numerische Ableitung der Sprungantwort
%h = diff(sprung); % Länge N-1
%h = [h, 0]; % auf Länge N erweitern
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS