Inverse laplace transform from system identification data
4 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Star Strider
2020-11-10
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 个评论
Star Strider
2020-11-11
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
2020-11-12
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!
.
更多回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Dynamic System Models 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!