Parameters identification with Levenberg-Marquardt algorithm

close all
clear all ;

A = dlmread('10kg.csv', ',', 0, 0 );
 
eps=A(:,2);
esp_dot=A(:,5);
esp_2dot=A(:,7);
P_log=A(:,3)';
 
L_0=0.343;  
g=9.81; 
dt=0.02;
Fdis=10*g;
m=Fdis/g;
 
x=eps'*L_0;
x_dot=esp_dot'*L_0;
x_2dot=esp_2dot'*L_0;

N=size(P_log');
M=40 ;
count=0;
 
W=1;

c_0=1;   
k_0=50;  
n_0=10;    
a1_0=-50;  
a0_0=50;


c=c_0;
k=k_0;
n=n_0;
a1=a1_0;
a0=a0_0;

b=0.001;
E_0=0;

for j=1:N
  r=P_log(j)-(m*x_2dot(j)+c*x_dot(j)+k*x(j).^(1/n)+Fdis)/(a1*x(j)+a0);
 
  if j==N
    E_0 = E_0 + W*r*r;
  else
    E_0 = E_0 + r*r;
  end
 
end

E_0;
 
pause(1) ;
 
E=E_0*2;

%while(E_0<E)
 
for i=1:M
 
  count = count + 1;
  all_log(:,count)=[c;k;n;a1;a0];
  
  if E_0<=0.0005
    break;
  end
 
  Dfcc=0;
  Dfck=0;
  Dfcn=0;
  Dfca1=0;
  Dfca0=0;
 
  Dfkc=0;
  Dfkk=0;
  Dfkn=0;
  Dfka1=0;
  Dfka0=0;
 
  Dfnc=0;
  Dfnk=0;
  Dfnn=0;
  Dfna1=0;
  Dfna0=0;
 
  Dfa1c=0;
  Dfa1k=0;
  Dfa1n=0;
  Dfa1a1=0;
  Dfa1a0=0;
 
  Dfa0c=0;
  Dfa0k=0;
  Dfa0n=0;
  Dfa0a1=0;
  Dfa0a0=0;
 
  DfcR=0;
  DfkR=0;
  DfnR=0;
  Dfa1R=0;
  Dfa0R=0;
 
  %E_0=0;
 
  for j=1:N
    
      
    if x(j)==0
      x(j)=1*10^(-300);
    end
        
    dfc=-x_dot(j)/(a1*x(j)+a0);
    dfk=-x(j)^(1/n)/(a1*x(j)+a0);
    dfn=(m*x_2dot(j)+c*x_dot(j)+k*x(j)^(1/n)+Fdis)/(a1*x(j)+a0)*1/n^2*log(x(j));
    dfa1=(m*x_2dot(j)+c*x_dot(j)+k*x(j)^(1/n)+Fdis)*x(j)/(a1*x(j)+a0)^2;
    dfa0=(m*x_2dot(j)+c*x_dot(j)+k*x(j)^(1/n)+Fdis)/(a1*x(j)+a0)^2;
  
    r=P_log(j)-(m*x_2dot(j)+c*x_dot(j)+k*x(j)^(1/n)+Fdis)/(a1*x(j)+a0);
    r_log(j)=r;
        
    E_0 = E_0 +r*r;
        
    Dfcc = Dfcc + dfc*dfc;
    Dfck = Dfck + dfc*dfk;
    Dfcn = Dfcn + dfc*dfn;
    Dfca1 = Dfca1 + dfc*dfa1;
    Dfca0 = Dfca0 + dfc*dfa0;
        
    Dfkc = Dfkc + dfk*dfc;
    Dfkk = Dfkk + dfk*dfk;
    Dfkn = Dfkn + dfk*dfn;
    Dfka1 = Dfka1 + dfk*dfa1;
    Dfka0 = Dfka0 + dfk*dfa0;
        
    Dfnc = Dfnc + dfn*dfc;
    Dfnk = Dfnk + dfn*dfk;
    Dfnn = Dfnn + dfn*dfn;
    Dfna1 = Dfna1 + dfn*dfa1;
    Dfna0 = Dfna0 + dfn*dfa0;
        
    Dfa1c = Dfa1c + dfa1*dfc;
    Dfa1k = Dfa1k + dfa1*dfk;
    Dfa1n = Dfa1n + dfa1*dfn;
    Dfa1a1 = Dfa1a1 + dfa1*dfa1;
    Dfa1a0 = Dfa1a0 + dfa1*dfa0;
        
    Dfa0c = Dfa0c + dfa0*dfc;
    Dfa0k = Dfa0k + dfa0*dfk;
    Dfa0n = Dfa0n + dfa0*dfn;
    Dfa0a1 = Dfa0a1 + dfa0*dfa1;
    Dfa0a0 = Dfa0a0 + dfa0*dfa0;
        
    DfcR = DfcR + dfc*r;
    DfkR = DfkR + dfk*r;
    DfnR = DfnR + dfn*r;
    Dfa1R = Dfa1R + dfa1*r;
    Dfa0R = Dfa0R + dfa0*r;
        
  end
        
  H=[Dfcc*(1+b) Dfck Dfcn Dfca1 Dfca0;...
    Dfkc Dfkk*(1+b) Dfkn Dfka1 Dfka0;...
    Dfnc Dfnk Dfnn*(1+b) Dfna1 Dfna0;...
    Dfa1c Dfa1k Dfa1n Dfa1a1*(1+b) Dfa1a0;...
    Dfa0c Dfa0k Dfa0n Dfa0a1 Dfa0a0*(1+b)];
  
  R=[-DfcR;-DfkR;-DfnR;-Dfa1R;-Dfa0R];
        
  delta=pinv(H)*R;
        
  ct=c+delta(1);
  kt=k+delta(2);
  nt=n+delta(3);
  a1t=a1+delta(4);
  a0t=a0+delta(5);
        
  E=0;
  for l=1:N
    r=P_log(l)-(m*x_2dot(l)+ct*x_dot(l)+kt*x(l)^(1/n)+Fdis)/(a1t*x(l)+a0t);
    
    if l==N
      E = E + W*r*r;
    else
      E = E + r*r;
    end
        
  end
  
  fprintf('Current fitting error after GN = %f\n',E) ;
  pause(1) ;
  
  ctt=ct;
  ktt=kt;
  ntt=nt;
  a1tt=a1t;
  a0tt=a0t;
        
  while(E_0<=E)
        
    b=10*b;
        
    H=[Dfcc*(1+b) Dfck Dfca1 Dfca0;...
      Dfkc Dfkk*(1+b) Dfka1 Dfka0;...
      Dfnc Dfnk Dfnn*(1+b) Dfna1 Dfna0;...
      Dfa1c Dfa1k Dfa1a1*(1+b) Dfa1a0;...
      Dfa0c Dfa0k Dfa0a1 Dfa0a0*(1+b)];
    
    R=[-DfcR;-DfkR;-DfnR;-Dfa1R;-Dfa0R];
        
    delta=0.2*inv(H)*R;
            
    ct=ctt+delta(1);
    kt=ktt+delta(2);
    nt=n+delta(3);
    a1t=a1tt+delta(4);
    a0t=a0tt+delta(5);
            
    count = count + 1;
    all_log(:,count)=[ct;kt;nt;a1t;a0t];
    count;
            
    E=0;
    for l=1:N
      r=P_log(l)-(m*x_2dot(l)+ct*x_dot(l)+kt*x(l)^(1/nt)+Fdis)/(a1t*x(l)+a0t);
            
      if l==N
        E = E + W*r*r;
      else
        E = E + r*r;
      end
      
    end
            
    fprintf('Current fitting error after Steepest descent = %f\n',E) ;
    
    if E_0>E || abs(delta(5))<0.00001
      break;
    end
            
  end   
        
  b=0.1*b;
  
  c=ct;
  k=kt;
  n=nt;
  a1=a1t;
  a0=a0t;
    
end 
 

 


 

c
k
n
a1
a0
 
E=0;
 
for j=1:N
  r=P_log(j)-(m*x_2dot(j)+c*x_dot(j)+k*x(j)^(1/n)+Fdis)/(a1*x(j)+a0);
  E = E + r*r;
end
 
E

i=0:count-1;
figure(1);
a=axes();
set(a,'fontsize',20) ;
 
subplot(2,3,1);
plot(i,all_log(1,:),'o-');
ylabel('c','fontsize',17);
xlabel('count');
grid on ;
 
subplot(2,3,2);
plot(i,all_log(2,:),'o-');
ylabel('k','fontsize',17);
xlabel('count');
grid on ;
 
subplot(2,3,3);
plot(i,all_log(3,:),'o-');
ylabel('a1','fontsize',17);
xlabel('count');
grid on ;
 
subplot(2,3,4);
plot(i,all_log(4,:),'o-');
ylabel('a0','fontsize',17);
xlabel('count');
grid on ;
 
subplot(2,3,5);
plot(i,all_log(5,:),'o-');
ylabel('n','fontsize',17);
xlabel('count');
grid on ;
 
%print -dpng '3.png' ;

Hello everyone, I can't identify the correct parameters, please help me to see where is wrong?

[Star Strider]

It appears that you want to do nonlinear parameter estimation. MATLAB has several funcitons for that, and have already resolved these issues. If you have a specific problem you would like help with (for example, estimating the parameters of a specific function), please describe the function, the parameters you want to estimate, and any constraints.

If you instead want help in understanding the Levenberg-Marquardt algorithm, please be specirfic as to what you want help with in that respect.

[XIAO LIU]

Thank you very much for your reply.

I have a specific problem would like help with (The specific conditions are shown in the figure above),

The function:(Fitting error)r=P-(m*x_2dot+c*x_dot+k*x^(1/n)+mg)/(a1*x+a0) Fitting criterion J as the picture shows.

The parameters I want to estimate:c k n a1 a0

Constraints:c＞0 k＞0 n＞0 a0?0(Maybe＞0) a1?0(Maybe＜0)

Thank you for your precious time.