Nutzung der gewonnenen Meßergebnisse zur Systemidentifikation des Einachsers mittels der Methode der kleinsten Quadrate - Teil 2
(EN google-translate)
(PL google-translate)
002identifikation.zip - Skripte und Daten zu nachfolgender Darstellung.
Identifikation
Die Daten liegen nun in den drei Matrizen datalsq1.mat, datalsq2.mat, datalsq3.mat mit folgenden Spalten vor:
Spalte Nr.: 1 2 3 4 5 6 Bedeutung: t/s | pwm | phi/rad | omega/(rad/s) | alfa/(rad/s^2) | dalfa/dt/(rad/s^3)
Code 0-1: Spalten der für die LSQ aufbereiteten Datenmatrizen
Wie zuvor dargestellt, erfolgt die Identifikation auf Grundlage eines PT2-Übertragungsgliedes bezüglich der Kippwinkelgeschwindigkeit ω:
Bild 0-1: Geschätzte Struktur der Systemmatrix der Regelstrecke.
$ \dot \alpha = p \cdot \omega + q \cdot \alpha + r \cdot pwm $
Formel 0-1: Algebraische Gleichung für LSQ
$ a= \frac {1}{r} & b=-q \cdot a & c=-p \cdot a $
Formel 0-2: Umrechnungen
Das entsprechende überbestimmte lineare Gleichungssystem zur Durchführung des LSQ mit Scilab ist dann:
$ B \cdot c = d $
Formel 0-3: Überbestimmtes LGS.
...mit:
$ B = \left[\begin{array}{cc} \omega _0 & \alpha _0 & pwm_0 \\ \omega _1 & \alpha _1 & pwm_1 \\ \omega _2 & \alpha _2 & pwm_2 \\ ... & ... & ...\end{array}\right] $
Formel 0-4: Verläufe der Zus \tan dsgrößen.
$ c = \left[\begin{array}{cc}p \\ q \\ r\end{array}\right] $
Formel 0-5: Unbekannte Parameter.
$ d = \left[\begin{array}{cc}\dot \alpha _0 \\ \dot \alpha _1 \\ \dot \alpha _2 \\ ...\end{array}\right] $
Formel 0-6: Inhomogenität.
Identifikation mittels Scilab
clear();
matrix = fscanfMat('datalsq1.mat');
zs = size(matrix);
zeilen = zs(1);
//omega alfa pwm
BB=[matrix(:,4),matrix(:,5),matrix(:,2)];
//pwm zeitverzögert in Wirksamkeit:
//BB=[matrix(:,4),matrix(:,5),[0;matrix(1:zeilen-1,2)]];
dd=matrix(:,6);
cc=lsq(BB,dd);
p=cc(1)
q=cc(2)
r=cc(3)
a=1/r
b=-q*a
c=-p*a
//Testweise das identifizierte System modellieren,
//mit den pwm-Werten aus der Identifikation beaufschlagen und
//mit dem gemessenen Verlauf von phi vergleichen:
A = [
0 1 0
0 0 1
0 -c/a -b/a
];
B = [ 0;0;1/a];
pwm = [matrix(:,2);matrix(zeilen,2)];
t = matrix(:,1)';
function f = rechteSeite(t,y)
// t_index = 1+round(t*100);
t_index = 1+round(t*100) +1; //auch hier verzögert wirksam
u=pwm(t_index);
f = A*y+B*u;
endfunction
y0 = [matrix(1,3),matrix(1,4),matrix(1,5)]';
t0 = 0;
y = ode(y0,t0,t,rechteSeite);
plot2d(t,y(1,:)); //Simulationsverlauf von phi
plot2d(t,matrix(:,3)); //gemessen
Code 0-2: Identifikation mittels LSQ und Verifikation, Skript ident1.sce.
Bild 0-2: Kippwinkel φ in Simulation und Messung: Simuliertes System schwingt auf.
Verbesserung: Das pwm-Signal ist aufgrund der Trägheit / Totzeit des Systems später wirksam, als es berechnet wird. Testweise wird ein Zeitschift um 0,01s für das pwm-Signal angenommen (Skript ident1b.sce):
BB=[matrix(:,4),matrix(:,5),[0;matrix(1:zeilen-1,2)]];
Code 0-3: Zeitversatz in der Wirksamkeit des PWM-Signals einführen.
Bild 0-3: Kippwinkel φ in Simulation und Messung: Simuliertes System driftet weg, schwingt aber nicht mehr auf.
Vermeidung unerlaubter Vorzeichen bei a bei der Identifikation durch Vorzeichenwechsel beim pwm-Signal (Skript ident1c.sce).
Bild 0-4: Kippwinkel φ in Simulation und Messung: Simuliertes System driftet weg, schwingt aber nicht mehr auf.
p=-2335.9637 q=-3.8490007 r=0.7723642 a=1.2947259 b=4.983401 c=3024.4327
Code 0-4: Identifizierte Parameter
Vergleich der Ergebnisse mit den anderen beiden Datensätzen
|
Bild 0-5: Drift mit datalsq2.mat, Skript ident2c.sce
p=-2396.7257 q=-5.2015687 r=0.7153613 a=1.3978951 b=7.2712475 c=3350.3711
Code 0-5: Identifizierte Parameter
|
Bild 0-6: Drift mit datalsq3.mat, Skript ident3c.sce
p=-2413.9817 q=-4.5329396 r=0.8706732 a=1.1485366 b=5.2062468 c=2772.5462
Code 0-6: Identifizierte Parameter
|
p=-2374.1533 q=-4.3571659 r=0.7734549 a=1.2929002 b=5.6333807 c=3069.5433
Code 0-7: Identifizierte Parameter