2. Anforderungen und Methoden für die Hardware-in-the-Loop Simulation zur Virtuellen Inbetriebnahme von Produktionssystemen

Heutige Produktionsanlagen sind hochkomplexe mechatronische Systeme. Die Wirkzusammenhänge zwischen den Anlagenkomponenten sind aufgrund des hohen Automatisierungsgrades derart komplex, dass deren Test und Inbetriebnahme eine technische Herausforderung darstellt. Zudem wird ein hohes Maß an Zuverlässigkeit und Fehlersicherheit vorausgesetzt. Daraus folgt ein hoher Qualitätsanspruch an die eingesetzte Steuerungssoftware und den Softwareentwicklungsprozess. Diesem kann nur durch den Einsatz moderner Entwicklungs- und Testmethoden Folge geleistet werden. Ohne solche modernen Methoden liegt die Gewissheit über die Funktion und das Zusammenspiel aller Automatisierungskomponenten erst bei der Inbetriebnahme der Anlage vor. Treten dann technische Probleme auf, bleiben zeitaufwendige und kostspielige Nacharbeiten meist nicht aus, insbesondere wenn die Anlage bereits beim Kunden steht. Dieses Risiko kann durch frühzeitiges Testen mit Hilfe moderner Simulationsmethoden minimiert werden.

Die Virtuelle Inbetriebnahme (VIBN) ist eine Simulationsmethode, bei der die reale Anlage ganz oder teilweise durch eine simulierte Anlage ersetzt und mit der realen Steuerungstechnik gekoppelt wird, mit dem Ziel das Steuerungssystem frühzeitig gegen das simulierte Verhalten der Produktionsanlage zu testen.

Anlagen- und Steuerungshersteller sehen in der VIBN das Potenzial die Entwicklungszeiten durch frühzeitige Tests zu reduzieren und das Risiko bei der realen Inbetriebnahme zu minimieren. Im Fokus stehen dabei die Einsparung von Kosten sowie die Steigerung der Wettbewerbsfähigkeit. Die VIBN hat folgende Vorteile:

Für die VIBN mit realen Steuerungssystemen eignet sich die aus dem mechatronischen Entwicklungsprozess bekannte Hardware-in-the-Loop Simulation (HiLS). Im Gegensatz zur Software-in-the-Loop Simulation (SiLS), bei der sowohl die Steuerung als auch die Anlage simuliert wird [1], wird bei der HiLS das reale Steuerungssystem inkl. Steuerungshardware verwendet. Dieses wird über den realen Feld- oder Antriebsbus mit einem Simulationsrechner verbunden, auf dem das Anlagenverhalten simuliert wird [2]. Der Vorteil dabei ist, dass aus Sicht der Steuerung exakt die gleichen Bedingungen wie im realen Betrieb gelten, was bildlich gesprochen einem „Umschalten“ von der realen Anlage auf eine virtuelle Anlage gleichkommt (Abb. 2.1).

Die HiLS stellt hohe Anforderungen an die Simulationsumgebung. Im Folgenden sind einige wichtige Anforderungen zusammengestellt.

Um das reale Bussystem einsetzen zu können, muss eine Hardwareschnittstelle am Simulationsrechner für die entsprechende Bushardware zur Verfügung stehen. Das Kommunikationsverhalten der Antriebs- bzw. Feldbus-Teilnehmers muss vollständig in der Simulation abgebildet sein, um einen störungsfreien Betrieb der Steuerung zu ermöglichen. Um die Daten der Steuerung vom Bussystem abgreifen zu können, muss die Simulationsumgebung über einen Bustreiber verfügen, der den Zugriff auf diese Daten in Echtzeit ermöglicht. Die Simulation mehrerer Feldbus- bzw. Antriebsbuskomponenten erfordert zudem, dass sowohl der Bustreiber als auch die Hardwareschnittstelle am Simulationsrechner den gleichzeitigen Betrieb mehrerer Busteilnehmer unterstützen.

Eine verlustfreie und zeitsynchrone Kommunikation zwischen Steuerung und den Feldbusteilnehmern ist im realen Betrieb zwingend notwendig. Dies muss auch für die Kommunikation zwischen Steuerung und Simulation gelten, da sonst das Zeitverhalten der virtuellen Anlage im Vergleich zur realen Anlage verzerrt wird und damit keinem realistischen Verhalten entspricht. Es muss stets gewährleistet sein, dass nach jedem Steuerungstakt ein Berechnungstakt der Simulation erfolgt. Die berechneten Ausgabewerte der Simulation müssen rechtzeitig an der Busschnittstelle zur Verfügung stehen. Um den geforderten Zeitdeterminismus zu gewährleisten, muss die Simulation echtzeitfähig sein. In der Simulationstechnik ist es üblich, die Echtzeitfähigkeit einer Simulation durch den Echtzeitfaktor (EZF) anzugeben. Der EZF gibt das Verhältnis zwischen der tatsächlich verbrauchten Verarbeitungszeit eines Simulationsschritts \({T}_{V}\) (Modellrechenzeit und Verarbeitungszeit des Betriebs- und Kommunikationssystems) und einer fest vorgegebenen Simulationszeitschrittweite \({h}_{SIM}\) an. Eine Simulation ist echtzeitfähig, wenn für jeden Simulationszeitschritt folgende Bedingung eingehalten wird:

Um das Systemverhalten von Produktionsmaschinen und -anlagen realistisch abbilden zu können, sind Simulationsmodelle erforderlich, die das dynamische Prozess- und Bewegungsverhalten der Anlage wiedergeben. Dieses wird im Wesentlichen durch den Produktionsprozess, die Maschinenmechanik und die Antriebsregelung bestimmt. Die daraus resultierenden Differentialgleichungen [3] müssen numerisch gelöst werden. Das Lösen von Differentialgleichungen in einer HiLS unter der Echtzeitbedingung (\(\text{EZF }\leq 1\)) erfordert numerische Rechenverfahren mit hoher Recheneffizienz und Robustheit. Dabei dürfen die durch das Rechenverfahren unumgänglichen numerischen Fehler unter keinen Umständen zu Fehlurteilen bei der Inbetriebnahme führen.

Abb. 2.2 zeigt beispielhaft den Aufbau eines Hardware-in-the-Loop Simulators.

In den folgenden Abschnitten werden die wichtigsten Komponenten und Funktionalitäten aus Abb. 2.2 erläutert.

Der Echtzeit-Rechenkern ist mit numerischen Lösungsverfahren zur Lösung der Modellgleichungen ausgestattet. Für eine HiLS ist eine zeitdeterministische und effiziente Berechnung erforderlich, um stets \(\text{EZF }\leq 1\) zu erreichen, was hohe Anforderungen an die Lösungsverfahren stellt. Um den Zeitdeterminismus gewährleisten zu können, muss der Rechenkern zudem auf einem echtzeitfähigen Simulationsrechner mit einem Echtzeit-Betriebssystem ausgeführt werden. Echtzeit-Betriebssysteme ermöglichen ein echtzeitfähiges Taskscheduling sowie eine echtzeitfähige Speicherverwaltung und sind für den Einsatz auf heutigen Standard-Rechnerarchitekturen (Standard-PC) auf dem Markt verfügbar.

Der Echtzeit-Rechenkern stellt Schnittstellen zum Austausch von Daten zur Verfügung. Ein zentraler Modellspeicher im globalen Adressraum (Shared Memory) ermöglicht den Datenaustausch zwischen Modell, Bustreiber im Echtzeit-Betriebssystem (Kernel-Space) und Benutzungsoberflächen im Standard-Betriebssystem (User-Space). Während im Kernel-Space die Synchronität der Zugriffe überwacht wird und eine maximale Zugriffszeit nicht überschritten werden darf, findet der Datenzugriff im User-Space komplett asynchron statt. Aufgrund der Asynchronität müssen die Daten zwischen dem Modellspeicher und den Anwendungen in FIFO-Datenpuffern (FIFO: First In First Out) zwischengespeichert werden, um einen Datenverlust zu vermeiden. So werden die berechneten Zustandsdaten (ZD: Innere Zustandsgrößen und modellierte Ausgangsgrößen zur Modellbeobachtung) im Kernel-Space nach jedem Simulatortakt aus dem Modellspeicher in die Datenpuffer zur Weiterverarbeitung im User-Space geschrieben. Umgekehrt werden die in den Datenpuffern zwischengespeicherten Eingangsdaten (ED: Modellierte Eingangsgrößen zur Modellbeeinflussung) aus dem User-Space vor jedem Simulatortakt auf den Modellspeicher im Kernel-Space geschrieben. Die Modelldaten (MD: Daten der Modellbeschreibung und -parametrierung) werden nur einmalig vor Beginn der Simulation zur Initialisierung des Modells und des Modellspeichers in den Kernel-Space übertragen.

Benutzungsschnittstellen (User-Interfaces) ermöglichen das Modellieren, Beeinflussen und Beobachten des Simulationsszenarios. Wie bereits erwähnt, können Daten aus dem Echtzeit-Rechenkern im Kernel-Space mit Anwendungen im User-Space ausgetauscht werden. Dies ermöglicht das Ankoppeln von Anwendungen zur Beeinflussung und Beobachtung der Simulation während der Laufzeit. Die Online-Beobachtung erlaubt dabei die Überwachung und Aufzeichnung aller berechneten Zustandsdaten (ZD) bspw. mit Hilfe von Zeitschrieben. Mittels Online-Beeinflussung können die Eingangsdaten (ED) ins Modell vom Benutzer manuell geschrieben werden. Dies ermöglicht die Einbindung von Bedienereingriffen, wie das Lösen einer Schutzeinrichtung oder das Betätigen eines Not-Aus-Schalters. Zudem können Störgrößen in das System eingebracht werden, wie Parameterschwankungen, Temperatureinflüsse oder Systemausfälle. Auf Basis der Benutzungsschnittstellen können externe Benutzungsoberflächen an den Simulator angeschlossen werden (bspw. über Ethernet mittels TCP/IP-Sockets), um einen individualisierten virtuellen Anlagenprüfstand für die Virtuelle Inbetriebnahme zu realisieren. Dies ermöglicht bspw. das Anbinden von 3D-Visualisierungswerkzeugen.

Bei einer HiLS werden zwischen Steuerung und Simulation Daten in Echtzeit über die realen Feld- und Antriebsbussysteme übertragen. Für einen verlustfreien Datenaustausch zwischen Steuerungs- und Simulationsrechner ist die Taktung und Zeitsynchronität des Simulators von Bedeutung. Im Folgenden wird zwischen zwei Arten der Simulatorkopplung unterschieden: der selbstgetakteten asynchronen Simulatorkopplung und der fremdgetakteten synchronen Simulatorkopplung.

Bei dieser Kopplungsart wird der Simulator über die eigene Systemuhr synchron zur realen Zeit getaktet. Dabei ist die Zeitbasis der Steuerung physikalisch unabhängig von der Zeitbasis des Simulators. Die Daten der Steuerung werden über das Bussystem an den Bustreiber übertragen (vgl. Abb. 2.2). Der Simulator tastet den Speicher des Bustreibers im Simulatortakt ab (Polling). Da der Bustreiber auf die Systemuhr der Steuerung synchronisiert ist und der Simulator auf die Systemuhr des Simulationsrechners, kommt es aufgrund der Asynchronität zwischen den beiden Systemuhren zu Zeitverschiebungen in der Datenübertragung. Solche Zeitverschiebungen können zu Fehlern in der Datenübertragung und somit im Modellverhalten führen, was Fehleinschätzungen bei der Virtuellen Inbetriebnahme zur Folge haben kann.

Der Vorteil dieser Kopplungsart ist die softwaretechnisch einfache Realisierbarkeit, da der Datenaustausch nur ein einfaches Abtasten von Speicherbereichen erforderlich macht. Viele HiL-Simulatoren mit geringen Anforderungen an die Synchronität zwischen Steuerung und Simulation beruhen auf diesem Prinzip.

Das Prinzip der fremdgetakteten Simulatorkopplung basiert darauf, dass Steuerung und Simulation auf Basis einer gemeinsamen Systemuhr auf die reale Zeit synchronisiert werden. Grundsätzlich wäre eine Synchronisierung zwischen Steuerung und Simulation gemäß den Prinzipien, die beim Network Time Protocol (NTP) oder beim Precision Time Protocol (PTP) in Computernetzwerken zum Einsatz kommen, möglich. Diese Protokolle verwenden als Transportprotokoll UDP (NTP auch TCP) und stellen spezifische Anforderungen an das eingesetzte Kommunikationssystem. Um zusätzliche Kommunikationsschnittstellen für die Synchronisierung der Rechner und einen daraus resultierenden Eingriff in das Steuerungssystem zu vermeiden, wird im Folgenden eine Methode vorgeschlagen, die zur Synchronisation den Bus-Interrupt ausnützt.

Ausgehend vom Busmaster in der Steuerung wird ein Bus-Interrupt auf dem Simulationsrechner ausgelöst, der wiederum den Bustreiber aufruft (vgl. Abb. 2.3). Der Bustreiber stellt die Kommunikation zwischen Bus und Simulation sicher und beauftragt das Betriebssystem die bereits wartende Simulationstask mittels Mechanismen der Interprozesssynchronisation (bspw. Semaphore) zu starten. In der Simulationstask wird das Anlagenmodell im Simulationsintervall [\({\sigma }_{n-1},{\sigma }_{n}\)] berechnet und die neu berechneten Zustandsdaten aktualisiert. Danach gibt die Simulationstask die Rechenzeit an das Betriebssystem bis zum nächsten Simulationsintervall [\({\sigma }_{n},{\sigma }_{n+1}\)] zurück.

Die Reaktionszeit zwischen Bus-Interrupt und Aktualisierung der Zustandsdaten setzt sich aus der Kommunikationsdauer T_COM, der Reaktionsdauer des Betriebssystems T_OS und der Simulationsdauer T_SIM zusammen. T_COM beinhaltet dabei die Bearbeitungszeit des Bustreibers, T_OS ist der Zeitraum, den das Betriebssystem benötigt, um die Simulationstask aufzurufen und T_SIM ist die benötigte Rechenzeit eines Simulationsintervalls. Um eine synchrone und verlustfreie Kommunikation sicherstellen zu können, muss die Dauer vom Bus-Interrupt bis zur Aktualisierung der Zustandsdaten immer kleiner oder gleich dem Bustakt h_BUS sein.

Da bei der fremdgetakteten Simulatorkopplung Steuerung und Simulation immer absolut synchron laufen, kann der Simulatortakt h_SIM gleich dem Steuerungs- bzw. Bustakt h_BUS gewählt werden. Somit kann die Bedingung aus Gl. (2.1) erweitert werden:

Bei Antriebsschnittstellen einer numerischen Bahnsteuerung (CNC) sind Bustakte von h_BUS \(=[1\dots 5]\,ms\) üblich, während Ablaufsteuerung (PLC) mit Bustakten h_BUS \(=[1\dots 200]\,ms\) betrieben werden. Die Einhaltung der Echtzeitbedingung nach Gl. (2.2) kann nur durch einen echtzeitfähigen Bustreiber, ein echtzeitfähiges Betriebssystem und eine echtzeitfähige Berechnung des Anlagenmodells erfüllt werden.

Diese Art der Simulatorkopplung kommt bei HiLS-Anwendung mit sehr hohen Anforderungen an die Synchronität zwischen Steuerung und Simulation zum Einsatz, wie bspw. bei der HiLS von NC-gesteuerten Servoantrieben.

Bei der Inbetriebnahme von hochdynamischen und stark belasteten Maschinenachsen, wie sie in einer modernen Produktion heute üblich sind, spielen die mechanischen Struktureigenschaften der Maschine in Kombination mit den Antriebsregelungen eine große Rolle. Nachgiebigkeiten und Eigenschwingungen steifer Strukturen lassen sich in guter Näherung durch lineare Modelle approximieren. Dabei kommen häufig lineare Finite-Elemente-Modelle (FE-Modelle) gekoppelt mit linearen Antriebs-Regler-Modellen zum Einsatz. Die schwachgedämpften und steifen metallischen Maschinenstrukturen weisen mehrere sehr unterschiedliche und zugleich sehr große Eigenschwingfrequenzen auf, während die angekoppelten Antriebs-Regler-Modelle oft sehr kleine elektrische Zeitkonstanten einbringen. Diese sogenannten steifen Systeme sind für die numerische Berechnung im Zeitbereich ungünstig, da große Eigenschwingfrequenzen und kleine Zeitkonstanten sehr kleine Zeitschrittweiten für eine robuste und ausreichend genaue Lösung erfordern und dies wiederum die Recheneffizienz vermindert.

Die Grundidee für die effiziente echtzeitfähige Berechnung linearer Systeme ist das Gesamtsystem in kleinere unabhängige Teilsysteme zu zerlegen, um jedes Teilsystem für sich einzeln lösen zu können. Dies ermöglicht neben einer höheren Recheneffizienz auch eine unabhängige numerische Betrachtung jedes Teilsystems, sodass Systemteile mit großen Eigenschwingfrequenzen numerisch anders behandelt werden können, wie die Systemteile mit kleinen Eigenschwingfrequenzen. Diesen Vorgang nennt man Systementkopplung. Abb. 2.4 zeigt beispielhaft die Systementkopplung eines FE-Modells von einem Gantry-Direktantrieb zur Vorbereitung für eine Echtzeitsimulation in einem HiL-Simulator.

Die Systemmatrizen des linearen FE-Modells werden aus der FE-Software exportiert und bspw. in Matlab/Simulink eingelesen und dort um lineare Antriebs-Regler-Modelle ergänzt. Nach Berechnung der Systementkopplung in Matlab werden die Matrizen des entkoppelten Systems dem HiL-Simulator zur Echtzeitsimulation (mittels numerischer Integration der linearen Teilsysteme) zur Verfügung gestellt [4, 5]. Im Folgenden wird die Systementkopplung linearer Systeme näher erläutert.

Die Modellfunktion eines linearen Systems mit konstanten Koeffizienten kann allgemein durch die lineare Zustandsraumdarstellungmit den konstanten Matrizen \(\textbf{A}\in {\mathbb{R}}^{w\times w}\) (Systemmatrix), \(\textbf{B}\in {\mathbb{R}}^{w\times v}\) (Eingangsmatrix), \(\textbf{C}\in {\mathbb{R}}^{r\times w}\) (Ausgangsmatrix) und \(\textbf{D}\in {\mathbb{R}}^{r\times v}\) (Durchgangsmatrix) angegeben werden. Dabei ist \(x\in {\mathbb{R}}^{w\times 1}\) der Vektor der Zustandsgrößen, \(y\in {\mathbb{R}}^{r\times 1}\) der Vektor der Ausgangsgrößen und \(u\in {\mathbb{R}}^{v\times 1}\) der Vektor der Eingangsgrößen.

Das Eigenschwingverhalten eines linearen Systems ist vollständig durch die Eigenwerte und Eigenvektoren (Moden) der Systemmatrix \(\textbf{A}\) beschrieben. Die Systemmatrix \(\textbf{A}\) ist eine quadratische Matrix, die mittels Ähnlichkeitstransformation in eine andere Darstellungsform transformiert werden kann, ohne dass sich deren Eigenwerte ändert. Geht man davon aus, dass eine reguläre zeitinvariante Transformationsmatrix \(\textbf{T}\) existiert und transformiert man damit die Zustandsgrößen und deren zeitliche Ableitungen

erhält man durch Einsetzen in Gl. (2.3) eine neue Systemdarstellung mit

Je nach Wahl der Transformationsmatrix \(\textbf{T}\) kann das System in verschiedene Darstellungsformen, in denen das System in unabhängige Teilsysteme entkoppelt ist, transformiert werden.

Der höchste Grad der Entkopplung ist erreicht, wenn die transformierte Systemmatrix \(\textbf{A}^{\prime}\) Diagonalform besitzt. Dabei entsprechen die Diagonalelemente den Eigenwerten \({\lambda }_{i}\) der ursprünglichen Systemmatrix \(\textbf{A}\). Die dazu notwendige Transformationsmatrix ist die Matrix der Eigenvektoren \({e}_{i}\) von \(\textbf{A}\) (Modalmatrix) [6, 7]:

Nach Transformation mit der Modalmatrix \({\textbf{T}}_{1}\) erhält man eine diagonalisierte Systemmatrix, unter der Voraussetzung, dass die Systemmatrix A ausschließlich einfache Eigenwerte λ₁≠λ₂≠⋯≠λ_w besitzt:

Treten komplexe Eigenwerte auf, ist eine weitere Transformation erforderlich, um die imaginären Anteile zu eliminieren. Die dazu notwendige Transformationsmatrix \({\textbf{T}}_{2}\) kann je nach Auftreten der komplexen Eigenwerte, wie in Abb. 2.5 beispielhaft illustriert, konstruiert werden.

Durch erneute Transformation der Diagonalmatrix mit der Transformationsmatrix \({\textbf{T}}_{2}\) erhält man die reelle Blockdiagonalmatrix des entkoppelten Systems:

Dabei ist \({\omega }_{i}\) die gedämpfte Eigenkreisfrequenz und \({\delta }_{i}\) die Abklingkonstante des Teilsystems \({\textbf{A}^{\prime}}_{i}\). Die transformierte Systemmatrix \(\textbf{A}^{\prime}\) ist vollständig entkoppelt und kann in unabhängige Teilsysteme \({\textbf{A}^{\prime}}_{i}\) zerlegt werden. Die maximale Größe eines Teilsystems ist dabei \({w}_{i}=2\).

Bei Systemen mit mehrfachen Eigenwerten \({\lambda }_{1}={\lambda }_{2}=\cdots ={\lambda }_{\mu }\) ist die Systemmatrix nicht zwingend vollständig diagonalisierbar. Es können sogenannte Jordanblöcke \({\textbf{J}}_{{\varvec{i}}}\) mit der Dimension \(\mu >1\) auftreten, die nicht entkoppelt werden können. Die daraus resultierende Jordanmatrix hat die Blockdiagonalform:

Die Koppelelemente (1-en) auf der Nebendiagonale, auch als Jordanbrücken bezeichnet, verhindern die weitere Zerlegung der Jordanblöcke. Treten mehrfache komplexe Eigenwerte auf, erfolgt die Transformation in eine reelle Jordanmatrix durch Konstruktion einer entsprechenden weiteren Transformationsmatrix, analog zur Vorgehensweise bei den einfachen Eigenwerten.

Die oben beschriebene Systementkopplung ermöglicht die Zerlegung des Systems in \(m\) kleinere Teilsysteme, die unabhängig numerisch gelöst werden können:

Zur Lösung der linearen Teilsysteme kann die zeitdiskrete Exponentiallösung herangezogen werden. Diese Methode ermöglicht ein äußerst robustes und sehr genaues numerisches Lösen der Teilsysteme. Die Differentialgleichung eines Teilsystems besitzt die Exponentiallösungmit den Anfangswerten \({\kappa }_{i,0}\) zum Zeitpunkt \({t}_{0}\). Das Integral in Gl. (2.12) lässt sich für beliebige Eingangsgrößen \(u\left(t\right)\) nicht analytisch lösen und muss durch eine numerische Näherungslösung ersetzt werden. Dazu werden die Eingangsgrößen über ein Abtast-Halteglied 0-ter Ordnung approximiert. So kann das Integral für stückweise konstante Eingangsgrößen \(u({t}_{n})\) im Zeitintervall \(\left[{t}_{n},{t}_{n+1}\right]\) gelöst werden und man erhält die Näherungslösungwobei \(\textbf{I}\) die Einheitsmatrix darstellt [6]. Durch Einsetzen in Gl. (2.12) erhält man schließlich die Lösung des Teilsystems in diesem Zeitintervall:

Aufgrund der Echtzeitanforderung und des damit verbundenen Zeitdeterminismus erfolgt die Diskretisierung in äquidistanten Zeitabständen mit der konstanten Zeitschrittweite \(h={t}_{n+1}-{t}_{n}\), die gleich oder ein ganzzahliger Teiler vom Simulatortakt \({h}_{SIM}\) sein muss.

Zusammenfassend lautet die Berechnungsvorschrift für die Berechnung eines linearen Teilsystems aus Gl. (2.11):mit den Verfahrensmatrizen

Diese Lösungsformel konvergiert immer unabhängig von der gewählten Schrittweite und der Steifigkeit des Systems. Der numerische Fehler resultiert lediglich aus der diskreten Abtastung der Eingangsgrößen.

Nach Lösen aller Teilsysteme für ein Zeitintervall, können die Ausgangsgrößen entsprechend Gl. (2.17) anhand der Zustandsgrößen für den aktuellen Zeitschritt berechnet werden:

Zum Vergleich soll der Berechnungsaufwand eines entkoppelten Systems mit dem zugrundeliegenden nicht entkoppelten System gegenübergestellt werden. Die Verfahrensmatrizen \({\varvec{\Phi}}\) und \({\varvec{\Omega}}\) der einzelnen Teilsysteme in Gl. (2.16) sind aufgrund der konstanten Matrizen \({\textbf{A}}_{i}^{\prime}\) und \({\textbf{B}}_{i}^{\prime}\) sowie der konstanten Zeitschrittweite \(h\) ebenfalls konstant und müssen nur einmalig in der Initialisierungsphase der Simulation berechnet werden. Somit beschränkt sich der Berechnungsaufwand in einem Simulationsschritt im Wesentlichen auf das Auswerten der Matrix-Vektor-Produkte in Gl. (2.15) und (2.17). Zur Quantifizierung werden die Additionen und Multiplikationen als Gleitpunktoperationen in FLOP (Floating Point Operations) gezählt. Anhand des folgenden vereinfachten Pseudocodes wird die erforderliche Anzahl an Gleitpunktoperationen für eine Matrix-Vektor-Multiplikation \(x=\textbf{A}\cdot b\) mit \(\textbf{A}\in {\mathbb{R}}^{n\times m}\), \(b\in {\mathbb{R}}^{m\times 1}\) abgeschätzt:

Für das Auswerten der Matrix-Vektor-Multiplikation sind demnach \(2mn\) FLOP erforderlich, wobei darauf hinzuweisen ist, dass aus Gründen der besseren Verständlichkeit auf mögliche Optimierung des Codes, um bspw. Additionen und Multiplikationen mit Null einzusparen, verzichtet wird. Überträgt man die Auswertung der Matrix-Vektor-Produkte entsprechend Gl. (2.15) und (2.17) auf das nichtentkoppelte System mit den Matrizen \(\textbf{A}\in {\mathbb{R}}^{w\times w}\), \(\textbf{B}\in {\mathbb{R}}^{w\times v}\), \(\textbf{C}\in {\mathbb{R}}^{r\times w}\) und \(\textbf{D}\in {\mathbb{R}}^{r\times v}\) werden pro Simulationsschritt

benötigt. Für die Bewertung des Rechenaufwands des entkoppelten Systems ist die Größe der Teilsysteme, die aus der Systementkopplung hervorgehen, ausschlaggebend. Die Dimension der Teilsysteme sei angegeben mit:

wobei \(\vartheta\) festlegt, ob es sich um einen reellen \((\vartheta =1)\) oder einen komplexen \((\vartheta =2)\) Eigenwert handelt und \(\mu\) die Dimension des entstehenden Jordanblockes bei einem mehrfachen Eigenwert angibt.

Die Anzahl der Rechenoperationen für das entkoppelte System ist die Summe der Matrix-Vektor-Produkte aller \(m\) Teilsysteme mit der Größe \({w}_{i}=p({\vartheta }_{i},{\mu }_{i})\). Die Auswertung der Berechnungsvorschrift nach Gl. (2.15) und (2.17) für alle Teilsysteme des entkoppelten Systems benötigt somit

Die Matrix-Vektor-Produkte für die Eingangs- und Ausgangsgrößen werden dabei analog zu dem nichtentkoppelten System berechnet, da die Eingangs-, Ausgangs- und Durchgangsmatrizen auch nach der Entkopplung i.A. vollbesetzte Matrizen sind.

Zur Veranschaulichung soll der Rechenaufwand eines nicht entkoppelten Systems mit einem entkoppelten System am Beispiel des FE-Modells in Abb. 2.4 betrachtet werden. Das Beispielmodell besitzt \(v=11\) Eingangsgrößen und \(r=17\) Ausgangsgrößen und keinen direkten Durchgriff der Eingangsgrößen auf die Ausgangsgrößen, wodurch \(\textbf{D}=\textbf{0}\) eine Nullmatrix wird und demzufolge nicht berücksichtigt werden muss. Damit ergibt sich für das nicht entkoppelte FE-Modell aus Gl. (2.18) in Abhängigkeit der \(w\) Zustandsgrößen

Zur Gegenüberstellung des Berechnungsaufwands soll das Gesamtmodell in verschiedenen Detaillierungsstufen mit \(w=\left\{\text{50,100,500,1000}\right\}\) Zustandsgrößen betrachtet werden. Das System besitzt ausschließlich einfache komplexe Polstellen und kann somit in gleich große Teilsysteme mit der Dimension \({w}_{i}=p\left(\text{2,1}\right)=2\) zerlegt werden. Damit ergibt sich für das entkoppelte FE-Modell aus Gl. (2.20)

Eine Gegenüberstellung der erforderlichen Gleitpunktoperationen für das nicht entkoppelte FE-Modell mit dem entkoppelten FE-Modell für verschiedene Systemgrößen ist in Tab. 2.1 zusammengefasst.

Die Anzahl der Zustandsgrößen \(w\) ist im Wesentlichen abhängig von der Detaillierung des FE-Modells. Zur Reduzierung der Modellgröße können Modellreduktionsverfahren zum Einsatz kommen, auf die jedoch in diesem Beitrag nicht weiter eingegangen werden soll. Eine Übersicht über gängige Modellreduktionsverfahren wird in [14] beschrieben.

Anhand des Beispiels ist gut zu erkennen, dass die Entkopplung und Zerlegung eines linearen dynamischen Systems in unabhängige Teilsysteme eine deutliche Reduzierung der Rechenoperationen pro Zeitschritt ermöglicht. Vor allem bei Systemen mit vielen Zustandsgrößen und vergleichsweise wenigen Ein- und Ausgangsgrößen macht sich ein hoher Entkopplungsgrad deutlich bemerkbar und ist häufig die Grundvoraussetzung für eine effiziente Simulation in Echtzeit. Zudem ist eine Entkopplung zwischen Systemteilen mit großen Eigenwerten und Systemteilen mit kleinen Eigenwerten gewährleistet und die Problematik bei der numerischen Behandlung steifer Systeme einfacher beherrschbar, da jedes Teilsystem unabhängig numerisch integriert werden kann.

Die Modellierung des dynamischen Prozess- und Bewegungsverhaltens von Produktionsmaschinen führt im Allgemeinen auf nichtlineare Differentialgleichungen für die keine allgemeine Methode zur Systementkopplung, wie im vorangegangenen Abschnitt beschrieben, existiert. Zur Lösung nichtlinearerer Differentialgleichungen wird üblicherweise auf klassische numerische Integrationsverfahren zurückgegriffen.

Jedes explizite gewöhnliche Differentialgleichungssystem kann in ein System von Differentialgleichungen 1. Ordnung überführt werden (Zustandsraumdarstellung):

Dabei ist \(x\left(t\right)\) der Vektor der Zustandsgrößen, \(u\left(t\right)\) der Vektor der Eingangsgrößen und \(f\) die Vektorfunktion der Modellgleichungen. Dieses Differentialgleichungssystem kann intervallweise im Zeitintervall \(\left[{t}_{n},{t}_{n+1}\right]\) numerisch gelöst werden mitwobei das Integral durch eine Näherungsformel approximiert wird:

Dabei ist \(Q\left(f\right)\) die approximierte Integralfläche von \(f\) (auch Quadraturformel genannt) und \(E\left(f\right)\) der Verfahrensfehler im Integrationsintervall \(\left[{t}_{n},{t}_{n+1}\right]\). Das Zeitintervall legt einen Integrationszeitschritt mit der Zeitschrittweite \(h={t}_{n+1}-{t}_{n}\) fest. Zur Veranschaulichung sei das implizite Trapezverfahren angegeben mit

Einen soliden Überblick über weitere Verfahren gibt [6] und eine sehr ausführliche theoretische Beschreibung wird in [8, 9, 10] gegeben.

Die Genauigkeit und Stabilität des numerischen Integrationsverfahrens sind vom Verfahren selbst und von den Eigenschaften des Systems, sprich von der Modellfunktion \(f\), abhängig. Um den Zusammenhang zwischen den Systemeigenschaften und der numerischen Integration herzustellen, wird die nichtlineare Modellfunktion \(f\) in einem beliebigen Betriebspunkt [\({x}_{B}\), \({u}_{B}\)] linearisiert. Man erhält die lineare Modellfunktion

wobei die Jacobi-Matrixdie Systemmatrix des linearisierten Systems im Betriebspunkt [\({x}_{B}\), \({u}_{B}\)] darstellt. Die Eigenwerte dieser Systemmatrix \({\lambda }_{i}\) beschreiben die dynamischen Eigenschaften im Betriebspunkt des Systems und bilden multipliziert mit der Zeitschrittweite \(h\) der numerischen Integration die Stabilitätsvariablen des Systems [8]

Abb. 2.6 zeigt exemplarisch den Einfluss einer Stabilitätsvariablen \(z\) auf die Genauigkeit und Stabilität numerischer Integrationsverfahren. In Abb. 2.6 (links) ist der Verlauf des numerischen Fehlers doppellogarithmisch über die Stabilitätsvariable \(z\) aufgetragen. Der numerische Fehler setzt sich aus dem Rundungsfehler des Rechensystems und dem Verfahrensfehler zusammen. Der Rundungsfehler dominiert bei sehr kleinen \(z\) nahe der Maschinenpräzision. Aufgrund der hohen Zahlendarstellungsgenauigkeit heutiger 64Bit-Rechensysteme ist dieser Fehlereinfluss vernachlässigbar klein und wird daher hier nicht weiter vertieft. Viel größeren Einfluss hat der Verfahrensfehler. Dieser hängt direkt von der Konsistenzordnung \(r\) des Integrationsverfahrens ab. Dabei gilt vereinfacht ausgedrückt: Je höher die Konsistenzordnung des Verfahrens, desto genauer die numerische Lösung.

Abb. 2.6 (rechts) zeigt, welchen Einfluss die Stabilitätsvariable \(z\) auf die Stabilität der numerischen Integration hat. Jedes numerische Integrationsverfahren definiert ein Stabilitätsgebiet \(G\) in der komplexen Zahlenebene. Schließt das Stabilitätsgebiet \(G\) des Verfahrens alle Stabilitätsvariablen \({z}_{i}=h{\lambda }_{i}\) des Systems ein, konvergiert die Lösung der numerischen Integration. Die Form des Stabilitätsgebietes \(G\) hängt von der Wahl des Integrationsverfahrens ab.

Da die Eigenwerte und damit auch die Stabilitätsvariablen nur für lineare Systeme definiert sind, lassen sich daraus resultierende Genauigkeits- und Stabilitätsaussagen nur begrenzt auf nichtlineare Systeme übertragen. Jedoch kann es für ein besseres qualitatives Verständnis der Zusammenhänge zwischen den Systemeigenschaften und dem numerischen Verfahren beitragen, wie folgende Überlegung zeigt. Linearisiert man ein nichtlineares System in jedem Zeitschritt, erzeugen die daraus hervorgehenden Stabilitätsvariablen jeweils eine zeitvariante Ortskurve \({z}_{i}\left({t}_{n}\right)\) in der komplexen Zahlenebene. Bei konstanter Zeitschrittweite \(h\) bedeutet dies, dass eine stabile Lösung nichtlinearer Differentialgleichungen nicht garantiert werden kann, da die Ortskurven aus dem Stabilitätsgebiet des numerischen Integrationsverfahrens heraus laufen können.

Es existieren zwar numerische Integrationsverfahren, deren Stabilitätsgebiete die komplette linke Hälfte der komplexen Zahlenebene abdecken und somit unabhängig von der Lage der Stabilitätsvariablen eine absolut stabile Integration ermöglichen (A-Stabilität), die Dahlquist-Barriere besagt jedoch, dass diese Eigenschaft nur implizite Verfahren haben können [10]. Bei impliziten Verfahren treten unbekannte Lösungen \(x_{n+1}\) und Funktionswert \(f_{n+1}\) in der Quadraturformel \(Q(f)\) auf (vgl. Gl. 2.26), was dazu führt, dass diese Verfahren, insbesondere bei großen Systemen, sehr rechenaufwendig und in einer Echtzeitsimulation nur begrenzt einsetzbar sind.

Um eine effiziente und dennoch robuste Integration zu ermöglichen, werden häufig die deutlich recheneffizienteren expliziten Verfahren, bei denen nur bekannte Lösungen \(x_{n-i}\) und Funktionswerte \(f_{n-i}\) in die Quadraturformel \(Q(f)\) eingehen (vgl. Gl. 2.26), mit einer Zeitschrittweitensteuerung ausgestattet. Je nach Genauigkeitsanforderungen wird die Zeitschrittweite \(h\) zur Laufzeit der Berechnung angepasst. Dabei wird die Zeitschrittweite \(h\) derart variiert, dass ein ausreichend geringer Verfahrensfehler bei gleichzeitig effizienter Berechnung erzielt wird. Dadurch wird auch erreicht, dass alle Ortskurven \({z}_{i}\left({t}_{n}\right)\) des Modells in dem Stabilitätsgebiet \(G\) des numerischen Verfahrens verlaufen, da die Zeitschrittweite \(h\) den Abstand der Stabilitätsvariablen zum Ursprung \(\left|{z}_{i}\left({t}_{n}\right)\right|=h\left|{\lambda }_{i}\left({t}_{n}\right)\right|\) beeinflusst, und somit durch Verkleinerung von h die Ortskurven an den Ursprung herangezogen werden. Nachteil der Zeitschrittweitensteuerung ist, dass die Anzahl der erforderlichen Integrationsschritte innerhalb eines festen Zeitintervalls variiert, wodurch eine zeitdeterministische Berechnung in einer Echtzeitsimulation nicht mehr gewährleistet werden kann.

Die vorangegangenen Überlegungen zeigen, dass es sehr sinnvoll und häufig auch notwendig sein kann, das Integrationsverfahren für den Einsatz in einer Echtzeitsimulation gezielt auszuwählen. Somit sind Integrationsverfahren mit möglichst hoher Recheneffizienz bei ausreichender Stabilität (große Stabilitätsgebiete) und Genauigkeit (hohe Konsistenzordnung) gesucht. Als sehr recheneffiziente Integrationsverfahren sind die sogenannten linearen Mehrschrittverfahren (MSV) bekannt. Diese Verfahren nutzen für die Berechnung der Integralfläche \(Q\left(f\right)\) zusätzliche zurückliegende (bereits berechnete) Lösungsschritte \({x}_{n-1},{x}_{n-2},\dots\) und Funktionswerte \({f}_{n-1},{f}_{n-2},\dots\), um eine höhere Konsistenzordnung zu erreichen, und benötigen dadurch weniger Auswertungen der meist rechenaufwendigen Modellfunktion \(f\) pro Integrationsschritt im Vergleich zu anderen Verfahrensklassen mit gleicher Konsistenzordnung. Da bei komplexen Modellen der Rechenaufwand für die Auswertung der Modellfunktion \(f\) den Gesamtrechenaufwand für das Verfahren dominiert, ist es erforderlich die Anzahl der Auswertungen der Modellfunktionen möglichst gering zu halten. Die allgemeine Berechnungsvorschrift der MSV lautet:

Durch die Wahl der Koeffizienten \({\alpha }_{j}\) und \({\beta }_{j}\) sowie der Schrittanzahl \(k\) können verschiedene Typen von MSV mit unterschiedlichen Eigenschaften konstruiert werden. Bekannte Vertreter der MSV sind die Adams-Verfahren und die BDF-Verfahren (BDF: Backward Differentiation Formulas). Abb. 2.7 zeigt beispielhaft die Stabilitätsgebiete der expliziten Adams-Bashforth-Verfahren AB(\(k\)) und der impliziten Adams-Moulton-Verfahren AM(\(k\)) für \(k=1-5\).

Vereinfacht kann man sagen, dass die Konsistenzordnung mit der Schrittanzahl steigt und somit Verfahren mit hoher Schrittanzahl genauer sind als Verfahren mit kleiner Schrittanzahl. Die AB-Verfahren sind explizite MSV (\({\beta }_{0}=0\)) und dadurch sehr recheneffizient, da die Modellfunktion \({f}_{n}\) zum Zeitpunkt \({t}_{n}\) explizit ausgewertet werden kann und nur einmal pro Integrationsschritt berechnet werden muss, unabhängig von der Schrittanzahl des Verfahrens. Allerdings weisen diese Verfahren nur kleine Stabilitätsgebiete auf. Die AM-Verfahren sind implizite MSV (\({\beta }_{0}\ne 0\)) und besitzen sehr große Stabilitätsgebiete, sind jedoch sehr viel rechenaufwendiger als die AB-Verfahren. Dies liegt daran, dass die Modellfunktion \({f}_{n+1}\) zum Zeitpunkt \({t}_{n}\) noch nicht bekannt sein kann und eine explizite Auswertung nicht möglich ist. Die implizite Auswertung erfordert die numerische Lösung eines nichtlinearen algebraischen Gleichungssystems, was zahlreiche Aufrufe der Modellfunktion mit sich bringt.

Um eine bessere numerische Stabilität und Genauigkeit bei gleichzeitig hoher Recheneffizienz zu erzielen, kann eine Kombination beider Verfahren anhand der sogenannten Prädiktor-Korrektor-Methode (PECE) eingesetzt werden. Dabei wird durch einen vorgeschalteten expliziten Integrationszeitschritt mit einem AB-Verfahren der Zustandsvektor \({\tilde{x }}_{n+1}\) prädiziert (Predict), der zur Auswertung der Modellfunktion (Evaluate) \(f\left({\tilde{x }}_{n+1}\right)\) für den nachgeschalteten impliziten Integrationsschritt mit einem AM-Verfahren zur Korrektur genutzt werden kann (Correct). Eine erneute Auswertung der Modellfunktion (Evaluate) mit dem korrigierten Zustandsvektor \({x}_{n+1}\) liefert eine gute Näherung von \({f}_{n+1}\). Die Berechnungsvorschrift eines Adams-Bashforth-Moulton-PECE-Integrationsschritts lautet somit:

Dadurch wird die explizite Auswertung einer impliziten Berechnungsvorschrift möglich, wodurch eine Verbesserung der numerischen Stabilität und Genauigkeit gegenüber einem gewöhnlichen expliziten MSV erzielt wird [8, 10]. Allerdings wird die Stabilität des ursprünglichen impliziten MSV nicht erreicht. Abb. 2.8 zeigt die Stabilitätsgebiete des Adams-Bashforth-Moulton-PECE-Verfahrens (kurz ABM-Verfahren).

Um nun die Verbesserung gegenüber den AB-Verfahren einordnen zu können, sollte beachtet werden, dass die Konsistenzordnung bei den ABM-Verfahren \(r=k+1\) entspricht. D. h., neben der um eins höheren Konsistenzordnung und der damit verbundenen höheren Genauigkeit gegenüber den AB(\(k\))-Verfahren, zeigen die ABM(\(k\))-Verfahren auch bessere Stabilitätseigenschaften. Stellt man nämlich die Verfahren mit gleicher Konsistenzordnung \(r\) gegenüber, d. h. AB(\(r\)) jeweils mit ABM(\(r-1\)) für \(r\ge 2\), so fallen die Stabilitätsgebiete bei den ABM-Verfahren und die damit verbundene höhere Robustheit deutlich größer aus (vgl. Abb. 2.7 links und Abb. 2.8).

Die ABM-Verfahren benötigen zwei Auswertungen der Modellfunktion \({f}\) pro PECE-Integrationsschritt, unabhängig von der Konsistenzordnung. Das heißt, der Rechenaufwand ist in etwa doppelt so hoch wie bei den AB-Verfahren, jedoch um ein Vielfaches geringer wie bei den impliziten AM-Verfahren.

Durch mehrfaches Iterieren über den Evaluate-Correct-Schritt kann die Genauigkeit zusätzlich verbessert werden, jedoch kostet jeder Schritt eine weitere Auswertung der Modellfunktion. Diese Verfahren heißen dann P(EC)^ME-Verfahren mit M zusätzlichen Evaluate-Correct-Schritten.

Auf Basis der in diesem Beitrag beschriebenen Methoden wurde in den Jahren 2002–2005 im Rahmen verschiedener Forschungsprojekte am Institut für Steuerungstechnik der Werkzeugmaschinen und Fertigungseinrichtungen (ISW) an der Universität Stuttgart ein HiL-Simulator realisiert. Der Simulator basierte auf der in Abb. 2.2 beschriebenen Architektur. Der Echtzeit-Rechenkern wurde in einer ersten Version für die Windows Echtzeiterweiterung RTX entwickelt und später auf das TwinCAT Echtzeitsystem der Fa. Beckhoff Automation GmbH & Co.KG portiert. Der Simulator konnte selbstgetaktet asynchron und fremdgetaktet synchron mit Steuerungssystemen über den Feldbus gekoppelt werden. Als Busschnittstellen wurde Sercos-Interface und Profibus realisiert. Damit konnten HiL-Simulationen für Servoantriebe mit Bustakten unter 1 ms durchgeführt werden. In den Echtzeit-Rechenkern wurden die in diesem Beitrag beschriebenen Systementkopplung für lineare FE-Modelle implementiert. Für nichtlineare Modelle wurden Integrationsverfahren auf Basis der Adams-Mehrschrittverfahren mit PECE-Methode umgesetzt.

Die erste VIBN auf Basis einer HiLS am ISW wurde an einer real verfügbaren Forschungsmaschine mit neuartiger Hexapod-Parallelkinematik durchgeführt (Abb. 2.9). Das Modell umfasste ein nichtlineares Mehrkörpermodell der Maschinenkinematik gekoppelt mit 6 linearen Antriebs- und Reglermodellen.

Die Buskommunikation für die Antriebe erfolgte über Sercos-Interface. Als NC-Steuerung kam Beckhoff TwinCAT CNC zum Einsatz.

Die Verifizierung des Modells erfolgte anhand der gemessenen Stromverläufe aller 6 Antriebe an der realen Maschine. Zur Analyse und Optimierung wurde ein Testprogramm auf der realen Steuerung abgefahren und die dabei berechneten Zustandsgrößen der virtuellen Maschine farblich über den Bahnverlauf visualisiert. Damit konnten Problemzonen hinsichtlich der Vorschubgeschwindigkeit oder auftretender Antriebsstörkräfte auf dem Werkstück einfach lokalisiert, gezielt analysiert und durch Parameteranpassungen und Programmänderungen in der Steuerung optimiert werden [11].

Auf Basis des Prototypen wurde im Jahr 2005 damit begonnen die kommerzielle HiLS-Plattform ISG-virtuos von der Industriellen Steuerungstechnik GmbH in enger Kooperation mit dem ISW zu entwickeln. ISG-virtuos kam in zahlreiche Forschungs- und Industrieprojekten zum Einsatz [4, 12, 13]. Heute ist ISG-virtuos ein etabliertes Werkzeug für die VIBN von komplexen Steuerungssystemen.

Hoher Wettbewerbsdruck und kurze Produktentwicklungszyklen erfordern den Einsatz immer leistungsfähigerer Produktionsanlagen. Mit wachsenden Anforderungen an die Produktivität steigt die Komplexität dieser Anlagen. Die Beherrschbarkeit dieser Komplexität kann nur durch den Einsatz moderner Entwicklungs- und Testmethoden erreicht werden. Treten zu einem späten Zeitpunkt in der Entwicklung oder erst bei der Inbetriebnahme einer Anlage Fehler auf, bleiben zeitaufwendige und kostspielige Nacharbeiten meist nicht aus.

Das Ersetzen der realen Anlage durch eine virtuelle Anlage ermöglicht das Testen und Inbetriebnehmen von Steuerungssystemen bereits in einer frühen Phase des Entwicklungsprozesses einer Produktionsanlage. Dieser Beitrag stellt Methoden für die Virtuelle Inbetriebnahme mithilfe einer HiLS vor. Dabei wird beschrieben, wie eine verlustfreie und zeitsynchrone Ankopplung der Simulation an ein reales Steuerungssystem über den realen Feld-/Antriebsbus möglich ist. Zudem werden in diesem Beitrag Methoden für die Echtzeitberechnung von linearen und nichtlinearen Modellen mit Blick auf die numerischen Lösungsverfahren vorgestellt.