Inverse laplace transform from system identification data

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.

Answer 1

[Star Strider]

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;

.

[martin Kostovcik]

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?

[Star Strider]

I would need your data (attach it either in the original format or as a .mat file to your original Question or to a Comment), and to see your code to figure out what the problem is.

I do not have enough information to proceed.

[martin Kostovcik]

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]

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]

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]

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.