lsqcurvefit for multiple variables optimization

Hi experts,

I have a set of data which is in the cartesian coordinates (xOy) but needs to be moved a distance (say xc, yc) and rotate an angle (theta) (still cartesian coordinates) to fitting the theoretical value calculated from ode45 function.

In previous code, I transform data by manual procedure and then use lsqcurvefit to fitting the data, with B and L are optimization variables.

lb=[0, 1];
ub=[2, 2.8];
p0=[0.5 2.7];
p=lsqcurvefit(@(p,y_exp) pendant_bubble_arc_d(p,y_exp),p0,y_exp,x_exp,lb,ub);

function x_model = pendant_bubble_arc_d(p, y_exp) 

global l_exp
B = p(1); 
L=p(2);
% initial conditions
y0= zeros(5,1);
sspan = [0 l_exp/L]; y_solver=[];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',5e-2);
% call the ODE solver
[~, y_solver] = ode45(@(t,y) pendant_bubble_arc_f(t,y,B), sspan, y0, options);

% model prediction 

x_model = interp1( y_solver(:,2)*L, y_solver(:,1)*L, y_exp,'spline',2); 

return;

The code run successfully and result is correct. However, my PI also needs xc, yc, theta are optimization values for lsqcurvefit function.

And this my idea for the problem:

lb=[0, 1, 0, 0, 0];
ub=[2, 2.8, 5, 5, 5];
p0=[0.5, 2, 2.5, 0.8, 1 ];

p= lsqcurvefit(@(p,y_exp) lsq_function(p,y_exp),p0,y_exp,x_exp,lb,ub);

function x_model = lsq_function(p, y_exp) 

global l_exp
B = p(1); 
L=p(2);
yc=p(4);
xc=p(3);
theta=p(5);
% initial conditions
y0= zeros(5,1);
sspan = [0 l_exp/L]; y_solver=[];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',5e-2);
% call the ODE solver
[~, y_solver] = ode45(@(t,y) YL_function(t,y,B), sspan, y0, options);

% transforming coordinates
R=[cosd(theta) sind(theta); -sind(theta) cosd(theta)]; 
y_trans(:,1)=y_solver(:,1)+xc;
y_trans(:,2)=y_solver(:,2)+yc;
y_rotate=R*y_trans';

% model prediction 

x_model = interp1(y_rotate(:,2)*L, y_rotate(:,1)*L, y_exp,'spline',2); 
return;

Although the code does not have any error, but it does not return the values I expect. The code just return value slightly change from initial guess (few percents). I just learning this and do not have sufficient knowledge about matlab. So if you have any ideas how I should modify the code, please let me know.

Thank you!

P/s: if the question is unclear, please pointed it and I will explain.

function f=YL_function(s,y,B)

r = y(1); % radial distance
h = y(2); % height
theta = y(3); % angle
V = y(4); % volume
A = y(5); % surface area

f(1) = cos(theta);
f(2) = sin(theta);
if s==0 
    f(3) = 1/B;
else
    f(3) = 2/B-h-sin(theta)./r;
end
f(4) = pi*r^.2*sin(theta);
f(5) = 2*pi*r;

f=f';

This is function YL_function (also pendant_bubble_arc_f in upper code).

Answer 1

[Torsten]

We need your "YL_function" to see what's happening.

I wonder whether your transformation is correct.

I would have thought

y_rotate = [cosd(theta) -sind(theta); sind(theta) cosd(theta)]*[y_solver(:,1),y_solver(:,2)].' + [xc;yc];

but maybe I'm mistaken.

[Nhat Nguyen]

I tried but it is the same.

I also add the YL_function to the post, please see the edited content.

Thank you!

[Torsten]

In your ODE function, you use theta in radians, but in your transformation matrix R, theta must be in degrees. Is that correct though ?

And what is l_exp ? It's not set in the code you supplied.

And your experimental data y_exp are missing to run the code.

I guess your transformation must read

y_rotate = [cos(theta) -sin(theta); sin(theta) cos(theta)]*[y_solver(:,1)-xc,y_solver(:,2)-yc].' + [xc;yc];

but I cannot tell without further information.

The above defines a rotation of (y(1),y(2)) around (xc,yc) by an angle of theta.

[Nhat Nguyen]

Sorry for the unclear, the theta in ODE is a different variable compared to theta inlsq_function.

I have attached the data file, please see it.

l_exp is arclength of the experiment curve xy, which about 5.41 for this set of data

[Torsten]

global l_exp
l_exp = 5.41;
x_exp = readmatrix('x.txt');
y_exp = readmatrix('y.txt');
hold on
plot(y_exp,x_exp)

lb=[0, 1, 0, 0, 0];
ub=[2, 2.8, 5, 5, 5];
p0=[0.5, 2, 2.5, 0.8, 1 ];
p= lsqcurvefit(@(p,y_exp) lsq_function(p,y_exp),p0,y_exp,x_exp,lb,ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x_model = lsq_function(p,y_exp);
plot(y_exp,x_model)
hold off

function x_model = lsq_function(p, y_exp) 
global l_exp
B = p(1); 
L = p(2);
yc = p(4);
xc = p(3);
theta = p(5);
% initial conditions
y0= zeros(5,1);
sspan = linspace(0,l_exp/L,numel(y_exp));
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',5e-2);
% call the ODE solver
[~, y_solver] = ode45(@(t,y) YL_function(t,y,B), sspan, y0, options);
% transforming coordinates
R=[cosd(theta) sind(theta); -sind(theta) cosd(theta)]; 
y_trans(:,1)=y_solver(:,1)+xc;
y_trans(:,2)=y_solver(:,2)+yc;
y_rotate=R*y_trans';
y_rotate = y_rotate.';
x_model = interp1(y_rotate(:,2)*L, y_rotate(:,1)*L, y_exp,'spline',2);
end

function f=YL_function(s,y,B)
r = y(1); % radial distance
h = y(2); % height
theta = y(3); % angle
V = y(4); % volume
A = y(5); % surface area
f(1) = cos(theta);
f(2) = sin(theta);
if s==0 
    f(3) = 1/B;
else
    f(3) = 2/B-h-sin(theta)./r;
end
f(4) = pi*r^.2*sin(theta);
f(5) = 2*pi*r;
f=f';
end

[Nhat Nguyen]

Amazing, thank you!

Btw can you explain why changing sspan is make this code working instead of old expression?

[Torsten]

Btw can you explain why changing sspan is make this code working instead of old expression?

This didn't make the code work. The main point was to transpose y_rotate:

y_rotate = y_rotate.';

[Nhat Nguyen]

Thank you very much!