Step response of transfer function different in ODE45 than in 'step' and 'lsim'

I want to simulate a double integrator system (a simple mass block) that is controlled with a PD controller for positioning.

I'm both trying to implement an ODE45 version, and comparing it to 'step' or 'lsim'. Later, I want to incorperate a nonlinear system and prefer to use the ODE45, so that is why I am trying this.

The result is not promising:

Furthermore, if I increase the Kp = 4, the feedback system is critically damped, and the discrepency is even larger!

This is the code

close all


m = 1;       % Mass
Kp = 1;
Kd = Kp;


x_desired = 1.0;     % Desired position
v_desired = 0.0;     % Desired velocity

% Initial conditions
x0 = 0.0;            % Initial position
v0 = 0.0;            % Initial velocity
initial_conditions = [x0; v0];

% Define the ODE function
odefun = @(t, y) [ ...
    y(2); 
    (Kp * (x_desired - y(1)) + Kd * (v_desired - y(2))) / m 
];

% Solve the ODE
[t, y] = ode45(odefun, [0 10], [0 0]);

% Plot the results
figure(1);
plot(t, y(:, 1), 'LineWidth', 2); % Position
xlabel('Time (s)');
ylabel('Position (x)');
title('Position vs Time');

sys = tf([1],[m 0 0]);
PDcontroller = tf([Kd Kp],[1]);
sysFB = feedback(sys*PDcontroller,1);


figure
[y,tout]=step(sysFB);

figure(1);
hold on
plot(tout,y)
legend(["ode 45","Step"])

 figure
 pzmap(sysFB)

 
 
% %% Compare a



 figure
 pzmap(sysFB)

Answer 1

[Paul]

Hi dikkemulle,

The compensation is different between the two cases.

For the ode45 implementation (ignoring v_desired for simplicity) we have

u = Kp*(xd - y1) - Kd*y2 = Kp*(xd - y1) - Kd*s*y1

For the LTI object implementation we have

u = (Kp + Kd*s)*(xd - y1)

In the latter case, the Kd*s term is applied to xd, which is not true the former.

Assuming the ode45 implemenation is desired, see below for how to achieve that using LTI objects.

If the LTI object implementation is desired, we'll have to do a bit more work to deal with differentiating the step input of xd.

m = 1;       % Mass
Kp = 1;
Kd = Kp;
x_desired = 1.0;     % Desired position
v_desired = 0.0;     % Desired velocity
% Initial conditions
x0 = 0.0;            % Initial position
v0 = 0.0;            % Initial velocity
initial_conditions = [x0; v0];
% Define the ODE function
odefun = @(t, y) [ ...
    y(2); 
    (Kp * (x_desired - y(1)) + Kd * (v_desired - y(2))) / m 
];
% Solve the ODE
[t, y] = ode45(odefun, [0 10], [0 0]);
% Plot the results
figure(1);
plot(t, y(:, 1), 'LineWidth', 2,'DisplayName','ode45'); % Position
xlabel('Time (s)');
ylabel('Position (x)');
title('Position vs Time');

sys = tf([1],[m 0 0]);

%PDcontroller = tf([Kd Kp],[1]);
%sysFB = feedback(sys*PDcontroller,1);

sys1 = feedback(sys,tf([Kd 0],1)); % Kd feedack only on the velocity
sysFB = feedback(sys1*Kp,1)

sysFB =
 
       1
  -----------
  s^2 + s + 1
 
Continuous-time transfer function.

[yout,tout] = step(sysFB);

hold on
plot(tout,yout,'DisplayName','step');
legend

[dikkemulle]

Thank you!

I'm trying to understand this in block diagram. The controller I assumed when using 'step'

and the controller I assumed in the ODE45 code:

Is this correct?

[Paul]

Almost. For the ODE45 case the equation for u is: u = P*e - D*y (the block diagram is missing the negative sign at the sum junction that forms u).