In einer Reihe von Blogposts werden wir uns mit linearen zeitinvarianten zeitkontinuierlichen Systemen beschäftigen. Dabei können wir diese wiederum in SISO (Single Input Single Output) oder MIMO (Multi Input Multi Output) Systemen unterscheiden. Die klassische Regelungstechnik, also zwischen 1900 und 1950, beschäftigte sich fast ausschließlich mit SISO Systemen im Frequenzbereich. Zwar hat sich seitdem viel getan, aber noch heute kann man die Grundlagen der Regelungstechnik in diesem Setting am besten verstehen, und dominiert eigentlich noch immer die Prozessautomatisierung.
Die drei gleichwertigen Darstellungen
Obwohl wir uns im Weiteren auf SISO Systeme beschränken, gelten die Ausführungen im Prinzip auch für MIMO Systeme. Für MIMO System treten aber neue Objekte auf, wie Übertragungsmatrizen, Gewichtungsmatrizen, …
Diese Zusammenhänge wollen wir nun skizzenhaft untersuchen ohne alle mathematische Aspekte genau zu beschreiben.
Impulsantwort
Wir wollen mit \(g(t)\) jene Funktion bezeichnen, die entsteht, wenn wir am Eingang einen Impuls (Delta-Funktion) \(u(t)=\delta(t)\) aufschalten. Eine Delta-Funktion ist weder physikalisch realisierbar, noch ist sie eigentlich eine Funktion, sondern eine Distribution. Uns soll das hier aber nicht weiter stören. Es gilt also,
\[ g(t) = y^{\delta}(t) \qquad \text{erzeugt durch} \qquad u(t) = \delta(t). \]
Die Funktion \(g(t)\) wird nun als Impulsantwort (auch Gewichtungsfunktion) bezeichnet. Für einen allgemeinen Eingang \(u(t)\) können wir mithilfe der Impulsantwort \(g(t)\) die allgemeine Systemantwort \(y(t)\) durch das Anwenden des Faltungsintegrals berechnen
\[ \begin{split} y(t) = \int_{0}^{t} g(t-\tau)u(t) \end{split}. \]
Übertragungsfunktion
Die Übertragungsfunktion
\[ G(s) = \mathcal{L} \left\{ g(t) \right\} = \int_{0^{-}}^{\infty} g(t) e^{-st} dt. \]
entsteht aus der (einseitigen) Laplace Transformation der Impulsantwort. Die Laplace Transformation erweist sich besonders nützlich, da aus Differenzialgleichungen in Zeitbereich, algebraische Gleichungen im Frequenzbereich entstehen, und Faltungsintegrale in einfache Produkte übergehen. So gilt im Laplace Bereich,
\[ Y(s) = G(s)\cdot U(s). \]
In Zukunft für uns besonders nützlich wird sich die Laplace Transformation der zeitlichen Ableitung einer Funktion erweisen:
\[ \begin{split} \mathcal{L} \left\{ \frac{d}{dt} f(t) \right\} &= \int_{0^{-}}^{\infty} \frac{d}{dt}f(t)e^{-st}dt \\ &=\left[ f(t)e^{-st} \right]_{t=0^{-}}^{t=\infty} - \int_{0^{-}}^{\infty}f(t)\left(-se^{-st}\right)dt \\ &= f(0^{-}) + s\mathcal{L}\left\{ f(t) \right\} \\ &= f(0^{-}) + sF(s) \end{split} \]
Auch für höhere Ableitungen kann die Laplace Transformation durchgeführt werden und wir erhalten unter der Annahme, dass alle Anfangsbedingungen Null sind (\(f(0^{-}) = 0, \dot{f}(0^{-}) = 0, \cdots\)), die Beziehungen
| Zeitbereich | Frequenzbereich |
|---|---|
| \(f(t)\) | \(F(s)\) |
| \(\frac{d}{dt}f(t)\) | \(sF(s)\) |
| \(\frac{d^2}{dt^2}f(t)\) | \(s^2F(s)\) |
| \(\vdots\) | \(\vdots\) |
Oft können physikalische Prozesse mit Differenzialgleichungen beschrieben werden, und liegen in der Form
\[ \begin{split} & y^{(n)}(t)a_n + y^{(n-1)}(t)a_{n-1} + \cdots + \dot{y}(t)a_1 + y(t)a_{0} \\ & \quad = u^{(n)}(t)b_n + u^{(n-1)}(t)b_{n-1} + \cdots + \dot{u}(t)b_1 + u(t)b_{0} \end{split} \]
in Textbüchern vor. Durch das Anwenden der obigen Beziehungen können wir die Differenzialgleichung in die algebraische Gleichung
\[ \begin{split} & s^nY(s)a_n + s^{n-1}Y(s)a_{n-1} + \cdots + sY(s)a_1 + Y(s)a_{0} \\ & \quad = s^nU(s)b_n + s^{n-1}U(s)b_{n-1} + \cdots + sU(s)b_1 + U(s)b_{0} \end{split} \]
überführen. Das Herausheben von \(Y(s)\) und \(U(s)\) erzeugt
\[ \begin{split} & (s^na_n + s^{n-1}a_{n-1} + \cdots + sa_1 + a_{0})Y(s) \\ & \quad = (s^nb_n + s^{n-1}b_{n-1} + \cdots + sb_1 + b_{0})U(s) \end{split} \]
und durch das Umstellen der Gleichung erhalten wir die Übertragungsfunktion
\[ G(s) = \frac{Y(s)}{U(s)} = \frac{s^nb_n + s^{n-1}b_{n-1} + \cdots + sb_1 + b_{0}}{s^na_n + s^{n-1}a_{n-1} + \cdots + sa_1 + a_{0}}. \]
Zustandsraum
Die Zustandsraumdarstellung kann mit
\[ \begin{split} \dot{x}(t) &= Ax(t) + bu(t) \\ y(t) &= c^Tx(t) + du(t) \end{split} \]
angegeben werden. Diese Darstellung ist seit den späten 1950er bzw. frühen 1960er Jahren gebräuchlich und wurde vor allem durch Rudolf Kalman in die Control-Community eingeführt. Kalman hat in dieser Darstellung viele Beiträge für die Regelungstheorie geliefert und ist oft Namensgeber von Methoden (z.B. LQR, Kalman-Filter, Ho-Kalman Identifikation, Kalman-Zerlegung, Steuerbarkeit und Beobachtbarkeit nach Kalman, …).
Wenn jedoch keine vollständige Zustandsmessung vorliegt, sondern nur das Eingangs-Ausgangsverhalten, dann ist die obige Darstellung keineswegs eindeutig. Eine beliebige invertierbare Transformation
\[ x = T \tilde{x}, \qquad \tilde{x} = T^{-1}x \]
erzeugt ein transformiertes Modell
\[ \begin{split} \dot{\tilde{x}}(t) &= T^{-1}AT\tilde{x}(t) + T^{-1}bu(t) \\ y(t) &= c^TT\tilde{x}(t) + du(t) \end{split} \]
welches jedoch dasselbe Eingangs-Ausgangsverhalten besitzt. Gewisse Transformationen haben sich für die Regelungstechnik als hilfreich erwiesen und den Darstellungen wurden Namen zugewiesen. Wichtige Darstellungen: Regelungsnormalform (1. Standardform), Beobachtbarkeitsnormalform (2 Standardform, Jordan-Darstellung, balancierte Darstellung, …
Als Beispiel wollen wir die Regelungsnormalform darstellen: Durch eine spezielle Matrix \(T\), (welche mithilfe der Erreichbarkeitsmatrix errechnet wird) kann die Regelungsnormalform mit
\[ \begin{split} \dot{x} &= A_cx + b_cu \\ y &= c_c^Tx + d_cu \end{split} \rightarrow \begin{cases} \frac{d}{dt} \underset{\dot{\tilde{x}}}{\underbrace{ \begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ \tilde{x}_{n-1} \\ \tilde{x}_n \end{bmatrix}}} = \underset{A_c = T^{-1}AT}{\underbrace{ \begin{bmatrix} 0 & 1 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \\ -a_0 & -a_1 & \dots & -a_{n-1} \end{bmatrix}}} \underset{\tilde{x}}{\underbrace{ \begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ \tilde{x}_{n-1} \\ \tilde{x}_n \end{bmatrix}}} + \underset{b_c = T^{-1}b}{\underbrace{\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}}}u, \\ y = \underset{c_c^T=c^TT}{\underbrace{\begin{bmatrix} \tilde{c}_0 & \tilde{c}_1 & \dots & \tilde{c}_{n-1} \end{bmatrix}}} \underset{\tilde{x}}{\underbrace{ \begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ \tilde{x}_{n-1} \\ \tilde{x}_n \end{bmatrix}}} + \underset{d_c}{\underbrace{d}} u \end{cases} \]
angegeben werden. Für eine detaillierte Herleitung sei auf die entsprechende Literatur verwiesen.
Verbindung zwischen den Darstellungen
Vom Zustandsraum zur Übertragungsfunktion
Mit der Beziehung \(\mathcal{L}(\dot{x}(t)) = sX(s)\) kann die Zustandsraumdarstellung
\[ \begin{split} \dot{x}(t) &= Ax(t) + bu(t) \\ y(t) &= c^Tx(t) + du(t) \end{split} \]
in den Laplacebereich
\[ \begin{split} sX(s) &= AX(s) + bU(s) \\ Y(s) &= c^TX(s) + dU(s) \end{split} \]
überführt werden. Wird die Zustandsgleichung nach dem Zustand
\[ (sI-A)X(s) = bU(s) \quad \rightarrow X(s) = (sI-A)^{-1}bU(s) \]
aufgelöst, und in die Ausgangsgleichung eingesetzt, erhalten wir die Übertragungsfunktion
\[ G(s) = \frac{Y(s)}{U(s)} = c^T(sI-A)^{-1}b + d. \]
Diese Überführung wird in der Praxis sehr oft verwendet. Auch der umgekehrte Weg ist von großer Relevanz und wird als Realisierungsproblem bezeichnet. Gegeben sei die Übertragungsfunktion
\[ Y(s) = \frac{b_0 + b_1s + \dots + b_{n-1}s^{n-1} + b_ns^n}{a_0 + a_1s + \dots + a_{n-1}s^{n-1} + 1s^n}U(s) \]
mit einem monischen Nenner-Polynom.
Eine Möglichkeit der Realisierung ist die Regelungsnormalform. Mit den Faktoren
\[ \tilde{b}_0 = (b_0 - b_n a_0), \quad \tilde{b}_1 = (b_1 - b_n a_1), \quad \cdots \quad \tilde{b}_{n-1} = (b_{n-1} - b_n a_{n-1}) \]
können wir die Regelungsnormalform
\[ \begin{split} \dot{x} &= A_cx + b_cu \\ y &= c_c^Tx + d_cu \end{split} \rightarrow \begin{cases} \frac{d}{dt} \underset{= x}{\underbrace{ \begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ _{n-1} \\ \tilde{x}_n \end{bmatrix}}} = \underset{A_c}{\underbrace{ \begin{bmatrix} 0 & 1 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \\ -a_0 & -a_1 & \dots & -a_{n-1} \end{bmatrix}}} \underset{= x}{\underbrace{\begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ \tilde{x}_{n-1} \\ \tilde{x}_n \end{bmatrix}}} + \underset{b_c}{\underbrace{\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}}}u, \\ y = \underset{c_c^T}{\underbrace{\begin{bmatrix} \tilde{b}_0 & \tilde{b}_1 & \dots & \tilde{b}_{n-1} \end{bmatrix}}} \underset{= x}{\underbrace{\begin{bmatrix} \tilde{x}_1 \\ \tilde{x}_2 \\ \vdots \\ \tilde{x}_{n-1} \\ \tilde{x}_n \end{bmatrix}}} + \underset{d_c}{\underbrace{b_n}} u \end{cases} \]
errechnen. Für \(b_n=0\) gilt \(\tilde{b}_i=b_i\) und \(d_c = 0\), somit wird die Überführung besonders leicht.
Es sei angemerkt, dass die Realisierung von der Darstellung der Übertragungsfunktion abhängig ist. Oft werden die Polynome der Übertragungsfunktion
\[ G(s) = \frac{b_0s^m + b_1s^{m-1} + \cdots + b_m}{1s^n + a_1s^{n-1} + \cdots + a_n} \]
absteigen geordnet, d.h. mit der höchsten Potenz zuerst, wobei manchmal (wie hier) die Koeffizienten umgekehrt angeordnet sein können. Dadurch erstehen alternative Realisierung, die sich jedoch qualitativ nicht unterscheiden.
Übertragungsfunktion, Markov-Parameter und Impulsantwort
Für \(d=0\) sind die Markov-Parameter eindeutig durch die Übertragungsfunktion bestimmt, da gilt
\[ \begin{split} G(s) &= c^T\left(sE-A\right)^{-1}b \\ &= \frac{1}{s}c^T\left(E-\frac{A}{s}\right)^{-1}b \\ &= \frac{1}{s}c^T\left(E+\frac{A}{s}+\frac{A^2}{s^2}+\cdots\right)b = \frac{1}{s}\underset{m_1}{\underbrace{c^Tb}} + \frac{1}{s^2}\underset{m_2}{\underbrace{c^TAb}} + \frac{1}{s^3}\underset{m_3}{\underbrace{c^TA^2b}} + \cdots \\ &= \sum_{k=1}^{\infty} m_ks^{-k}, \end{split} \]
wobei die Markov-Parameter durch \(m_k = c^TA^{k-1}b, \quad k=1,2,\cdots\) gegeben sind.
Auch in der Impulsantwort
\[ g(t) = c^T\Phi(t)b = c^T\left( E+At+A^2\frac{t^2}{2!}+\cdots \right)b = \underset{m_1}{\underbrace{c^Tb}} + \underset{m_2}{\underbrace{c^TAb}}t + \underset{m_3}{\underbrace{c^TA^2b}}\frac{t^2}{2!} + \cdots \]
treten die Markov-Parameter auf, und kommen durch die zeitlichen Ableitungen an der Stelle \(t=0\)
\[ \frac{d^{k-1}}{dt^{k-1}}g(t)\left|_{t=0} \right. = c^TA^{k-1}b = m_k \]
exakt zum Vorschein.
Markov-Parameter, Hankel-Matrix zum Zustandsraum
Wenn die Impulsantwort \(g(t)\) geben ist, können die Markov-Parameter errechnet werden. Sind die Markov-Parameter bekannt, so konstruiert sich die Hankel-Matrix wie folgt
\[ \begin{split} H &= \begin{bmatrix} m_1 & m_2 & \cdots & m_n \\ m_2 & m_3 & \cdots & m_{n+1} \\ \vdots & \vdots & \ddots & \vdots \\ m_n & m_{n+1} & \cdots & m_{2n-1} \end{bmatrix} \\ &= \begin{bmatrix} c^Tb & c^TAb & \cdots & c^TA^{n-1}b \\ c^TAb & c^TA^2b & \cdots & c^TA^nb \\ \vdots & \vdots & \ddots & \vdots \\ c^TA^{n-1}b & c^TA^nb & \cdots & c^TA^{2n-2}b \end{bmatrix}. \end{split} \]
Die Hankel-Matrix
\[ H = \begin{bmatrix} c^T \\ c^TA \\ \vdots \\ c^TA^{n-1} \\ \end{bmatrix} \begin{bmatrix} b & Ab & \cdots A^{n-1}b \\ \end{bmatrix} = \mathcal{O}\mathcal{R} \]
lässt sich in die Beobachtbarkeitsmatrix \(\mathcal{O}(A,B)\) und in die Erreichbarkeitsmatrix \(\mathcal{R}(A,c^T)\) zerlegen. Daraus können die Matrizen \(A\), \(b\) und \(c^T\) selektiert werden.
Alle drei Darstellungen bieten Methoden der Modellierung, Identifizierung, Systemanalyse und Reglersynthese.
SISO Systeme
Viele klassische Konzepte und Reglerentwurfsverfahren können (effektiv) nur für SISO Systeme angewendet werden. SISO Systeme besitzen nur einen Eingang und einen Ausgang, womit sich sogenannte einschleifige Regelkreise ergeben.
Die meisten PID-basierten Autotuning-Verfahren und Faustformeln wurden für SISO Systeme entwickelt und sind im SISO Fall einfach anzuwenden. Die PID-Einstellregeln nach Ziegler-Nichols haben eine große Bekanntheit erlangt. Es gibt aber hunderte alternative Verfahren, wobei die konkrete Anzahl von der Zählweise abhängt. Die meisten Verfahren basieren auf einfachen Ersatzmodellen realer Strecken. Die praktische Bedeutung der PID Regelung darf auf keinen Fall unterschätzt werden. Noch immer gibt es Industriezweige, wo die eingesetzten Regler sich wie folgt verteilen:
- ca. 90% sind PI-Regler
- ca. 5% sind PID-Regler
- ca. 5% sind fortschrittliche Regler, insbesondere MPC-Variationen
Diese Verteilung ist insbesondere in der Prozessautomatisierung anzutreffen. PID Regler werden aber auch in der modernen Raumfahrt (z.B. CHEOPS Weltraumteleskop) eingesetzt.
Der Frequenzbereich erlaubt elegante Entwurfsverfahren. Verschiedene Frequenzkennlinienverfahren stellen sich bei etwas Übung als erstaunlich effektiv heraus und sind sehr gut nachvollziehbar. Vor allem Robustheitseigenschaften lassen sich dort mithilfe grafischer Tools wie dem Bode-Diagramm oder der Nyquist-Ortskurve gut darstellen. Auch robuste Synthesen wie \(\mathcal{H}_2, \mathcal{H}_{\infty}, ...\) sind für SISO Systeme mithilfe geeigneter Programmpakete relativ einfach möglich.
Bevor wir uns aber dem Reglerdesign widmen, müssen wir Modelle anlegen. Wir starten mit Übertragungsfunktionen.
Referenzen
- Feedback Systems: An Introduction for Scientists and Engineers, Second Edition, 2020/2022 (Karl J. Aström, Richard M. Murray)
- Multivariable Feedback Control: Analysis and Design, 2005 (Sigurd Skogestad, Ian Postlethwaite)
- Data-Driven Science and Engineering: Machine Learning, Dynamical Systems and Control, 2019 (Steven L. Brunton, J. Nathan Kutz)