Discrete Derivative block (Simulink) behaves differently than expected

I'm discretizing a known-working continuous system model that includes x state outputs from the simulink Varying State Space block, but the discrete derivative block is not behaving as expected.

The SS model is a basic 2nd order plant, and I'm just giving a step input for this post for debugging/illustration purposes:

zpk(plant) =

From input "u" to output "theta":

3.1978e+05

-----------------------

(s^2 + 360s + 3.198e05)

The rate transition block is to 100 kHz (1e-5 s)

The State Space block outputs are:

y: position, and

x: x1 is position, and x2 is velocity.

I'm comparing the x (state vector) outputs to a discrete derivative on y, since the discrete derivative of position should be velocity.

I expected that they would match, but I'm seeing the discrete derivative behave very differently than expected: see red step-wise line, which is "y discDeriv", and does not seem to touch either sloped line (ie much delayed).

The other outputs "x2 @rate" and "dy/dt @ rate" match, but "y discDeriv" does not.

Why is Simulink Discrete Derivative not outputting the same result as the other output paths, specifically the x2 vel state ("x2@rate") output?

What is the Discrete Derivative block actually doing? The documentation seemed sparse, and what was there, seemed (to my interpretation) to support that it should behave the same as the other methods i tried.

Answer 1

[Paul]

Hi John,

How did you get the yellow colors on the rate transition blocks on the signals coming out of those? If you turned on Sample Time Colors, I wonder why the output of the Discrete Derivative block isn't red.

Anyway, I'm pretty sure what's happening is that the solver is running a at higher rate than 1e-5. The discrete derivative block is updating at a sample time of 1e-5 secs, and at each sample time computes

y_discderiv[n] = (y@rate[n] - y@rate[n-1])/1e-5, where the samples of y@rate are taken at 1e-5.

But the continuous signal y is updated faster (is my assertion), and the output of the Derivative block is essentially

dy/dt = (y(t) - y(tprev))/(t - tprev)

where t is current sim time and tprev is the time of the last major step. If t - tprev < 1e-5 (as I expect), then dy/dt won't match y_discderiv (at the 1e-5 sample points) because they are computed based on different quantities.

If you want them to match, you can change your solver settings to fixed-step with a step size of 1e-5. But there is no guarantee that a step size of 1e-5 is small enough to accurately simulate the continuous plant dynamics (you'd have to check that), keeping in mind that whatever solver settings you have are already using a smaller time step.

[Paul]

Here's an example

Create the model

bdclose('all');
hsys = new_system('derivtest1');
hconstant = add_block('simulink/Sources/Constant','derivtest1/Step');

hplant = add_block('simulink/Continuous/Transfer Fcn','derivtest1/Plant');
set_param(hplant,'Numerator','3.1978e+05');
set_param(hplant,'Denominator','[1 , 360 , 3.198e05]');
add_line(hsys,get_param(hconstant,'PortHandles').Outport, ...
              get_param(hplant,   'PortHandles').Inport);


hderiv = add_block('simulink/Continuous/Derivative','derivtest1/deriv');
add_line(hsys,get_param(hplant,'PortHandles').Outport, ...
              get_param(hderiv,'PortHandles').Inport);

hout1 = add_block('simulink/Sinks/To Workspace','derivtest1/out1');
set_param(hout1,'VariableName','dydt');
add_line(hsys,get_param(hderiv,'PortHandles').Outport, ...
              get_param(hout1, 'PortHandles').Inport);

hrate1 = add_block('simulink/Signal Attributes/Rate Transition','derivtest1/Rate1');
set_param(hrate1,'OutPortSampleTime','1e-5');
add_line(hsys,get_param(hplant,'PortHandles').Outport, ...
              get_param(hrate1,'PortHandles').Inport);

hdiscderiv = add_block('simulink/Discrete/Discrete Derivative','derivtest1/DiscDeriv');
add_line(hsys,get_param(hrate1,    'PortHandles').Outport, ...
              get_param(hdiscderiv,'PortHandles').Inport);

hout2 = add_block('simulink/Sinks/To Workspace','derivtest1/out2');
set_param(hout2,'VariableName','yDiscDeriv');
add_line(hsys,get_param(hdiscderiv,'PortHandles').Outport, ...
              get_param(hout2, 'PortHandles').Inport);


hrate2 = add_block('simulink/Signal Attributes/Rate Transition','derivtest1/Rate2');
set_param(hrate2,'OutPortSampleTime','1e-5');
add_line(hsys,get_param(hderiv,'PortHandles').Outport, ...
              get_param(hrate2,'PortHandles').Inport);


hout3 = add_block('simulink/Sinks/To Workspace','derivtest1/out3');
set_param(hout3,'VariableName','dydtAtRate');
add_line(hsys,get_param(hrate2,'PortHandles').Outport, ...
              get_param(hout3,'PortHandles').Inport);

Ensure use of variable step solver with a max step less than the 1e-5 Sample Time for the discrete derivative

cset = getActiveConfigSet('derivtest1');
set_param(cset,'StopTime','1e-4');
set_param(cset,'Solver','VariableStepAuto');
set_param(cset,'MaxStep','0.2e-5');
out = sim('derivtest1',cset);
figure
plot(out.dydt,'-x','DisplayName','dydt'),
hold on
plot(out.dydtAtRate,'-o','DisplayName','dydtAtRate');
plot(out.yDiscDeriv,'-s','DisplayName','yDiscDeriv');
legend

We see that dydt is computed at a much higher rate than yDiscDeriv. Basically, computing the derivative at a high rate then sampling the derivative at a low rate is NOT the same as sampling the position at a low rate and then computing derivative.

Change the solver to fixed step with the same step size used to compute the discrete derivative

set_param(cset,'Solver','FixedStepAuto');
set_param(cset,'FixedStep','1e-5')
out = sim('derivtest1',cset);
figure
plot(out.dydt,'-x','DisplayName','dydt'),
hold on
plot(out.dydtAtRate,'-o','DisplayName','dydtAtRate');
plot(out.yDiscDeriv,'-s','DisplayName','yDiscDeriv');
legend

bdclose('all')

The results match, but the accuracy of continuous states might suffer because the integration step size is increased by at least a factor of 1e-5/0.2e-5 = 5. And there might be other reasons to desire a variable step solver depending on other aspects of the model.

If you really want to roll your own estimate of a derivative with a variable step solver, you could use a Memory block to delay the position by one major time step, subtract that from the current position, then divided that result by the differene between the current clock time and the previous clock time from another Memory block. Of course, the result won't, in general, have a fixed sample time, which might be what you're looking for if you're going to use the estimated derivative in discrete-time control loop.

Also, dydt is computed with a difference formula applied to position, so it won't exactly match the velocity state from the plant that is computed via numerical integration of the differential equation.

[John]

@Paul "How did you get the yellow colors on the rate transition blocks on the signals coming out of those? If you turned on Sample Time Colors, I wonder why the output of the Discrete Derivative block isn't red."

Yeah, it was sample time colors -- and good question on why Disc Deriv isn't red.

"Anyway, I'm pretty sure what's happening is that the solver is running a at higher rate than 1e-5. The discrete derivative block is updating at a sample time of 1e-5 secs, and at each sample time computes

y_discderiv[n] = (y@rate[n] - y@rate[n-1])/1e-5, where the samples of y@rate are taken at 1e-5.

But the continuous signal y is updated faster (is my assertion), and the output of the Derivative block is essentially

dy/dt = (y(t) - y(tprev))/(t - tprev)

where t is current sim time and tprev is the time of the last major step. If t - tprev < 1e-5 (as I expect), then dy/dt won't match y_discderiv (at the 1e-5 sample points) because they are computed based on different quantities."

Not sure i'm following... Even when the solver rate is very high (>1e-5), the Disc Deriv block is pulling from a series of high-rate points that are the same value, and it will still have a 0 deriv value until a step change occurs.

As we can see, the (t-tprev) gives the same magnitude step on the Rate Transition changes -- just shifted by half the amplitude lower for the very first step.

If the solver were the issue (i'm not sure, just raising the question), wouldn't the denominator change and the step not be the correct magnitude for all steps? Why would only the first step be half the magnitude?

"We see that dydt is computed at a much higher rate than yDiscDeriv. Basically, computing the derivative at a high rate then sampling the derivative at a low rate is NOT the same as sampling the position at a low rate and then computing derivative."

Good point -- agreed. But that scenario is specifically what disc deriv is for, i thought? Ie lower sampling, then compute the derivative.

I do see how increasing the rate before the Disc Deriv causes "y discDeriv" to approach dy/dt, as expected.

It's still not clear to me why it's shifted by half the amplitude relative to x2@rate, since the input to the disc deriv ("y@rate") is at the same value as y, ie the inputs to Disc Deriv seems like they should be the same as the inputs to dy/dt. Since dy/dt is then sampled at the same rate as Disc Deriv, I thought they'd match.

But maybe I'm still missing your fundamental point...i'll read back through your post and think on it some more.

[Paul]

I'm not hung up on what the amplitude of the error is.

Still curious why your Sample Time Colors isn't marking that block and its output as red. That is really weird. You might want to turn on the Sample Time text annotations and legend as well (if you haven't already).

Anyway, from the 'x' symbols on my first plot we can see tha that the simulation is stepping at Tsim = 2e-6, at least for the time 0 <= t <= 1e-5.

Let's see what happens at t = 1e-5, i.e., at the first sample (after t = 0) where yDiscDeriv is computed. At t0 = 1e-5, we have

dydt(t0) = ( y(t0) - y(t0 - Tsim) ) / Tsim

yDiscDeriv(t0) = ( y(t0) - y(0) ) / 1e-5

Those two quantities are not the same, with the former generally being a better approximation to ydot. The same will hold true at t = 2*t0, t = 3*t0, etc.

[John]

"Still curious why your Sample Time Colors isn't marking that block and its output as red. That is really weird. You might want to turn on the Sample Time text annotations and legend as well (if you haven't already)."

Huh: the timing colors question is even more odd when i show the timing legend: it's blank. I also turned on the Timing Text information overlay, but that's also not appearing. So, I don't know why it's not showing what's expected... In other models this has always worked.

I'll hold off on replying to the rest, since maybe this impacts how the model responds and what it's showing...

[Paul]

The timing text is there, that's the "Cont" and "D2" symbols. The legend, as strange as it looks, indicates there's another discrete signal in the model at Sample Time D1 = 1e-7. Where's that in the model?

[John]

Ah, that's the fixed-step solver time. It goes away if i change back to variable-step.

Not sure why the Disc Deriv block isn't yellow, but it shows D1 so it looks okay (perhaps just a bug)