Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

On The Solution of Matrix Ricatti Equation ODE

  1. Feb 14, 2013 #1
    I've become a little confused about why noone cares to actually explicitly solve the Matrix Ricatti Differential Equation (RDE) of the form:

    $$ -{\dot{P}} = Q + PA^T + A^TP + PBB^TP $$ where [itex] BB^T, Q, P [/itex] are a positive-definite matrices, and [itex] A, BB^T, Q, P \in \mathbb{R}^{n \times n} [/itex]

    This equation pops up all the time in controls. A solution to this ode also appears to be avoided for some reason, even though deriving it (in my opinion) is nontrivial. Neither professors nor textbooks seem to address solving this other than to comment away the need to discuss a solution by saying "Often you need to solve this numerically". I asked about this and still felt curious but rather than asking more questions and just looking stupid, I wanted to fill what seems to be a void with my solution. It was also a good exercise for me.

    Note: In the following derivation, [itex] P,R,X,J [/itex] are time-varying matrices, yet I will only explicitly show a time argument when this needs to be emphasized.

    Step 1: Introduce change of coordinates.

    $$ R(t) = P(t) - R_o $$
    where [itex] R_o [/itex] is constant and must be positive-definite (PD). Solving for [itex] P [/itex] and substituting this into the RDE and collecting terms results in
    $$ -{\dot{R}} = (Q + R_oA + A^TR_o - R_oBB^TR_o) + R(A-BB^TR_o)+(A-BB^TR_o)^TR - RBB^TR $$
    and select [itex] R_o [/itex] to satisfy
    $$0=Q + R_oA + A^TR_o - R_oBB^TR_o $$
    to get rid of the constant, resulting in a new RDE below
    $$ -{\dot{R}} = R(A-BB^TR_o)+(A-BB^TR_o)^TR - RBB^TR $$
    Notice that the constraint resulting in the equation above is actually the Algebraic RE (no derivative term) and the solution to this will be discussed later, and remember now that [itex] R_o [/itex] is no longer arbitrary. In fact, for my application (I believe) I may not need to solve for this matrix due to something really really really convenient, but this post should solve the RDE given the solution to this algebraic

    Step 2: Decompose [itex] R(t) [/itex]

    define matrices $$ X(t), J(t) $$ as both N x N matrices such that they satisfy
    $$R(t) = X^T(t)J(t)X(t)$$
    and note that while [itex] J(t) [/itex] is positive-definite (PD), [itex] X(t) [/itex] is not. Plugging this into the RDE results in (this will get messy)
    $$ -{\dot{R}} = -{\dot{X}}^TJX - X^T{\dot{J}}X - X^TJ{\dot{X}} = X^TJX(A-BB^TR_o)-(A-BB^TR_o)^TX^TJX - X^TJXBB^TX^TJX $$

    Step 3: Undetermined Coefficients

    At this point you do "pattern matching" or "undetermined coefficients". You get two equations from this (two identifical equations for [itex] X(t) [/itex]). The first equation is:
    $${\dot{X}} = X(BB^TR_o-A)$$
    and the second is
    $${\dot{J}} = JXBB^TX^TJ$$

    Step 4: A Solution.

    The first equation in step 3 represents the matrix ODE for linear time-invariant systems.
    $${\dot{X}} = X(BB^TR_o-A)$$
    Inverting the system and integrating yields
    $${\dot{X}}(t)X^{-1}(t) = (BB^TR_o-A)$$
    $$\int_{t}^{T}{\dot{X}}(s)X^{-1}(s)ds = \int_{t}^{T}(BB^TR_o-A)ds$$
    $$log(X(T-t)) = (BB^TR_o-A)(T-t)$$
    and exponentiating...
    $$X(T-t) = e^{(BB^TR_o-A)(T-t)}$$

    The second equation in step 3 is a bit more complex and requires some more manipulation and thought, and the previously solved solution [itex] X(t) [/itex]. Inverse [itex] J [/itex] on both sides to get
    $$J^{-1}{\dot{J}}J^{-1} = XBB^TX^T$$
    Integrating (the odd looking limits again come from my particular application, optimal control).
    $$\int_{t}^{T}J(s)^{-1}{\dot{J(s)}}J(s)^{-1}ds = \int_{t}^{T}X(s)BB^TX^T(s)ds$$
    This integral looks nasty, but on the left hand side is actually setup for chain rule of PD matrices, on the right hand side is a controllability grammian for an LTI system with the system matrix [itex] BB^TR_o-A [/itex]. But keep in mind controllability grammians have to be full rank (i.e., invertible) for all time in order for the related system to be controllable. To be sure that it plays the role of the controllability grammian as it does in linear systems, trying to solve for [itex] J(t) [/itex]
    $$J^{-1}(T-t) = \int_{t}^{T}X(s)BB^TX^T(s)ds$$
    Inversing gives
    $$J(T-t) = \left( \int_{t}^{T}X(s)BB^TX^T(s)ds \right) ^{-1}$$
    and THAT is why the controllability grammian needs to be invertible for all time. If the controllability grammian is not full rank, we cannot control the system, and thus we cannot invert this matrix to create our solution to the system.
    Combining [itex] X(T-t), J(T-t) [/itex] we have
    $$R(t) = X^T(T-t)J(T-t)X(T-t) = e^{(BB^TR_o-A)^Tt}\left( \int_{t}^{T}X(s)BB^TX^T(s)ds \right) ^{-1}e^{(BB^TR_o-A)t}$$
    And yet we still need to add the constant term to find the solution in the original coordinate system
    $$P(t) = R(t) + R_o = X^T(T-t)J(T-t)X(T-t) + R_o = e^{(BB^TR_o-A)^T(T-t)}\left( \int_{t}^{-1}X(s)BB^TX^T(s)ds \right) ^{-1}e^{(BB^TR_o-A)(T-t)} + R_o$$

    I've reduced the problem of finding [itex] P(t) [/itex] for
    $$ -{\dot{P}} = Q + PA + A^TP - PBB^TP$$
    to finding [itex] R_o [/itex] that solves
    $$0=Q + R_oA + A^TR_o - R_oBB^TR_o $$

    I left out some details I think, partly because I'm still learning them for myself, and some things may only make sense here if you're using this for optimal control like me. Regardless, I think I've left the meat and potatoes here. Let me know what you think!
    Last edited: Feb 14, 2013
  2. jcsd
  3. Feb 15, 2013 #2
    I put some effort into this and if you want to print this out for w/e reason or save this, I've made a pdf file same post here

    https://dl.dropbox.com/u/25124237/Blog/jimseebooks/Matrix%20RDE%20Solution%20140213.pdf [Broken]
    Last edited by a moderator: May 6, 2017
  4. Feb 16, 2013 #3


    User Avatar
    Gold Member

    Well, in your first step it should be:

    [tex] +R_0 BB^T R_0[/tex] and not with a minus sign, also the rest of RHS should be with plus signs, cause you plug [tex] P=R+R_0 [/tex].

    I've noticed some more misprints, such as you forgot Transpose in some of A's.

    A piece of advice don't hurry, be precise and slow.
  5. Feb 17, 2013 #4
    Actually it looks like the I wrote the original equation wrong. It should be

    $$ -{\dot{P}} = Q + PA^T + A^TP - PBB^TP $$

    Last edited: Feb 17, 2013
  6. Feb 18, 2013 #5


    User Avatar
    Gold Member

    OK, I spotted another mistake (or so I think).

    The integral

    Shoule be [tex] \log(X(T))-\log(X(t))=\log (X(T)/X(t))[/tex] and not [tex] \log X(T-t)[/tex]
  7. Feb 18, 2013 #6
    Yes I believe you're right.

    Aside from these notation mistakes, do you believe the final answer to be correct?
  8. Feb 18, 2013 #7


    User Avatar
    Gold Member

    I'll be reading it tomorrow again, by the weekend I hope I'll have time to respond to you.

    I just spotted these mistakes first, btw I think it's best to read your own stuff once more and be sure that every step is justified, and write your justification down for every equality.

    I myself have exams the next two days and next monday, not sure why I am being so altruistic... :-)
  9. Feb 18, 2013 #8
    I understand. I have related issues.
  10. Feb 18, 2013 #9
    I haven't studied your derivation, but I believe writing the differential Riccati system as a algebraic Riccati system is called Chandrasekhar decomposition. The problem has now moved to solving the algebraic system, which is still difficult to solve in general.
    There is also the Bernoulli substitution method, which transforms the Riccati ode system to a (larger) system of linear ode's (Hamilton ode's).

    A large number of books have been written on the subject of solving matrix riccati equations, they might give you some more ideas.
  11. Feb 20, 2013 #10
    Thanks for the ideas!

    I'll see if I can get my hands on some of those books.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook