Output of a transfer function with input

Hello there.

I have a transfer function:

M=2;
b=14;
q=20;
num=[1];
G=tf(num,[M b q])

G =
 
          1
  -----------------
  2 s^2 + 14 s + 20
 
Continuous-time transfer function.

And a input: 2*(heaviside(t)-heaviside(t-2))

Applying laplace: 2/s - (2*exp(-2*s))/s

How do I find out the tf output using the input I just mentioned?

I think I might use lsim, but I have no idea how to.

Thanks

Answer 1

[Paul]

The HeavisideAtOrigin should be 1 when using lsim

M=2;
b=14;
q=20;
num=[1];
G=tf(num,[M b q]);

syms t
F = 2*(heaviside(t)-heaviside(t-2));
T = linspace(0, 3);
sympref('HeavisideAtOrigin', 0);
u = double(subs(F, t, T));
y1 = lsim(G, u, T);

sympref('HeavisideAtOrigin', 1);
u = double(subs(F, t, T));
y2 = lsim(G, u, T);

% exact answer
syms s
y3 = ilaplace(1/(M*s^2 + b*s + q)*(2/s - 2*exp(-2*s)/s),s,t);
y3 = double(subs(y3,t,T)).';

plot(T,y1,T,y2,T,y3,'o')
legend('HeavisdeOrigin = 0','HeavisdeOrigin = 1','exact','Location','South')

[Sam Chak]

Hi Paul,

Is this analysis true for all LTI systems? Or sometimes it has to be assigned to 0, or 0.5? Or maybe it works due to the algorithm of the Heaviside function, specifically only in MATLAB?

I ponder on the issue because the default value of the Heaviside function at the origin is mathematically defined at 0.5. In reality, the signal actually jumps from 0 to 1 at the origin, maybe with a slight delay.

[Paul]

Hi Sam,

In continuous-time, the actual value of heaviside(0) doesn't matter because t = 0 is an infinitesimally small point. I don't know if that's the correct mathematical term, but the idea is that it's so small that it doesn't matter. The only thing that matters is that h(t) = 1 for t>0. That's why these integrals are all the same

syms t
int(heaviside(t),t,0,1)

ans = 

1

sympref('HeavisideAtOrigin',0);
int(heaviside(t),t,0,1)

ans = 

1

sympref('HeavisideAtOrigin',100);
int(heaviside(t),t,0,1)

ans = 

1

And the Laplace transforms would be the same as well.

However, lsim discretizes the continuous-time LTI system and then propagates the approximate solution in discrete-time. In that case, the value of the input signal at t = 0 makes a difference. Consider, for example, what happens when lsim uses the zoh approximation (it doesn't always do this), where the discrete-time propagation reflects an assumption that the input is constant over the interval between samples. In that case, the best assumption is that the input = 1 for 0 <= t < dt, because that's the value of heaviside(t) over that interval, except possibly at that infinitessimally small point at t = 0. So we want heaviside(0) = 1.

When simulating a system defined as a discrete-time system (as opposed to using a disrcetized approximation to approximate a continuous-time response), the input u[0] would be whatever is defined for the problem. However, the standard definiton (that I'm aware of) in discrete-time control systems is that the discrete-time unit step is defined as u[n] =1 for n >= 0. So again, if using heaviside for the discrete-time unit step, it's still correct that heaviside(0) = 1.

[Paul]

Another factor to consider is that the effect of the 2 second delay is to introduce a discontinuous, step change in the input at t = 2. If the T vector doesn't hit T = 2 excactly, then lsim will step right over the discontinuity. If the system is strictly proper, than that discontinuity is smoothed out by the system dynamics, and so might not be much of a problem if the system dynamics are slow. But if the system dynamics are fast, that might not be the case. In the limit if the system is proper, i.e., with direct feedthrough, then there should be a discontinuity in the output at T = 2. If the time vector doesn't hit T = 2 exactly, then there will be inaccuracy in the lsim solution for T>2. For example

M=2;
b=14;
q=20;

Change num so that G is proper

num=[.01 0 1];
G=tf(num,[M b q]);
syms t
F = 2*(heaviside(t)-heaviside(t-2));

Change T so that any(T == 2) is false

T = linspace(0, 2.99);
sympref('HeavisideAtOrigin', 0);
u = double(subs(F, t, T));
y1 = lsim(G, u, T);
sympref('HeavisideAtOrigin', 1);
u = double(subs(F, t, T));
y2 = lsim(G, u, T);
% exact answer
syms s
y3 = ilaplace((0.01*s^2 + 1)/(M*s^2 + b*s + q)*(2/s - 2*exp(-2*s)/s),s,t);
y3 = double(subs(y3,t,T)).';
plot(T,y1,T,y2,T,y3,'o')
legend('HeavisdeOrigin = 0','HeavisdeOrigin = 1','exact','Location','South')

The example clearly shows the problem with heaviside(0) = 0 at T = 0. And you can see, the lsim output, even for heaviside(0) = 1, doesn't match the exact solution for T>2. That's because lsim didn't catch the change in the input until

T(find(T>=2,1))

ans = 2.0235

which is why the lsim output lags the exact solution.

[Walter Roberson]

[I moved Paul's material to an Answer of its own, as he made significant improvements and commetary to my solution, and think he deserves to be awarded the aclaim.]

[Paul]

@Walter Roberson

Thanks for the kind words.

[Paul]

We can also examine some other approaches.

Let's consider the case with the proper transfer function so we can see the discontinuity. And we'll make the delay of the second step a not-nice value.

M = 2;
b = 14;
q = 20;
num = [.03 0 1];
G = tf(num,[M b q]);
delay = 1.96;

First, the exact symbolic solution

syms s t
ysym = ilaplace((0.03*s^2 + 1)/(M*s^2 + b*s + q)*(2/s - 2*exp(-delay*s)/s),s,t);

One option that comes to mind would be as follows

tempsys = G - tf(G.num{:},G.den{:},'InputDelay',delay);
[y1,t1] = step(tempsys,5); % allow the CST to pick the time vector
y1 = 2*y1; % gain of 2 on the input
figure
fplot(ysym,[0 5],'-o')
hold on
plot(t1,y1),grid

At a a distance, that seems pretty good. Zoom in

copyobj(gca,figure)

xlim([1.85 2.1])

We see that the autogenerated time vector from lsim does not include a point at the delay time (I was surprised by this), but lsim does seem to be computing a very accurate solution on both sides of the discontinuity.

tempsys has an internal delay

hasInternalDelay(tempsys)

ans = logical
   1

I think I've seen somewhere in the doc that lsim handles systems with internal delays differently than systems without internal delays. I wouldn't be surprised if the solution was computed accurately at the delay time, even if the returned solution doesn't include that point. But that's just speculation on my part.

If one were so inclined to really catch the effect of the step at the delay time, we can use the properties of superposition and time invariance.

First, define a time vector that has a point exactly (as close as we can get) at the delay time

t2 = linspace(0,delay,150);
t2 = [t2, (delay+t2(2)):t2(2):5];

Compute the step response with the gain

ystep = 2*step(G,t2);

They delayed step response is

ystepdelayed = [zeros(149,1) ; ystep(1:(numel(t2)-149))];

The total response is the difference between the step and the delayed step.

y2 = ystep - ystepdelayed;

Plot and zoom

figure
fplot(ysym,[0 5],'-o')
hold on
plot(t2,y2),grid
copyobj(gca,figure)

xlim([1.85 2.1])

Now y2 has a point at the delay time that matches the sym response at t = delay+.

We can manually add another point into y2 and t2 for the value at t = delay- to get a clean step at the delay

t3 = [t2(1:149) , delay, t2(150:end)];
y3 = [y2(1:149) ; ystep(150); y2(150:end)];

figure
fplot(ysym,[0 5],'-o')
hold on
plot(t3,y3),grid
copyobj(gca,figure)

xlim([1.85 2.1])

Answer 2

[Walter Roberson]

M=2;
b=14;
q=20;
num=[1];
G=tf(num,[M b q])

G =
 
          1
  -----------------
  2 s^2 + 14 s + 20
 
Continuous-time transfer function.

syms t
F = 2*(heaviside(t)-heaviside(t-2));
T = linspace(0, 3);
sympref('HeavisideAtOrigin', 0);
u = double(subs(F, t, T));

plot(T, u);
hold on

lsim(G, u, T); xlim([-.1 3.1]); ylim([-0.1 2.1]);