Combine lsqcurvefit with fsolve

Hi,

my goal is to optimize the length L0, L1, L2 and L3 of the bars of the following mechanism so that the theta2 angle follows the trend of my input data:

This is the mechanism:

And this is the data that I want to fit:

theta1 = linspace(0,pi,n); 
theta2 = [0.53,0.40,0.32,0.28,0.27,0.27,0.28,0.29,0.31,0.34,0.36,0.40,0.43,0.47,0.50,0.55,0.59,0.64,0.69,0.75]; %(data to fit)

which they have this trend:

My strategy was to mainly use the lsqcurvefit function and have the kinematic closure equation solved recursively using the fsolve function:

clc;
clear all;

n = 20;
theta1 = linspace(0,pi,n); 
theta2 = [0.53,0.40,0.32,0.28,0.27,0.27,0.28,0.29,0.31,0.34,0.36,0.40,0.43,0.47,0.50,0.55,0.59,0.64,0.69,0.75]; %(data to fit)
L0 = [0,0,0,0]; % Initial lenght of the bars, in the order: L0,L1,L2 and L3

for i = 1:n
    y0 = [0.5,1.1]; % Initial guess of theta2 and theta3
    
    options = optimset('display', 'off'); 

    y(i,:) = fsolve(@(y)fourbar(y,theta1(i)),y0,options);
end

%Objective: find the L0, L1, L2 and L3 optimized to make the function...
...interpolating the data (theta1,theta2)
L = lsqcurvefit(@fourbar,L0,theta1,theta2);

function Phi = fourbar(y,L,theta1)

    u1 = [cos(theta1), sin(theta1)];
    u2 = [cos(y(1)), sin(y(1))];
    u3 = [cos(y(2)), sin(y(2))];

    Phi = L(1)*u1 + L(2)*u2 - L(3)*u3 - [L(4), 0];
end

But unfortunately it doesn't work and I don't know how to solve the problem.

[Matt J]

It doesn't look like you have enough constraints to determine four L(i). If I multiply all L(i) by a common positive scalar, the angles don't change.

[Gregory Cottone]

Thank you for the answer.

If I set L0 to a determinate value, how can I optimize the others?

Thank you.

[Matt J]

I'm still not sure there are enough constraints. For example L1=L2=0, L3=L0 is always a solution.

Also, do you not have theta3 data as well? That would reduce the problem to a set of linear equations.

[Gregory Cottone]

Unfortunately I don't have any theta3 data to fit in my problem.

Answer 1

[Matt J]

I would just use lsqnonlin.

theta2 = [0.53,0.40,0.32,0.28,0.27,0.27,0.28,0.29,0.31,0.34,0.36,0.40,0.43,0.47,0.50,0.55,0.59,0.64,0.69,0.75]; %(data to fit)
theta1 = linspace(0,pi,numel(theta2)); 
L0=1;

opts=optimoptions(@lsqnonlin,'StepTolerance',1e-12,'OptimalityTolerance',1e-12,'FunctionTolerance',1e-12);
[L,resn,res]=lsqnonlin(@(L) resid(L,L0,theta1,theta2),  [0.5,0.5,1]*L0, [0,0,0],[1,inf inf]*L0);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

L

L = 1×3

0.5047 1.4891 1.0920

function fval=resid(L,L0,theta1,theta2)

  dX=L(1)*cos(theta1(:))+L(2)*cos(theta2(:));
  dY=L(1)*sin(theta1(:))+L(2)*sin(theta2(:));
  
  fval=(L0-dX).^2 + dY.^2-L(3).^2;
  
end

[Gregory Cottone]

Really thank you!

It seems work!

At this point my only question is what you have done in the function "resid"?

How you started with:

function Phi = fourbar(y,L,theta1)
    u1 = [cos(theta1), sin(theta1)];
    u2 = [cos(y(1)), sin(y(1))];
    u3 = [cos(y(2)), sin(y(2))];
    Phi = L(1)*u1 + L(2)*u2 - L(3)*u3 - [L(4), 0];
end

and come up with:

function fval=resid(L,L0,theta1,theta2)
  dX=L(1)*cos(theta1(:))+L(2)*cos(theta2(:));
  dY=L(1)*sin(theta1(:))+L(2)*sin(theta2(:));
  fval=(L0-dX).^2 + dY.^2-L(3).^2;
end

Thank you!

[Matt J]

I did not start with fourbar (I don't understand it). I started with your diagram. Basically, the components of u3 are,

u3=[dX-L0,dY]

and the equations we want to solve are

norm(L3)^2=norm(u3)^2

[Gregory Cottone]

Now I understand! Thank you Matt.