import numpy as np
import scipy.linalg as linalg
import scipy.signal as signal
import control as ctDie Regler-Evaluierung (policy evaluation) ist ein zentraler Schritt für viele Methoden aus dem bestärkenden Lernen (reinforcement learning) und der Regler-Suche (policy search). In diesem Notebook wollen wir die Regler-Evaluierung für ein LQR Problem mit zufälligen Anfangszuständen untersuchen. Aufbauend auf diesen Ergebnissen werden wir im nächsten Blogeintrag policy search-Methoden besprechen.
“Spur Trick” und Erwartungswert
In vielen wissenschaftlichen Aufsätzen werden die Spur einer Matrix und der Erwartungswert gemeinsam verwendet, ohne jedoch diese Schritte genau zu beschreiben. Dem geübten Studieren sind diese Beziehungen oft bekannt, Anfänger stehen jedoch diesen Umformungen oft ratlos gegenüber.
Spur Trick für Quadratische Formen
Starten wir mit dem Matrixprodukt \(x^TMx\), wobei \(x\) ein \(n \times 1\)-Spaltenvektor und \(M\) eine \(n \times n\) Matrix ist. Das Matrixprodukt
\[ x^TMx \]
ergibt dann einen skalaren Wert, oder anderes ausgedrückt, eine \(1 \times 1\)-Matrix. Für die Spur einer Matrixmultiplikation
\[ tr(ABC) = tr(CAB) = tr(BCA) \]
gilt die Invarianz der Spur unter zyklischen Vertauschungen. Da es sich bei \(x^TMx\) um einen skalaren Wert handelt, können wir die Spur \(tr()\) anwenden, ohne den Wert zu verändern. Weiters können wir die zyklische Vertauschung durchführen und erhalten
\[ \begin{split} x^TMx &= tr(x^TMx) \\ &= tr(xx^TM) \\ &= tr(Mxx^T). \end{split} \]
Erwartungswert / Spur
Nun wollen wir annehmen, dass es sich bei \(x\) um eine Realisierung eines Zufallsvektors \(X\) handelt. Es soll der Erwartungswert
\[ \mathbb{E}_{x \sim X}(x^TMx) \]
ausgewertet werden. Dabei erweist sich der oben angeführt Spur Trick als hilfreich:
\[ \begin{split} \mathbb{E}(x^TMx) &= \mathbb{E}(tr(x^TMx)) \\ &= \mathbb{E}(tr(xx^TM)) = tr(\mathbb{E}(xx^T) M ) \\ &= \mathbb{E}(tr(Mxx^T)) = tr(M \mathbb{E}(xx^T)) \end{split} \]
Der Erwartungswert \(\mathbb{E}(xx^T)\) kann durch
\[ \mathbb{E}(xx^T) = Cov(x,x) -\mathbb{E}(x)\mathbb{E}(x)^T \]
ausgedrückt werden, wobei \(Cov(x,x)\) die Kovarianz und \(\mathbb{E}(x)\) der Erwartungswert des Zufallsvektors ist.
Wenn wir nun annehmen, dass \(x\) aus einer Normalverteilung \(\mathbb{N}(0,\Sigma_0)\) gezogen wird, dann ergibt sich für die Auswertung des Erwartungswerts
\[ \mathbb{E}(x^TMx) = tr(M \Sigma_0). \]
Diese Beziehung wird sich im Folgenden als sehr nützlich erweisen.
Regler-Suche / Optimierung ( engl. policy search / optimization)
Die Methoden der Regler-Suche (engl. policy search aka policy optimization) basieren auf der Idee, eine Reglerstruktur mit Parametern vorzugeben. Optimiert wird dann direkt über diese Reglerparameter. Als Reglerstukturen kommen z.B. PID-, Youla-, Matrix- oder Neuronale-Parametrierung in Frage. Uns stehen eine Vielzahl von Methoden für diese Art der Reglersuche zur Verfügung, sowohl modellbasierte wie modellfreie Spielarten.
LQR / \(\mathcal{H}_2\) Zustandsrückführung
Für Regelungstechniker ist es wie immer hilfreich mit einem Problem zu starten, welches wir sehr gut verstehen. Daher bietet sich wieder das LQR Problem an. Die lineare Kontrolltheorie hat sich mittlerweile zu einem Benschmarkproblem für das bestärkende Lernen entwickelt. Hier konnten in den letzten Jahren viele theoretische Einsichten gewonnen werden.
Gegeben sei ein LQR Problem, wobei wir jedoch die Regelstruktur mit \(u_k=-Kx_k\) vorgeben. Zusätzlich sei für den Anfangszustand \(x_0\) eine Gauß-Verteilung \(\mathbb{N}(0,\Sigma_0)\) angenommen. Nun lautet das LQR Problem:
\[ C(K) = \underset{K}{\min} \mathbb{E} \left( \sum_{k=0}^{\infty} {x_k^TQx_k + u_k^TRu_k} \right) \]
\[ \begin{split} \text{subject to} \quad & x_{k+1} = A x_k + B u_k \\ & u = -K x_k \\ & x_0 \sim \mathbb{N}(0,\Sigma_0) \end{split} \]
Dieses Problem wird auch als optimales \(\mathcal{H}_2\)-Kontrollproblem bezeichnet, da die zu erwartenden Kosten, hier gleich der \(\mathcal{H}_2\)-Norm des linearen Systems vom Störeingang zum virtuellen Ausgang sind (hier \(z = \begin{bmatrix} Q^{1/2}x \\ R^{1/2}u \end{bmatrix}\)). Der Störeingang ist hier nur ein Impuls, der die Anfangsbedingungen erzeugt.
Wir wollen im Weitern zwei Methoden für die Regler-Evaluierung besprechen und in Python implementieren.
Regler-Evaluierung über die Bellman-Gleichung
Die Regler-Evaluierung kann nun mit Hilfe der Lyapunov Gleichung (Bellman-Gleichung)
\[ P_K = (A-BK)^TP_K(A-BK) + Q + K^TRK \]
durchgeführt werden. Es gilt die Beziehung
\[ C(K) = \mathbb{E}_{x_0 \sim \mathbb{N}(0,\Sigma_0)} x_0^T P_K x_0, \]
welche durch das Anwenden der obigen Regeln, durch
\[ C(K) = tr(P\Sigma_0) \]
ausgewertet wird.
Algorithmus 1
- Gegebei sei
\[ A,B,Q,R,K,\Sigma_0 \]
- Löse die Lyapunov-Gleichung (Bellman-Gleichung)
\[ P_K = (A-BK)^TP_K(A-BK) + Q + K^TRK \]
- Regler-Evaluierung
\[ C(K) = tr(P\Sigma_0) \]
Regler-Evaluierung über die Zustands-Korrelationsmatrix
Wir setzen hier den Regler \(u_k = -Kx_k\) in \(x_k^TQx_k+u_k^TRu_k\) ein, und erhalten durch die Eigenschaften der Matrixspur:
\[ x_k^TQx_k+x_k^TK^TRKx_k = x_k^T(Q+K^TRK)x_k = tr((Q+K^TRK)x_kx_k^T) \]
Nun können wir durch die Linearität der Spur, der Summe und des Erwartungswert die Kostenfunktion umformen, und erhalten
\[ C(K) = \mathbb{E} \left( \sum_{k=0}^{\infty}\left[ tr((Q+K^TRK)x_kx_k^T) \right] \right) = tr\left((Q+K^TRK) \sum_{k=0}^{\infty} \mathbb{E}\left[x_kx_k^T \right]\right). \]
Für den zeitlichen Verlauf des Zustands \(x_k\), gilt der Zusammenhang
\[ x_k = (A-BK)^k x_0, \]
welcher in den Erwartungswert \(\mathbb{E}(x_kx_k^T)\) eingesetzt, die Gleichung
\[ \mathbb{E}(x_kx_k^T) = (A-BK)^k \mathbb{E}(x_0x_0^T) ((A-BK)^T)^k = (A-BK)^k \Sigma_0 ((A-BK)^T)^k \]
erzeugt. Rückeingesetzt in die Summe, wird die Zustands-Korrelationsmatrix
\[ \Sigma_K = \sum_{k=0}^{\infty} \mathbb{E}(x_kx_k^T) = \sum_{k=0}^{\infty} (A-BK)^k \Sigma_0 ((A-BK)^T)^k \]
erzeugt. Diese Matrix kann für einen stabilisierenden Regler \(K\) durch die Lyapunov-Gleichung
\[ \Sigma_K = (A-BK)\Sigma_K(A-BK)^T + \Sigma_0 \]
gelöst werden (Vergleiche dazu die Gramsche-Matrix controllability gramian). Die Regler-Evaluierung lautet nun
\[ C(K) = tr((Q+K^TRK)\Sigma_K). \]
Algorithmus 2
- Gegebei sei
\[ A,B,Q,R,K,\Sigma_0 \]
- Löse die Lyapunov-Gleichung (Zustands-Korrelationsmatrix)
\[ \Sigma_K = (A-BK)\Sigma_K(A-BK)^T + \Sigma_0 \]
- Regler-Evaluierung
\[ C(K) = tr((Q+K^TRK)\Sigma_K) \]
Code
sys_d = ct.drss(2) # create a stable time discrete system
A = sys_d.A
B = sys_d.B
print(np.abs(linalg.eigvals(A)) < 1.) # check if system is stable[ True True]Q = np.eye(2)*5
R = np.eye(1)*10def policy_evaluation_P(A,B,Q,R,K,S):
""" Policy Evaluation via Lyapunov Equation (Bellman Equation)
(A-B K)^T P_K (A-B K) - P_K + Q + K^TRK = 0
Returns: Value Functions (Matrix)
"""
Ac = A-B@K
Qc = Q+K.T@R@K
P_K = linalg.solve_discrete_lyapunov(Ac.T,Qc)
C_K = np.trace(P_K@S)
return C_K, P_Kdef policy_evaluation_S(A,B,Q,R,K,S):
""" Policy Evaluation via Lyapunov Equation ()
(A-B K) S_K (A-B K)^T - S_K + S_0 = 0
Returns: Value Functions (Matrix)
"""
Ac = A-B@K
S_K = linalg.solve_discrete_lyapunov(Ac,S)
C_K = np.trace((Q+K.T@R@K)@S_K)
return C_K, S_K# a stable policy is needed, design by pole-placement
Ac_eig = np.sort(np.array([0.5,0.6]))
fsf1 = signal.place_poles(A, B, Ac_eig)
K = fsf1.gain_matrix
# covarinace of initial state
Sigma_0 = np.eye(2)*5
# check eigenvalues of closed-loop
print(np.abs(Ac_eig))
print(np.sort(np.abs(linalg.eigvals(A-B@K))))[0.5 0.6]
[0.5 0.6]# algo 1
C_K_algo1, P_K = policy_evaluation_P(A,B,Q,R,K,Sigma_0)
# algo 2
C_K_algo2, Sigma_K = policy_evaluation_S(A,B,Q,R,K,Sigma_0)print(C_K_algo1)
print(C_K_algo2)
print(np.allclose(C_K_algo1, C_K_algo2))837.4845833988895
837.4845833988902
TrueWir sehen also, dass beide Algorithmen zum gleichen Ergebnis führen.
Fazit
Wir haben hier zwei Methoden der Regler-Evaluierung besprochen, welche wir im nächsten Blogeintrag für eine Regler-Suche verwenden werden. Diese Methoden sind derzeit im Reinforcement Learning und im Robot Learning sehr beliebt.