Curve-fitting by Levenberg–Marquardt method??

I am new to Matlab, but I want to write a matlab code to estimate the values of spring and damper coefficients to fit a function to data curve. I have a simultaneous ode equation for acceleration of a planar rigid body(mass=m) connected to a pair of spring and damper in parallel in x-axis and another pair in y-axis directions as follows.

M*D2X(t)=-Cx*DX(t)-Kx*X(t)+w*DY(t)

M*D2Y(t)=-Cy*DY(t)-Ky*Y(t)+w*DX(t)

where: M=mass C=damper coefficient K=spring coefficient w =yaw rate (constant value) In addition to the functions above, I have x-axis and y-axis direction acceleration data from experimental measurement. So, how can I estimate four parameters Kx ,Cx , Ky and Cy?? I appreciate any hint or help!!

Answer 1

[Matt J]

So X,Y are positions or accelerations? If they are positions and you have acceleration data, D2X and D2Y you can integration that data numerically to obtain X,Y, DX,DY. Then, your equations become linear and you can use total lienar least squares to solve them.

[Dave Yadeta]

X and Y are positions. I have D2X and D2Y acceleration data, but i needed my model eqn to fit to the acceleration data curve. So I need to estimate the coefficient valuses C_x,K_x,C_y and K_y. I tried the following in matlab.

Model eqn;

function d2x = mod1(t,x,dx) % Model implementation

global alpha beta delta gamma;

d2x = [0;0]; dx=[9.72;0]; d2x(1) = -alpha * dx(1)-beta*x(1)+ 2*dx(2); d2x(2) = -delta* dx(2) -gamma*x(2)-2*dx(1);

Expected acceleration func:

function y = accex(a,x) % To return expected value of acceleration in x-axis %direction given the a set of time points and coefficients. % a-list of parameter values.

global alpha beta delta gamma;

alpha = a(1); beta = a(2); delta =a(3); gamma =a(4);

[t,y] = ode15s('mod1',[0 x(end)],[0,0]); y = interp1(t,y(:,1),x); options = statset('Display','iter');

data upload:

load Filtered.mat % acceleration data

time1 = Time(3001:4800,1); % Time series

acce1 = Xcfc180(3001:4800,1);% acceleration data in X-axis direction

acce2 =Ycfc180(3001:4800,1);% acceleration data in Y-axis direction

output coefficent values:

options = statset('Display', 'iter'); format long [values] = nlinfit(time1, acce1, 'accex', [1250 166000 3492 10000], options)

what is wrong here?

[Matt J]

I'm not sure that anything is wrong, and you haven't mentioned any symptoms. However, I wonder if you need anthing that complicated. Your differential equations are linear, after all, yet you seem to be attempting a very complicated non-linear method of solution. If you have measurements of D2X and D2Y, you can compute DX,DY,X,and Y as follows

   DX=trapz(time,DX2);
   DY=trapz(time,DY2);

   X=trapz(time,DX);
   Y=trapz(time,DY);

Now you construct linear equations in your unknown parameters

 Ax=[DX(:), X(:)];   %equation Ax*CxKx=bx
 bx=[w*DY(:)-M*D2X(:)];

 Ay=[DY(:), Y(:)];   %equation Ay*CyKy=by
 by=[w*DX(:)-M*D2Y(:)];

If your acceleration data is sufficiently non-noisy it may be enough to do a straight least squares linear solve

 CxKx = Ax\bx;
   Cx=CxKx(1); Kx=CxKx(2);

 CyKy = Ay\by;
      Cy=CyKy(1); Ky=CyKy(2);

A more robust solution would probably be to use Total Least Squares

[U,S,V] = svd([Ax,bx],0);
 CxKx = -V(1:2,3)/V(3,3);

and similarly for Ky and Cy.

[Dave Yadeta]

I have not noticed such a short path before I saw your answer. I get nice estimate now. Thank you very much!