First off, let me say this subject is something I am not familiar with, and only started working on it a few days ago, so sorry in advance if I say something really stupid.
I have been trying to find a numerical routine that allows me to solve the matrix equation $ A^{T}L + LA + L(B_{1} + B_{2})L + Q = 0$, where $Q$ and the $B$'s are symmetric and real valued.
I have tried converting it to the standard CARE, which is $A^{T}L + LA - L(BR^{-1}B^{T})L + Q = 0$, by eigen decomposition and schur decomposition. The decomposition works well, but using a numerical routine equivalent to python's solve_continuous_are doesn't return symmetric $L$ as a solution, because the matrix $R^{-1}$ that is obtained is negative, so the decomposition into the Hamiltonian matrix doesn't quite work out I guess.
I have tried finding iterative methods, but I only found one that solves a system of the form $ A^{T}L + LA + -L(BB^{T} - CC^{T})L + QQ^{T} = 0$.
I haven't tried implementing this one yet, since my original $B_{1}$ and $B_{2}$ are real and symmetric, so I cannot decompose them to obtain real $B$ and $C$, for obvious reasons (I guess?).
Maybe the solution to this is actually quite simpler than what I'm thinking and I'm just burned out, but this has had me stuck for the last 3 days, I have been working tens of hours on this and cannot find any solution, so any input would be immensely appreciated.
The CARE can be used to solve for the full state feedback gain of the LQR problem, which solves (omitting cross terms for now)
$$ \begin{array}{rl} \min & \int_0^\infty x^\top Q\, x + u^\top R\, u\, dt \\ s.t. & \dot{x} = A\,x + B\, u \end{array} $$
where the pair $(A,B)$ has to be stabilizable, $Q\geq0$ (positive semi definite) and $R>0$ (positive definite). There might actually be a solution even if $(A,B)$ is not stabilizable, namely when $Q$ does not 'observe' any of the unstable modes, but normally one is interested in finding a controller that stabilizes the entire system so this should never be the case. The reason why $Q$ can't have negative eigenvalues is because then a solution that minimizes the integral would have states that blow up to infinity. Similarly the same holds why $R$ can't have negative eigenvalues. Basically any negative eigenvalue for $Q$ or $R$ would result in trying to find the minimum of a quadratic function with a saddle point. Lastly $R$ can't have eigenvalues of zero, because that would allow some parts of $u$ to blow up to infinity in order to force $x$ to zero as fast as possible (instantly) without penalty (in the discrete case it is allowed, which would result into a deadbeat controller).
The LQR problem can then indeed be formulated as the following CARE
$$ A^\top L + L\,A - L\,B\,R^{-1} B^\top L + Q = 0. $$
It can be noted that the solution for $P$ also allows you to evaluate the cost function integral given the initial conditions of $x$ when the LQR controller is used, namely
$$ \int_0^\infty x^\top Q\, x + u^\top R\, u\, dt = x(0)^\top L\,x(0). $$
The constraints on $Q$ and $R$ can also be translated to your CARE
$$ A^\top L + L\,A - L (B_1 + B_2) L + Q = 0. $$
Namely when $B_1 + B_2 \geq 0$ and $Q \geq 0$ then the eigenvalues of $L$ should at least all be non-negative and finite. The 'input' might still go to infinity, such that at least some of the modes of $x$ go immediately to zero. But this would imply that the cost function should evaluate to a finite value such that for any non zero initial condition ones gets $x(0)^\top L\,x(0) \geq 0$, so $L$ could be positive semi definite.
However when either $B_1 + B_2 < 0$ or $Q < 0$ you would get the saddle point case again. This would mean that $x$ would go off to infinity with $x(0)^\top L\,x(0) = -\infty$, so $L$ would have eigenvalues at $-\infty$. And I believe CARE solvers use iterative methods, but iterative methods have an hard time converging to $-\infty$.
So if you have $B_1 + B_2 < 0$ ($R^{-1} < 0$) there shouldn't be a sensible solution or you might have gotten a minus sign wrong.