I have the following system of coupled Lyapunov equations over the time horizont $0 \leq t \leq T$:
$$
\begin{align}
\dot{Q} &= AQ + QA^\intercal – BB^\intercal\\
\dot{P} &= AP + PA^\intercal + BB^\intercal\\
\Sigma_0^{-1} &= Q(0)^{-1} + P(0)^{-1}\\
\Sigma_T^{-1} &= Q(T)^{-1} + P(T)^{-1},
\end{align}
$$
where $Q(t),P(t)$ are invertible for all $t\in[0,T]$ (but not necessarily sign definite), $\Sigma_0,\Sigma_T > 0$. I know that the general solution to the Lypanunov equation $\dot{\Pi} = \Pi A + A^\intercal \Pi + Q(t)$ for some matrix $\Pi(t)$ with $\Pi(0) = \Pi_0$ is
$$
\Pi(t) = \Phi(t,0)\Pi_0\Phi(t,0)^\intercal + \int_{0}^{T}\Phi(t,\tau)Q(t)\Phi(t,\tau)^\intercal dt
$$
The problem is that this only holds for $\Pi(0) = \Pi_0$, but in my case I have coupled, nonlinear, boundary constraints. I realize this is a BVP, and I guess it is possible to solve with numerical techniques such as BVP4c in MATLAB, but I tried this, and it couldn't come up with a solution. Any thoughts are appreciated!
Best Answer
Instead of solving for the solution for $Q$ and $P$ one could solve for $Q^{-1}$ and $P^{-1}$. This would make your boundary constraints linear and thus easier to solve. To avoid having to keep denoting the inverse I define new variables as $\mathcal{Q} = Q^{-1}$ and $\mathcal{P} = P^{-1}$. For the dynamics of these new variable one can use that $X\,X^{-1} = I$, so $d/dt\,(X\,X^{-1}) = \dot{X}\,X^{-1} + X\,\dot{X}^{-1} = 0$ which solved for $\dot{X}^{-1}$ gives $\dot{X}^{-1} = -X^{-1}\,\dot{X}\,X^{-1}$. Using this and the initial differential equations gives
\begin{align} \dot{\mathcal{Q}} &= -\mathcal{Q}\,A - A^\top \mathcal{Q} + \mathcal{Q}\,B\,B^\top \mathcal{Q}, \\ \dot{\mathcal{P}} &= -\mathcal{P}\,A - A^\top \mathcal{P} - \mathcal{P}\,B\,B^\top \mathcal{P}. \end{align}
It can be noted that the original differential equations are linear in $Q$ and $P$, but the ones in $\mathcal{Q}$ and $\mathcal{P}$ are not. So there is a trade-off between easier boundary conditions and more difficult dynamics.
Instead of making your boundary constraints easier to solve you could also use an analytical solution to the differential equations. Namely, by solving for stationary solutions and vectorizing both differential equation can be transformed into an autonomous linear time invariant system. The stationary solutions are denoted with $\bar{Q}$ and $\bar{P}$, which can be found by solving $A\,\bar{Q} + \bar{Q}\,A^\top = B\,B^\top$ and $A\,\bar{P} + \bar{P}\,A^\top = -B\,B^\top$ respectively. Those equations can be seen as Sylvester equations, which have one unique solution if the spectra of $A$ and $-A^\top$ are disjoint, which should be true if $A$ is non-singular. It can be noted that $\bar{Q} = -\bar{P}$ due to the difference of minus sign in front of the $B\,B^\top$ term. By defining intermediate variables $\hat{Q} = Q - \bar{Q}$ and $\hat{P} = P - \bar{P}$, their dynamics can be shown to be
\begin{align} \dot{\hat{Q}} &= \hat{Q}\,A + A^\top \hat{Q}, \\ \dot{\hat{P}} &= \hat{P}\,A + A^\top \hat{P}. \end{align}
By using vectorization and the Kronecker product those differential equations can also be written as
\begin{align} \dot{x}_1 &= M\,x_1, \\ \dot{x}_2 &= M\,x_2, \end{align}
with $x_1 = \text{vec}(\hat{Q})$, $x_2 = \text{vec}(\hat{P})$ and $M = I \otimes A + A \otimes I$. Their solutions can be written as $x_i(t) = e^{M\,t}\,x_i(0)$, with $e^{X}$ the matrix exponential of $X$. Therefore, the values of $Q(T)$ and $P(T)$ can be found with
\begin{align} Q(T) &= \text{vec}^{-1}\!\left(e^{M\,T}\,\text{vec}(Q(0) - \bar{Q})\right) + \bar{Q}, \\ P(T) &= \text{vec}^{-1}\!\left(e^{M\,T}\,\text{vec}(P(0) - \bar{P})\right) + \bar{P}. \end{align}
It can be noted that such analytical solution also exists for the dynamics of $\mathcal{Q}$ and $\mathcal{P}$, but in that case one would have to use $\hat{\mathcal{Q}}^{-1} = \mathcal{Q} - \bar{\mathcal{Q}}$ and $\hat{\mathcal{P}}^{-1} = \mathcal{P} - \bar{\mathcal{P}}$. However, that would make the boundary constraints nonlinear again which defeats to purpose of that coordinate transformation.