fplot and laplace transform

i have a system given in state space representation ( matrices A,B,C,D in the following code) and i need to plot the steady state output given the input vector u = [sin(t) 0]; i have written the following code:

A = [-4 1 0; 1 -2 -2; 0 2 -2];
B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];
s = sym('s');
G = C * inv(s*eye(size(A,1)) - A) * B + D;
t = sym('t');
u = [sin(t) 0];
U = laplace(u1);
Y = G*U1';
y = ilaplace(Y1);
fplot(y(1))
fplot(y(2))

however, Matlab doesn't print anything and gives me a warning; moreover, if i try to print the expression of y, Matlab gives me an "implicit" expression in the sense that it is not a function of the symbolic variable t, but something like ilaplace(function(s)); can someone help me?

[Paul]

Hi Luca,

The posted code has the line

U = laplace(u1)

but u1 isn't defined?

The system will only have a steady state output if it is stable. So stability should always be checked first.

Assuming the system is stable, are you supposed to plot only the steady state output? Or are you supposed to plot the full output in response to the input until the output reaches close to steady state?

If the former, there's really no need to go through the inverse Laplace transform if you've learned about the theory of steady state output to sinusoidal inputs. I guess it depends on the intent of the problem .....

[Luca Virecci Fana]

yes, you are right: these three lines

U = laplace(u1);
Y = G*U1';
y = ilaplace(Y1);

should be replaced with these three lines

U = laplace(u);
Y = G*U';
y = ilaplace(Y);

(i have copied here the code incorrectly);

anyway i am explicitely requested to compute the steady state response using the laplace and inverse laplace transform

[Paul]

i am explicitely requested to compute the steady state response using the laplace and inverse laplace transform

Does that mean that you have to compute Laplace transform of the output, then compute the inverse Laplace transform, then extract the steady state output, in closed form? If so, it will be difficult, if even feasible, to find the closed form, inverse Laplace transform of a fifth order response. Are you sure that's what you're supposed to do?

A = [-4 1 0; 1 -2 -2; 0 2 -2];
B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];
syms s t
G = C * inv(s*eye(size(A,1)) - A) * B + D;
u = [sin(t); 0];
U = laplace(u);
Y = simplify(G*U)

Y = 

y = ilaplace(Y)

y = 

Matlab can't find the full, closed form solution. But, if you look carefully you can extract the steady state solution if that's all you need.

[Luca Virecci Fana]

My exact assignment is:

"Calculate the output of the system (in steady-state) and corresponding amplitude amplification for the following input signal, using Laplace and Inverse Laplace Transforms: u(t) = [sin(t) 0]"

(the system is the one given with matrices A, B, C, D in the code)

Is it not correct to laplace transform the input, then multiply it by G to obtain the laplace transform of the output and then inverse laplace transform the output?

[Paul]

Yes, that would be correct, if it were feasible. The part of the problem statement about using Inverse Laplace Transforms is the part that's troubling. All you really need to solve this problem is G(s). You don't need the Laplace transform of U and you don't need the inverse Laplace transform of Y. So it's still not clear to me what the problem statement is really expecting you to do.

[Luca Virecci Fana]

I know that I could use the frequency response theorem to compute the steady state output and I should just compute the magnitude and the phase of the components of G at frequency 1 to do that, but this is the assignment

Answer 1

[Star Strider]

I am not certain what the problem is. The inversion appears to be correct, and the result is what I would expect it to be. Since Laplace transforms are defined for , it is necessary to define that as the second argument to fplot.

After the inversion, ‘y’ is a function of ‘t’ as expected, and the plot appears to be appropriate —

A = [-4 1 0; 1 -2 -2; 0 2 -2];
B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];
s = sym('s');
G = C * inv(s*eye(size(A,1)) - A) * B + D;
t = sym('t');
u = [sin(t) 0];
U = laplace(u.');
Y = G*U;
y = ilaplace(Y)

y = 

figure
fplot(y, [0 10])
grid
xlabel('t')
ylabel('System Output')
legend('y_1','y_2', 'Location','best')

Is this the result you want?

.

[Luca Virecci Fana]

why have you put u.' as argument of laplace and not just u? maybe that is my error

[Star Strider]

The input to a state space system needs to be a column vector (or matrix of column vectors depending on the system) or a scalar. Defining ‘u’ as a column vector as the input makes it conformable for matrix multiplication with ‘B’ and ‘D’. I assumed that is actually what you want it to be, so I transposed it.

.

[Luca Virecci Fana]

anyway, than you very much. it works!

then i should compute the amplitude amplification given by the system. using the definition, i have tried this:

u_2norm = sqrt(int(sum((abs(u)).^2),-Inf,Inf));

y_2norm = sqrt(int(sum((abs(y)).^2),-Inf,Inf));

ampl = y_2norm / u_2norm;

but it doesn't work; i think both u_2norm and y_2norm should be infinite, but i expect their ratio is not; anyway, if i try to print both of them, Matlab prints an expression and not Inf; do you know how to solve this? maybe is there a Matlab function to compute the amplitude amplification automatically?

[Luca Virecci Fana]

i have read the reply you have delete. you have obtained zero as amplification, but it sounds me weird; could it be because of the calculation of y_2norm?

[Star Strider]

My pleasure!

The ‘u_2norm’ is indeed infinite. In order to get a value for ‘y_2norm’, it is necessary to simplify it, and that is taking forever. I am doing that offline (giving it a 500 iteration limit, usually enough to produce a reasonable result), and will post the result when it finishes.

EDIT — (8 Apr 2023 at 17:37)

The simplification on MATLAB Online does not give a significantly simplified result for ‘y_2norm’ , however it calculates:

ampl =

0

so I suspect that the numeric value of ‘y_2norm’ is also infinite. I am runninig it offline as well on this computer, and (after 26 minutes) without calculating an actual value for ‘y_2norm’. (the result is also a symbolic expression) returned the same result for ‘ampl’.

If I get any further reults, I will post them as edits to this Comment.

.

[Luca Virecci Fana]

do you know why, if i try to evaluate y at a certain time (t=3 for example), using

subs(y, t, 3)

matlab returns an expression and not a numerical value?

[Star Strider]

The expression is being evaluated, however you need to use the vpa function to see the numeric result —

A = [-4 1 0; 1 -2 -2; 0 2 -2];
B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];
s = sym('s');
G = C * inv(s*eye(size(A,1)) - A) * B + D;
t = sym('t');
u = [sin(t) 0];
U = laplace(u.');
Y = G*U;
y = ilaplace(Y)                                         % Function

y = 

y_3 = subs(y, t, 3)                                     % Symbolic Solution

y_3 = 

y_3 = simplify(y_3, 500)                                % Simplified Symbolic Solution

y_3 = 

y_3 = vpa(y_3)                                          % Numeric Symbolic Solution

y_3 = 

format long
y_3 = double(y_3)                                       % Double Precision Solution (Can Be Used Outside Of The Symbolic Math Toolbox)

y_3 = 2×1

0.282757162323995 0.096723511065013

In other instances, the symbolic result is because MATLAB cannot find a numeric solution for a specific result, and so returns the symbolic expression. Sometimes, the simplify function can do this (with perhaps 500 or more permitted iterations), however other times, not (as in this instance).

.

[Walter Roberson]

If you need the full symbolic expression without the sum of roots in it, then see my comment to Paul's answer, where I link to code I posted before that will expand the expression. But the result will be messy, and will contain lots of terms like sqrt(sin(3)*(1-sqrt(5))^2*1i)/(sqrt(cos(7)^2+3i)) that very few people will be able to intuitively understand the purpose of.

Answer 2

[Paul]

Based on the comment thread in the question, I think the best you can do symbolically would be some something like this:

A = [-4 1 0; 1 -2 -2; 0 2 -2];
eig(A)

ans = 
  -4.2483 + 0.0000i
  -1.8759 + 1.8822i
  -1.8759 - 1.8822i

B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];
syms s t
G = C * inv(s*eye(size(A,1)) - A) * B + D;
u = [sin(t); 0];
U = laplace(u);
Y = simplify(G*U)

Y = 

y = ilaplace(Y)

y = 

If we look carefully at the two elements of y we see that each has terms in sin(t) and cos(t) and then a bunch of other stuff. That other stuff comes from the impulse response of the plant, which all decays to zero becasue the plant is stable (all the eigenvalues of A have negative real parts). We can see this more clearly by

vpa(y,5)

ans = 

Every exponential term will decay to zero in steady state. Hence, the steady state output are the terms that do not decay to zero, i.e.,

yss = [67/44*sin(t) - 3/44*cos(t); 15/22*sin(t)]

yss = 

Verify the answer based on mag and phase of G

G1 = subs(G,s,1j)

G1 = 

M = abs(G1)

M = 

P = angle(G1)

P = 

simplify(yss(1) - M(1,1)*sin(t + P(1,1)))

ans = 

0

[Walter Roberson]

At some point, I posted code to locate sum of root() and expand those, so it would be possible to get a version of y without the summation or root. It would be a bit ugly though.

https://www.mathworks.com/matlabcentral/answers/890822-how-can-i-integrate-a-function-which-has-three-variables-but-has-to-integrate-it-with-respect-to-onl#answer_759472

[Paul]

Would that be different than vpa(y) ?

[Walter Roberson]

The code should be able to get the full symbolic version of y.

Which, to be sure, would probably be mostly unreadable. The meaning of the full symbolic expression of roots of a cubic gets ugly. It is only useful for completionists, or in rare cases for proving complicated identities.

[Luca Virecci Fana]

how did you get yss? just looking at y expression and neglecting terms containing exponential terms?

[Paul]

In short, yes. Those terms can be neglected once you realize where they come from and that they will all decay to zero. Therefore, the steady state solution is composed of the sin(t) and cos(t) terms that don't decay to zero.

[Paul]

Because we know the system is stable, the steady state solution can be obtained thusly.

A = [-4 1 0; 1 -2 -2; 0 2 -2];
B = [1 1; 0 1; 1 0];
C = [1 0 1; 1 1 2];
D = [1 0; 0 0];

syms s t
u = [sin(t); 0];
U = laplace(u)

U = 

G = C * inv(s*eye(size(A,1)) - A) * B + D;
G = simplify(G);

yss = subs(G(:,1),s,1i)/(1i + 1i)*ilaplace(1/(s - 1i)) + subs(G(:,1),s,-1i)/(-1i - 1i)*ilaplace(1/(s + 1i));
yss = rewrite(yss,'sincos')

yss = 

Though I don't think it's needed to answer this specific question, here's an approach to get the full, symbolic form of y. It takes advantage of the fact that the eigenvalues of A are distinct and can be computed symbolically.

e = eig(sym(A))

e = 

Y = simplify(G*U)

Y = 

e = [e ; sym(1i) ; sym(-1i)];
y1 = 0;
y2 = 0;
num1 = numden(Y(1));
num2 = numden(Y(2));
for ii = 1:5
    y1 = y1 + subs(num1/prod(s-e(e~=e(ii))),s,e(ii))*ilaplace(1/(s-e(ii)));
    y2 = y2 + subs(num2/prod(s-e(e~=e(ii))),s,e(ii))*ilaplace(1/(s-e(ii)));
end
y1

y1 = 

y2

y2 = 

Verify that [y1;y2] == y as computed from ilaplace

y = ilaplace(Y);
figure
subplot(211);
fplot(real([y1-y(1) , y2 - y(2)]),[0 20])
subplot(212);
fplot(imag([y1-y(1) , y2 - y(2)]),[0 20])