Optimierung des LQR-Steuerungssystems

Einführung


Auf Habr wurden mehrere Artikel [1,2,3] veröffentlicht, die sich direkt oder indirekt auf dieses Thema beziehen. In diesem Zusammenhang ist die Veröffentlichung [1] mit dem Titel „Mathematik an den Fingern: Linear-Quadratischer Regler“ zu übersehen, die im Volksmund das Funktionsprinzip des optimalen LQR-Reglers erklärt.

Ich wollte dieses Thema fortsetzen, nachdem ich die praktische Anwendung der dynamischen Optimierungsmethode untersucht hatte, aber bereits an einem konkreten Beispiel mit Python. Zunächst einige Worte zur Terminologie und Methode der dynamischen Optimierung.

Optimierungsmethoden werden in statische und dynamische unterteilt. Das Steuerobjekt befindet sich unter dem Einfluss verschiedener externer und interner Faktoren in einem kontinuierlichen Bewegungszustand. Daher wird die Auswertung des Regelergebnisses für die Regelzeit T gegeben, und dies ist die Aufgabe der dynamischen Optimierung.

Mit Methoden der dynamischen Optimierung werden Probleme im Zusammenhang mit der Verteilung begrenzter Ressourcen über einen bestimmten Zeitraum gelöst und die Zielfunktion als integrale Funktion geschrieben.

Der mathematische Apparat zur Lösung solcher Probleme sind Variationsmethoden: klassische Variationsrechnung, das Prinzip des maximalen L.S. Pontryagin und dynamische Programmierung R. Bellman.

Die Analyse und Synthese von Steuerungssystemen erfolgt im Zeit-, Frequenz- und Zustandsraum. Die Analyse und Synthese von Steuerungssystemen im Zustandsraum wird in den Lehrplan aufgenommen. Die in den Schulungsunterlagen unter Verwendung des SQR-Controllers vorgestellten Techniken sind jedoch für Matlab konzipiert und enthalten keine praktischen Analysebeispiele.

Der Zweck dieser Veröffentlichung ist es, die Methoden zur Analyse und Synthese linearer Steuerungssysteme im Zustandsraum anhand des bekannten Beispiels zur Optimierung des Steuerungssystems eines elektrischen Antriebs und eines Flugzeugs mithilfe der Programmiersprache Python zu betrachten.

Mathematische Begründung der dynamischen Optimierungsmethode


Zur Optimierung werden die Kriterien Amplitude (MO), Symmetrie (CO) und Kompromissoptimum (KO) verwendet.

Bei der Lösung von Optimierungsproblemen im Zustandsraum wird ein lineares stationäres System durch Vektormatrixgleichungen gegeben:

 dotx= fracdxdt=A cdotx+B cdotu;y=C cdotx, (1)

Integrale Kriterien für minimalen Steuerenergieverbrauch und maximale Geschwindigkeit werden von den Funktionalen festgelegt:

Jx= int infty0(xTQx+uTRu+2xTNu)dt rightarrowmin, (2)

Ju= int infty0(yTQy+uTRu+2yTNu)dt rightarrowmin. (3)

Das Steuergesetz u liegt in Form einer linearen Rückkopplung auf die Zustandsvariablen x oder auf die Ausgangsvariablen y vor:

u= pmKx,u= pmKy

Eine solche Kontrolle minimiert die quadratischen Qualitätskriterien (2), (3). In den Beziehungen (1) ÷ (3) sind Q und R symmetrische positive definitive Matrizen der Dimension [n × n] bzw. [m × m]; K ist eine Matrix konstanter Dimensionskoeffizienten [m × n], deren Werte nicht begrenzt sind. Wenn der Eingabeparameter N weggelassen wird, wird angenommen, dass er Null ist.

Die Lösung für dieses Problem, das als Problem der linearen integralen quadratischen Optimierung (LQR-Optimierung) bezeichnet wird, im Zustandsraum wird durch den Ausdruck bestimmt:

u=R1BTPx

wobei die Matrix P die Ricatti-Gleichung erfüllen muss:

ATP+PA+PBR1BTP+Q=0

Die Kriterien (2), (3) sind ebenfalls widersprüchlich, da die Implementierung des ersten Terms eine Quelle unendlich hoher Leistung und für den zweiten eine Quelle unendlich niedriger Leistung erfordert. Dies kann durch den folgenden Ausdruck erklärt werden:

Jx= int infty0xTQxdt ,

Welches ist die Norm Mod(x) Vektor x, d.h. ein Maß für seine Schwingung im Regelungsprozess und nimmt daher bei schnellen Transienten mit weniger Schwingung kleinere Werte an, und der Ausdruck:

Ju= int infty0uTRudt

ist ein Maß für die zur Steuerung verwendete Energiemenge und eine Strafe für die Energiekosten des Systems.

Die Gewichtsmatrizen Q, R und N hängen von den Bedingungen der entsprechenden Koordinaten ab. Wenn ein Element dieser Matrizen gleich Null ist, unterliegt die entsprechende Koordinate keinen Einschränkungen. In der Praxis erfolgt die Auswahl der Werte der Matrizen Q, R, N nach der Methode der Expertenschätzung, des Versuchs, des Fehlers und hängt von der Erfahrung und dem Wissen des Konstrukteurs ab.

Um solche Probleme zu lösen, wurden die folgenden MATLAB-Operatoren verwendet:
 beginbmatrixR,S,E endbmatrix=lqr(A,B,Q,R,N) und  beginbmatrixR,S,E endbmatrix=lqry(Ps,Q,R,N) die die Funktionen (2), (3) durch den Zustandsvektor x oder durch den Ausgangsvektor y minimieren.

Verwaltungsobjektmodell Ps=ss(A,B,C,D). Das Ergebnis der Berechnung ist die Matrix K der optimalen Rückkopplungen in Bezug auf die Zustandsvariablen x, die Lösung der Ricatti-Gleichung S und die Eigenwerte E = eiq (A-BK) des Regelungssystems.

Funktionskomponenten:

Jx=x0TP1x0,Ju=x0TP2x0

wobei x0 der Vektor der Anfangsbedingungen ist; P1 und P2 - unbekannte Matrizen, die eine Lösung der Lyapunov-Matrixgleichungen sind. Sie werden mithilfe von Funktionen gefunden. P1=lyap(NN,Q) und P2=lyap(NN,KTRK) , NN=(A+BK)T

Merkmale der Implementierung der dynamischen Optimierungsmethode mit Python


Obwohl die Python Control Systems Library [4] folgende Funktionen hat: control.lqr, control.lyap, ist die Verwendung von control.lqr nur möglich, wenn das Slycot-Modul installiert ist, was ein Problem darstellt [5]. Bei Verwendung der Lyap-Funktion im Kontext einer Aufgabe führt die Optimierung zu einer control.exception.ControlArgument: Q muss ein symmetrischer Matrixfehler sein [6].

Um die Funktionen lqr () und lyap () zu implementieren, habe ich scipy.linalg verwendet:

from numpy import* from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E P1=solve_lyapunov(NN,Ct*Q*C) 

Übrigens sollten Sie die Kontrolle nicht ganz aufgeben, da die Funktionen step (), pol (), ss (), tf (), feedback (), acker (), place () und andere gut funktionieren.

Ein Beispiel für die LQR-Optimierung in einem elektrischen Antrieb

(Ein Beispiel stammt aus der Publikation [7])

Betrachten Sie die Synthese eines linear-quadratischen Reglers, der das Kriterium (2) für ein im Zustandsraum durch Matrizen definiertes Steuerobjekt erfüllt:

A = \ begin {bmatrix} & -100 & 0 & 0 \\ & 143.678 & -16.667 & -195.402 \\ & 0 & 1.046 & 0 \ end {bmatrix}; B = \ begin {bmatrix} 2300 \\ 0 \\ 0 \\\ end {bmatrix}; C = \ begin {bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \ end {bmatrix}; D = 0



Folgendes wird als Zustandsvariable betrachtet: x1 - Wandlerspannung V; x2 - Motorstrom A; x3 - Winkelgeschwindigkeit c1 . Dies ist das TP-DPT-System mit HB: Motor R nom = 30 kW; Unom = 220 V; Inom = 147 A ;;  omegaω 0 = 169 c1 ;;  omegaω max = 187 c1 ;; Nennwiderstandsmoment Mnom = 150 N * m; Vielzahl des Anlaufstroms = 2; Thyristorwandler: Unom = 230 V; Uy = 10 B; Inom = 300 A; Vielzahl von kurzzeitigen Überströmen = 1,2.

Bei der Lösung des Problems akzeptieren wir die Diagonale der Matrix Q. Als Ergebnis der Modellierung wurde festgestellt, dass die Minimalwerte der Matrixelemente R = 84 und Q=[[0,01,0,0],[0,0,01,0],[0,0,0,01]]. In diesem Fall wird ein monotoner Übergangsprozess der Winkelgeschwindigkeit des Motors beobachtet (Abb. 1). Bei R = 840 Q=[[0,01,0,0],[0,0,01,0],[0,0,0,01]] Der Übergangsprozess (Abb. 2) erfüllt das MO-Kriterium. Die Berechnung der Matrizen P1 und P2 erfolgte bei x0 = [220 147 162].

Auflistung des Programms (Abb. 1).
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[84]] Q=matrix([[0.01, 0, 0],[0, 0.01, 0],[0, 0, 0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0, 0.01,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



Abb. 1

Auflistung des Programms (Abb. 2).
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[840]] Q=matrix([[0.01, 0, 0],[0, 0.01, 0],[0, 0, 0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0, 0.01,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



Abb. 2
Durch geeignete Auswahl der Matrizen R und Q ist es möglich, den Anlaufstrom des Motors auf akzeptable Werte von (2–2,5) In (Abb. 3) zu reduzieren. Zum Beispiel mit R = 840 und
Q = ([[[0,01,0,0]], [0,0,88,0], [0,0,0,01]], sein Wert beträgt 292 A und die Übergangsprozesszeit unter diesen Bedingungen beträgt 1,57 s.

Auflistung des Programms (Abb. 3).
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[840]] Q=matrix([[0.01,0,0],[0,0.88,0],[0,0,0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0,0.88,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



Abb.3

In allen betrachteten Fällen sind die Rückkopplungen von Spannung und Strom negativ und die Geschwindigkeit positiv, was gemäß den Stabilitätsanforderungen unerwünscht ist. Darüber hinaus ist das synthetisierte System für die Aufgabe astatisch und für die Last statisch. Daher betrachten wir die Synthese eines PI-Reglers im Zustandsraum mit einer zusätzlichen Zustandsvariablen x4 - Übertragungskoeffizient des Integrators.

Wir stellen die Ausgangsinformationen in Form von Matrizen dar:

A = \ begin {bmatrix} & -100 & 0 & 0 & 0 \\ & 143.678 & -16.667 & -195.402 & 0 \\ & 0 & 1.046 & 0 & 0 \\ & 0 & 0 & 1 & 0 \ end {bmatrix}; B = \ begin {bmatrix} 2300 \\ 0 \\ 0 \\ 0 \\\ end {bmatrix}; C = Auge (4)); D = 0

Aufgabentransienten, die dem MO-Kriterium entsprechen, wurden bei R = 100, q11 = q22 = q33 = 0,001 und q44 = 200 erhalten. Fig. 4 zeigt die Transienten der Zustandsvariablen, die den Astatismus des Systems nach Aufgabe und nach Last bestätigen.

Programmliste (Abb. 4).
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0,0],[143.678, -16.667,-195.402,0],[0,1.046,0,0] ,[0,0,1,0]]); B=matrix([[2300], [0], [0],[0]]); C=matrix([[1, 0, 0,0],[0, 1, 0,0],[0, 0, 1,0],[0, 0, 0,1]]); D=matrix([[0],[0],[0],[0]]); #    R=[[100]]; Q=matrix([[0.001, 0, 0,0],[0, 0.001, 0,0],[0, 0, 0.01,0],[0, 0,0,200]]); #  LQR-  (K,S,E)=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162],[0]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure(num=None, figsize=(9, 7), dpi=120, facecolor='w', edgecolor='k') plt.suptitle('      LQR -',size=14) ax1 = fig.add_subplot(411) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('out(1)') ax1.grid(True) ax2 = fig.add_subplot(412) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('out(2)') ax2.grid(True) ax3 = fig.add_subplot(413) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('out(3)') ax3.grid(True) ax4 = fig.add_subplot(414) y,x=step(w[3,0]) ax4.plot(x,y,"b") plt.ylabel('out(4)') ax4.grid(True) plt.xlabel(', .') plt.show() 



Abb. 4
Um die Matrix K zu bestimmen, hat die Python-Steuerungssystembibliothek zwei Funktionen: K = acker (A, B, s) und K = Ort (A, B, s), wobei s der Vektor ist - die Zeile der gewünschten Pole der Übertragungsfunktion des Regelungssystems. Der erste Befehl kann nur für Systeme mit einem Eingang in u für n ≤ 5 verwendet werden. Die zweite hat keine solchen Einschränkungen, aber die Vielzahl der Pole sollte den Rang der Matrix B nicht überschreiten. Ein Beispiel für die Verwendung des Operators acker () ist in der Auflistung und in (Abb. 5) dargestellt.

Programmliste (Abb. 5).
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt i=(-1)**0.5 A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) p=[[-9.71+14.97*i -9.71-14.97*i -15.39 -99.72]]; k=acker(A,B,p) H=AB*k; Wss=ss(H,B,C,D); step(Wss) w=tf(Wss) fig = plt.figure() plt.suptitle('   acker()') ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") ax3.grid(True) plt.xlabel(', .') plt.show() 




Abb. 5

Fazit


Die Implementierung der in der Veröffentlichung angegebenen LQR-Optimierung des elektrischen Antriebs unter Berücksichtigung der Verwendung neuer Funktionen lqr () und lyap () ermöglicht es uns, die Verwendung des Lizenzprogramms MATLAB aufzugeben, wenn wir den entsprechenden Abschnitt der Steuerungstheorie studieren.

Ein Beispiel für die LQR-Optimierung in einem Flugzeug (LA)

(Ein Beispiel stammt aus der Veröffentlichung von Astrom und Mruray, Kapitel 5 [8] und wurde vom Autor des Artikels fertiggestellt, um die Funktion lqr () zu ersetzen und die Terminologie auf allgemein akzeptiert zu reduzieren.)

Der theoretische Teil, eine kurze Theorie, die LQR-Optimierung, wurde bereits oben betrachtet. Fahren wir also mit der Code-Analyse mit den entsprechenden Kommentaren fort:

Erforderliche Bibliotheken und LQR-Controller-Funktion.

 from scipy.linalg import* # scipy   lqr from numpy import * # numpy    from matplotlib.pyplot import * #     from control.matlab import * #  control..ss()  control..step() #     lqr() m = 4; #    . J = 0.0475; #        ***2 #́-       #  . r = 0.25; #     . g = 9.8; #     /**2 c = 0.05; #   () def lqr(A,B,Q,R): X =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*X)) E= eig(AB*K)[0] return X, K, E 

Anfangsbedingungen und Grundmatrizen A, B, C, D Modelle.

 xe = [0, 0, 0, 0, 0, 0];#   x   () ue = [0, m*g]; #   #   A = matrix( [[ 0, 0, 0, 1, 0, 0], [ 0, 0, 0, 0, 1, 0], [ 0, 0, 0, 0, 0, 1], [ 0, 0, (-ue[0]*sin(xe[2]) - ue[1]*cos(xe[2]))/m, -c/m, 0, 0], [ 0, 0, (ue[0]*cos(xe[2]) - ue[1]*sin(xe[2]))/m, 0, -c/m, 0], [ 0, 0, 0, 0, 0, 0 ]]) #   B = matrix( [[0, 0], [0, 0], [0, 0], [cos(xe[2])/m, -sin(xe[2])/m], [sin(xe[2])/m, cos(xe[2])/m], [r/J, 0]]) #   C = matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]) D = matrix([[0, 0], [0, 0]]) 

Wir bauen Ein- und Ausgänge, die schrittweisen Änderungen der xy-Position entsprechen. Die Vektoren xd und yd entsprechen dem gewünschten stationären Zustand im System. Die Matrizen Cx und Cy sind die entsprechenden Ausgaben des Modells. Um die Dynamik des Systems zu beschreiben, verwenden wir ein System von Vektormatrixgleichungen:

\ left \ {\ begin {matrix} xdot = Ax + Bu => xdot = (A-BK) x + xd \\ u = -K (x-xd) \\ y = Cx \\ \ end {matrix} \ richtig.



Die Dynamik einer geschlossenen Schleife kann mit der Funktion step () für K \ cdot xd und K \ cdot xd als Eingabevektoren modelliert werden, wobei:

 xd = matrix([[1], [0], [0], [0], [0], [0]]); yd = matrix([[0], [1], [0], [0], [0], [0]]); 

Die aktuelle Steuerungsbibliothek 0.8.1 unterstützt nur einen Teil der SISO Matlab-Funktionen zum Erstellen von Rückkopplungsreglern. Um Daten vom Regler zu lesen, müssen daher Indexvektoren lat (), alt () für die laterale -x- und vertikale -y-Dynamik erstellt werden.

 lat = (0,2,3,5); alt = (1,4); #    Ax = (A[lat, :])[:, lat]; Bx = B[lat, 0]; Cx = C[0, lat]; Dx = D[0, 0]; Ay = (A[alt, :])[:, alt]; By = B[alt, 1]; Cy = C[1, alt]; Dy = D[1, 1]; #      K Qx1 = diag([1, 1, 1, 1, 1, 1]); Qu1a = diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1a) K1a = matrix(K); #      x H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx); (Yx, Tx) = step(H1ax, T=linspace(0,10,100)); #     y H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy); (Yy, Ty) = step(H1ay, T=linspace(0,10,100)); 

Die Grafiken werden nacheinander angezeigt, eine für jede Aufgabe.

 figure(); title("   x  y") plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); plot([0, 10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 

Grafik:



Der Einfluss der Amplitude der Eingabeaktionen
 #     Qu1a = diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1a) K1a = matrix(K); H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx); Qu1b = (40**2)*diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1b) K1b = matrix(K); H1bx = ss(Ax - Bx*K1b[0,lat], Bx*K1b[0,lat]*xd[lat,:],Cx, Dx); Qu1c = (200**2)*diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1c) K1c = matrix(K); H1cx = ss(Ax - Bx*K1c[0,lat], Bx*K1c[0,lat]*xd[lat,:],Cx, Dx); figure(); [Y1, T1] = step(H1ax, T=linspace(0,10,100)); [Y2, T2] = step(H1bx, T=linspace(0,10,100)); [Y3, T3] = step(H1cx, T=linspace(0,10,100)); plot(T1.T, Y1.T, 'r-'); plot(T2.T, Y2.T, 'b-'); plot(T3.T, Y3.T, 'g-'); plot([0 ,10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); title("   ") legend(('1', '$1,6\cdot 10^{3}$','$4\cdot 10^{4}$'), loc='lower right'); text(5.3, 0.4, 'rho'); ylabel(' '); xlabel(''); grid() show() 


Grafik:



Einschwingverhalten
 #  Qx    Qx2 = CT * C; Qu2 = 0.1 * diag([1, 1]); X,K,E =lqr(A, B, Qx2, Qu2) K2 = matrix(K); H2x = ss(Ax - Bx*K2[0,lat], Bx*K2[0,lat]*xd[lat,:], Cx, Dx); H2y = ss(Ay - By*K2[1,alt], By*K2[1,alt]*yd[alt,:], Cy, Dy); figure(); [Y2x, T2x] = step(H2x, T=linspace(0,10,100)); [Y2y, T2y] = step(H2y, T=linspace(0,10,100)); plot(T2x.T, Y2x.T, T2y.T, Y2y.T) plot([0, 10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); title(" ") ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 


Grafik:



Physikalisch fundiertes Einschwingverhalten
 #    #       #   1 .     10 . Qx3 = diag([100, 10, 2*pi/5, 0, 0, 0]); Qu3 = 0.1 * diag([1, 10]); X,K,E =lqr(A, B, Qx3, Qu3) K3 = matrix(K); H3x = ss(Ax - Bx*K3[0,lat], Bx*K3[0,lat]*xd[lat,:], Cx, Dx); H3y = ss(Ay - By*K3[1,alt], By*K3[1,alt]*yd[alt,:], Cy, Dy); figure(); [Y3x, T3x] = step(H3x, T=linspace(0,10,100)); [Y3y, T3y] = step(H3y, T=linspace(0,10,100)); plot(T3x.T, Y3x.T, T3y.T, Y3y.T) axis([0, 10, -0.1, 1.4]); title("   ") ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 


Grafik:



Fazit:


Die Implementierung der LQR-Optimierung in dem in der Veröffentlichung vorgestellten Flugzeug unter Berücksichtigung der Verwendung der neuen Funktion lqr () ermöglicht es uns, die Verwendung des Lizenzprogramms MATLAB beim Studium des entsprechenden Abschnitts der Steuerungstheorie aufzugeben.

Referenzen


1. Mathematik an den Fingern: linear-quadratischer Regler.

2. Der Zustandsraum bei den Problemen des Entwurfs optimaler Steuerungssysteme.

3. Verwenden der Python Control Systems Library zum Entwerfen automatischer Steuerungssysteme.

4. Python Control Systems Library.

5. Python - Das Slycot-Modul kann nicht in Spyder (RuntimeError & ImportError) importiert werden.

6. Python Control Systems Library.

7. Kriterien für eine optimale Steuerung und LQR-Optimierung in einem elektrischen Antrieb.

Source: https://habr.com/ru/post/de441582/


All Articles