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

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])


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.
• Thanks, I did not know about this method. Calling H the Hamiltonian matrix, the diagonalization can be performed in matlab with [V,D] = eig(H), where V is the matrix of eigenvectors and D the diagonal matrix of eigenvalues. However, also in this case the values computed are very long. – sound wave Mar 23 '20 at 18:21
• Using V and D, how can I construct the matrices $U_1$ and $U_2$ cited in the page you linked? In particular, see i.imgur.com/WbfKw7I.png – sound wave Mar 23 '20 at 22:31
• $\begin{bmatrix}U_1 \\ U_2\end{bmatrix}\in\mathbb{C}^{2n\times n}$ ($n=2$ in your case) is obtained by concatenating horizontally the $n$ eigenvectors that are associated to eigenvalues in the left half-plane. In (succinct) Matlab, U = V(:, real(diag(D)) < 0). If you have eigenvalues with real part zero, then it is trickier (but that shouldn't happen in your case). – Federico Poloni Mar 24 '20 at 7:40
• If the eigenvalues are complex, $U$ (as we constructed it) is complex. If the eigenvalues are real, $U$ is real. That wiki page is assuming that $H$ may have complex eigenvalues, but not imaginary eigenvalues (that's two different things). – Federico Poloni Mar 25 '20 at 7:30