help code with lsqnonlin

3 views (last 30 days)
franklin momo
franklin momo on 1 Jun 2022
Commented: Bjorn Gustavsson on 19 Jun 2022
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
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.

Sign in to comment.

Answers (2)

Bjorn Gustavsson
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
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
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);

Sign in to comment.


franklin momo
franklin momo on 19 Jun 2022
greeting
thank for your reply. i know that there are input values thats will not give an answers, i just want to skip those input values and go to the next. i am try to input some thing like 1000 values and put them one by one will be easy so i write a code and want that the code skip all input values and continu running to the next one. so i will have many input values and many answers. we can [0 0 0 0 0 0 0 0] as answers for the input values that does not have answers.
below the code that i have written
i am trying to check the code that you wrote if it can help me.
% 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];
Rmin = [1 10 20 30 40 80 160 320 640 1250 2500 5000];
Rmax = [10 20 30 40 80 160 320 640 1250 2500 5000 10000];
Nreel = [100 200 30 20 40 60 200 180 60 20 5 1];
N = 916;
a1=0.3; a11=0.34;
b1=8.5; b11=11.5;
c1=1.1; c11=1.2;
a2=0.59; a22=0.69;
b2=119; b22=121;
c2=1.1; c22=1.9;
b3=2258; b33=2261;
c3=1; c33=1.2;
a=(a1+a11)/2; b=(b1+b11)/2; c=(c1+c11)/2; d=(a2+a22)/2; e=(b2+b22)/2;
f=(c2+c22)/2; g=(b3+b33)/2;h=(c3+c33)/2;
x1(1)=a1; r1(1)=b1; t1(1)=c1; x2(1)=a2; r2(1)=b2; t2(1)=c2; r3(1)=b3; t3(1)=c3;
for i=1:2
x1(i+1)=x1(i)+a;
r1(i+1)=r1(i)+b;
t1(i+1)=t1(i)+c;
x2(i+1)=x2(i)+d;
r2(i+1)=r2(i)+e;
t2(i+1)=t2(i)+f;
r3(i+1)=r3(i)+g;
t3(i+1)=t3(i)+h;
end
for h=1:3
for i=1:3
for j=1:3
for k=1:3
for l=1:3
for m=1:3
for n=1:3
for o=1:3
x0 = [x1(o) r1(n) t1(m) x2(l) r2(k) t2(j) r3(i) t3(h)]
lb = [0.2 8 1 0.55 117 1 2257 1];
ub = [0.4 12 1.5 0.7 122 2 2262 1.5];
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
[x,resnorm,residual,exitflag,output] = lsqnonlin(@(x)paramter(x,N,Rmin,Rmax,Nreel),x0,[],[],options)
% x3(i) = 1-x1(o)-x2(i)
end
end
end
end
end
end
end
end
function [Nest1] = paramter(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))./log10(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
and the coomiplation give
>> ext4
x0 =
1.0e+03 *
0.0003 0.0085 0.0011 0.0006 0.1190 0.0011 2.2580 0.0010
Local minimum possible.
lsqnonlin stopped because the final change in the sum of squares relative to
its initial value is less than the default value of the function tolerance.
<stopping criteria details>
x =
1.0e+03 *
0.0002 0.0109 0.0010 0.0007 0.3017 0.0023 2.2580 0.0010
resnorm =
1.6480e+04
residual =
Columns 1 through 11
-99.9898 21.9651 -28.7112 -17.0779 -9.3213 52.7751 9.8833 18.1139 33.1513 21.9651 -1.9510
Column 12
-0.8086
exitflag =
3
output =
struct with fields:
iterations: 20
funcCount: 190
stepsize: 0.0097
cgiterations: []
firstorderopt: 3.5155
algorithm: 'levenberg-marquardt'
message: 'Local minimum possible.↵↵lsqnonlin stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.↵↵Stopping criteria details:↵↵Optimization stopped because the relative sum of squares (r) is changingby less than options.FunctionTolerance = 1.000000e-06.↵↵Optimization Metric Optionsrelative change r = 3.73e-07 FunctionTolerance = 1e-06 (default)'
x0 =
1.0e+03 *
0.0006 0.0085 0.0011 0.0006 0.1190 0.0011 2.2580 0.0010
Error using erf
Input must be real and full.
Error in ext4>paramter (line 64)
Nest1 =
((N.*x(1)./2).*(erf((sqrt(2)/2).*(log(Rmax./x(2))/log(x(3))))-erf((sqrt(2)./2).*(log(Rmin./x(2))./log10(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);
Error in ext4>@(x)paramter(x,N,Rmin,Rmax,Nreel)
Error in levenbergMarquardt (line 212)
trialCostFun = feval(funfcn{3},reshape(trialX,xShape),varargin{:});
Error in lsqncommon (line 174)
levenbergMarquardt(funfcn,xC,flags.verbosity,options,defaultopt,initVals.F,initVals.J, ...
Error in lsqnonlin (line 240)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,caller,...
Error in ext4 (line 50)
[x,resnorm,residual,exitflag,output] =
lsqnonlin(@(x)paramter(x,N,Rmin,Rmax,Nreel),x0,[],[],options)
>>
we can see the code start running and stop when the input value do not give an anwers. i wwant to skip those input values by give them the answers [0 0 0 0 0 0 0 0] and continu to the next value

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!