Continuous vs discrete time

Hello,

I am trying to design EKF in MATLAB, where my system dynamics (transition model) has this form --> (1), ​which is continuous in time, but my measurement model is discrete. So in order to convert Eq (1) to a discrete domain, I tried the ODE45 solver, and I expected the result in this form --> (2). But I have been told that the MATLAB ODE45 solver will produce only continuous results, i.e. integrating continuous-time domain differential equations using an ODE45 solver gives results in continuous-time, not in discrete, is that true? So for this reason I used forward Euler (first order RK). Which one is the correct one to transform equation (1) into Eq (2) since EKF in MATLAB expects everything in discrete?

Thanks in advance

Best

Hashir

[Paul]

Hi Hashir,

The function ode45 does not convert a differential equation into the form of a discrete-time, approximating difference equation. Rather, ode45 solves for an approximate numerical solution to the differential equation, which can be used in an EKF implementation.

Regarding "... since EKF in MATLAB expects everything in discrete?" Does this refer to a specific EKF function that ships with Matlab? If so, can you provided a link to the relevant doc page?

[VAHA]

Hi @Paul Paul,

Regarding this statement : "since EKF in MATLAB expects everything in discrete?" what I meant is, in MATLAB the ExtendedKalmanFilter expects the vehicle model and measuement model in discrete-time. Exactly the ode45 approximates the differential equation numerically. So is there any method to convert Eq (1) into (2) completely?

Answer 1

[Star Strider]

It is difficult to follow what you are doing as you describe it. It might be best to use the Control System Toolbox functions if you have access to them, and the extendedKalmanFilter function (introduced in R2016b, improved in R2020b) might be particularly helpful. (I have limited experience with Kalman filters in general, so this reference is likely the best I can do.)

You can force the MATLAB ODE solvers to produce output at specific times by giving it a vector of (ideally a constant sampling interval) times using the linspace function.

[VAHA]

Thanks Strider. Actually if I want to use extendedKalmanFilter function, all my vehicle dynamics and measurement model must be discrete in time. But as you can see from my question, my dynamic model is continuous in time, even if I use constant sampling time, say T=0.1second and odesolver evaluates the function during constant time intervals using linspace say between 0 to 10 seconds with 100 steps, will I get the results in discrete time? Because in general we use odesolver to integrate the differential equations but then results will be in continuous time not in discrete (subject to correction). But here my concern, how do I transform my Eq (1) to Eq (2) using any MATLAB default function?

[Star Strider]

There are several methods of transforming continuous to discrete time. Computers by definition use discrete time, so you will need to sample your data appropriately to use the discrete methods that are available.

[VAHA]

Actually, that is my question. I know a lot of numerical integration methods. I tried forward Euler and RK45 methods, of course, they will give discrete results. But ODE45 solver which is based on RK45 gives values at the given time, and it looks like discrete, but they are continuous. That is how I interpreted. Even if I integrated my Eq (1) using ODE45 solver, will I get Eq (2)? That is all about my doubt.

[Star Strider]

If you have the individual time instants in a vector, and solve the differential equations at those time instants using the MATLAB ODE integrators, you should get the result you want. The ODE integrators iwll interpolate the results of their integrations to the time instants your code specifies. See the ode45 documentation for tspan for details (those should be common to other solvers as well).

[Torsten]

There will always be a discrete formula behind an ODE solver for advancing the solution from time t_k to time t_(k+1). Your formula (2) e.g. is valid for one-step methods like RK, Euler etc.. (For multistep methods, x_(k+1) will not only depend on x_k, but on solutions that are more than one time step apart.) But these update formulae can change from step to step and will in general look very complicated. This is necessary for that the discrete advancement in time comes close to the continuous solution.

So in short: You will get equation (2) from equation (1) by using ODE45, but the f_k change with t_k and x_k and can look very complicated.

[Paul]

Hashir,

I don't think you need to worry about the form of the update equation. For an EKF with a continous-time system and discrete measurements, you just need to be able to propagate the state and error covariance from one time step to the next using whatever way you want. You can do that using Euler integeration, or ode45, or whatever.

In fact, if you track through the first example on the doc page that @Star Strider linked, you'll see that the state equation in vdpStateFcn.m is actually doing Euler integration of a continuous time differential equation, but I don't see why that couldn't be replaced with a call to ode45. (though I don't know how that will impact the estimation of the Jacobians in extendedKalmanFilter)

dbtype vdpStateFcn.m 25:33

25    % Euler integration of continuous-time dynamics x'=f(x) with sample time dt
26    dt = 0.05; % [s] Sample time
27    x = x + vdpStateFcnContinuous(x)*dt;
28    end
29    
30    function dxdt = vdpStateFcnContinuous(x)
31    %vdpStateFcnContinuous Evaluate the van der Pol ODEs for mu = 1
32    dxdt = [x(2); (1-x(1)^2)*x(2)-x(1)];
33    end

One thing to consider is that for an EKF for a continuous-time system it's perhaps more correct to use this form of the EKF equations where the error covariance is propagated with the state via a differential equation as should be the case for a continuous-time system. However, extendedKalmanFilter seems use a discrete-time implementation as shown here.

[VAHA]

Hi @Paul Paul,

Thank you for your response. Actually, you are right, the Continuous-Discrete EKF takes a continuous-time vehicle dynamic model and discrete time measurement model. During the prediction phase, the error covariance is predicted by integrating the Algebraic Riccati equation. In this EKF, during the prediction phase, the state and the error covariance can be integrated forward in time using ode45, and I did it using RK45 integration. But my question is, can I use discrete EKf if I have a vehicle model in continuous and a measurement model in discrete? Currently, I am using Forward Euler (for discrete EKF with continuous time dynamic model and discrete measurement) to discretize my vehicle's dynamic model as , and the state transition matrix, , where is the jacobian of the dynamic model ), but it is too inaccurate. So I am skeptical about the solution obtained from the ODE45 solver, that is, if I use the ODE45 solver can I get the discrete one? Because in that case, I may also need to discretize the Jacobian also as I did with forward Euler.

[Paul]

I was following everything, for the most part, until I got to the part about the state transition matrix. The EKF on the linked page, either the continous-discrete or the discrete-discrete, does not use a state transition matrix.

Here's an example that shows using ode45 for the state prediction is, at least, feasible using extendedKalmanFilter for a continuous system with discrete meausrements. However, I'm still not clear on how the extendedKalmanFilter estimates the Jacobian for the state equation. Also, the extendedKalmanFilter is propagating the error covariance for a discrete time system as opposed to integrating the differential equation for the error covariance. Consequently, we have to be careful about the definiton of ProcessNoise, which is a spectral density for continuous time systems and has to be converted to a covariance for discrete-time propagation, if you want ProcessNoise to represent a physical process as opposed to just a user-defined parameter. Anyway, the code below shows how to use ode45 for the state propagation, w/o worrying about any other issues.

Model the truth system and plot the solution

vdpode = @(t,x) [x(2); (1-x(1)^2)*x(2)-x(1)];
x0 = [2;0];

sol = ode45(vdpode,[0 5],x0,odeset('MaxStep',0.05));

figure
plot(sol.x,sol.y)
hold on

Define EKF measurement function

measurementfunction = @(x) x(1);

Define the state prediction function using Euler integration and the EKF with default parameters.

state1 = @(x,dt) x + vdpode(0,x)*dt;
ekf1 = extendedKalmanFilter(state1,measurementfunction,x0);

Execute the EKF. Use the truth data as measurements. Plot the result.

dt = 0.05;
tvec = 0:dt:5;
xhat1 = nan(2,numel(tvec));

for ii = 1:numel(tvec)
    meas = deval(sol,tvec(ii),1);
    xhat1(:,ii) = correct(ekf1,meas);
    predict(ekf1,dt);
end

plot(tvec,xhat1,'o'),grid

Define state prediction function using ode45 and the EKF with default parameters.

% the vdp ode is not dependent on time, so use [0 dt] as the tspan.
state2 = @(x,dt) ode45(vdpode,[0 dt],x).y(:,end);
ekf2 = extendedKalmanFilter(state2,measurementfunction,x0);

Execute and plot

xhat2 = nan(2,numel(tvec));

for ii = 1:numel(tvec)
    meas = deval(sol,tvec(ii),1);
    xhat2(:,ii) = correct(ekf2,meas);
    predict(ekf2,dt);
end

plot(tvec,xhat2,'x'),grid

[VAHA]

THank you for all your valuable information