help code with lsqnonlin
4 views (last 30 days)
Show older comments
greenting
please i want to know where i have made an error of this code. the code try to estimate 8 unkwons values of x for a function with lsqnonlin.
this is the code
1-
2-
3- format long
4- Rmin = [0.001 0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0];
5- Rmax = [0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0 10.0];
6- Nreel = [100 200 30 20 40 60 200 180 60 20 5 1];
7- N = 916;
8-
9- Nest = ((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3))))-erf((sqrt(2)./2).*(log(Rmin./x(2)). /log(x(3)))))+(N.*x(4)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(5))./log(x(6))))-erf((sqrt(2)/2).*(log(Rmin./x(5))./log(x(6)))))+(N.*(1-x(1)-x(4))./2).*(erf((sqrt(2)/2).*(log(Rmax./x(7))./log(x(8))))-erf((sqrt(2)/2).*(log(Rmin./x(7))./log(x(8)))))-Nreel);
10-
11- % x0 = [3.,3.,3.,3.,3.,3.,3.,3.];
12- x0 = [0.11;1;1.5;0.86;1;2;1;2];
13- lb = zeros(8,1);
14- lb(:) = -inf;
15- ub = zeros(8,1);
16- ub(:) = inf;
17-
19- % options = optimoptions('lsqnonlin','Display','iter');
20- % options.Algorithm = 'levenberg-marquardt';
21- % 'trust-region-reflective'
22- options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
24- [x,resnorm,residual,exitflag,output] = lsqnonlin(@(x)Nest,x0,[],[],options)
25-
but when is run the code that's the answer
extimation
Undefined function or variable 'x'.
Error in extimation (line 9)
Nest =
((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3))))-erf((sqrt(2)./2).*(log(Rmin./x(2))./log(x(3)))))+(N.*x(4)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(5))./log(x(6))))-erf((sqrt(2)/2).*(log(Rmin./x(5))./log(x(6)))))+(N.*(1-x(1)-x(4))./2).*(erf((sqrt(2)/2).*(log(Rmax./x(7))./log(x(8))))-erf((sqrt(2)/2).*(log(Rmin./x(7))./log(x(8)))))-Nreel);
>>
1 Comment
Matt J
on 1 Jun 2022
I recommend not including line numbers when you post your code. It makes it harder for responders to copy/paste it.
Answers (2)
Bjorn Gustavsson
on 1 Jun 2022
Edited: Bjorn Gustavsson
on 1 Jun 2022
At the very least you'll have to make Nest into a function - otherwise you have nothing to optimize. Perhaps something like this:
Rmin = [0.001 0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0];
Rmax = [0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0 10.0];
Nreel = [100 200 30 20 40 60 200 180 60 20 5 1];
N = 916;
% Not that this is the function-definition: (using @(x) creates a function-handle to a
% function that depends on x and takes all other variables in the
% definition as they are at the time of function-handle-creation
Nest = @(x) ((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3))))-erf((sqrt(2)./2).*(log(Rmin./x(2)). /log(x(3)))))+(N.*x(4)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(5))./log(x(6))))-erf((sqrt(2)/2).*(log(Rmin./x(5))./log(x(6)))))+(N.*(1-x(1)-x(4))./2).*(erf((sqrt(2)/2).*(log(Rmax./x(7))./log(x(8))))-erf((sqrt(2)/2).*(log(Rmin./x(7))./log(x(8)))))-Nreel);
% x0 = [3.,3.,3.,3.,3.,3.,3.,3.];
x0 = [0.11;1;1.5;0.86;1;2;1;2];
lb = zeros(8,1);
lb(:) = -inf;
ub = zeros(8,1);
ub(:) = inf;
% options = optimoptions('lsqnonlin','Display','iter');
% options.Algorithm = 'levenberg-marquardt';
% 'trust-region-reflective'
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
[x,resnorm,residual,exitflag,output] = lsqnonlin(@(x) Nest(x),x0,[],[],options);
% Above we generate a dynamic function that depends on x. This to make the
% slightly more general procedure easier to get to:
Nest2 = @(x,Rmax,Rmin,Nreel) ((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3))))-erf((sqrt(2)./2).*(log(Rmin./x(2)). /log(x(3)))))+(N.*x(4)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(5))./log(x(6))))-erf((sqrt(2)/2).*(log(Rmin./x(5))./log(x(6)))))+(N.*(1-x(1)-x(4))./2).*(erf((sqrt(2)/2).*(log(Rmax./x(7))./log(x(8))))-erf((sqrt(2)/2).*(log(Rmin./x(7))./log(x(8)))))-Nreel);
% Nest2 is now a function that takes 4 arguments. This makes it possible to
% re-use this function for multiple different Rmax, Rmin or Nreel
% parameters and fit for the correspongin optimal x:
[x,resnorm,residual,exitflag,output] = lsqnonlin(@(x) Nest2(x,Rmax,Rmin,Nreel),x0,[],[],options);
HTH
8 Comments
Bjorn Gustavsson
on 17 Jun 2022
Edited: Bjorn Gustavsson
on 19 Jun 2022
That sounds like "a very dubious idea", because that would be to dynamically change/modify the function you're trying to optimize during the optimization. My suggestion is that you instead constrain the search to the region of parameters that gives sensible outputs (the way I understand your problem and the error is that you have some physics/principled based reasons that all the calls to log should give real results). To get towards that I've rewritten your function as:
function [Nest1] = paramters(x,N,Rmin,Rmax,Nreel)
Nest1 = ((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3)))) - ...
erf((sqrt(2)./2).*(log(Rmin./x(2))./log(x(3))))) + ...
(N.*x(4)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(5))./log(x(6)))) - ...
erf((sqrt(2)/2).*(log(Rmin./x(5))./log(x(6))))) + ...
(N.*(1-x(1)-x(4))./2).*( ...
erf((sqrt(2)/2).*(log(Rmax./x(7))./log(x(8)))) - ...
erf((sqrt(2)/2).*(log(Rmin./x(7))./log(x(8)))))-Nreel);
end
This gives you one erf-call per line, and makes it easier to isolate the possible real-log-violators. If I get it those should be: log(Rmax./x(2)), log(x(3)), log(Rmin./x(2)), log(Rmax./x(5)), log(x(6)), log(Rmin./x(5)), log(Rmax./x(7)), log(x(8)), and log(Rmin./x(7)). Since all of Rmin and Rmax are positive this means that all of x([2 3 5 6 7 8]) has to be positive. Therefore you could (perhaps: should) enforce lower bound on those components of x. Perhaps something like:
Rmin = [0.001 0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0];
Rmax = [0.01 0.02 0.03 0.04 0.08 0.16 0.32 0.64 1.25 2.5 5.0 10.0];
Nreel = [100 200 30 20 40 60 200 180 60 20 5 1];
N = 916;
x0 = [0.11;1;1.5;0.86;1;2;1;2];
lb(:) = -inf(8,1);
ub(:) = inf(8,1);
lb([2 3 5 6 7 8]) = 0;
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
[x,resnorm,residual,exitflag,output] = lsqnonlin(@(x)paramter(x,N,Rmin,Rmax,Nreel),...
x0,lb,ub,...
options);
Since the first and fourth elements of lb are -inf they will remain unconstrained, while the values you allow for the other components are constrained to be positive.
HTH
Bjorn Gustavsson
on 19 Jun 2022
The lower bounds should be just larger than zero to be sure that we never get there (from the help:
X = lsqnonlin(FUN,X0,LB,UB) defines a set of lower and upper bounds on
the design variables, X, so that the solution is in the range LB <= X
<= UB.)
So change:
lb([2 3 5 6 7 8]) = 0;
to something like:
lb([2 3 5 6 7 8]) = eps(1);
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!