Piece-wise non-linear fitting with fminsearch - issue with quadratic function

Hello everybody.

I have a series of experimental data (tabella_raw.txt, column #4) for which I want to find the best interpolation function. If you have a look at the experimental data named 'qse' you can see that the best option is a piece-wise interpolation: three functions, in particular, the first part quadratic (from 0 to 1 h), as well as the second (from 1 to 2), and the third is hyperbolic (from 2 to end). The function must be continuous. I tried to find at first some suitable expressions for the three parts and then, by adopting fminsearch that minimizes the error between the two curves, and I got some kind of interpolation curve (the red one). In blue you can see the experimental data.

However, there is something wrong or better to say, something I don't like. At 1 and 2 (yellow dots) I want to have the maximum and the minimum values of the function. Is there any possibile way to set up some paramters contraints to obtain something like the second picture?

I attach the matlab code, as well.

For any doubts, please write me.

Thanks in advance.

Maja

[Mathieu NOE]

just for my own fun, I wanted to see if we could put a first order model (so exponential functions) to fit your data (you say "interpolation" but you meant model fitting probably)

so this is it , and obviously it coud be improved (usind an optimizer) , especially on the 3rd segement

code FYI

clc
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                    INTERPOLAZIONE CURVE DINAMICO                        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


data=readmatrix('tabella_raw.txt');

tsi=data(:,1); tse=data(:,2); qsi=data(:,3); qse=data(:,4); 
ti=data(:,5); to=data(:,6);

%__________________________________________________________________________
%                      INTERPOLAZIONE CURVE
%__________________________________________________________________________

% usare la scala x ogni ora per fitting corretto
n  = 500;
ind = (1:n);
x=(ind/60)'; 
dT = 10; % delta T impulso (10°)
qse = qse(ind);

%% flusso termico esterno qse
b1 = -2; % exponant term
% first segment
[m1,i1] = max(qse);
y1 = m1*(1-exp(b1*x));
y1(i1:end) = NaN;

% second segment
[m2,i2] = min(qse);
[m3,i3] = max(y1);
b2 = -2; % exponant term

y2 = m3 - (m1-m2)*(1-exp(b2*(x-x(i1))));
y2(1:i1) = NaN;
y2(i2:end) = NaN;

% third segment
b3 = -3; % exponant term
y3 = m2*exp(b3*(x-x(i2)));
y3(1:i2) = NaN;

plot(x,qse)
hold on
plot(x,y1) 
plot(x,y2) 
plot(x,y3)
hold off

[Maja]

Thanks a lot for the reply.

I'll try to implement exponential functions, especially for the first and the second part which seem perfect as fitting curves for the experimental data (and yes, I meant fitting, not interpolating function!!).

[Mathieu NOE]

hello again

still playing with it. Found a better fit function for the 3rd segment

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                    INTERPOLAZIONE CURVE DINAMICO                        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data=readmatrix('tabella_raw.txt');

tsi=data(:,1); tse=data(:,2); qsi=data(:,3); qse=data(:,4); 
ti=data(:,5); to=data(:,6);

%__________________________________________________________________________
%                      INTERPOLAZIONE CURVE
%__________________________________________________________________________

% usare la scala x ogni ora per fitting corretto
n  = 500;
ind = (1:n);
x=(ind/60)'; 
dT = 10; % delta T impulso (10°)
qse = qse(ind);

%% flusso termico esterno qse
b1 = -2.1; % exponant term
% first segment
[m1,i1] = max(qse);
y1 = m1*(1-exp(b1*x));
y1(i1:end) = NaN;

% second segment
[m2,i2] = min(qse);
[m3,i3] = max(y1);
b2 = -2.4; % exponant term
xx = x-x(i1);
y2 = m3 - (m1-m2)*(1-exp(b2*xx));
y2(1:i1-1) = NaN;
y2(i2:end) = NaN;

% third segment
b3 = -1; % exponant term
b4 = -1; % exponant term
xx = x-x(i2);
y3 = m2*exp(b3*xx + b4*xx.^0.25);
y3(1:i2-1) = NaN;

plot(x,qse)
hold on
plot(x,y1) 
plot(x,y2) 
plot(x,y3)
hold off

[Maja]

Super, it is even better! Thanks a lot.

How would you do in order to have a continous function at x = 1 h and x = 2 h ?

[Maja]

clc

clear all

close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %

% INTERPOLAZIONE CURVE DINAMICO %

% %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data=readmatrix('tabella_raw.txt');

tsi=data(:,1); tse=data(:,2); qsi=data(:,3); qse=data(:,4);

ti=data(:,5); to=data(:,6);

%__________________________________________________________________________

% INTERPOLAZIONE CURVE

%__________________________________________________________________________

% usare la scala x ogni ora per fitting corretto

x=((0:1:length(tsi(:,1))-1)/60)';

dT = 10; % delta T impulso (10°)

%% flusso termico interno qsi

output1=stima_parametri_int(x,qsi); % altezza, curvatura, ..., ultimo parametro è offset

fun = @(params, x) params(1)*exp(-params(2)*(log(x/params(3)).^2))+params(4);

plot(x,qsi)

hold on

plot(x,fun(output1,x))

hold off

% non entra in gioco offset nel calcolo della U

U_I_int=exp(1/(output1(2)*4))*pi^(1/2)*output1(1)*output1(3)/(output1(2)^(1/2))/dT

U_I_int_a=(sum(fun(output1,x))-output1(4)*length(x(:,1)))/60/dT

U_I_int_exp=(sum(qsi)-output1(4)*length(qsi))/60/dT

% %% integrale

% x2 = linspace(0, 49, 50);

% p_a=0;

% p_ex=0;

%

% for i = 1: length(x)

% int_a(i)= p_a + fun(output1,x(i))/60-output1(4)/60;

% int_ex(i)= p_ex + qsi(i)/60-output1(4)/60;

% p_a= int_a(i);

% p_ex=int_ex(i);

% end

% int_a_int = int_a;

% int_ex_int = int_ex;

%

% %mostro nel grafico punti ogni ora

% for i = 2: 50

% y_a(i)=fun(output1,x((i-1)*60));

% y_ex(i)=qsi((i-1)*60);

% end

%

%

% p_a=0;

% p_ex=0;

% for i = 2: 50

% int_punti_ex(i)=p_ex + y_ex(i)-output1(4);

% int_punti_a(i)=p_a + y_a(i)-output1(4);

%

% p_ex=int_punti_ex(i);

% p_a=int_punti_a(i);

% end

%

%

% plot(x, int_a_int, LineWidth=1, Color='b')

% hold on

% plot(x, int_ex_int, LineWidth=1, Color='r')

% hold on

% plot(x2, int_punti_ex, '--o', Color='r', MarkerSize=5)

% hold on

% plot(x2, int_punti_a, '--o', Color='b', MarkerSize=5)

% hold off

%% flusso termico esterno qse

output2=stima_parametri_est(x,qse); % m1 b1 offset m2 b2 b3 b4 %

fun2 = @(params, x) (params(1)*(1-exp(params(2).*x))+params(3)).*(x>=0 & x<=1)+...

(params(1)*(1-exp(params(2)))+params(3)-(params(1)-params(4))*(1-exp(params(5).*(x-1)))).*(x>1 & x<=2)+...

(params(4)*exp(params(6).*(x-2)+params(7).*(x-2).^0.25)+params(3)).*(x>2);

x2 = linspace(0, 49, 50);

for i = 2: 50

y_ex(i)=qse((i-1)*60);

y_a(i)=fun2(output2,x((i-1)*60));

end

plot(x,qse)

hold on

plot(x,fun2(output2,x))

% hold on

% plot(x2,y_ex)

% hold on

% plot(x2,y_a)

hold off

xlim([0, 24]);

U_I_est_a=(sum(fun2(output2,x))-output2(3)*length(x(:,1)))/60/dT

U_I_est_exp=(sum(qse)-output2(3)*length(qse))/60/dT

% %% integrale

% p_a=0;

% p_ex=0;

%

% for i = 1: length(x)

% int_a(i)= p_a + fun2(output2,x(i))/60-output2(3)/60;

% int_ex(i)= p_ex + qse(i)/60-output2(3)/60;

% p_a= int_a(i);

% p_ex=int_ex(i);

% end

% int_a_int = int_a;

% int_ex_int = int_ex;

%

% %mostro nel grafico punti ogni ora

% for i = 2: 50

% y_a(i)=fun2(output2,x((i-1)*60));

% y_ex(i)=qse((i-1)*60);

% end

%

%

% p_a=0;

% p_ex=0;

% for i = 2: 50

% int_punti_ex(i)=p_ex + y_ex(i)-output2(3);

% int_punti_a(i)=p_a + y_a(i)-output2(3);

%

% p_ex=int_punti_ex(i);

% p_a=int_punti_a(i);

% end

%

%

% plot(x, int_a_int, LineWidth=1, Color='b')

% hold on

% plot(x, int_ex_int, LineWidth=1, Color='r')

% hold on

% plot(x2, int_punti_ex, '--o', Color='r', MarkerSize=5)

% hold on

% plot(x2, int_punti_a, '--o', Color='b', MarkerSize=5)

% hold off

%% Integrale considerando la somma delle funzioni

funzx=fun2(output2,x)-output2(3);

funzy=fun(output1,x)-output1(4);

U_somma=sum(funzx+funzy)/60/dT/2;

[Maja]

function [optimal_params]=stima_parametri_int(x,funzy)

fun = @(params, x) params(1)*exp(-params(2)*(log(x/params(3)).^2))+params(4);

error_function = @(params) sum((fun(params, x) - funzy).^2);

initial_guess = [1, 1, 1, 1];

optimal_params = fminsearch(error_function, initial_guess);

end

[Maja]

function [optimal_params]=stima_parametri_est(x,funzy)

fun = @(params, x) (params(1)*(1-exp(params(2).*x))+params(3)).*(x>=0 & x<=1)+...

(params(1)*(1-exp(params(2)))+params(3)-(params(1)-params(4))*(1-exp(params(5).*(x-1)))).*(x>1 & x<=2)+...

(params(4)*exp(params(6).*(x-2)+params(7).*(x-2).^0.25)+params(3)).*(x>2);

error_function = @(params) sum((fun(params, x) - funzy).^2);

initial_guess = [8.6 -2.1 -0.035 -4.7 -2.4 -1 -1];

optimal_params = fminsearch(error_function, initial_guess);

end

[Mathieu NOE]

hello

just regarding your comment about having one single continuous fitting function

this is what I could do , simply by using a sigmoid function to apply like a weight window (could be done with some kind of step function)

clc
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                    INTERPOLAZIONE CURVE DINAMICO                        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data=readmatrix('tabella_raw.txt');

tsi=data(:,1); tse=data(:,2); qsi=data(:,3); qse=data(:,4); 
ti=data(:,5); to=data(:,6);

%__________________________________________________________________________
%                      INTERPOLAZIONE CURVE
%__________________________________________________________________________

% usare la scala x ogni ora per fitting corretto
n  = 500;
ind = (1:n);
x=(ind/60)'; 
dT = 10; % delta T impulso (10°)
qse = qse(ind);

%% flusso termico esterno qse
% sigmoid a parameter
a_sig = -1000; % a large negative value so the sigmoid is like a steped window (sharp transition)
% first segment
[m1,i1] = max(qse);
b1 = -2.1; % exponant term
y1 = m1*(1-exp(b1*x));
x1 = x(i1);
s1 = mysigmoid(x,x1,a_sig);
y1 = y1.*s1;

% second segment
[m2,i2] = min(qse);
[m3,i3] = max(y1);
b2 = -2.4; % exponant term
y2 = m3 - (m1-m2)*(1-exp(b2*(x-x1)));
x2 = x(i2);
s2 = mysigmoid(x,x2,a_sig);
y2 = y2.*(1-s1).*s2;


% third segment
b3 = -1; % exponant term
b4 = -1; % exponant term
y3 = m2*exp(b3*(x-x2) + b4*(x-x2).^0.25);
y3 = y3.*(1-s2);

y = y1+y2+y3;

plot(x,qse)
hold on
plot(x,y) 
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = mysigmoid(x,c,a)
% sigmoid evaluates a simple sigmoid function along x: 
% 
%         ______1________
%     y =        -a(x-c)
%          1 + e^
% 
%% Syntax 
% 
% y = sigmoid(x)
% y = sigmoid(x,c)
% y = sigmoid(x,c,a)
% 
 
y = 1./(1 + exp(-a.*(x-c)));

end

[Maja]

Super! Thanks again.

[Mathieu NOE]

tada !

now the full monty with fminsearch at the rescue

see answer section

Answer 1

[Mathieu NOE]

as announced above in the teaser , please find now the full code with fminsearch for the fine tuning

I tooked a while to make it work as I wanted but I am happy to have achieved this (in parallel to my official work today :))

clc
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                    INTERPOLAZIONE CURVE DINAMICO                        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data=readmatrix('tabella_raw.txt');

tsi=data(:,1); tse=data(:,2); qsi=data(:,3); qse=data(:,4); 
ti=data(:,5); to=data(:,6);

%__________________________________________________________________________
%                      INTERPOLAZIONE CURVE
%__________________________________________________________________________

% usare la scala x ogni ora per fitting corretto
n  = 500;
ind = (1:n);
x=(ind/60)'; 
dT = 10; % delta T impulso (10°)
qse = qse(ind);

%% flusso termico esterno qse

%% step1 : best expo fit of first segment to compute exactly t1 with best accuracy
[m1,i1] = max(qse); % initial guess (i1) where is t1 , but this is not yet accurate enough
[k, yInf, y0, ~] = fitExponential(x(1:i1), qse(1:i1));
yFit            = yInf + (y0-yInf) * exp(-k*(x-x(1)));
% refined t1 / x1 value (important for the step2 good performance
ind = find(yFit>qse);
ind = ind(ind>i1);
ind = ind(1);
x1 = x(ind); % or t1 if you prefer

%% step2 : do the general equation fit
% sigmoid a parameter (fixed paramter)
a_sig = -1000; % a large negative value so the sigmoid is like a steped window (sharp transition); -1000 is fine

% compute some vlaues we need to supply as IC for fminsearch
[m1,i1] = max(qse);
s1 = mysigmoid(x,x1,a_sig);
[m2,i2] = min(qse);
x2 = x(i2); % or t2 if you prefer
s2 = mysigmoid(x,x2,a_sig);

% curve fit using fminsearch
smooth_factor = 0.1; % 1 is too large , 0.1 is good
b1 = -2; % exponant term
b2 = -2; % exponant term
b3 = -1; % exponant term
b4 = -1; % exponant term

f = @(a,b,c,d,e,f,g,h,l,m,x) a*(1-exp(b*x)).*s1 + (c - d*(1-exp(e*(abs(x-x1)).^m))).*(1-s1).*s2 + (f*exp(g*(x-x2) + b4*(abs(x-x2)).^l)).*(1-s2);
obj_fun = @(p) norm(f(p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8),p(9),p(10),x)-qse) + smooth_factor*sum(abs(gradient(f(p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8),p(9),p(10),x))));
IC = [m1 b1 m1 (m1-m2) b2 m2 b3 b4 0.25 0.75];
sol = fminsearch(obj_fun, IC);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
d_sol = sol(4);
e_sol = sol(5);
f_sol = sol(6);
g_sol = sol(7);
h_sol = sol(8);
l_sol = sol(9);
m_sol = sol(10);
yfit = f(a_sol,b_sol,c_sol,d_sol,e_sol,f_sol,g_sol,h_sol,l_sol,m_sol, x);

plot(x,qse)
hold on
plot(x,yfit) 
hold off
legend('data','fminsearch fit');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [k, yInf, y0, yFit] = fitExponential(x, y)

% FITEXPONENTIAL fits a time series to a single exponential curve. 
% [k, yInf, y0] = fitExponential(x, y)
%
% The fitted curve reads: yFit = yInf + (y0-yInf) * exp(-k*(x-x0)).
% Here yInf is the fitted steady state value, y0 is the fitted initial
% value, and k is the fitted rate constant for the decay. Least mean square
% fit is used in the estimation of the parameters.
% 
% Outputs:
% * k: Relaxation rate
% * yInf: Final steady state
% * y0: Initial state
% * yFit: Fitted time series
% 
% Last modified on 06/26/2012
% by Jing Chen

% improve accuracy by subtracting large baseline
yBase           = y(1);
y               = y - y(1); 

fh_objective    = @(param) norm(param(2)+(param(3)-param(2))*exp(-param(1)*(x-x(1))) - y, 2);

initGuess(1)    = -(y(2)-y(1))/(x(2)-x(1))/(y(1)-y(end));
initGuess(2)    = y(end);
initGuess(3)    = y(1);
param           = fminsearch(fh_objective,initGuess);

k               = param(1);
yInf            = param(2) + yBase;
y0              = param(3) + yBase;

yFit            = yInf + (y0-yInf) * exp(-k*(x-x(1)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = mysigmoid(x,c,a)
% sigmoid evaluates a simple sigmoid function along x: 
% 
%         ______1________
%     y =        -a(x-c)
%          1 + e^
% 
%% Syntax 
% 
% y = sigmoid(x)
% y = sigmoid(x,c)
% y = sigmoid(x,c,a)
% 
 
y = 1./(1 + exp(-a.*(x-c)));

end