1
$\begingroup$

Given the Algebratic Riccati Equation (ARE) $$A^T X + XA + XRX + Q = 0$$ where $A,R,Q \in \mathbb R^{n \times n}$, we are interested in the matrix $X$ that solves this equation. If we define the $2n \times 2n$ Hamiltonian matrix $$H = \begin{bmatrix} A & R \\ -Q & -A^T\end{bmatrix}$$ and let $X = X_2 X_1^{-1}$, then $H$ satisfies the equation $$H\begin{bmatrix} X_1 \\ X_2\end{bmatrix} = \begin{bmatrix} X_1 \\ X_2\end{bmatrix} \Lambda$$ which represents an eigenvalue problem. That is, if $v_i \in \mathbb R^{2n}$ is the $i$th column of $\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$, and $\lambda_i$ is the $i$th corresponding eigenvalue, then $H$ satisfies the equation $$Hv_i = \lambda_i v_i$$ I'm aware that there are better ways of solving this ARE, but I'm interested in solving this eigenvalue problem to obtain $X_1$ and $X_2$, and then $X$. However, when computing the eigenvalues and eigenvectors of $H$ using, for example, NumPy, I obtain a $2n \times 2n$ matrix of eigenvectors. I'm not sure how to relate this matrix to $\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$, which is $2n \times n$.

$\endgroup$
4
  • $\begingroup$ The algebraic Riccati usually admits several solutions, and most likely you are interested in finding the so-called stabilizing solution. That is the eigenvalues of the closed-loop matrix $A-BB^T X$ all have negative real-part. If this is the case, You use the eigenvectors of the Hamiltonian corresponding to the eigenvalues with negative real-part as well. Stack all these eigenvectors together and you get the desired $2n\times n$ matrix. $\endgroup$
    – DerZwirbel
    Nov 11, 2022 at 10:09
  • $\begingroup$ @DerZwirbel ahh I see. If I’m not necessarily interested in a stabilizing solution, does this mean that any combination of $n$ eigenvectors of $H$, out of a total of $2n$, can be used to construct $\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$, and then from this obtain $X$? $\endgroup$
    – mhdadk
    Nov 11, 2022 at 11:24
  • $\begingroup$ Roughly Speaking Yes. But It might be the case that $X_1$ is singular for some combinations. But anyway, If $X_1$ is nonsingular, then $X_2 X_1^{-1}$ solves the algebraic Riccati equation. $\endgroup$
    – DerZwirbel
    Nov 11, 2022 at 12:34
  • $\begingroup$ @DerZwirbel Please post your contribution as an answer; comments are only meant to suggest improvements to posts. $\endgroup$ Nov 11, 2022 at 15:04

1 Answer 1

0
$\begingroup$

For reference, here is some Python code implementing the suggestions by @DerZwirbel. The following ARE is solved in this case: $$ XA + A^T X - XBR^{-1}B^T X + Q = 0 $$

import numpy as np
import scipy.linalg

# generate random test matrices (Q and R need to be positive semi-definite and positive definite respectively)
A = 3 * np.random.randn(2,2)
B = 3 * np.random.randn(2,2)
Q = 3 * np.random.randn(2,2)
Q = Q.T @ Q
R = 3 * np.random.randn(2,2)
R = R.T @ R

# compute the Hamiltonian matrix, its eigenvalues, and its eigenvectors
H = np.block([[A,-B @ np.linalg.solve(R,B.T)],[-Q,-A.T]])
eigs,eigvs = np.linalg.eig(H)

# find out which eigenvalues are negative
stable_idxs = np.where(np.real(eigs)<0)[0]

# split the array of eigenvectors into two separate blocks to create the X_1 and X_2 matrices,
# such that X = X_2 * X_1^{-1}
H_eigvs_top,H_eigvs_bottom = np.split(eigvs,2,axis=0)
X1 = H_eigvs_top[:,stable_idxs]
X2 = H_eigvs_bottom[:,stable_idxs]
X_hat = X2 @ np.linalg.inv(X1)

# as a sanity check, compute the solution using SciPy's implementation and compare
# this solution to the one we obtained manually
X = scipy.linalg.solve_continuous_are(A,B,Q,R,balanced=True)
np.testing.assert_allclose(actual=X_hat,desired=X,err_msg="X and X_hat do not match.")

# evaluate the value of the ARE using X and X_hat (both arrays should be close to 0)
ARE = A.T @ X + X @ A - X @ B @ np.linalg.solve(R,B.T) @ X + Q
ARE_hat = A.T @ X_hat + X_hat @ A - X_hat @ B @ np.linalg.solve(R,B.T) @ X_hat + Q

print(f"Eigenvalues of Hamiltonian matrix:\n{eigs}\n")
print(f"X1:\n{X1}\n")
print(f"X2:\n{X2}\n")
print(f"Ground truth solution X:\n{X}\n")
print(f"Solution X using manual method:\n{X_hat}\n")
print(f"Ground truth value of ARE (should be close to 0):\n{ARE}\n")
print(f"Value of ARE using manual method (should be close to 0):\n{ARE_hat}\n")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.