11

The essence of my question is the following: I have a system of two ODEs. One has an initial-value constraint and the other has a final-value constraint. This can be thought of as a single system with an initial-value constraint on some variables and a final-value constraint on the others.

Here are the details:

I am trying to use a continuous-time finite-horizon LQR controller to drive a linear dynamical system. I'd like to continue using the Python ecosystem.

The system is in the form $\dot{x}(t) = Ax(t) + Bu(t)$, subject to $x(0)=x_0$

The LQR solution generates a matrix $K(t)$ such that the optimal control input u(t), linear in $x(t)$, is $u(t) = K(t)x(t)$.

where $K(t) = R^{-1} B^T P(t)$

and $P(t)$ is the solution to a continuous time Riccati differential equation (note that this $P(t)$ is a matrix)

$\dot{P}(t) = -A^T P(t) - P(t) A + P(t) B R^{-1} B^T P(t) + Q$ subject to $P(t_f) = Q$

$A$, $B$, $x_0$, $Q$, $Qf$, $R$, $t_f$ are all given.

In English: you have some dynamical system that starts in state $x_0$. The LQR controller generates a feedback matrix to use between time $0$ and $t_f$ ($t_f$ is commonly called the time-horizon of the problem)

Note that the two ODEs are coupled only in one direction -- the solution to $P(t)$ does not depend on $x(t)$. Therefore one way to solve the problem is to reverse the Riccati equation in order to turn the final-value problem into an initial-value problem and find a numerical solution between time $0$ and $t_f$ using a standard ODE integrator. I can then use this numerical solution to find $x(t)$. This concerns me because the numerical ODE solver for x(t) will not necessarily sample the ODE at the same times as the times in the numerical solution to $P(t). Maybe there is some clever way to enforce this.

The other way I foresee of solving the problem is to solve the sytem together, but I don't know how to deal with the mix of initial-value and final-value constraints. Are these problems computationally heavy to solve? Can I do it in SciPy / Python?

Gus
  • 251
  • 2
  • 9

4 Answers4

8

I disagree with the other answers. Since you only have one-way coupling between the reverse- and forward-time problems, it will be much more efficient to solve them in sequence, as you first propose. You simply need a solution $P(t)$ that can be evaluated at any time $t \in [0,t_f]$.

You can do this by interpolating between the output values. I recommend that you use a Runge-Kutta method that supports dense output. For instance, scipy.integrate.ode.dopri5 is based on such a method. So you should be able to specify very finely spaced output times without forcing the integrator to take very small steps (assuming that the scipy interface to it is implemented correctly).

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
  • Yes, this is indeed simpler. You can use any method to generate $P(t)$ and $P'(t)$ at whatever points the integrator finds necessary, and interpolate in between with a cubic Hermite spline if O(h^4) accuracy is enough. – Arnold Neumaier Mar 28 '12 at 12:39
6

This is called a two-point boundary value problem and is well studied.

The shooting method is very simple to program but may be extremely unstable numerically.

The standard way to solve these problems is using a multiple shooting approach and solving the corresponding nonlinear system of equations by a standard nonlinear solver. For a list of solvers for nonlinear systems of equations, see, e.g.,
http://www.mat.univie.ac.at/~neum/glopt/software_l.html#nonlin

You take as variables the states on a regular grid in time (usually no very fine grid is needed), and as equations the boundary conditions and the mappings that map the time t variables to the time t+h variables. This gives as many equations as variables. You only need to provide the routines for evaluating this mapping for a given configuration of states on the grid, and the nonlinear solver does everything else. (Perhaps you need multiple starting points if your initial guesses are poor.)

Wikipedia http://en.wikipedia.org/wiki/Direct_multiple_shooting_method has a useful description of the process, if the above description isn't detailed enough for you. The book by Stoer/Bulirsch cited there gives complete details.

Arnold Neumaier
  • 11,318
  • 20
  • 47
5

I don't know how to do it in Python, but the keyword you want to look for in the literature is the "shooting method". That's the name of a method that solves problems that have both initial and final value constraints.

Wolfgang Bangerth
  • 55,373
  • 59
  • 119
1

AUTO can solve two point BVPs and has a python interface and is relatively easy to install. http://www.ma.hw.ac.uk/~gabriel/auto07/node6.html.

If you go the route of wanting to solve P(t) first and feed it to the other ODE as an input, then an efficient way to set that up is using PyDSTool. PyDSTool is very easy to install on any platform, see http://pydstool.sf.net. It will, by default, only use linear interpolation though for your previously computed solution (so compute that at fine time resolution). However, you can force PyDSTool to step to exactly the desired time points even with an adaptive integrator (although that might be inefficient and lead to inaccuracies). But with small enough max time steps, the linear interpolation and a fast integrator (Dopri is built in) for the second system means you'll be fine for "regular" systems like this.

RHC
  • 211
  • 1
  • 5