ghorn / dynobud

your dynamic optimization buddy
GNU Lesser General Public License v3.0
27 stars 5 forks source link

estimation support #17

Open peddie opened 10 years ago

peddie commented 10 years ago

I'd use it more if there were a high-level trajectory estimation with irregular sensor timings, akin to the OCP one.

ghorn commented 10 years ago

Thanks for the suggestion. The thing that makes this non-trivial is that the irregular sensor timings make it NOT akin to the OCP one.

If you help me figure out the type you need, I can try to think of the best way to implement it.

-- x: differential states
-- z: algebraic variables
-- u: controls
-- p: parameters
-- r: dae residual
-- h: sensor measurement
-- w: state disturbance
data EstimationPhase x z u p r h w =
  EstimationPhase {
             -- | the system dynamics of the stage: @f(x'(t), x(t), z(t), u(t), p(t), w(t), t)@
             estDae :: forall a. Floating a => x a -> x a -> z a -> u a -> p a -> w a -> a -> r a

             -- | the sensor measurement function @h(x(t),t)@
           , estSensorFun :: forall a. Floating a => x a -> a -> h a

             -- | list of (time, data) sensor measurement tuples, time may be irregular
           , estSensorData :: [(Double, h Double)]

             -- | Lagrange objective for disturbances @J_{L,w}(w(t))@
           , estDisturbanceObjective :: forall a. Floating a => w a -> a

             -- | objective function for sensor mismatch @J_h(h(t_k), h(t))@
           , estSensorObjective :: forall a. Floating a => h a -> h a -> a
           }

The above type takes in system dynamics (which allows for disturbances), a measurement function, real measurement data and two objective functions. The total objective function is minimizing disturbances (estDisturbanceObjective) while also minimizing sensor mismatch (estSensorObjective).

This leads to a few questions

I'm gonna ask Moritz to weigh in on this one.

@peddie does this type fit your problem? Is there a more restricted type which fits?

rienq commented 10 years ago

Assuming the irregular timings would be defined online, I believe "integration methods with continuous output" is the correct answer to this.. ;-)

ghorn commented 10 years ago

This is an offline problem, but I don't care if you know the answer :). I do know how to calculate the state and expected sensor functions at an arbitrary time, but I don't know how to weigh sensor errors on an arbitrary grid.

Consider sensor measurements taken at time [1, 5, 6, 7, 10, 100] seconds. If you want to minimize "overall sensor error", I don't think you weigh them all the same. Maybe there is no better option...

ghorn commented 10 years ago

*I don't care that it is an offline problem if you know the answer. Not "I don't care if you know the answer or not" :p

rienq commented 10 years ago

Hmm, unless one measurement is more accurate than another one according to you, I don't see a reason to use a different weighting.. right?

ghorn commented 10 years ago

re-posted from email thread:

Consider a multiple shooting system with 101 timesteps of 1 second, and 7 sensor data h_k on a regular grid at these time points: [1, 10, 11, 12, 13, 14, 70] seconds. You now have 100 "disturbance controls" w_k to solve for. I think the objective should be:

Sum_0^99{ Norm2(w_k, Q)^2/100 } + Sum_[1,10,11,12,13,14,70]{ Norm2( h_k - h(x_k), R )^2/7 }

It seems to me that this is the Right Way statistically, but I would also expect the w_k residual to be larger near the clustered sensor measurements because of mis-modeling. My question is: is the best objective under the circumstances?

ghorn commented 10 years ago

from the email thread:

ok, the weighting.

first, i see no reason to divide the R by 7. each measurement adds info, so 8 are better than 7, so no division is necessary here.

for state noise, you should think of it like a least squares integral, this gives the right weighting when different timesteps are used.

how to weight state noise vs measurement noise is a hard problem, even if you have measured controls and an idea how big their errors are. the problem is that the noise is never white noise. (for both state noise and measurement noise).

one way to consistently deal with this issue is to create a linear filter that gives the right noise when its input is white noise. then you use this filter output as your state or measurement noise in the problem.

white state noise --> filter W--> w white measurement noise --> filter V --> v

but we would need to think longer about this - and in practice these autocorrelations are not known either.

best regards,moritz

peddie commented 10 years ago

Oh, I just understood what you're doing with the w_k vector. Sorry for the noise about that.

how is the sensor mismatch objective added up?

Can you explain more? I don't know what you mean.

Should the estimation horizon be from the min sensor time to max sensor time, or user-defined?

User-defined, please :)

I will have to look at that type in more detail.