Solve a 2-dimensional optimal control problem via Riccati nonlinear equation

Question

Consider the 2-dimensional optimal control problem of the LQR kind $$ \min_u \int_0^\infty (x^T Q x + u^TRu) \, dt \quad\text{such that}\quad \begin{cases}\dot x(t) = Ax(t)+Bu(t) \\ x(0) = \begin{pmatrix}1\\-1\end{pmatrix}\end{cases} $$ with $x=\begin{pmatrix}x_1 \\ x_2\end{pmatrix}$, $u=\begin{pmatrix}u_1\\u_2\end{pmatrix}$, $Q=B=I$ (identity), $R=\gamma I$, $A=\begin{pmatrix}-\alpha & \alpha\\\beta & -\beta\end{pmatrix}$, $\alpha>\beta>0$ and $\gamma>0$.

Solve it using the stationary Riccati equation $$ 0 = Q+A^TS+SA-SBR^{-1}BS,\quad \text{with}\quad S=\begin{pmatrix}s_1 & s_2 \\ s_2 & s_3\end{pmatrix}. $$

In order to find $x$ and $u$ we have to combine the equation for $\dot x$ with the equation for the optimal control $u = -R^{-1}B^TSx = -\frac1\gamma Sx$. So, firstly, we have to find the expression for $S$.

The Riccati equation is simplified as $$ 0 = I+A^TS+SA-\frac1\gamma S^2 $$ which is equivalent to the following nonlinear system $$ \begin{cases} -\dfrac{s_1^2}{\gamma}-2\alpha s_1 -\dfrac{s_2^2}{\gamma}+2\beta s_2+1=0\\ \alpha s_1 - \alpha s_2 - \beta s_2 + \beta s_3 - \dfrac{s_1s_2}{\gamma}-\dfrac{s_2s_3}{\gamma}=0\\ -\dfrac{s_2^2}{\gamma}+2\alpha s_2 - \dfrac{s_3^2}{\gamma} - 2\beta s_3+1 = 0 \end{cases} $$

How to solve such a non-linear system to find expressions for $s_1, s_2$ and $s_3$? I think there is a fast way to solve either the system or directly the Riccati equation in matrix form, but I don't know how.

The solutions provided by Matlab, using the code below, are so long that it let me think is not the correct way to solve the problem

syms x y z a b g
eqn1 = 0 == -x^2/g-2*a*x-y^2/g+2*b*y+1;
eqn2 = 0 == a*x-a*y-b*y+b*z-x*y/g-y*z/g;
eqn3 = 0 == -y^2/g+2*a*y-z^2/g-2*b*z+1;
[x,y,z] = solve([eqn1, eqn2, eqn3], [x, y, z])

Answer 1

The solution of a Riccati equation can be found by determining the eigenvectors of the Hamiltonian; see e.g. here on Wikipedia. This is your best hope for a closed-form symbolic solution, in my view. In your case, the Hamiltonian is [ \begin{bmatrix} -\alpha & \alpha & -\gamma^{-1} & 0\\ \beta & -\beta & 0 & -\gamma^{-1}\\ -1 & 0 & \alpha & -\beta\\ 0 & -1 & -\alpha & \beta \end{bmatrix}. ] I don't see immediately how to diagonalize it, but maybe someone here has a more trained eye than me.