LabVIEW-Simulationsprogramme und Fallstudien zur modellprädiktiven Regelung

Krämer, Wolfgang Working Paper LabVIEW-Simulationsprogramme und Fallstudien zur modellprädiktiven Regelung Arbeitsberichte - Working Papers, No. 32 Provided in Cooperation with: Technische Hochschule Ingolstadt (THI) Suggested Citation: Krämer, Wolfgang (2014) : LabVIEW-Simulationsprogramme und Fallstudien zur modellprädiktiven Regelung, Arbeitsberichte - Working Papers, No. 32, Technische Hochschule Ingolstadt (THI), Ingolstadt, https://nbn-resolving.de/urn:nbn:de:bvb:573-5936 This Version is available at: https://hdl.handle.net/10419/202584 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc-nd/3.0/de/ Heft Nr. 32 aus der Reihe „ Arbeitsberichte – Working Papers“ ISSN 1612 -6483 Ingolstadt, im November 2014 Working Paper Abstract Es werden LabVIEW -Simulationsprogramme für das regelungstechnische Praktikum beschrieben, die intensiven Gebrauch von MATLAB -Skript-Nodes machen. Ein verfahrenstechnischer Prozess und ein Antriebssystem werden modellprädiktiv geregelt. Die Auswirkungen von Begrenzungen, die Teil der Regelungsaufgabe si nd, werden untersucht. Prof. Dr. Wolfgang Krämer LabVIEW -Simulationsprogramme und Fallstudien zur modell - prädiktiven Regelung 1 Vorwort When all is said and all is done, usually more is said than done.1 In diesem Bericht stelle ich die von mir während meines Forschungsfreisemesters an der University of Wisconsin – Madison zum Zweck meiner Weiterbildung durchgeführten Arbeiten und erzielten Ergebnisse vor. Meinen Fachkollegen an der THI bin ich zu tiefem Dank dafür verpflichtet, dass sie während meines Forschungsfreisemesters meine Lehrveranstaltungen übernommen und mir damit die Forschungsarbeiten in Madison ermöglicht haben. Ferner danke ich der Fakultät für Maschinenbau und der Hochschulleitung für die Gewährung des Freisemesters. Mein besonderer Dank gilt meinem Kollegen Prof. James B. Rawlings am Department of Chemical and Biological Engineering der University of Wisconsin – Madison für die freundliche Aufnahme in seine Forschungsgruppe, seine großzügige fachliche Unterstützung und Offenheit. Die Diskussionen mit ihm waren stets sehr interessant und hilfreich. Außerdem danke ich den Doktoranden (graduate students) Michael Risbeck und Nishith Patel für ihre Vorarbeiten, auf denen ich gut aufbauen konnte, und ihre wertvolle Unterstützung während meines Aufenthalts. Wie im Folgenden ausgeführt, erstellte ich Simulationsmodelle mit Hilfe der weit verbreiteten Systemdesignsoftware LabVIEW und konnte mich dabei in dieses komplexe Programmiersystem einarbeiten. Damit kann ich dieses zukünftig im Messtechnikpraktikum an der THI einsetzen. Außerdem führte ich Fallstudien zur modellprädiktiven Regelung durch, wobei ich die Programmiersprache MATLAB verwendete. Neben dem Kennenlernen dieser modernen Regelungsmethode konnte ich dabei meine Erfahrungen mit MATLAB vertiefen und neue Funktionen dieser weltweit eingesetzten Software kennenlernen. Somit war das Forschungssemester für meine zukünftige Arbeit an der THI äußerst lohnend. Ingolstadt, im November 2014, Wolfgang Krämer 1 Verfasser unbekannt. 2 Inhalt Vorwort ............................................................................................................... 1 1 Einleitung ......................................................................................................... 3 2 LabVIEW-Simulationsprogramme mit MATLAB-Skript-Nodes ......................... 4 2.1 Rechnerkonfiguration und Softwarewerkzeuge ......................................... 4 2.2 Beschreibung eines Simulationsprogramms ............................................. 7 3 Fallstudien zur modellprädiktiven Regelung (MPC) ....................................... 13 3.1 Einführung in lineare MPC ...................................................................... 13 3.2 Regelung eines Rührkesselreaktors ....................................................... 17 3.2 Regelung eines elastischen Antriebs ...................................................... 22 3.3.1 Lageregelung mit Servointegrierer ................................................... 23 3.3.2 Lageregelung mit Störgrößenschätzung ........................................... 27 3.3.3 Vergleich .......................................................................................... 31 3.4 Fazit ........................................................................................................ 31 Literaturhinweise .............................................................................................. 32 Anhang A: MATLAB-Funktionen für die Wassersäulen-Simulation .................. 33 Anhang B: Parameters for Standpipe Simulation ............................................. 35 3 1 Einleitung Dieser Bericht besteht aus zwei voneinander ziemlich unabhängigen Teilen: • Im nächsten Kapitel wird über die Entwicklung von EchtzeitSimulationsprogrammen für den Einsatz im regelungstechnischen Praktikum am Department of Chemical and Biological Engineering der University of Wisconsin - Madison berichtet. Dabei kommen die bekannte Systemdesignsoftware LabVIEW und die Programmiersprache MATLAB zum Einsatz. Insbesondere werden eine Programmstruktur und Regeln für die intensive Benutzung von MATLAB-Skript-Nodes für die Echtzeitberechnungen vorgestellt. • In Kapitel 3 werden Fallstudien zur modellprädiktiven Regelung (Model Predictive Control, MPC), einem modernen Verfahren der Regelungstechnik, durchgeführt. Als Regelstrecken werden je ein typisches Beispiel aus dem Maschinenbau bzw. der Antriebstechnik und aus der Verfahrenstechnik verwendet. Es wird insbesondere untersucht, wie sich die Möglichkeit, Begrenzungen bei der Formulierung der Regelungsaufgabe direkt zu berücksichtigen, auf das Regelverhalten auswirken. Diese Möglichkeit ist ein Spezifikum von MPC. 4 2 LabVIEW-Simulationsprogramme mit MATLAB-Skript-Nodes Um die zunehmende Zahl von Bachelorstudierenden (undergraduate students) am Department of Chemical and Biological Engineering der University of Wisconsin - Madison zu bewältigen, muss das Labor für das regelungstechnische Praktikum ausgebaut werden. Aus Raumund Kostengründen sollen jedoch keine weiteren Versuche aufgebaut werden, sondern es sollen Simulationen erstellt werden, die das Verhalten der realen Versuchsaufbauten möglichst realistisch nachbilden. Für folgende Versuchsaufbauten, die in [1] beschrieben sind, werden Simulationen erstellt: • Zwei-Tank-System bestehend aus zwei übereinander angeordneten, hintereinander geschalteten Behältern. Dieses Experiment wird für die Untersuchung von dynamischen Systemantworten verwendet. • Zwei hintereinander geschaltete Rührkesselreaktoren, in denen Malachitgrün mit Natronlauge reagiert. Der Versuchsaufbau dient der Anwendung von Identifikationsverfahren. • Hoher Wasserbehälter (Wassersäule) für die Entwicklung einer Füllstandregelung. • Kontinuierlicher beheizbarer Rührkessel für kombinierte Temperaturund Füllstandregelung (Mehrgrößenregelung). 2.1 Rechnerkonfiguration und Softwarewerkzeuge Bei den vorhandenen Praktikumsversuchen erfolgt die Datenerfassung, Steuerung und Regelung mit Hilfe von Rechnern, die LabVIEW-Programme ausführen. LabVIEW ist ein weit verbreitetes Programmiersystem der Firma National Instruments für Personal Computer zur Datenerfassung, Steuerung und Regelung, s. z.B. [2], [3]. LabVIEW-Programme werden als virtuelle Instrumente (Virtual Instruments, VIs) bezeichnet. Sie bestehen aus einer Benutzeroberfläche, Frontpanel genannt, und einem Blockschaltbild, in dem die Programmfunktionen auf grafische Weise in Form von Signalflussplänen programmiert werden. Bei den Praktikumsversuchen erfolgt die Erfassung von analogen Messsignalen und die Ausgabe von analogen Stellund Steuersignalen an die Versuchsaufbauten über ein Interface, das A/Dund D/A-Wandler enthält und mit dem 5 Rechner über eine serielle Schnittstelle kommuniziert. Bild 2.1 zeigt die messund regelungstechnische Konfiguration der vorhandenen Versuche. Bild 2.1: Messund regelungstechnische Konfiguration der vorhandenen Laborversuche Um eine hohe Realitätsnähe der neuen Versuche zu erreichen, bei denen das Verhalten des Versuchsaufbaus simuliert wird, sollen für die Datenerfassung, Steuerung und Regelung die gleichen VIs wie bei den realen Versuchsaufbauten verwendet werden. Außerdem sollen analoge Messund Stellsignale über das Interface eingelesen bzw. ausgegeben werden, um etwaige Effekte der A/Dund D/A-Wandlungen und der analogen Signalübertragung zu berücksichtigen. Damit bietet es sich an, einen zweiten Rechner für die Simulation des Versuchsaufbaus sowie ein zweites Interface zu verwenden. Die Konfiguration der neuen Versuche zeigt Bild 2.2. Bild 2.2: Rechnerkonfiguration der neuen Versuche Wie in Bild 2.2 angegeben, werden auch die Simulationsprogramme als LabVIEW-VIs erstellt, um die Anzahl der verwendeten Softwarewerkzeuge in Grenzen zu halten. Wie oben erwähnt, werden die VIs in Form von Blockschaltbzw. Signalflussbildern grafisch programmiert. Blockschaltbilder werden bei der Programmierung von komplexeren Systemen, wie sie die Modelle der Versuchsaufbauten darstellen, schnell unübersichtlich. Deshalb wird von der Möglichkeit Gebrauch gemacht, in LabVIEW sog. MATLAB-Skript-Nodes zu verwenden, in denen MATLAB-Skripte programmiert werden können. Diese Kombination von LabVIEW und MATLAB hat folgende Vorteile: • Ansprechende Benutzeroberflächen können in Form von LabVIEWFrontpanels einfach erstellt werden. • Es ergeben sich übersichtliche LabVIEW-Blockschaltbilder. Rechner mit LabVIEW-VI (Datenerfassung und Regelung) Interface (A/Du. D/AWandler) Realer Versuchsaufbau (Regelstrecke) Rechner mit LabVIEW-VI (Datenerfassung und Regelung ) Interface (A/Du. D/AWandler) Interface (A/Du. D/AWandler) Rechner mit LabVIEW-VI (Simulation des Versuchsaufbaus) 6 • Bewährte und wohlbekannte MATLAB-Funktionen, z.B. für numerische Integration, können verwendet werden • Messund Stellsignale werden auf beiden Rechnern mit Hilfe entsprechender LabVIEW-Funktionen und der Interfaces auf einheitliche Weise eingelesen bzw. ausgegeben. Außerdem gibt es in der Forschungsgruppe von Prof. Rawlings viel Knowhow zur MATLAB-Programmierung. Bei der Erstellung der Simulationsprogramme hat es sich als sinnvoll erwiesen, folgende Regeln zu beachten: • Für die Initialisierung und für die Durchführung der Simulation wird jeweils ein MATLAB-Skript-Node verwendet. • Berechnungen werden so weit wie möglich in den MATLAB-Skript-Nodes durchgeführt. Dies ergibt übersichtliche LabVIEW-Blockschaltbilder. • Die Simulationsberechnungen werden in einer MATLAB-Funktion (mFile) zusammengefasst. Im entsprechenden Skript-Node wird nur diese Funktion aufgerufen. Die MATLAB-Funktion kann damit zunächst ohne die Verwendung von LabVIEW entwickelt und getestet werden. • Eingangsgrößen der MATLAB-Funktion werden in Arrays zusammengefasst. Die MATLAB-Funktion liefert ihre Ausgangsgrößen zusammengefasst in Clustern, die bei Bedarf im LabVIEW-Blockschaltbild in ihre Komponenten zerlegt werden. Dies führt zu einem kompakten Aufruf der MATLAB-Funktion. Die letzten beiden Regeln erhöhen auch die Zuverlässigkeit der Simulationsprogramme. Wenn umfangreiche Berechnungen direkt in Skript-Nodes durchgeführt werden, treten immer wieder unverständliche Laufzeitfehler auf; ob dies ein generelles Problem oder auf die Programminstallation zurückzuführen ist, ist unklar. Abschließend sei noch bemerkt: • Die Simulationsberechnungen müssen selbstverständlich in Echtzeit durchgeführt werden. • Bei den Simulationen wird Messund Prozessrauschen modelliert, um eine möglichst hohe Realitätsnähe zu erreichen. 13 3 Fallstudien zur modellprädiktiven Regelung (MPC) 3.1 Einführung in lineare MPC Die modellprädiktive Regelung (Model Predictive Control, MPC) ist ein modernes Verfahren der Regelungstechnik, das insbesondere in der chemischen Industrie erfolgreich eingesetzt wird. Es wird im Folgenden kurz vorgestellt, soweit es für die Fallstudien erforderlich ist; die Darstellung folgt [4]. Bild 3.1 zeigt schematisch eine vollständige modellprädiktive Regelung mit dem Regler (regulator), der Zustandsschätzung (estimator), der Sollwertvorgabe (target selector) sowie der Regelstrecke (plant). Es handelt sich um eine zeitdiskrete Regelung. Bild 3.1: Modellprädiktive Regelung aus [4] Der Regler Ausgangspunkt für den Entwurf eines linearen MPC-Reglers ist die lineare zeitdiskrete Beschreibung der Regelstrecke in Zustandsraumbeschreibung )()()1( kuBkxAkx ⋅ + ⋅ = + , k = 0, 1, 2, … (3.1) und )()( kxCky ⋅ = (3.2) 14 mit der Anfangsbedingung x(0) = x 0 . Dabei sind wie üblich n Rx∈ der Zustandsvektor, m Ru∈ der Stellgrößenvektor, p Ry ∈ der Ausgangsgrößenvektor oder Vektor der gemessenen Größen, nn R A × ∈ die Systemmatrix, mn R B × ∈ die Eingangsmatrix und np RC × ∈ die Ausgangsmatrix. Aufgabe des Reglers ist, den Anfangszustand x(0) in den Endzustand null zu überführen, so dass das quadratische Kriterium [ ] )()'( 2 1 )()'()()'( 2 1 1 0 NxPNxkRukukQxkxV f N k ++= ∑ − = (3.3) mit den symmetrischen Gewichtsmatrizen Q ≥ 0 , R > 0 und P f ≥ 0 minimal wird; das Matrizenpaar (A, Q) muss stabilisierbar sein. N ist der zeitdiskrete Vorhersagebzw. Prädiktionshorizont. Zusätzlich sollen die Zustandsund Stellgrößenbegrenzungen 1,,0,)()( − = ≤ + NkdkGxkDu K (3.4) erfüllt werden. Durch diese Begrenzungen wird die Regelungsaufgabe nichtlinear. Die Minimierung des Kriteriums (3.3), so dass die Bedingungen (3.1) und (3.4) erfüllt werden, stellt ein sog. quadratisches Programm dar, das in jedem Abtastzeitschritt gelöst wird. Ohne Begrenzungen (3.4) ergibt die Lösung des quadratischen Programms eine lineare Zustandsrückführung (optimale Rückführung) )()( kxKku ⋅ − = (3.5) mit der Rückführmatrix ABRBBK )1('))1('( 1 Π+Π= − . (3.6) Dabei wird die Matrix Π (1) als Lösung der rückwärts-Riccati-Iteration AkBRBkBBkAAkAQk )('))('()(')(')1( 1 Π+ΠΠ−Π+=−Π − , (3.7) k = N, N-1, …, 2 mit der Endbedingung f PN =Π )( bestimmt. Für einen unendlichen Horizont, N → ∞, ergibt dies den wohlbekannten Riccati-Regler bzw. „Linear Quadratic Regulator (LQR)“. Man erhält den LQR auch, indem man die Matrix Π als Lösung der diskreten algebraischen Riccati-Gl. 0')'('' 1 =+Π+ΠΠ−Π−Π − QABRBBBAAA (3.8) bestimmt und anstatt Π(1) in Gl. (3.6) einsetzt. Dieselbe Rückführung ergibt sich bei der oben beschriebenen Iteration unabhängig vom Horizont N , wenn man als Gewichtsmatrix für die Endwerte der Zustandsgrößen P f = Π wählt. Dies wird 15 in [5] empfohlen, um die bekannten guten Eigenschaften des LQR zu erhalten, wenn keine Begrenzungen wirksam sind. Das quadratische Programm mit Begrenzungen (3.4) wird numerisch gelöst. Die Zustandsschätzung In der Regel sind die n Zustandsgrößen x unbekannt; nur die p < n Ausgangsgrößen y werden gemessen. Für die Rückführung (3.5) bzw. die Lösung des quadratischen Programms werden deshalb die Zustandsgrößen mit Hilfe eines linearen stationären Kalmanfilters bzw. eines „Linear Quadratic Estimators (LQE)“ geschätzt. Dabei wird davon ausgegangen, dass auf die Regelstrecke Prozessund Messstörungen w bzw. v wirken, so dass sie durch die Gln. )()()()1( kwkuBkxAkx + ⋅ + ⋅ = + , (3.9) )()()( kvkxCky + ⋅ = (3.10) beschrieben wird. Die Störungen werden als mittelwertfreies weißes Rauschen mit den Kovarianzmatrizen w QkwkwE =⋅ ))'()(( , v RkvkvE =⋅ ))'()(( (3.11) modelliert; w und v seien unkorreliert. Die Schätzung erfolgt in zwei Schritten. Im Korrekturschritt werden die Schätzwerte )( ˆ kx der Zustandsgrößen x zum Zeitpunkt k berechnet: ))( ˆ )(()( ˆ )( ˆkxCkyLkxkx −− −+= . (3.12) Dann erfolgt die Prädiktion der Schätzwerte für den nächsten Zeitpunkt k +1: )()( ˆ )1( ˆkBukxAkx +=+ − (3.13) Die Korrekturmatrix L ist 1 )'(' − += RCPCPCL , (3.14) mit der Lösung P der diskreten Riccati-Gl. 0')'('' 1 =++−− − QCPARCPCAPCPAPA . (3.15) Die Sollwertvorgabe Wie oben erwähnt hat der Regler die Aufgabe, die Zustandsgrößen x und nach Gl. (3.2) auch die Ausgangsgrößen y zu null zu machen. Falls die Sollwerte y sp ≠ 0 sind, werden Werte x s und u s für die Zustandsbzw. Stellgrößen berechnet, die im stationären Zustand y = y sp ergeben. Da in den Fallstudien nur die 16 Störgrößenausregelung betrachtet wird, d.h. y sp = 0 ist, wird auf diese Berechnung hier nicht eingegangen. Der Regler arbeitet dann mit den Abweichungen s xxx −= ~ und s uuu −= ~ (3.16) der Zustandsbzw. Stellgrößen anstatt mit x und u . Die stationären Zustandsund Stellgrößen x s bzw. u s hängen auch von eventuell auftretenden konstanten Störungen ab, wie im Folgenden erläutert wird. Berücksichtigung von Störgrößen Konstante nicht gemessene Störgrößen sollen zu keinen bleibenden Regelabweichungen führen. Solche Störgrößen werden mit zeitdiskreten Integratoren der Form )()()1( kwkdkd d +=+ , (3.17) modelliert, die durch weißes Rauschen w d angeregt werden. Es wird ein Störmodell der Ordnung n d angesetzt, d.h. d n Rd ∈ und d n Rw∈ . Die Gln. (3.17) werden mit den Zustandsund Ausgangsgln. für die Zustandsschätzung (3.9) bzw. (3.10) kombiniert, so dass sich die erweiterte Zustandsraumbeschreibung wu B d x I BA d x d +       +             =       + 00 , (3.18) [ ] v d x CCy d +       = (3.19) ergibt. Entsprechend diesem Vorgehen werden also – im Unterschied zu einer Störgrößenschätzung [6] – immer alle Zustandsund Störgrößen geschätzt. Damit das System detektierbar ist, müssen die Matrizen Bd und Cd so gewählt werden, dass d d d nn CC BAI +=      −− Rang (3.20) ist. Daraus ergibt sich, dass pn d ≤ (3.21) sein muss, d.h. die Anzahl der modellierten Störgrößen darf die Zahl der Messgrößen nicht überschreiten. Es ist offensichtlich, dass mit m Stellgrößen nicht mehr als m Sollwerte korrekt eingestellt werden können. Sollte m < p sein, werden aus den p Messgrößen nc ≤ m Regelgrößen 17 )()( kHCykr = , (3.22) c n R r ∈ ausgewählt, die stationär genau auf ihre Sollwerte rsp eingestellt werden sollen. Hierfür ist es nach [4] notwendig, dass p Störgrößen geschätzt werden, d.h. nd = p ist. Dabei wird keine physikalisch richtige Modellierung der Störungen vorausgesetzt; nur die Bedingung (3.20) muss erfüllt sein. Dies ist günstig, wenn die Wirkungsweise der Störungen unbekannt ist. Können jedoch die Störungen physikalisch richtig modelliert werden, ist auch mit Störmodellen niedrigerer Ordnung nd < p eine stationär genaue Störgrößenausregelung möglich. Dafür müssen in der erweiterten Zustandsraumbeschreibung (3.18) und (3.19) die richtigen Störeingangsund -ausgangsmatrizen Bd bzw. Cd verwendet werden. Mit einem Kalman-Filter für das erweiterte System (3.18), (3.19) können Schätzwerte x ˆ und d ˆ für die Zustandsund Störgrößen ermittelt werden. Die Werte d ˆ werden als Schätzung der stationären Störgrößen betrachtet, d.h. dd s ˆ ˆ = . Dann treten für n c = m mit den stationären Zustandsund Stellgrößen         −      −− =       − sd sd s s dHCr dB HC BAI u xˆ ˆ 0 sp 1 (3.23) keine bleibenden Abweichungen der Regelgrößen r auf. Mit den so bestimmten stationären Werten xs und us werden die Abweichungen (3.16) berechnet. 3.2 Regelung eines Rührkesselreaktors Als erste Fallstudie wird die modellprädiktive lineare Regelung eines kontinuierlichen Rührkesselreaktors nach [4] betrachtet. Beschreibung der Regelstrecke Eine schematische Darstellung des Reaktors, in dem eine irreversible Reaktion A → B erster Ordnung in einer Flüssigphase abläuft, zeigt Bild 3.2. Sein Verhalten wird durch folgendes nichtlineare Differentialgleichungssystem beschrieben:       − − − =RT E ck hr ccF dt dc exp )( 0 2 00 π , )( 2 exp )( 0 2 00 TT Cr U RT E ck C H hr TTF dt dT c pp −+      −∆− + − = ρρ π , 2 0 r FF dt dh π − =. 18 Bild 3.2: Schematische Darstellung des Rührkesselreaktors aus [4] Die Regelgrößen sind der Füllstand h und die molare Konzentration c des Reaktanden A. Die Reaktortemperatur T ist eine zusätzliche Zustandsgröße. Die Stellgrößen sind die Kühlmitteltemperatur T c und der Abflussvolumenstrom F. Der Zuflussstrom F 0 wirkt als Störung. Die Parameter sind in Tabelle 3.1 angegeben. Tabelle 3.1: Parameter des Rührkesselreaktors aus [4] Die Zustandsgrößen haben im stationären Zustand die Werte c s = 0,87778 kmol/l, T s = 324,5 K, h s = 0,659 m; die Stellgrößen sind T c,s = 300K, F s = 0,1 m 3 /min. Alle Zustandsgrößen werden gemessen. 19 Für die Reglerberechnungen wird ein linearisiertes Modell )()()()1( kdbkBukAxkx d ++=+ , )()( kCxky = für die Abweichungen           − − − == s s s hh TT cc yx ,       − − = s scc FF TT u , , s FFd ,00 −= bestimmt. Bei einer Abtastzeit von T 0 = 1 min sind die Modellmatrizen           − −− = 100 44,253279,0703,9 00728,000338,0,268.0 A ,           = 100 010 001 C ,           − − = 637,60 91,97297,1 1655,000537,0 B ,          − = 637,6 64,69 1175,0 d b . Als Gewichtsmatrizen werden ( ) 222 /1,/1,/1Diag sss hTcQ =, ( ) 22 /1,1/Diag sc,s FTR = gewählt. Regelung ohne Begrenzungen Für die Regelung wird ein linearer MPC nach den Gln. (3.5) – (3.7) mit den angegebenen Gewichtsmatrizen verwendet; es werden keine Begrenzungen berücksichtigt. Zunächst wird der Einfluss des Vorhersagehorizonts untersucht, indem die Eigenwerte des geschlossenen Regelkreises für verschiedene Horizonte N berechnet werden. Die Gewichtsmatrix der Endwerte wird der Einfachheit halber Pf = Q gesetzt; diese Festlegung wird in der Praxis häufig verwendet. Sowohl die Regelstrecke als auch der geschlossene Regelkreis haben einen reellen Eigenwert und ein konjugiert komplexes Eigenwertpaar. Bild 3.2.a zeigt alle Eigenwerte in der komplexen Zahlenebene, Bild 3.2b zeigt den reellen Eigenwert in Abhängigkeit von N . Man erkennt, dass die Eigenwerte sehr schnell gegen ihre Grenzwerte für N →∞ konvergieren, so dass schon mit kleinen Horizonten gutes Regelkreisverhalten zu erwarten ist. 20 -0.5 0 0.5 1 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Eigenwerte in der reellen Zahlenebene Realteil Imaginärteil 0 2 4 6 8 10 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 Reeller Eigenwert Horizont N Eigenwert Bild 3.2: Eigenwerte des geschlossenen Regelkreises in Abhängigkeit des Horizonts N . (×: Eigenwerte der Regelstrecke; ○: Eigenwerte für N →∞, blau: Einheitskreis) Zur Untersuchung des Regelverhaltens wird eine Simulation aus [4] nachvollzogen. Dabei wird für die Reglerberechnungen eine in der Forschungsgruppe neu geschriebene MATLAB-Funktion lmpc_matlab verwendet und getestet, die für die Lösung des oben beschriebenen quadratischen Programms die Funktion quadprog aus der Optimization Toolbox verwendet. Als äußere Erregung des Regelkreises wird von einer nicht gemessenen Störung in Form einer Erhöhung des Zuflusses F0 um zehn Prozent zum Zeitpunkt t = 10 min ausgegangen. Der Horizont des MPC ist N = 5 , als Gewichtsmatrix für den Endzustand wird die Lösung der stationären Riccati-Gl. (3.8), d.h Pf = Π, verwendet. Ein Störmodell (3.17) der Ordnung nd = 3 dient der Vermeidung bleibender Regelabweichungen. Damit ist ein Zustandsschätzer erforderlich, obwohl alle Zustandsgrößen der Regelstrecke gemessen werden. Als Störeinund -ausgangsmatrix werden für das erweiterte Modell (3.18)           − = 637,600 91,9700 1655,000 d B bzw.           = 010 000 001 d C gewählt. Die Kovarianzmatrizen sind )1,,,,,(Diag εεεεε = w Q und ( ) 2227 ,,Diag105,2 sssv hTcR − ⋅= mit ε = 10-8 ; weitere Details s. [4]. Damit ergibt sich das gleiche Simulationsergebnis wie in [4], womit bestätigt ist, dass die neue Funktion lmpc_matlab für Regelaufgaben ohne Begrenzungen korrekt arbeitet. Nun wird die Störung auf 25 % von F0 erhöht. Alle Parameter der Regelung bleiben unverändert. Es wird verglichen, wie sich Regelkreise mit dem nichtli- 21 nearen und mit dem linearisierten Modell der Regelstrecke verhalten. Bild 3.3 zeigt das Ergebnis der Simulationen, links den Verlauf der Zustandsgrößen und rechts die Stellgrößenverläufe. Man erkennt, dass die Nichtlinearität der Regelstrecke einen deutlichen Einfluss auf die Konzentration c sowie auf die Stellgröße Kühlmitteltemperatur Tc hat. Bei den anderen Größen gibt es geringere Unterschiede zwischen der nichtlinearen und der linearen Simulation. In beiden Fällen haben die Regelgrößen Konzentration c und Füllstand h keine bleibenden Abweichungen, d.h. die Störung wird vollständig ausgeregelt, da die Bedingung (3.20) erfüllt und nd = p = 3 ist. 0 10 20 30 40 50 0.84 0.86 0.88 c (kmol/m 3 ) 0 10 20 30 40 50 324 326 328 330 T (K) 0 10 20 30 40 50 0.7 0.8 0.9 1 h (m) Time nichtlinear linearisiert 0 10 20 30 40 50 298 299 300 301 302 Tc (K) 0 10 20 30 40 50 0.1 0.12 0.14 0.16 F (m 3 /min) Time nichtlinear linearisiert Bild 3.3: Verhalten der Regelung ohne Begrenzung bei einer Störung des Zuflusses um 25 % mit nichtlinearem und linearisiertem Streckenmodell, Zeit in min Regelung mit Begrenzungen Die Regelgrößen c und h weichen während des in Bild 3.3 gezeigten Einschwingvorgangs zeitweise erheblich von ihren stationären bzw. Sollwerten ab. Um diese Abweichungen zu verringern, wird versucht, Zustandsgrößenbegrenzungen bei den Reglerberechnungen zu berücksichtigen. Hierfür werden die Ungleichungen (3.4) verwendet; die Funktion lmpc_matlab wird etwas erweitert. Die Simulationen erfolgen mit dem nichtlinearen Regelstreckenmodell. Zunächst wird eine Begrenzung c ≥ cmin = 0,86 kmol/m3 für die Konzentration vorgegeben. Dies hat keinerlei Auswirkungen auf das Regelkreisverhalten. Der Grund dafür ist, dass der Regler mit den geschätzten Zustandsgrößen arbeitet 22 und die Schätzwert c ˆ der Konzentration deutlich über dem wahren Wert liegt. Auch größere Werte cmin < cs führen zu keinem anderen Verhalten. Sodann wird versucht, die starke Überhöhung des Füllstands durch Vorgabe einer Begrenzung h ≤ hmax = 0,7 m zu verringern. Dadurch verändert sich der Maximalwert des Füllstands jedoch nicht. Der maximale Füllstand wird zur Zeit t = 11 min erreicht, d.h. unmittelbar nach dem Eintreten der Störung; er ist die erste Regelabweichung, die auftritt. Dieser Wert kann vom Regler nicht beeinflusst werden, da der Regler erst auf diese Abweichung reagiert. Für t > 11 min ergeben sich mit der Begrenzung allerdings niedrigere Füllstände. Um den maximalen Füllstand zu reduzieren, muss die Abtastzeit verringert werden. Die Abweichungen der Stellgrößen von ihren stationären Werten sind gering, so dass eine Stellgrößenbegrenzung nicht sinnvoll erscheint. 3.2 Regelung eines elastischen Antriebs Die Lageregelung elastischer Antriebe stellt eine typische Aufgabenstellung der Antriebstechnik dar [6], [7]. Bild 3.4 zeigt einen Zweimassenschwinger als einfachstes Modell eines ungedämpften elastischen Antriebs. Bild 3.4: Modell eines elastischen Antriebs aus [7] Bei den Positionen von Motor und Last y1 bzw. y2 kann es sich um translatorische oder rotatorische Freiheitsgrade handeln. Dementsprechend sind m1 und m2 die Massen oder Trägheitsmomente von Motor bzw. Last. Die Feder mit der Steifigkeit c beschreibt eine elastische Verbindung von Motor und Last, z.B. durch ein Getriebe. Die Stellgröße u stellt ein Antriebsmoment oder eine Antriebskraft dar, die Störgröße d ist ein Lastmoment oder eine Lastkraft. Das Verhalten des Antriebs wird durch die Bewegungsgln. dbubKyyM dbb , +=+ && mit der regulären Massenund der singulären Steifigkeitsmatrix       = 2 1 0 0 m m M bzw.       − − =cc cc K d 29 0 5 10 15 20 25 30 -3 -2 -1 0 1 2 3 t[s] Positionen y Motor Last 0 5 10 15 20 25 30 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 t[s] Stellgröße u Bild 3.9: Störsprungantwort der Regelung ohne Begrenzungen Regelung mit Begrenzungen Zunächst wird wieder die Stellgröße betragsmäßig auf umax = 1,2 begrenzt. Bild 3.10 zeigt das sich damit ergebende Regelverhalten. Die Begrenzung führt zu keiner Vergrößerung der Einbruchtiefe und nur zu einer geringfügig verlängerten Ausregelzeit. Das Regelverhalten ist besser als mit dem Servointegtierer mit Stellgrößenbegrenzung, vgl. Bild 3.7. 0 5 10 15 20 25 30 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 t [s] Positionen y Motor Last 0 5 10 15 20 25 30 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 t [s] Stellgröße u Bild 3.10: Störsprungantwort der Regelung mit │ u │ ≤ 1,2 30 Es ist möglich, die Stellgröße noch stärker zu begrenzen, ohne dass der Regelkreis instabil wird. Dies ist aber aus praktischer Sicht irrelevant, da sich dann ein viel zu schwach gedämpftes Regelkreisverhalten ergibt. Nun wird die Regelgröße auf │ r │ ≤ rmax mit rmax = 2 begrenzt. Wie Bild 3.11 zeigt, ist dies möglich, führt aber natürlich zu größeren Ausschlägen der Stellgröße und zu heftigeren Motorbewegungen. Die Regelgröße kann noch stärker begrenzt werden, womit selbstverständlich auch die negativen Auswirkungen zunehmen. So ergibt sich z.B. für rmax = 0,5 betragsmäßig eine maximale Stellgröße │ u │ max ≈ 20 und eine maximale Motorposition │ y1 │ max ≈ 4 . Realistischer Weise kann die Regelgröße also nur in Maßen begrenzt werden. 0 5 10 15 20 25 30 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 t [s] Positionen y Motor Last 0 5 10 15 20 25 30 -2.5 -2 -1.5 -1 -0.5 0 0.5 t [s] Stellgröße u Bild 3.11: Störsprungantwort der Regelung mit │ r │ ≤ 2 Es ist auch möglich, sowohl die Regelals auch die Stellgröße zu begrenzen. So ergibt sich z.B. für rmax = 2 , umax = 1,5 oder für rmax = 2,25 , umax = 1,2 jeweils ordentliches Regelverhalten; die Begrenzungen werden eingehalten. Wird allerdings ein Wert der genannten Kombinationen von rmax und umax weiter verringert, wird das quadratische Programm (3.1)-(3.4) unlösbar. Da in der Praxis die maximalen Störgrößen häufig nicht genau bekannt sind, sollte man bei der gleichzeitigen Begrenzung der Stellund der Regelgröße sehr vorsichtig sein, um dieses Problem zu vermeiden. 31 3.3.3 Vergleich Vergleicht man die Regelungen mit Servointegrierer und mit Störgrößenschätzung, so scheint für die Regelung des elastischen Antriebs die Störgrößenschätzung vorteilhaft zu sein. Stellgrößenbeschränkungen führen zu einer geringeren Beeinträchtigung des Regelverhaltens, und die Stellgröße kann stärker beschränkt werden, ohne dass der Regelkreis instabil wird. Außerdem kann mit der Störgrößenschätzung auch die Regelgröße beschränkt werden. Selbst eine gleichzeitige Beschränkung der Regelund der Stellgröße ist möglich, wobei jedoch wie oben erwähnt vorsichtig vorgegangen werden sollte. 3.4 Fazit Lineare MPC ist eine attraktive Methode der Regelungstechnik. Ihre Besonderheit besteht darin, dass Begrenzungen von Stellund Zustandsgrößen explizit vorgegeben und berücksichtigt werden können. Wie die Fallstudien zeigen, ist jedoch bei der Formulierung von Begrenzungen, insbesondere von Zustandsgrößenbegrenzungen Vorsicht geboten, da zu restriktive Vorgaben dazu führen können, dass der Regelkreis instabil wird oder dass das quadratische Programm, mit welchem die Stellgrößen berechnet werden, unlösbar wird. Unter Umständen können Zustandsgrößenbeschränkungen unwirksam sein. Die Fallstudien zeigen ebenfalls, dass das Konvergenzverhalten der RiccatiIteration (3.7) stark unterschiedlich sein kann. Dieses Verhalten kann Hinweise auf die Wahl des Vorhersagehorizonts N geben. Auf neuere Weiterentwicklungen von MPC, insbesondere nichtlineare MPC [8], verteilte MPC [9], hybride MPC und ökonomische MPC [10], bei der sich das Optimierungskriterium nicht auf technische sondern auf wirtschaftliche Größen, wie z.B. den Ertrag eine Chemieanlage, bezieht, wird in diesem Bericht nicht eingegangen. 32 Literaturhinweise [1] Rawlings, J. B.; M. D. Graham; W. H. Ray: Process Dynamics and Control Laboratory Manual CBE 470, Second Edition. Department of Chemical and Biological Engineering; University of Wisconsin – Madison, 1999. [2] Systemdesignsoftware NI LabVIEW. www.ni.com/labview/d, aufgerufen am 10.10.2014. [3] Georgi, W.; E. Metin,.: Einführung in LabVIEW, 5. Auflage. München, 2012. [4] Rawlings, J. B.; D. Q. Maine: Model Predictive Control: Theory and Design. Nob Hill Publishing, Madison, WI, 2009. [5] Rawlings, J. B.; D. Q. Maine: Postface to “Model Predictive Control: Theory and Design”. Nob Hill Publishing, Madison, WI, 2009. [6] Juen, G.; Lageregelung elastischer Antriebe dargestellt an einem Radioteleskop. VDI-Verlag, Düsseldorf, 1987. [7] Krämer, W.: Lageregelung elastischer Antriebe durch Ausgangsrückführungen. VDI-Verlag, Düsseldorf, 1991. [8] Angeli, D.; J. B. Rawlings: Receding Horizon Cost Estimation and Control for Nonlinear Plants. 8th IFAC Symposium on Nonlinear Control Systems (NOLCOS), Bologna, Italy, September 2010. [9] Brett, T. S.; S. J. Wright; J. B. Rawlings: Cooperative Distributed Model Predictive Control for Nonlinear Systems. J. Proc. Cont. 21, 2011, S. 698-704. [10] Rawlings, J. B.; D. Angeli; C. Bates: Fundamentals of Economic Model Predictive Control. IEEE Conference on Decision and Control (CDC), S. 3851-3861, Maui, HI, December 2012. 33 Anhang A: MATLAB-Funktionen für die WassersäulenSimulation Funktion Standpipe_calc function outputs = Standpipe_calc(inputs, outputs, dt, param, … options) % Inputs v_Valve = inputs(1); % Valve voltage normOut = inputs(2); % Normal outflow valve open / closed distOut = inputs(3); % Disturbance outflow valve open / closed much_noise = inputs(4); % Much measurement noise due to splashing % Steady state inflow rate qIn_ss = max([0, (param.vIntercept + param.vSlope * v_Valve) * … .qIn_Max]); % Noise of inflow rate if qIn_ss == 0 noise_q = 0; else noise_q = param.stdev_q * randn(1); end % Cd*A for outflow param.Cda = (normOut + distOut) * param.CdA; % Initial conditions qIn = outputs(1); h = outputs(2); % Solve DEs [t,x] = ode45(@(t, x) Standpipe_DE(t, x, qIn_ss, noise_q, param), ... [0 dt/2 dt], [qIn, h], options); % Determine Outputs outputs(1) = x(3,1); % outputs(1) = qIn if x(3,2) >= param.hMax outputs(2) = param.hMax; % outputs(2) = h outputs(3) = 1; % outputs(3) = overflow else outputs(2) = x(3,2); outputs(3) = 0; end outputs(4) = param.sIntercept + param.sSlope*outputs(2) … + (param.stdev_v_lo + much_noise*(param.stdev_v_hi … - param.stdev_v_lo))*randn(1); % senor voltage end % of function Standpipe_calc 34 Funktion Standpipe_DE function dxdt = Standpipe_DE( t, x, qIn_ss, noise_q, param) % %Calculation of derivatives for standpipe simulation % % x(1): Inflow rate (qIn) % x(2): Level in standpipe (h) dxdt = zeros(2,1); % column vector % Outflow rate qOut = param.Cda * sqrt(2 * param.g * max([x(2), 0])); % Derivatives dxdt(1) = (qIn_ss - x(1)) / param.T1_Valve; dxdt(2) = (max([x(1)+noise_q, 0]) - qOut)/param.A; end 35 Anhang B: Parameters for Standpipe Simulation Remark: For all calculations in the model, the units cm and s are used. For display, the quantities are converted to US-units, if necessary. The diameter of the standpipe d = 4.5 cm was measured at the lab setup. Therefore, the cross sectional area of the standpipe is A = 15.5 cm2 . According to the manual, the height of the standpipe is hMax = 150 cm . Also according to the manual, there are the following correlations for the valve and sensor voltages: • h = 0 cm → VSensor = 3 V • h = 150 cm → VSensor = 10 V • VValve = 1 V → Valve completely closed • VValve = 5 V → Valve completely open Therefore, the following parameters hold: • for the sensor: sIntercept = 3 V , sSlope = 0.0467 V/cm • for the relative opening of the valve vIntercept = -0.25 , vSlope = 0.25 / V The product Cd · AV of the flow rate coefficient and the cross sectional area of the outflow valve is determined from the lab report of X. Teng, I. Tonner, and W. Han from April 17, 2013. From this, s cm 1285.0 2= A gAC Vd , c.f. p. 7. Therefore, Cd · AV = 0.045 cm2 . This value is used for both, the nominal outflow and the step disturbance valve, as the valves are equal. The little height difference of the valves is neglected. (It could easily be accounted for in the simulation, if desired). In afore mentioned report, the following regression results are given: • Level at steady state: h = 93 cm/V · VValve - 203 cm (figure 1) • Flow rate at steady state: q = 0.072 gpm/V · VSensor + 0.0929 gpm (figure 3) • For the sensor voltage, as above: VSensor = 0.0467 V/cm · h + 3 V From these data, the maximum inflow rate, if the valve is fully open, is calculated to be qin,Max = 75.1 cm3/s . The time lag of the valve is assumed to be T1Valve = 0.5 s . For the noises, the following standard deviations were chosen: • sq = 1 cm3/s for the inflow rate • sV,lo = 0.1 V for the sensor voltage in normal operation • sV,hi = 0.5 V for the sensor voltage with inflow splashing, c.f. section 10 of part D in the manual. Heft Nr. 32 aus der Reihe „ Arbeitsberichte – Working Papers“ ISSN 1612 -6483 Ingolstadt, im November 2014 Working Paper Impressum Herausgeber Der Präsident der Technischen Hochschule Ingolstadt Esplanade 10, 85049 Ingolstadt Telefon: +49 841 9348-0 Fax: +49 841 9348-2000 E -Mail: [email protected] Druck Hausdruck Die Beiträge aus der Reihe „Arbeitsberichte – Working Papers“ erscheinen in unregel mäßigen Abständen. Alle Rechte, insbesondere das Recht der Vervielfältigung und Verbreitung sowie der Übersetzung vorbehalten. Nachdruck, auch auszugsweise, ist gegen Quellenangabe gestattet, Belegexemplar erbeten. Internet Alle Themen aus der Reihe „Arbeitsberichte – Working Papers“, können Sie unter der Adresse www.thi.de nachlesen. Prof. Dr. Wolfgang Krämer LabVIEW -Simulationsprogramme und Fallstudien zur modell - prädiktiven Regelung