imaginary value from solution nonlinier least square

im trying to fitting data from reflectance measurement with lavenberg marquardt algorithm , ive got this code LMFnlsq in file exchange but i dont understand how to use it in my case

here is the equation

R = (0.6/(a + b))*(sqrt(3*a*(a + b))+(1/c))*(exp(-sqrt(3*a*(a + b ))*c))/(c^2)

a and b is parameter that im trying to solve

R is known vector from measurement reflectance

c is known vector from distance in measurement

im using the code like this

r = [0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.6]';

data = [0.787357544 0.63199586 0.393109753 0.271890049 0.179074593 0.133498607 0.112568249 0.05274273]';

res = @(x) ((0.6./(x(1)+x(2))).*(sqrt(3.*x(1).*(x(1)+x(2)))+(1./r)).* (exp(-sqrt(3*x(1).*(x(1)+x(2)))*r))./(r.^2)) - data;

x0 = [0.1,10];

[x,ssq,cnt] = LMFnlsq(res,x0)

hold on

plot(r,data)

plot(r,res(x)+ data,'y'), grid

the problem is it sure did the iteration but the final value is in imaginary solution like this

x =

   0.0528 - 0.2730i
   1.1469 + 0.6248i

ssq =

    0.0091

cnt =

    29

anything wrong with the code?

Answer 1

[Miroslav Balda]

Let us reformulate the problem:

data = (q + 1./r).*exp(-r*q).*(0.6/(p(1)*r.^2)),

where

q = sqrt(p(2));

The problem may be solved without any troubles by introducing a penalty function

(q<=0)*w

into vector of residuals. Thus, the main script could be

% Surpan.m      Application of LMFnlsq
%%%%%%%%%%%    12. 3. 2013
% Answers: Imaginary value from solution nonlinear least square
clc
clear all
close all
%
global data r p0 w
%
data = [0.779543143 0.624771688 0.387776261 0.266986051 0.174732672 ...
      0.130313689 0.110540421 0.051412959]';  %   column vector
r    = [0.8 0.9 1 1.1 1.2 1.3 1.4 1.6]';        %   column vector
w    = 10;
%
while 1
% function inp from   www.mathworks.com/matlabcentral/fileexchange/9033
  w = inp('weight', w);       %   weight of penalty function
  if w<=0, break, end
  p0 = rand(2,1);
  [p,ssq,cnt] = LMFnlsq(@resid,ones(2,1),'Display',-100,'MaxIter',1e3);
  p = p.*p0;
  a = p(2)/(3*p(1))
  b = p(1)-a
  res = resid(p./p0);
  plot(r,data,'o-',  r,res(1:end-1)+data,'-or')
end

The function for evaluating residuals has the form:

function res = resid(p)   %   Residuals of Surpans' equations
%%%%%%%%%%%%%%%%%%%%%%%%%
global data r w p0
%    
p = p.*p0;              %   p(1) = (a+b)
q = sqrt(p(2));         %   q = sqrt(3*a*(a+b))
%                        
res = [(q + 1./r).*exp(-r*q).*(0.6./(p(1)*r.^2)) - data
       (q>=0)*w         %   penalty function
      ];

The nonlinear fit is made by the function LMFnlsq, which may be found on

www.mathworks.com/matlabcentral/fileexchange/17534

Answer 2

[Matt J]

You have square roots in your function, which will produce complex values when their argument goes negative. You must use another algorithm, one that lets you constrain 3*a*(a + b))+(1/c) to a strictly positive lower bound.

[Surpan]

is that wrong if im using the code in res above with real?

res = @x real(....)

but the parameter a and b do have constraint a in 0.09-0.8 b in 7-15 but i dont know how to use it in the code

[Matt J]

Adding real(...) isn't a great solution because then your function would be non-differentiable where 3.*x(1).*(x(1)+x(2)))+(1./r) = 0.

If you have the Optimization Toolbox, there are lots of alternative solvers that would let you incorporate constraints.

[Surpan]

yeah i think so

but i must solve it with lavenberg marquardt algorithm so i search in FEX and got that LMFnlsq from prof misrolav

[Matt J]

but i must solve it with lavenberg marquardt algorithm

Meaning it's an assignment? Then get clarification from the one who assigned it to you. LM is not applicable to constrained problems.