Inverse laplace transform from system identification data

Question

martin Kostovcik 2020-11-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/643400-inverse-laplace-transform-from-system-identification-data

评论： Star Strider 2020-11-12

model.txt

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-11-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/643400-inverse-laplace-transform-from-system-identification-data#answer_540315

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

martin Kostovcik 2020-11-11

在 MATLAB Online 中打开

But I still have some troubles with this projekt because when I try that I still don't get correct number on output

but this output is from system identification toolbox and its look good but (not sure). (imput data pics.1)

then I estimate the model

(setup of identification toolbox pics.2)

(plot of model pics.3)

the Result of estimate :

But when I try estimate inverse laplace transformation by :

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)

I always get wrong output.

I set time to 50 sec. and output have to be around 160 degrees celsius.

>> syms s
>> F = (0.0003038)/(s^2 + 0.0122*s + 0.0001576 )
 
F =
 
2802060424796481/(9223372036854775808*(s^2 + (61*s)/5000 + 5814413732033251/36893488147419103232))
 
>> f = vpa(simplify(ilaplace(F), 'Steps',500))
 
f =
 
0.027688062224839641393809837346696*exp(-0.0061*t)*sin(0.01097223769337868513904354310076*t)
 
>> f_fcn = matlabFunction(f)
f_fcn =
  function_handle with value:
    @(t)exp(t.*-6.1e-3).*sin(t.*1.097223769337869e-2).*2.768806222483964e-2
>> t = 50
t =
    50
>> 0.027688062224839641393809837346696*exp(-0.0061*t)*sin(0.01097223769337868513904354310076*t)
ans =
    0.0106

Ans have to be around 160 but output is still wrong.

But I don't know why its like this can anyone help me?

martin Kostovcik 2020-11-11

在 MATLAB Online 中打开

model.txt

Oke I try to explain you all what you need

I have my data from PLC system output (exemple data) I attach the model.txt there are the data as first ist a time secound its a temperature

load model.txt
time = model(:,1);
temp = model(:,2);

then I have to imput data into a system identification toolbox (I attach PrintScreen you what steps I following)

Then I select the estimate --> and select the Transfer funcion models

I set that like :

Then press the estimate then I get the result as transfer function :

then I try to solve this transfer funcion (my last try was by dsolve)

f = dsolve('D2y+0.0122*Dy+0.0001576*y = 0.0003038','Dy(1)=1,y(102.47)=102.47')

then I alway try to test that like:

syms t 
model= double(subs(f,t,50))

the result of this steps are :

f = dsolve('D2y+0.0122*Dy+0.0001576*y = 0.0003038','Dy(1)=1,y(102.47)=102.47')
%Warning: Support of character vectors and strings will be removed in a future release. Use sym objects to
%define differential equations instead. 
%> In dsolve (line 126) 
 
f =
 
(exp(-(61*t)/10000)*(197000000*cos((10247*12039^(1/2))/1000000)*exp(61/10000)*sin((12039^(1/2)*t)/10000) - 197000000*sin((10247*12039^(1/2))/1000000)*exp(61/10000)*cos((12039^(1/2)*t)/10000) + 120821724*cos(12039^(1/2)/10000)*exp(625067/1000000)*sin((12039^(1/2)*t)/10000) - 120821724*sin(12039^(1/2)/10000)*exp(625067/1000000)*cos((12039^(1/2)*t)/10000) + 2316475*cos(12039^(1/2)/10000)*sin((10247*12039^(1/2))/1000000)*exp((61*t)/10000) - 2316475*cos((10247*12039^(1/2))/1000000)*sin(12039^(1/2)/10000)*exp((61*t)/10000) + 37975*12039^(1/2)*cos(12039^(1/2)/10000)*cos((10247*12039^(1/2))/1000000)*exp((61*t)/10000) + 37975*12039^(1/2)*sin(12039^(1/2)/10000)*sin((10247*12039^(1/2))/1000000)*exp((61*t)/10000) + 1980684*12039^(1/2)*cos(12039^(1/2)/10000)*exp(625067/1000000)*cos((12039^(1/2)*t)/10000) + 1980684*12039^(1/2)*sin(12039^(1/2)/10000)*exp(625067/1000000)*sin((12039^(1/2)*t)/10000)))/(19700*(61*cos(12039^(1/2)/10000)*sin((10247*12039^(1/2))/1000000) - 61*cos((10247*12039^(1/2))/1000000)*sin(12039^(1/2)/10000) + 12039^(1/2)*cos(12039^(1/2)/10000)*cos((10247*12039^(1/2))/1000000) + 12039^(1/2)*sin(12039^(1/2)/10000)*sin((10247*12039^(1/2))/1000000)))
 

then next step is a test my model we compare the model with the measured values

(from model.txt) by the time of 50 secound for exemple :

>> syms t 
>> model= double(subs(f,t,50))
model =
  131.1831

But result in 50 secound must be 158.33 degrees celsius but was 131.1831 degrees celsius

I try to understand why it is like that first I think its must be something like

amplification of data but when I try to find some constant for amplification but its always different number when I change a Time of model

>> model= double(subs(f,t,50))

where 158.33 is required temperature and model return the model temperature to estimate the Constant

    resultConstant = 158.33/model
 resultConstant = 1.3049

But for t = 10 sec its

  >> model= double(subs(f,t,10))
  resultConstant = 117.73/model
  resultConstant = 0.9703

but this is to the dead end.

Now I need the estimate math. function of the model and when I solve this some how I apply that for a next 8 model beacuse I have to do estimate models for different spectrum of temperatures for exemple this one is for 100-200 degrees celsius but only for a active part of heating.

Forward thanks you for your help on this project .

Star Strider 2020-11-11

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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:

Success!

.

请先登录，再进行评论。

Answer 2

martin Kostovcik 2020-11-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/643400-inverse-laplace-transform-from-system-identification-data#answer_542293

Hi, thank you for your time spent on this project

I appreciate it when I find solution I will publish that. I hope I can prove it.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Inverse laplace transform from system identification data

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Inverse laplace transform from system identification data

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论