Quote:Original post by swiftcoder
I am very surprised how well my trial/error steering behaviour works, so I think for this particular project I will forgo the analytical approach. Still very interested though, as I would like to build a more complex simulation in the future.

Wow, that was quick! Yeah, I've been pretty much assuming that you have the problem solved to your satisfaction by your steering behavior and that this is all academic -- so I was doing this mainly to satisfy my own interest (and for anybody else who reads this, since people don't talk about optimal control much in this forum but it seems useful).

Also, while you were replying, I edited my previous post; there's a new paragraph at the end that starts "The only thing that I've glossed over here..." While I've actually used the algorithm I've described to solve a problem with a single goal state [it was just how to stabilize a simple harmonic oscillator in minimum time with force inputs of plus or minus 1 (it's more complicated than I'd have thought!)], I haven't done this for a terminal manifold (what we have here); that's why I had to add that last paragraph.

Anyway, I'll post back with the details later. I don't think the costate equations will be that bad; there's a nice general solution; I just don't remember it off the top of my head.

Yegods, I'm an idiot. Almost all the information I needed was staring me in the face on the Wikipedia page. I don't have the transversality condition on me (sorry, I forgot to check my notebook at the office. Tomorrow!), but we can do the rest...

We want to minimize the total time this takes, so the cost is just the integral of the constant function equal to 1 for all time (this integrates to t). That is, L=1. And recall, (defining X=[x,y,vx, vy, θ]T),

f(X) = (d/dt) X

is defined by

(d/dt)x = vx
(d/dt)y = vy
(d/dt)vx = a cos(θ)
(d/dt)vy = a sin(θ)
(d/dt)θ = ω

so the Hamiltonian is, (where λ=[λ1, λ2,..., λ5]T),

H = L + λT f(X)
= 1 + λ1 vx + λ2 vy + λ3 a cos(θ) + λ4 a sin(θ) + λ5 ω

where

(d/dt)λ = -dH/dX

or, elementwise,

(d/dt)λ1 = -dH/dx = 0
(d/dt)λ2 = -dH/dy = 0
(d/dt)λ3 = -dH/dvx = λ1
(d/dt)λ4 = -dH/dvy = λ2
(d/dt)λ5 = -dH/dθ = -a λ3 sin(θ) + λ4 a cos(θ)

and of course we also know the differential equation for the state, so together with the costate equations (above), we have a coupled system of 2*5=10 scalar differential equations (or, a 10-d vector differential equation), which we can simulate backwards from a goal state (and any admissable final costate, at that goal state) through time.

(Also, remember that at any instant, we choose the (a, ω) that minimizes the Hamiltonian.)

This gives us almost everything we need. All that's left is to know which final costates go with which goal states; that's given by the transversality condition, and I'll post back with that later.

AWESOME! This stuff is coming back. (I should really have remembered this stuff better, but if I'm relearning it now, I'll be happy with that!)

[Edited by - Emergent on October 27, 2008 7:43:30 PM]

Ok, FINALLY:

The transversality condition with free end times is,

(L + dφ/dT + λTf)t=T = 0

where the only thing you haven't seen before is φ, the terminal cost (which in this case we do not have, so φ=0). And since also, in our case, L=1, we have,

1 + (λTf)t=T = 0

or -- and if I could put this in a big bold box, I would,

(λT f)t=T = -1 .

This is the pure min-time transversality condition. In the case of our asteroids spaceship,

(λ1 vx + λ2 vy + λ3 a cos(θ) + λ4 a sin(θ) + λ5 ω)t=T = -1

w00t!


[Edited by - Emergent on October 28, 2008 3:54:50 PM]

Oh! And since at time t=T, we require x=y=vx=vy=0, this reduces to,

(λ3 a cos(θ) + λ4 a sin(θ) + λ5 ω)t=T = -1

This characterizes the optimal final costates.

But finally, there's one thing to note, which is a restatement of this: Since f and L are not time-varying, the Hamiltonian H is constant along optimal trajectories. Since from the above, and the definition of the Hamiltonian, we see that Ht=T=0, and we know that H is constant, then we know that, in fact, H=0 along optimal trajectories.
[EDIT: The struck statement above is only true when dH/du=0. This is what occurs when we minimize H with respect to u when u is continuous and unconstrained. But in our case, u can only take one of a few discrete values, and there is actually no guarantee that the Hamiltonian is zero along optimal trajectories -- just minimized.]

Ah hah! Wow; how had I forgotten that?

So, FINALLY, this gives us an algorithm (sort of... I've left out handling trajectories that intersect each other. When this happens, you want to discard the slower trajectory).

- For every final u=(a, ω) // There are six    - Solve H=0 for λfinal // Numerically, perhaps (e.g., w. Newton-Raphson)    - For each λfinal found above [there are multiple (infinite, actually) solutions], do:        - Using this λfinal and your goal Xfinal, simulate          backwards in time for, say, 10 seconds.  At every timestep of this simulation,          use whichever u minimizes the Hamiltonian.  And record the (X, u) pairs in          some data structure as outlined before.

So, I've left out how to deal with trajectories that intersect each other, but hopefully this gives the general idea.

[Edited by - Emergent on November 3, 2008 6:30:54 PM]

Quote:Original post by swiftcoder
Sweet, that actually works out a fair bit simpler than I expected.

:-)

Also, FYI, I just corrected a few mistakes in the algorithm above. And, since I was kind of fuzzy about how to iterate over the values of λfinal that meet the transversality condition, let me clear that up. In this case, here's one way:

Since,

λ3 a cos(θ) + λ4 a sin(θ) + λ5 ω = -1

then let's express (λ3, λ4) in polar coordinates (A, δ), where,

A = a sqrt(λ32 + λ42)
δ = atan2(λ4, λ3)

so the transversality condition becomes,

A cos(θ + δ) + λ5 ω = -1

or

cos(θ + δ) = -(λ5 ω + 1)/A

and changing coordinates again, defining α=-1/A,

cos(θ + δ) = α λ5 ω + α


There are only two solutions in (0, 2 π) for δ. They are given by,

δ1 = cos-1(α λ5 ω + α)
and
δ2 = δ1 + π

so what we can do is grid over (α, λ5), and for each point solve for δ.

[EDIT: I made this algorithm up. There are better, more well-established techniques. One in particular is known as Value Iteration. This discussion is continued in swiftcoder and my replies to this newer thread. (I.e.: The following post is pretty much obsolete now).]

A completely different approach

Here, I just made this up, but see if it makes sense. It doesn't use the calculus of variations at all, and is very intuitive. The idea depends on two different discretizations:

1. There are 'particles' in state space.
2. The state space is partitioned into small 'cells;' e.g., a uniform grid of small cubes.

Then, the idea is,

Initialization: Start with a list of particles containing feasible final states. These should be a good 'sample' of the terminal manifold. If there is a single goal state, then this is a single particle. Also, maintain a spatial structure of cells (e.g., a simple grid). Each cell has a corresponding 'cost;' this is initialized to infinity.

Then, iterate the following:

Each particle "splits" into N particles, where N is the number of possible control inputs, and each of these particles is simulated (as x_dot = f(x,u) with the corresponding constant u) until it leaves the cell that it began in at the beginning of this step. When it reaches a new cell, if its cost is lower than that currently stored on the cell, overwrite the label, and otherwise destroy the particle.

Do this until the list of particles is empty.

Now, at the end, each cell contains an approximate value for the Bellman Value Function, V(X).

[Edited by - Emergent on February 23, 2009 11:43:10 PM]

If you haven't already you should look at a game called Gravity Well. They do exactly what you are talking about. I think they have source code out there somewhere.

For those of us a little less math savvy could you explain, in simpler terms, the formulas you're using? I'm working on a little home project and you described perfectly a problem I was coming up too. I would love to see how you resolved it.

Thank you in advance.