Inverse laplace transform from system identification data
4 views (last 30 days)
Show older comments
martin Kostovcik
on 10 Nov 2020
Commented: Star Strider
on 12 Nov 2020
Hello guys, I really need help with my project how can I get Inverse laplace transform of my transfer function, but for the first: I use system identification toolbox when I put my data into app then i get transfer function of my process(furnace)
Its looks like :
1.204 (+/- 1.7) s + 0.002519 (+/- 0.003312)
----------------------------------------------------
s^2 + 0.1662 (+/- 0.2069) s + 0.008011 (+/- 0.01117)
then I try do Inverse laplace transform because I need math. model of temperature but when I use the function "ilaplace()"
like :
sym s;
f = (1.204*s + 0.002519) / (s^2 + 0.1662*s + 0.008011)
ilaplace(f)
then I have to get math. funcion of temperature but its not correct bacause when I put my time into function I don´t got correct output.
I always get some different numbers like I have get, and I still stucked on this step.
(301*exp(-(831*t)/10000)*(cos((124455849802476899811^(1/2)*t)/335544320000) - (17570055355847101387741*124455849802476899811^(1/2)*sin((124455849802476899811^(1/2)*t)/335544320000))/80447337606977714844798893948928))/250
I send a exemple data(time,temp.) for a test if you want but I think I do all correct can anyone help my with this?
I dont know if I do something wrong or I just follow wrong steps.
0 Comments
Accepted Answer
Star Strider
on 10 Nov 2020
First, simplify the result, use vpa to eliminate the fractions in the symbolic result, then use matlabFunction to produce an anonymous function version:
syms s
F = (1.204*s + 0.002519) / (s^2 + 0.1662*s + 0.008011)
f = vpa(simplify(ilaplace(F), 'Steps',500))
f_fcn = matlabFunction(f)
producing:
f =
1.204*exp(-0.0831*t)*(cos(0.033247405913845380174138784260994*t) - 2.4365151229809367326763139899363*sin(0.033247405913845380174138784260994*t))
f_fcn =
function_handle with value:
@(t)exp(t.*-8.31e-2).*(cos(t.*3.324740591384538e-2)-sin(t.*3.324740591384538e-2).*2.436515122980937).*1.204
or (with a bit of editing):
f_fcn = @(t)exp(t.*-8.31e-2).*(cos(t.*3.324740591384538e-2)-sin(t.*3.324740591384538e-2).*2.436515122980937).*1.204;
.
5 Comments
Star Strider
on 11 Nov 2020
I have been working on this for a few hours, with several different transfer function and state-space models, and I cannot get it to work, even though the compare function output says it should, and gives a good fit. (The compare, sim, and lsim functions filter the input data using the transfer function coefficients, at least as I understand the documentation.)
The plot of the imaginary part of the transfer function definitely implies (to me at least) that the transfer function has 2 poles and 2 zeros, and compare indicates a good fit to the data.
I am stopping here for now. If I come up with a new approach, I will try it, otherwise I will delete my Answer in a few hours. I have no idea why the inverse Laplace transform of the transfer function is not working.
The code I used:
model = load('model.txt');
time = model(:,1);
temp = model(:,2);
u = time;
figure
plot(time, temp)
grid
Ts = mean(diff(time));
Fs = 1/Ts;
Fn = Fs/2;
L = numel(time);
mean_temp = mean(temp);
FTtemp = fft(temp)/L;
FTu = fft(u)/L;
FTtf = FTtemp./FTu;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTtf(Iv))*2)
grid
title('Fourier Transform of the Transfer Function')
figure
plot(Fv, imag(FTtf(Iv)))
grid
title('Im\{FTtf\}')
idtemp = iddata(temp, u, Ts);
tf_sysz = tfest(idtemp, 4, 4, 'Ts',Ts) % Discrete Transfer Fucntion
tf_syss = tfest(idtemp, 2, 2) % Continuous Transfer Function
ss_syss = ss(tf_syss) % Continuous State-Space
figure
compare(idtemp, tf_sysz)
figure
compare(idtemp, tf_syss)
figure
pzmap(tf_syss)
ylim([-1 1]*0.001)
syms s t
F = vpa(poly2sym(tf_syss.Numerator, s), 5) / vpa(poly2sym(tf_syss.Denominator, s), 5)
F = vpa(partfrac(F), 5)
f = vpa(simplify(ilaplace(F,s,t), 'Steps',500), 5)
% f = subs(f,{dirac(t)},{heaviside(t)})
f_fcn = matlabFunction(f)
Out = f_fcn(50)
figure
plot(time, f_fcn(time))
grid
.
Star Strider
on 12 Nov 2020
I just now figured it out!
The input is a ramp function:
u(t) = k*t
the Laplace transform of it being:
U(s) = k/s^2
and the initial condition is:
f(0) = temp(1)
the Laplace transform of that being:
F(0) = temp(1)/s
multiplying the identified Laplace transform by ‘U(s)’ and adding the initial condition:
syms s t
F = vpa(poly2sym(tf_syss.Numerator, s), 5) / vpa(poly2sym(tf_syss.Denominator, s), 5)
F = F * 1/s^2 + temp(1)/s
F = vpa(F, 5)
f = vpa(simplify(ilaplace(F,s,t), 'Steps',500), 5)
f_fcn = matlabFunction(f)
Out = f_fcn(50)
figure
plot(time, f_fcn(time))
grid
xlabel('Time (s)')
ylabel('Temperature (°C)')
produces (with vpa):
F =
102.47/s + (103.35*s^2 + 51.099*s + 0.27354)/(s^2*(s^2 + 17.931*s + 0.92082))
and its inverse:
f =
0.29706*t - 5.6367*exp(-17.88*t) - 44.071*exp(-0.051501*t) + 152.18
producing:
Out =
163.6751
and:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/411618/image.png)
Success!
.
More Answers (1)
See Also
Categories
Find more on Dynamic System Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!