Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • Adam.Schrey/lecture-tutorials
1 result
Select Git revision
Loading items
Show changes
Commits on Source (62)
Show whitespace changes
Inline Side-by-side
Showing
with 14 additions and 13899 deletions
Index.ipynb
View file @ 9415dc13
%% Cell type:markdown id: tags:
# Systemtheorie 2 und JupyterHub
## Wie greife ich auf JupyterHub zu?
Loggen Sie sich unter <a href="https://jupyter.rwth-aachen.de/hub/spawn?profile=sys2"> https://jupyter.rwth-aachen.de/hub/spawn?profile=sys2 </a> mit Ihrem RWTH-Account ein.
## Wie nutze ich JupyterHub?
Auf der Weboberfläche sollten auf der linken Seite verschiedene Ordner angezeigt werden. Unter dem Pfad "lecture-tutorials/notebooks" sind Dateien (Notebooks) zum Thema "Systemtheorie 2".
Auf der Weboberfläche sollten auf der linken Seite verschiedene Ordner angezeigt werden. Unter dem Pfad "sys2/sys2-jupyter-notebooks/lecture_examples" sowie "sys2/sys2-jupyter-notebooks/exam_examples" sind Dateien (Notebooks) zum Thema "Systemtheorie 2".
In diesen Notebooks können Beispiele in Form von Code-Blöcken ausgeführt werden. Meist geben diese einen statischen Graphen aus. Jedoch gibt es auch oft interaktive Elemente, bei denen Werte über einen Schieberegler oder Knopf verändert werden können.
In diesen Notebooks können Beispiele in Form von Code-Blöcken ausgeführt werden. Meist geben diese einen statischen Graphen aus. Jedoch gibt es auch interaktive Elemente, bei denen Werte über einen Schieberegler oder Knopf verändert werden können.
Die Code-Blöcke können im Notebook geändert werden und diese Änderung kann unter dem Unterpunkt "Git" auch wieder rückgängig gemacht werden.
"V01_zeitdiskrete_systeme.ipynb" bis "V10_kalmanfilter_demonstration.ipynb" sind Notebooks zu den Themen der Vorlesung. In "Sandkasten.ipynb" kann man beliebige Funktionen untersuchen, eigene Simulationen erstellen und Verschiedenes ausprobieren.
"In der Template-Datei "sys2/sys2-jupyter-notebooks/template.ipynb" kann man zudem beliebige Funktionen untersuchen, eigene Simulationen (mittels Python 3 Syntax) erstellen und Verschiedenes ausprobieren.
### Ausführen der Code-Blöcke in den Folien:
<img src="introduction_pictures/run.png"/>
### Beispiel eines interaktiven Elementes:
<img src="introduction_pictures/interact.png"/>
### Änderungen rückgängig machen:
<img src="introduction_pictures/git.png"/>
<img src="introduction_pictures/discard_changes.png"/>
## Python und Matlab
Ältere Notebooks, welche Matlab-Code verwenden, können unter dem Pfad "lecture-tutorials/notebooks/alt" gefunden werden.Aber auch andere Notebooks können Matlab-Code beinhalten.Wenn man zwischen verschiedenen Notebooks mit Python-Code und Matlab-Code wechselt, muss man den Kernel womöglichmanuell ändern (Octave Kernel für Matlab):
## Python vs. Matlab
Einige der Jupyter-Notebooks können Matlab-Code beinhalten. Wenn man zwischen verschiedenen Notebooks mit Python-Code und Matlab-Code wechselt, muss man den Kernel gegebenenfalls manuell umstellen (Octave Kernel für Matlab; Python 3 (ipykernel) für Python).
### Ändern des Kernels:
<img src="introduction_pictures/switch_kernel.png"/>
## Fehlerbehebung
Sollte es zu Fehlern kommen oder das System in eine Endlosschleife geraten, so kann man unter dem Menü-Punkt "Kernel" >> "Restart Kernel" die aktuellen Berechnungen und gespeicherten Variablen zurücksetzen.
Sollten die Notebooks oder zugehöriger Python-Code fehlerhaft sein, können Sie dies gerne melden.
Auch können Sie das Systemtheorie 2 Jupyter auf "Werkseinstellungen" zurücksetzen (Achtung: alle selbst getätigten Änderungen gehen dabei verloren!) indem Sie unter dem Menü-Punkt "Git" >> "Open Git Repository in Terminal" den Befehl <code>git fetch origin && git reset --hard origin/master</code> eingeben.
Sollten die Notebooks oder zugehöriger Code trotzdem nicht richtig funktionieren, so können Sie dies jederzeit gerne an "acs-teaching-sys2@eonerc.rwth-aachen.de" melden.
......
interactive-slides/Vorlesung 01 - Zeitdiskrete Systeme.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 02 - Die z-Transformation.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 04 - Regelung diskreter Systeme.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 06 - Normalformen in Zustandsraumdarstellung.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 07 - Steuerbarkeit und Beobachtbarkeit .ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 08 - Regelung im Zustandsraum.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 09 - Zustandsrekonstruktion mittels Beobachter.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
interactive-slides/Vorlesung 10 - Kalman Filter.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
notebooks/alt/V01_zeitdiskrete_systeme.ipynb deleted 100644 → 0
View file @ 68879533
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V1-Zeitdiskrete Systeme </span>
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b>Wie kann sich die Abtastfrequenz auf eine Zahlenfolge auswirken?</b>
- Wir nehmen folgende Abtastfrequenz: $f_{a}=10.000 Hz$
- Hier gilt: $T=\frac{1}{f_{a}}$ , $K_{max}=20.000$ , $k=[1...K_{max}]$ , $t=Tk$.
- Wähle einen Ton mit 440 Hz. Es gilt:$\omega=2\pi440$
- Die Zahlenfolge wird definiert mit:$y=sin(\omega.t)$
<b> Codes: </b>
<br><br>
Folgender Matlab Befehl erzeugt die Wiedergabe des Tons <code>y</code> mit der Frequenz <code>fa</code>:
%% Cell type:code id: tags:
``` octave
% Abstast Frequenz
fa = 10000;
T = 1/fa;
maxk = 20000;
k=1:maxk;
t=k.*T;
% A Tone 440 Hz
om = 2*pi*440;
y = sin(om*t);
```
%% Cell type:code id: tags:
``` octave
soundsc(y,fa)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Ändern Sie nun die Abtastfrequenz wie folgt:
%% Cell type:code id: tags:
``` octave
soundsc(y,fa/2)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Der Ton wurde <span style='color:red'> fehlerhaft </span> übersetzt und weicht von dem Originalton ab!
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #2 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b> Abtastung und Quantisierung</b>
- Die Frequenzen für das Signal 1 und 2 werden gewählt: $\omega_{1}=1$, $\omega_{2}=1+2\pi$
- Die Abtastkreisfrequenz und Abtastzeit entsprechen:$\omega_{T}=2\pi$, $T=\frac{2\pi}{\omega_T}$
- Des Weiteren wird für die Folge bestimmt:$K_{max}=200$, $k=[1...K_{max}]$ , $t=Tk$.
Im Folgenden erzeugen wir zwei abgetastete Signale <code>y1</code> und <code>y2</code>:
%% Cell type:code id: tags:
``` octave
% Frequenz Signal 1
om1 = 1;
% Frequenz Signal 2
om2 = 1+2*pi;
% Abtast Frequenz
om_T = 2*pi;
% Abtast Zeit
T = 2*pi/om_T;
maxk = 200;
k=1:maxk;
t=k.*T;
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Das Ergebnis kann mit folgendem Matlab Befehl geplottet werden:
%% Cell type:code id: tags:
``` octave
%Signal 1
y1 = sin(om1*k*T);
%Signal 2
y2 = sin(om2*k*T);
plot(t,y1,t,y2)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Durch Abtastung erhält man nur eine unvollständige Information über das betrachtete
Signal. Es gibt sehr viele Funktionen die durch die Abtastwerte verlaufen können.
<br>
<span style='color:Gray'>Beispiel </span> : Es werden zwei sinusförmige Signale betrachtet:
<br>
$y_1(t)=sin(\omega_1t)$ , $y_2(t)=sin(\omega_2t)$
die auf die Abtastfolgen führen:
<br> $y_1(k)=sin(k\omega_1t)$ , $y_2(k)=sin(k\omega_2t)$
Die Abtastwerte für beide Funktionen sind an allen Abtastpunkten gleich wenn
$sin(k\omega_1T)=sin(k\omega_2T)$
$(k=0, 1,...)$ gilt, was in
$\omega_1=\omega_2+l\frac{2\pi}{kT}$
umgeformt werden kann.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
- $\omega_1=\omega_2+l\frac{2\pi}{kT}$
($l$ ganzzahlig)
<br><br>
- $\omega_1=\omega_2+\omega{T}$
($ \omega_T$ Abtastkreisfrequenz)
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Erfüllen zwei Frequenzen die obere Beziehung, so können die von der Sinusschwingung
mit der Kreisfrequenz $\omega_1$ erzeugten Abtastwerte auch durch eine Sinusschwingung mit der
Kreisfrequenz $\omega_2$ erzeugt werden.
<span style='color:OrangeRed'>Das heißt Signale mit diesen beiden Frequenzen können hinter dem Abtaster nicht mehr
unterschieden werden! </span>
notebooks/alt/V02_die_z-transformation.ipynb deleted 100644 → 0
View file @ 68879533
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V2 - Die z-Transformation</span>
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Gegeben ist die Funktion im Zähler mit: $8z^{-1}$
Des Weiteren ist die Nennerfunktion gegeben mit:$1-3z^{-1}+2z^{-2}$
%% Cell type:code id: tags:
``` octave
N = [8 0]
D = [1 -3 2]
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Da <code>N</code> und <code>D</code> Vektorpolynome sind , entspricht eine Entfaltung(.eng.: deconvolution) einer einfachen
Division. Die Matlab Funktion hierfür lautet:
%% Cell type:code id: tags:
``` octave
[fo R] = deconv(N,D)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Bei erstmaliger Anwendung der Funktion erhalten wir den ersten Koeffizienten der Folge <code>fo</code> und den
Rest ( eng.: remainder) <code>R</code>.
Der Rest <code>R</code> wird mit $z^{-1}$ multipliziert und weiter dividiert
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f2 R] = deconv(R,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f3 R] = deconv(R,D)
```
%% Cell type:markdown id: tags:
### Die Potenzreihenentwicklung
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<br><br>
<b> Mathematische Berechnungen: </b>
<br><br>
Es ist $f(k)$ für $k=0, 1, 2,...$ zu bestimmen, wenn
$F_z(z)= \frac{8z^{-1}}{(z-1)(z-2)}$
egeben ist.
<span style='color:OrangeRed'> Lösung:</span>
Durch Umschreiben folgt:
$F_z(z)= \frac{8z^{-1}}{1-3z^{-1}+2z^{-2}}$
Durch Division erhält man:
$F_z(z)=8z^{-1}+24z^{-2}+56z^{-3}+120z^{-4}+...$
Aus dieser Beziehung können direkt die Werte der Zahlenfolge $f(k)$ abgelesen werden:
$f(0)=0$, $f(1)=8$, $f(2)=24$, $f(3)=56$, $f(4)=120$,...
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #2 </span>
%% Cell type:code id: tags:
``` octave
N = [8 2 3 2 1];
D = [1 0 0 0 0];
```
%% Cell type:code id: tags:
``` octave
[fo R] = deconv(N,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f1 R] = deconv(R,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f2 R] = deconv(R,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f3 R] = deconv(R,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f4 R] = deconv(R,D)
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f5 R] = deconv(R,D)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
- Die Koeffizienten der Folge <code>f</code> entsprechen den Elementen des Zählervektors <code>N</code>.
- <code>R</code> wird lediglich = $ 0$.
- Gezeigt wurde ein <span style='color:OrangeRed'>„Finite Impulse Response“ </span>(kurz: <span style='color:OrangeRed'>FIR </span>) Filter.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #3 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Es gelten weiterhin die Werte aus dem <span style='color:Gray'>Beispiel #1 </span>:
$num$=$8z^{-1}$
$den$=$1-3z^{-1}+2z^{-2}$
Wir bestimmen außerdem die Anzahl der Abtastpunkte:
$npoint$ = $20$
Wie zuvor, wird wiederholt die Entfaltung, also Division durchgeführt. Der Code wird in
einer Funktion „ longdiv “ zusammengefasst.
%% Cell type:code id: tags:
``` octave
function y = longdiv(num,den,npoint)
[y(1) R] = deconv(num,den);
for i=2:npoint
num = conv(R,[1 0]);
[f R] = deconv(num,den);
y(i) = f(length(f));
end
end
```
%% Cell type:code id: tags:
``` octave
num = [8 0];
den = [1 -3 2];
npoint = 20;
y = longdiv(num, den, npoint);
plot(y)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Das System ist instabil.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #4 </span>
%% Cell type:code id: tags:
``` octave
pkg load control
Ts = 0.1;
p = -1;
npoint = 20;
y1 = longdiv([1 0],[1 -exp(Ts*p)],npoint);
t =(1:npoint)*Ts-Ts;
```
%% Cell type:code id: tags:
``` octave
yc = impulse(tf(1,[1 1]),t);
plot(t,yc,'r',t,y1,'*')
```
%% Cell type:code id: tags:
``` octave
yh = longdiv(1-exp(Ts*p),[1 -exp(Ts*p)],npoint);
plot(t,yh,'+')
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #5 </span>
%% Cell type:code id: tags:
``` octave
function [G Ph] = BodeZOH(Ts,ommax,dom)
k=0;
kmax = round(ommax/dom);
for k=1:kmax
om(k) = k*dom;
Renum = 1 -cos(-om(k)*Ts);
Imnum = -sin(-om(k)*Ts);
Gnum = sqrt(Renum*Renum+Imnum*Imnum);
G(k) = Gnum/(om(k)*Ts);
Ph(k) = atan2(Imnum,Renum)-pi/2;
end
end
Ts = 0.1;
ommax = 100;
dom = 0.1;
om=[0.1:0.1:100];
[G Ph] = BodeZOH(Ts,ommax,dom);
subplot(2,1,1), semilogx(om,20*log10(G));
subplot(2,1,2), semilogx(om,Ph*90/pi);
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Amplitude ist im z Bereich kleiner als 20dB.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #6 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
- Wir stellen drei Verfahren vor, wie in Matlab die Transformation von kontinuierlichen Systemen zu diskreten Systemen erfolgen kann.
- Dabei wird die Matlabfunktion<code>c2d</code> verwendet
- Gegeben sei das System:
$G(s)=\frac{1}{s+1}$ und $T_s=0,1$
%% Cell type:code id: tags:
``` octave
pkg load control
```
%% Cell type:code id: tags:
``` octave
disp('Abtastzeit')
Ts = 0.1
disp('Kontinuierliche Systeme')
G = tf(1,[1 1])
disp('Beispiel 1: ZOH')
Gh = c2d(G,Ts,'zoh')
disp('Beispiel 2: Impulse')
Gz = c2d(G,Ts,'impulse')
disp('Bespiel 3: Verfahren G(s)/s')
Gi = tf([1],[1 0])
Gs = G*Gi
Gsz = c2d(Gs,Ts,'impulse')
Gh = tf([1 -1], [1 0],Ts)
Gf = Gh*Gsz/Ts;
minreal(Gf)
```
%% Cell type:code id: tags:
``` octave
% Set the Octave Engine to run the simulatio
% Simulation Parameters
pkg load control
addpath("./Octsim");
% Start time
tini = 0;
% End time
tfinal = 5;
% Time Step
dt = 0.001;
% Number of data flows in the schematic
nflows = 5;
% Sampling time for discrete time
Ts = 0.1;
% Instance of the simulation schematic
sc = Schema(tini,tfinal,dt,nflows);
% List of components
c{1} = SinusoidalSignalSource(1,1,2*pi,0);
c{2} = TransferFunction(1,2,1,[1,1]); %G(s)
c{3} = ZOH(1,3,Ts);
c{4} = TransferFunction(3,4,1,[1,1]); %Example 1 - ZOH at the input
c{5} = DTTransferFunction(1,5,[0.09516],[1 -exp(-0.1)],Ts); %Example 2
c{6} = DTTransferFunction(1,6,[0.1 0],[1 -exp(-0.1)],Ts) ;%Example 3
sc.AddListComponents(c);
% Run the schematic and plot
out = sc.Run([ 2 3 4 5 6]);
%plot(out(1,:),out(2,:),out(1,:),out(3,:),out(1,:),out(4,:));
plot(out(1,:),out(2,:),out(1,:),out(4,:),out(1,:),out(5,:), out(1,:),out(6,:));
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #7 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Gegeben ist : $G(z)=\frac{z}{z-e^{-\tau}}$, und $U(z)=\frac{z}{z-1}$ mit $\tau=0,1$
Und nehmen wir an: $Y(z)=G(z)U(z)$
%% Cell type:code id: tags:
``` octave
pkg load control
```
%% Cell type:code id: tags:
``` octave
Ts = 0.1;
Gz = tf([1 0],[1 -exp(-0.1)],Ts);
%Impuls: longdiv(num,den,100)
y1=longdiv(cell2mat(Gz.num),cell2mat(Gz.den),100);
plot(y1);
```
%% Cell type:code id: tags:
``` octave
Uz = tf([1 0],[1 -1],Ts);
Yz = Gz*Uz;
y2=longdiv(cell2mat(Yz.num),cell2mat(Yz.den),100);
plot(y2);
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Übertragungsfunktion lautet stattdessen:
<br><br>
$G(z)$=$\frac{z(z-e^{-\tau})-(z-1)}{(z-1)(z-e^{-\tau})}$
%% Cell type:code id: tags:
``` octave
Gz2 = tf([1-exp(-0.1) 0],[1 -1-exp(-0.1) exp(-0.1)],Ts);
y3=longdiv(cell2mat(Gz2.num),cell2mat(Gz2.den),100);
plot(y3);
```
%% Cell type:code id: tags:
``` octave
```
notebooks/alt/V04_regelung_diskreter_systeme.ipynb deleted 100644 → 0
View file @ 68879533
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V4 - Regelalgorithmen für die digitale Regelung</span>
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b> Reglerauslegung im Frequenzbereich - Zeitkontinuierlicher PID </b>
Gegeben ist die Regelstrecke:
<br><br>
$G_S(s)= \frac{p_1p_2}{s^{2}-(p_1+p_2)s+p_1p_2}$
<b> Codes: </b>
<br><br>
In Matlab definieren wir:
%% Cell type:code id: tags:
``` octave
pkg load control
```
%% Cell type:code id: tags:
``` octave
p1 = -5;
p2 = -8;
Num = [p1*p2]
Den = [1 -(p1+p2) p1*p2]
Gs = tf(Num,Den)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Wir wählen für die Bandbreite $\omega_0=20$ rad /s.
Die Parameter des <span style='color:goldenrod'> PID Reglers </span> lassen sich
wie folgt berechnen:
%% Cell type:code id: tags:
``` octave
% PID
omo = 20
Kd = omo/(p1*p2)
Kp = Kd*(-(p1+p2))
Ki = Kd*p1*p2
Gr = tf([Kd Kp Ki],[1 0]);
```
%% Cell type:code id: tags:
``` octave
bode(Gs)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Das Bodediagramm beschreibt ein Verzögerungselement zweiter Ordnung (PT2). Der Phasenwinkel
verläuft von 0 bis -180 Grad. Die Regelstrecke ist stabil.
%% Cell type:code id: tags:
``` octave
Go = Gs*Gr
bode(Go)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Der Phasenwinkel bleibt konstant und beträgt -90° (der offene Regelkreis weist integrierendes Verhalten auf).
Die Phasenreserve ist 90°. Der geschlossene Regekreis hat somit ein stabiles und sehr robustes Verhalten.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #2 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b>Reglerauslegung im Frequenzbereich Zeitkontinuierlicher PI</b>
<br><br>
In diesem Beispiel wird anstelle eines PID Reglers, ein PI Regler verwendet. Es wird
weiterhin eine Regelstrecke mit zwei Polstellen betrachtet. Wir stellen ein mögliches
Verfahren vor, einen PI Regler auszulegen.
<br><br>
Ausgangspunkt: Wir geben die Bandbreite $\omega_0$ und die gewünschte Phasenreserve $\varphi_m$
vor.
<br><br>
Gesucht wird der Regler in der Form:
<br>
$R(s) = K_P + \frac{K_I}{s}$
<br><br>
Verfahren: Bei der Bandbreite $\omega_0$ soll die Phasenreserve $\varphi_m$ sein. Die Bandbreite ist als
Durchtrittsfrequenz zu verstehen. Die Durchtrittsfrequenz ist als die Frequenz definiert,
bei der für den Betrag des offenen Regelkreises gilt (mehr dazu in Systemtheorie I):
$|G_o(j\omega_0)|$ = 1
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Für den offenen Regelkreis gilt:
<br><br>
$|R(j\omega_0)(G(j\omega_0)| = 1e^{j(-\pi +\varphi_m)}$
<br><br>
Einsetzen des Frequenzgangs für den PI Regler:
<br><br>
$(K_P + \frac{K_I}{j\omega_0})|G(j\omega_0)|e^{j\varphi_G} = e^{j(-\pi +\varphi_m)}$
<br><br>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Für den PI Regler gilt dann:
<br><br>
$(K_P + \frac{K_I}{j\omega_0})$ = $\frac{e^{j(-\pi +\varphi_m - \varphi_G)}}{|G(j\omega_0)|}$
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Aus dem Koeffizientenvergleich folgen die Parameter und :
<br><br>
$K_P$=$\frac{cos(-\phi + \varphi_m -\varphi_G)}{|G(j\omega_0)|}$
<br><br>
$K_I$=$\frac{-\omega_0sin(-\phi + \varphi_m -\varphi_G)\omega_0}{|G(j\omega_0)|}$
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
In Matlab definieren wir zunächst die Regelstrecke (Definition bleibt gleich)
%% Cell type:code id: tags:
``` octave
p1 = -5
p2 = -8;
Num = [p1*p2]
Den = [1 -(p1+p2) p1*p2]
Gs = tf(Num,Den)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Wir wählen für die Bandbreite $\omega_0=20$ rad /s . Die Parameter des <span style='color:goldenrod'> PI Reglers </span> lassen sich
wie folgt bestimmen:
%% Cell type:code id: tags:
``` octave
% PID
omo = 10
Fi = 40*pi/180
[Go Fo] = bode(Gs,omo)
Fo = Fo*pi/180;
Kd = 0;
Kp = cos(-pi+Fi-Fo)/Go
Ki = -sin(-pi+Fi-Fo)*omo/Go
Gr = tf([Kd Kp Ki],[1 0]);
Go = Gs*Gr
```
%% Cell type:code id: tags:
``` octave
bode(Gs)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Das Bodediagramm beschreibt ein Verzögerungselement zweiter Ordnung. Der Phasenwinkel verläuft
von 0 bis -180 Grad. Die Regelstrecke ist stabil.
%% Cell type:code id: tags:
``` octave
bode(Go)
```
%% Cell type:code id: tags:
``` octave
nyquist(Go)
```
%% Cell type:code id: tags:
``` octave
margin(Go)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Es liegt eine Phasenreserve von 40 Grad bei der Frequenz $\omega_0=20$ rad/s vor.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #3 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b>Reglerauslegung im Frequenzbereich - Zeitdiskreter PI</b>
<br><br>
Es stellt sich nun die Frage, ob wir das Verfahren aus dem zeitkontinuierlichen Fall auch
im zeitdiskreten Bereich nutzen können?
<br><br>
Betrachtet wird wiederum die Regelstrecke:
<br><br>
$G_S(s)= \frac{p_1p_2}{s^{2}-(p_1+p_2)s+p_1p_2}$
<br><br>
<b> Codes: </b>
<br><br>
In Matlab definieren wir:
%% Cell type:code id: tags:
``` octave
p1 = -5;
p2 = -8;
Num = [p1*p2]
Den = [1 -(p1+p2) p1*p2]
Gs = tf(Num,Den)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Regelstrecke muss zunächst diskretisiert werden. Dafür verwenden wir in Matlab
folgende Befehle:
%% Cell type:code id: tags:
``` octave
Ts = 0.1;
Gsz = c2d(Gs,Ts,'zoh');
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Für den zeitdiskreten Regler verwenden wir zunächst die gleichen Parameter, die wir für
den zeitkontinuierlichen PI Regler bestimmen haben.
<br><br>
Für den zeitkontinuierlichen Fall wurden folgende Parameter bestimmt:
%% Cell type:code id: tags:
``` octave
% PID
omo = 10
Fi = 46*pi/180
[Go Fo] = bode(Gs,omo)
Fo = Fo*pi/180;
Kd = 0;
Kp = cos(-pi+Fi-Fo)/Go
Ki = -sin(-pi+Fi-Fo)*omo/Go
Gr = tf([Kd Kp Ki],[1 0]);
Go = Gs*Gr
```
%% Cell type:code id: tags:
``` octave
bode(Go)
```
%% Cell type:code id: tags:
``` octave
nyquist(Go)
```
%% Cell type:code id: tags:
``` octave
margin(Go)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Für den<span style='color:goldenrod'> zeitdiskreten Regler </span> verwenden wir zunächst die gleichen Parameter, die wir für
den <span style='color:goldenrod'> zeitkontinuierlichen PI Regler </span> bestimmen haben.
%% Cell type:code id: tags:
``` octave
Kr = Kp
Ti = Kr/Ki
Zn = 1/(1+Ts/Ti)
Numrz = [Kr/Zn -Kr]
Denrz = [1 -1]
Grz= tf(Numrz,Denrz,Ts)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Nun wird der zeitdiskrete offene Regelkreis im Frequenzbereich analysiert.
%% Cell type:code id: tags:
``` octave
Goz = Gsz*Grz
margin(Goz)
[Gzo Fzo] = bode(Goz,omo)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Wenn wir die Parameter des zeitkontinuierlichen Reglers übernehmen, beträgt die
Phasenreserve des zeitdiskreten Regelkreises 18 Grad (180°-162°=18°) anstatt der
gewünschten 40 Grad.
<br><br>
Durch die Diskretisierung wird die Phasenreserve kleiner als im zeitkontinuierlichen
Bereich
<br><br>
Unter bestimmten Bedingungen kann der Phasenverlust näherungsweise ausgerechnet
werden.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Der Phasenverlust beträgt im zeitdiskreten Fall:
<br><br>
$\Delta\varphi = \frac{\omega_0T_s}{2}$
<br><br>
In unserem Beispiel wäre der Phasenverlust:
<br><br>
$\Delta\varphi =\frac{10rad/s 0.1s}{2} \frac{180}{\pi}= 28,6479$
<br><br>
Damit der Phasenverlust bei derselben Bandbreite nicht zu groß wird , müssen wir die
Abtastzeit verkleinern . Wenn wir als Abtastzeit$T_2 = 0.02s$ wählen:
<br><br>
$\Delta\varphi =\frac{10rad/s 0.1s}{2} \frac{180}{\pi} = 5,7296$
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<span style='color:orange'> Vorgehen </span>: Wir möchten, dass der zeitdiskrete Regelkreis ebenfalls eine Phasenreserve
von 40 Grad aufweist.
<br><br>
Wir berechnen $K_P$ und $K_I$ für den zeitkontinuierlichen Fall. Anstatt$K_P$ und $K_I$ für eine
Phasenreserve von 40 Grad zu bestimmen, berechnen wir diese Parameter nun für eine
Phasenreserve von $40 + \Delta\phi$ .
<br><br>
Durch die Diskretisierung entsteht ein Phasenverlust von genau $\Delta\phi$ .
<br><br>
Schlussendlich erhält man einen zeitdiskreten Regler, der sicherstellt, dass der zeitdiskrete
Regelkreis eine Phasenreserve von 40 Grad aufweist (wie gewünscht)!
<br><br>
<span style='color:cyan'> Vorteil </span> dieses Verfahren: Bei ausreichend kleiner Abtastzeit kann der Reglerentwurf wie
im zeitkontinuierlichen Fall durchgeführt werden und der resultierende Regler als
zeitdiskreter Regler verwendet werden.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
In Matlab wählen wir zunächst anstatt der Abtastzeit 0.1s die Abtastzeit 0.02s aus. Damit
wird sichergestellt, dass der Phasenverlust zwischen dem zeitkontinuierlichen und dem
zeitdiskreten Fall nicht zu groß wird (der Phasenverlust wird dadurch von 28°für 0.1s auf
6° für 0.02s reduziert).
<b> Ihre Aufgabe: </b>
<br><br>
Wählen Sie die Phasenreserve größer !
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #4 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Gegeben ist die Regelstrecke
<br><br>
$G_S(s)= \frac{p_1p_2}{s^{2}-(p_1+p_2)s+p_1p_2}$
<br><br>
Dieses System 2. Ordnung wird in Matlab wie folgt definiert:
%% Cell type:code id: tags:
``` octave
p1 = -5
p2 = -8;
Num = [p1*p2]
Den = [1 -(p1+p2) p1*p2]
Gs = tf(Num,Den)
bode(Gs)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Regelstrecke wird diskretisiert weil nur für zeitdiskrete Systeme ein Deadbeat Regler
ausgelegt werden kann:
%% Cell type:code id: tags:
``` octave
Ts = 0.02
Gsz = c2d(Gs,Ts,'zoh')
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Wir definieren das gewünschte Verhalten des geschlossenen Regelkreises durch die
Eingabe des gewünschten Übertragungsverhaltens <code>Kw</code>
%% Cell type:code id: tags:
``` octave
Kw = tf(1,[1 0 0],Ts);
Grz = Kw/(Gsz*(1-Kw))
Gl = Gsz*Grz;h
Numz = cell2mat(Grz.num)
Denz = cell2mat(Grz.den)
```
%% Cell type:code id: tags:
``` octave
p1 = -5;
p2 = -8;
Num = [p1*p2];
Den = [1 -(p1+p2) p1*p2];
```
%% Cell type:code id: tags:
``` octave
Grz = tf([1 -1.757 0.7711 0 0], [0.00734 0.006731 -0.00734 0.006731 0 0], 0.02)
Numrz = cell2mat(Grz.num)
Denrz = cell2mat(Grz.den)
```
%% Cell type:code id: tags:
``` octave
% Set the Octave Engine to run the simulatio
% Simulation Parameters
pkg load control
addpath("./Octsim");
% Start time
tini = 0;
% End time
tfinal = 3;
% Time Step
dt = 0.001;
% Number of data flows in the schematic
nflows = 4;
% Sampling time for discrete time
Ts = 0.02;
% Instance of the simulation schematic
sc = Schema(tini,tfinal,dt,nflows);
% List of components
c{1} = StepSource(1,0,1,0.1);
c{2} = Sum(1,4,2,1,-1);
c{3} = DTTransferFunction(2,3,Numrz,Denrz,Ts);
c{4} = TransferFunction(3,4,Num,Den);
sc.AddListComponents(c);
% Run the schematic and plot
out = sc.Run([2 3 4]);
plot(out(1,:),out(3,:));
```
%% Cell type:code id: tags:
``` octave
plot(out(1,:),out(4,:));
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #5 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Für die gleiche Regelstrecke werden wir nun einen Regler mit
endlicher, aber nicht minimaler, Einstellzeit ( Deadbeat ) einsetzen.
<br><br>
Die Übertragungsfunktion dieses Reglers ergibt sich aus der Anforderung an das
Wunschübertragungsverhalten.
<br><br>
Das Sollübertragunsverhalten wird durch die Übertragungsfunktion <code>Kw</code> für den
geschlossenen Regelkreis definiert, die nun eine höhere Ordnung aufweist .
%% Cell type:code id: tags:
``` octave
Kw = tf([0.2 0.4 0.6 0.8 1],[3 0 0 0 0 0],Ts)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Der Regler wird wie im Fall vom Deadbeat Regler wie folgt definiert:
%% Cell type:code id: tags:
``` octave
Grz = Kw/(Gsz*(1-Kw))
Gl = Gsz*Grz;
Numz = cell2mat(Grz.num);
Denz = cell2mat(Grz.den);
```
%% Cell type:code id: tags:
``` octave
p1 = -5;
p2 = -8;
Num = [p1*p2]
Den = [1 -(p1+p2) p1*p2]
Gs = tf(Num,Den)
bode(Gs)
Ts = 0.02
Gsz = c2d(Gs,Ts,'zoh')
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Regelgröße erreicht die Führungsgröße
nicht in zwei Abtastschritten (minimaler
Anzahl an Abtastschritten), sondern es
werden mehr Schritte benötigt.
<br><br>
<span style='color:cyan'> Vorteil </span> der sich aus so einer Reglerauslegung
ergibt
≡ der Stellgrößenverlauf hat nicht so höhe Sprünge
wie in Fall einer Deadbeat Reglers
notebooks/alt/V06_normalformen_in_zustandsraumdarstellung.ipynb deleted 100644 → 0
View file @ 68879533
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V6 - Normalformen in Zustandsraumdarstellung</span>
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Wir arbeiten zunächst mit der folgenden Zustandsraumdarstellung:
<br> $A$ = $\left[ \begin{array}{rrrr}
-1 & 0 \\
1 & -2 \\
\end{array}\right] $
<br><br> $B$ = $\left[ \begin{array}{rrrr}
1 \\
0 \\
\end{array}\right] $
<br><br> $C$ = $\left[ \begin{array}{rrrr}
0 & 1 \\
\end{array}\right] $
<br><br> $D$ = $0$ .
%% Cell type:code id: tags:
``` octave
A = [-1 0;
1 -2]
B = [1; 0]
C = [0 1]
D = 0
```
%% Output
A =
-1 0
1 -2
B =
1
0
C =
0 1
D = 0
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Tranformationsmatrix <code>V</code> ist gegeben:
%% Cell type:code id: tags:
``` octave
V = [1 0;
1 1]
```
%% Output
V =
1 0
1 1
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Damit sieht die inverse Transformationsmatrix <code>V1</code> wie folgt aus:
%% Cell type:code id: tags:
``` octave
V1= [1 0;
-1 1]
```
%% Output
V1 =
1 0
-1 1
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Jetzt transfomieren wir die Matrizen auf Diagonalform:
<br> $A_{neu}$=$V^{-1}AV$ = $\left[ \begin{array}{rrrr}
-1 & 0 \\
0 & -2 \\
\end{array}\right] $
%% Cell type:code id: tags:
``` octave
A_neu = [-1 0;
0 -2]
```
%% Output
A_neu =
-1 0
0 -2
%% Cell type:markdown id: tags:
<br><br> $B_{neu}$=$V^{-1}B$ = $\left[ \begin{array}{rrrr}
1 \\
-1 \\
\end{array}\right] $
%% Cell type:code id: tags:
``` octave
B_neu =[1;
-1]
```
%% Output
B_neu =
1
-1
%% Cell type:markdown id: tags:
<br><br> $C_{neu}$=$CV$ = $\left[ \begin{array}{rrrr}
1 & 1 \\
\end{array}\right] $
%% Cell type:code id: tags:
``` octave
C_neu=[1 1]
```
%% Output
C_neu =
1 1
%% Cell type:code id: tags:
``` octave
Adc = V1*A*V
Bd = V1*B
Cd = C*V
```
%% Output
Adc =
-1 0
0 -2
Bd =
1
-1
Cd =
1 1
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<b>Frage</b>: Überprüfen Sie, ob das transformierte System eine diagonale Normalform hat?
%% Cell type:code id: tags:
``` octave
if Adc == A_neu
display("Transformationsmatrix erfolgreich verifiziert.")
end
```
%% Output
Transformationsmatrix erfolgreich verifiziert.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
In den folgenden Codeabschnitten bilden wir das nachfolgende Simulinkmodell nach.
Gegeben sind zwei Zustandsraumdarstellungen samt Ausganggrößen:
%% Cell type:markdown id: tags:
![rsz_1v6_normalformen_in_zustandsraumdarstellung-42.png](attachment:rsz_1v6_normalformen_in_zustandsraumdarstellung-42.png)
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Zunächst untersuchen wir das Ausgangsignal von der „originalen“
Zustandsraumdarstellung :
%% Cell type:code id: tags:
``` octave
% Set the Octave Engine to run the simulation
% Simulation Parameters
pkg load control
addpath("./Octsim");
% Start time
tini = 0;
% End time
tfinal = 7;
% Time Step
dt =0.01;
% Number of data flows in the schematic
nflows = 4;
xo = 0;
% Instance of the simulation schematic
sc = Schema(tini,tfinal,dt,nflows);
% List of components
c{1} = StepSource(1,0,1,1);
c{2} = StateSpace(1,2,A,B,C,D,xo);
c{3} = StateSpace(1,3,Adc,Bd,[0 1],0,xo);
c{4} = StateSpace(1,4,Adc,Bd,[1 0],0,xo);
c{5} = StateSpace(1,5,Adc,Bd,Cd,0,xo);
sc.AddListComponents(c);
% Run the schematic and plot
out = sc.Run([2 3 4 5]);
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Unter „Ausgang 1“ stellen wir das Ausgangsignal von der „originalen“
Zustandsraumdarstellung dar:
%% Cell type:code id: tags:
``` octave
%Ausgang 1
plot(out(1,:),out(2,:));
```
%% Output
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Unter „Ausgang 2“ können wir das Ausgangsignal von der „modifizierten“
Zustandsraumdarstellung in Normalform untersuchen:
%% Cell type:code id: tags:
``` octave
% Ausgang 2
plot(out(1,:),out(5,:));
```
%% Output
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Scope 2:
%% Cell type:code id: tags:
``` octave
%Scope 2:
plot(out(1,:),out(3,:),out(1,:),out(4,:));
```
%% Output
%% Cell type:code id: tags:
``` octave
```
notebooks/alt/V08_regelung_im_zustandsraum.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
notebooks/alt/V09_zustandsrekonstruktion_mittels_beobachter.ipynb deleted 100644 → 0
View file @ 68879533
Source diff could not be displayed: it is too large. Options to address this: view the blob.
notebooks/alt/V10_kalman_filter.ipynb deleted 100644 → 0
View file @ 68879533
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V10 - Kalman Filter</span>
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
SFUNTMPL General M-file S-function template:
<br><br> With M-file S-functions, you can define you own ordinary differential
equations (ODEs), discrete system equations, and/or just about
any type of algorithm to be used within a Simulink block diagram.
The general form of an M-File S-function syntax is:
<code> [SYS,X0,STR,TS] = SFUNC(T,X,U,FLAG,P1,...,Pn) </code>.
What is returned by <code>SFUNC</code> at a given point in time, <code>T</code>, depends on the
value of the <code>FLAG</code>, the current state vector,<code> X</code>, and the current
input vector, <code>U</code>.
%% Cell type:markdown id: tags:
| Flag |Result | Description|
| :- | -: | :-: |
| 0 | [SIZES,X0,STR,TS] | Initialization, return system sizes in SYS, initial state in X0, state ordering strings in STR, and sample times in TS.|
| 1 | DX | Return continuous state derivatives in SYS.|
| 2 | DS | Update discrete states SYS = X(n+1)
| 3 | Y | Return outputs in SYS.
| 4 | TNEXT | Return next time hit for variable step sample time in SYS.
| 5 | | Reserved for future (root finding).|
| 9 | [] | Termination, perform any cleanup SYS=[].
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Optional parameters, P1,...,Pn can be provided to the S-function and used during any <code>FLAG</code> operation.
When <code>SFUNC</code> is called with <code>FLAG</code> = 0, the following information
should be returned:
<code>SYS(1)</code> = Number of continuous states.
<br> <code>SYS(2)</code> = Number of discrete states.
<br><code>SYS(3)</code> = Number of outputs.
<br><code>SYS(4)</code> = Number of inputs.
Any of the first four elements in SYS can be specified
as -1 indicating that they are dynamically sized. The
actual length for all other flags will be equal to the
length of the input, <code>U</code>.
<br> <code>SYS(5)</code> = Reserved for root finding. Must be zero.
<br><code>SYS(6)</code> = Direct feedthrough flag (1=yes, 0=no). The s-function
has direct feedthrough if <code>U</code> is used during the <code>FLAG</code>=3
call. Setting this to 0 is akin to making a promise that
<code>U</code> will not be used during <code>FLAG</code>=3. If you break the promise
then unpredictable results will occur.
<br><code>SYS(7)</code> = Number of sample times. This is the number of rows in <code>TS</code>.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
The state vectors, <code>X</code> and <code>X0</code> consists of continuous states followed
by discrete states.
<br> <code>X0</code> = Initial state conditions or [] if no states.
<br><code>STR</code> = State ordering strings which is generally specified as [].
<br><code>TS</code> = An m-by-2 matrix containing the sample time. (period, offset) information. Where m = number of sampletimes. The ordering of the sample times must be
<br> <code>TS</code> =<br> [0 0, : Continuous sample time.
<br>0 1, : Continuous, but fixed in minor step sample time.
<br> PERIOD OFFSET, : Discrete sample time where
<br> PERIOD greater than 0 & OFFSET smaller than PERIOD.
<br> -2 0 ]; : Variable step discrete sample time
<br>where <code>FLAG</code>=4 is used to get time of next hit.
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
There can be more than one sample time providing
they are ordered such that they are monotonically
increasing. Only the needed sample times should be
specified in <code>TS</code>. When specifying than one
sample time, you must check for sample hits explicitly by
seeing if
<code>abs(round((T-OFFSET)/PERIOD) - (T-OFFSET)/PERIOD)</code>
is within a specified tolerance, generally 1e-8. This
tolerance is dependent upon your model's sampling times
and simulation time.
You can also specify that the sample time of the S-function
is inherited from the driving block. For functions which
change during minor steps, this is done by
specifying <code>SYS(7)</code> = 1 and <code>TS</code> = [-1 0]. For functions which
are held during minor steps, this is done by specifying
<code>SYS(7)</code> = 1 and <code>TS</code> = [-1 1].
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts] = kalman1(t,x,u,flag,sigma,Ts,xo,ro)
switch flag,
% The following outlines the general structure of an S-function.
% Initialization:
case 0,
[sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro);
% Derivatives:
case 1,
sys=mdlDerivatives(t,x,u);
% Update:
case 2,
sys=mdlUpdate(t,x,u,sigma);
% Outputs:
case 3,
sys=mdlOutputs(t,x,u);
% GetTimeOfNextVarHit:
case 4,
sys=mdlGetTimeOfNextVarHit(t,x,u,Ts);
% Terminate:
case 9,
sys=mdlTerminate(t,x,u);
% Unexpected flags:
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlInitializeSizes</code>
returns the sizes, initial conditions, and sample times for the S-function.
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro)
% Call simsizes for a sizes structure, fill it in and convert it to a sizes array.
% Note that in this example, the values are hard coded. This is not a recommended practice
% as the characteristics of the block are typically defined by the S-function parameters.
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 2;
sizes.NumOutputs = 1;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
% initialize the initial conditions:
x0 = [xo ro];
% str is always an empty matrix:
str = [];
% initialize the array of sample times:
ts = [Ts 0];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlDerivatives</code>
returns the derivatives for the continuous states.
%% Cell type:code id: tags:
``` python
function sys=mdlDerivatives(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlUpdate</code>
handles discrete state updates, sample time hits, and major time step
requirements.
%% Cell type:code id: tags:
``` python
function sys=mdlUpdate(t,x,u,sigma)
sys(2) = sigma*sigma*x(2)/(sigma*sigma+x(2));
K = sys(2)/(sigma*sigma);
sys(1) = x(1)+K*(u(1)-x(1));
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlOutputs</code>
returns the block outputs.
%% Cell type:code id: tags:
``` python
function sys=mdlOutputs(t,x,u)
sys(1) = x(1);
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code> mdlGetTimeOfNextVarHit</code>
returns the time of the next hit for this block. Note that the result is
absolute time. Note that this function is only used when you specify a
variable discrete-time sample time [-2 0] in the sample time array in
<code>mdlInitializeSizes</code>.
%% Cell type:code id: tags:
``` python
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code> mdlTerminate</code>
performs any end of simulation tasks.
%% Cell type:code id: tags:
``` python
function sys=mdlTerminate(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #2 </span>
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts] = kalman2(t,x,u,flag,sigma,Ts,xo,ro)
switch flag,
% The following outlines the general structure of an S-function.
% Initialization:
case 0,
[sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro);
% Derivatives:
case 1,
sys=mdlDerivatives(t,x,u);
% Update:
case 2,
sys=mdlUpdate(t,x,u,sigma);
% Outputs:
case 3,
sys=mdlOutputs(t,x,u);
% GetTimeOfNextVarHit
case 4,
sys=mdlGetTimeOfNextVarHit(t,x,u,Ts);
% Terminate:
case 9,
sys=mdlTerminate(t,x,u);
% Unexpected flags:
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlInitializeSizes</code>
returns the sizes, initial conditions, and sample times for the S-function.
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro)
% Call simsizes for a sizes structure, fill it in and convert it to a sizes array.
% Note that in this example, the values are hard coded. This is not a recommended practice
% as the characteristics of the block are typically defined by the S-function parameters.
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 2;
sizes.NumOutputs = 1;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
% initialize the initial conditions:
x0 = [xo ro];
% str is always an empty matrix:
str = [];
% initialize the array of sample times:
ts = [Ts 0];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlDerivatives</code>
returns the derivatives for the continuous states.
%% Cell type:code id: tags:
``` python
function sys=mdlDerivatives(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlUpdate</code>
handles discrete state updates, sample time hits, and major time step
requirements.
%% Cell type:code id: tags:
``` python
function sys=mdlUpdate(t,x,u,sigma)
sys(2) = sigma*sigma*x(2)/(sigma*sigma+x(2));
K = sys(2)/(sigma*sigma);
sys(1) = x(1)+K*(u(1)-x(1));
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlOutputs</code>
Return the block outputs.
%% Cell type:code id: tags:
``` python
function sys=mdlOutputs(t,x,u)
sys(1) = x(1);
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code> mdlGetTimeOfNextVarHit</code>
returns the time of the next hit for this block. Note that the result is absolute time. Note that this function is only used when you specify a variable discrete-time sample time [-2 0] in the sample time array in <code>mdlInitializeSizes</code>.
%% Cell type:code id: tags:
``` python
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlTerminate</code>
performs any end of simulation tasks.
%% Cell type:code id: tags:
``` python
function sys=mdlTerminate(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #3 </span>
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts] = kalman3(t,x,u,flag,Ts,xo,ro,r1,r2)
switch flag,
% The following outlines the general structure of an S-function.
% Initialization:
case 0,
[sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro);
% Derivatives:
case 1,
sys=mdlDerivatives(t,x,u);
% Update:
case 2,
sys=mdlUpdate(t,x,u,r1,r2);
% Outputs:
case 3,
sys=mdlOutputs(t,x,u);
% GetTimeOfNextVarHit
case 4,
sys=mdlGetTimeOfNextVarHit(t,x,u,Ts);
% Terminate:
case 9,
sys=mdlTerminate(t,x,u);
% Unexpected flags:
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlInitializeSizes</code>
returns the sizes, initial conditions, and sample times for the S-function.
%% Cell type:code id: tags:
``` python
sizes = simsizes;
```
%% Cell type:code id: tags:
``` python
function [sys,x0,str,ts]=mdlInitializeSizes(Ts,xo,ro)
% Call simsizes for a sizes structure, fill it in and convert it to a sizes array.
% Note that in this example, the values are hard coded. This is not a recommended
% practice as the characteristics of the block are typically defined by the S-function parameters.
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 2;
sizes.NumOutputs = 1;
sizes.NumInputs = 2;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
% initialize the initial conditions:
x0 = [xo ro];
% str is always an empty matrix:
str = [];
% initialize the array of sample times:
ts = [Ts 0];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlDerivatives</code>
returns the derivatives for the continuous states.
%% Cell type:code id: tags:
``` python
function sys=mdlDerivatives(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code> mdlUpdate</code>
handles discrete state updates, sample time hits, and major time step
requirements.
%% Cell type:code id: tags:
``` python
function sys=mdlUpdate(t,x,u,r1,r2)
sys(2) = r2*(r1+0.25*x(2))/(r2+r1+0.25*x(2));
K = sys(2)/(r2);
sys(1) = 0.5*x(1)+u(2)+K*(u(1)-(0.25*x(1)+u(2)));
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlOutputs</code>
returns the block outputs.
%% Cell type:code id: tags:
``` python
function sys=mdlOutputs(t,x,u)
sys(1) = x(1);
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code> mdlGetTimeOfNextVarHit</code>
returns the time of the next hit for this block. Note that the result is
absolute time. Note that this function is only used when you specify a
variable discrete-time sample time [-2 0] in the sample time array in
<code>mdlInitializeSizes</code>.
%% Cell type:code id: tags:
``` python
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sys = [];
end
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<code>mdlTerminate</code>
performs any end of simulation tasks.
%% Cell type:code id: tags:
``` python
function sys=mdlTerminate(t,x,u)
sys = [];
end
```
requirements.txt
View file @ 9415dc13
scilab_kernel
numpy
scipy
control
matplotlib
PyGnuplot
scilab-kernel==0.9.10
numpy==1.19.5
scipy==1.6.1
control==0.9.2
matplotlib==3.6.0
PyGnuplot==0.12.3