The issue of optimization (minimization) of the average relative error between experimental and calculated data

hello

I want to share the difficulties that I faced. Can someone help

problem statement:

there is a 'x' column where the values ​​of the independent variable are written and there is a 'y' column where the experimental values ​​of the dependent variable are written.

approximation model is considered:

y_calculate=A*x^B+C,

and based on this model, an objective function is created, which is equal to the average value of the relative deviations between y and y_calculate:

error_function = mean(abs(y - y_calculate)) / y)=mean(abs(y - =A*x^B+C)) / y);

Our goal is to select parameters A,B,C in such a way that 'error_function' takes the value of the global minimum.

I calculated the optimal values ​​of A, B, C and got:

A = 85.5880, B = -0.0460, C = 4.8824,

at which error function value for optimized parameters: 0.0285.

but I know in advance the specific values ​​of A, B, C:

A = 1005.6335852931, B = -1.59745963925582, C = 73.54149744754400,

at which error function value for specific parameters: 0.002680472178434,

which is much better than with optimization

Below is the code with visualization, which confirms the above.

clear
close all
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];

% Initial guesses for parameters A, B, C
initial_guess = [1, 1, 1];

% Error function
error_function = @(params) mean(abs(y - (params(1) * x.^params(2) + params(3))) ./ y);

% Optimization of parameters
optimized_params = fminsearch(error_function, initial_guess);

% Results of optimization
A_optimized = optimized_params(1);
B_optimized = optimized_params(2);
C_optimized = optimized_params(3);

% Calculation of the fitted function for optimized parameters
y_calculate_optimized = A_optimized * x.^B_optimized + C_optimized;

% Calculate and display the error function value for optimized parameters
value_error_optimized = error_function(optimized_params);
fprintf('Optimized parameters:\nA = %.4f\nB = %.4f\nC = %.4f\n', A_optimized, B_optimized, C_optimized);
fprintf(' error function value for optimized parameters: %.4f\n', value_error_optimized);

% Other specific parameters A, B, C
A_specific = 1005.63358529310;
B_specific = -1.59745963925582;
C_specific = 73.541497447544;

% Calculation of the fitted function for specific parameters
y_calculate_specific = A_specific * x.^B_specific + C_specific;

% Calculate and display the error function value for specific parameters
value_error_specific = error_function([A_specific, B_specific, C_specific]);
fprintf('Specific parameters:\nA = %.10f\nB = %.14f\nC = %.14f\n', A_specific, B_specific, C_specific);
fprintf(' error function value for specific parameters: %.4f\n', value_error_specific);

% Visualization
figure;
plot(x, y, 'bo-', 'DisplayName', 'Experimental data');
hold on;
plot(x, y_calculate_optimized, 'r--', 'DisplayName', 'Fitted model (Optimized)');
plot(x, y_calculate_specific, 'g-.', 'DisplayName', 'Fitted model (Specific)');
xlabel('x');
ylabel('y');
legend('Location', 'best');
title('Approximation of experimental data');
grid on;

Obviously, my optimization code does not lead to a global minimum of the objective function, since there is a better approximation for specific values ​​of A,B,C. Maybe this is caused by a random selection of the initial values ​​of the parameters A=1, B=1, c=1 and therefore my code is stuck in a local minimum?

who can write a code that will select the A,B,C parameters so as to achieve the global minimum of the target function 'error_function', for any initial iteration data of the variables A,B,C. Thoughts for testing: the value of the target function 'error_function' should not be worse (that is, more) than 0.002680472178434, which is obtained with the specific value of A,B,C: A = 1005.6335852931, B = -1.59745963925582, C = 73.54149744754400

Answer 1

[Matt J]

If you download minL1lin from

https://www.mathworks.com/matlabcentral/fileexchange/52795-constrained-minimum-l1-norm-solutions-of-linear-equations

then you can do,

% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];

% Initial guess for parameter B
initial_guess = 1;

% Error function
error_function = @(B) objFun(B, x,y);

% Optimization of parameters
B = fminsearch(error_function, initial_guess);

[fval,ABC]=error_function(B);

[A,B,C]=deal(ABC{:})

A = 1.0460e+03

B = -1.6184

C = 73.5535

fval

fval = 0.0025

function [fval,p]=objFun(theta, x,y)
arguments
  theta; x (:,1); y(:,1);
end

 B=theta;

 d=ones(size(y));
 Q=[x.^B, d]./y; 
 
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 p={A, B, C};

end

[roborrr]

Matt J, thank you for your response.

I ran your code in my matlab and saw the following error:

"Error: Line: 26 Column: 2

Function definitions in a script must appear at the end of the file.

Move all statements after the "objFun" function definition to before the first local function definition."

As I found out the reason of this error is that my version of matlab R(2018b) does not support the following construction:

arguments
  B; x (:,1); y(:,1);
end

this construction was introduced starting with version R(2019b)

please rework this part of the code so that it works in my version of matlab - R2018b.

[Matt J]

Just do,

function [fval,p]=objFun(theta, x,y)

 x=x(:); y=y(:); B=theta;

 d=ones(size(y));
 Q=[x.^B, d]./y; 
 
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 p={A, B, C};
end

[roborrr]

Dear Matt J, I checked the code you suggested and got a stunning result for the model

y=A*x^B+C

However, I did not understand much about how it works.

My goal is to check other models on the same experimental data x,y and choose the best model from them. I tried to do it myself, following your example, but I did not succeed and gave up.

Here are the models that I want to check:

Model 1

y_calculate = eta_inf + (eta_0 - eta_inf) / (1 + (lambda * x)^m);

In this model, two options need to be considered:

1.1- all parameters need to be optimized - eta_0 , eta_inf , lambda , m

1.2 - we assume that m=1.5 is known and we optimize only the parameters - eta_0, eta_inf, lambda

Model 2

y_calculate = eta_inf + (eta_0 - eta_inf) / (1 + (lambda * x)^a)^((1 - n) / a);

in this model, as in model 1, two cases are considered:

2.1. - all parameters need to be optimized - eta_0 , eta_inf , lambda , n,a

2.2 - we assume that parameters eta_0 , eta_inf , lambda , n is known and we optimize only the parameter a.

[Matt J]

Dear Matt J, I checked the code you suggested and got a stunning result for the model

Can I assume "stunning" means that it did what you want? If so, please Accept-click the answer to indicate that it resolved your question.

I tried to do it myself, following your example, but I did not succeed and gave up.

To see what went wrong, I would need to see what you tried.

all parameters need to be optimized - eta_inf , eta_inf , lambda , m

Here and elsewhere you have listed eta_inf twice, so I'm not sure what the model unknowns are.

[roborrr]

Can I assume "stunning" means that it did what you want? If so, please Accept-click the answer to indicate that it resolved your question.

Yes, "stunning" does mean that you did the initial task perfectly, which I did not even expect, for which I am very grateful. I simply forgot to click the "Accept" button, for which I apologize.

Here and elsewhere you have listed eta_inf twice, so I'm not sure what the model unknowns are.

I made a mistake in the letter, where eta_inf is repeated twice, while in place of the first eta_inf, there should be eta_0. I immediately edited the letter, apparently you read the old version of the letter. Please read the new version of the previous letter.

To see what went wrong, I would need to see what you tried.

I didn't try anything because I couldn't figure out how it all worked. How do you solve nonlinear problems with 'minL1lin', which solves the linear problem C*x=d. I tried to figure it out, but I couldn't. Is it possible to write similar codes, as you did for the model y=A*x^B+C, for the models I presented in the previous letter?

[Matt J]

Your objective function is, in general,

mean(    abs( y_calculate./y  -1 )     );

Originally, you had

y_calculate./y = A*x.^B./y  +  A./y
               = [x.^B, d] * [A;C]
               = Q * [A;C]

where we have defined the matrices,

 d=ones(size(y));
 Q=[x.^B, d]./y; 

Now take one of your new models. It can be put in the same form:

 y_calculate./y =  A* 1./(1 + (lambda * x).^m)./y  + C./y ;
                =  Q *[A;C]

but notice that Q has the new form,

      Q=[ 1./(1 + (lambda * x).^m), d]./y;        

So,

function [fval,p]=objFun(theta, x,y)

  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);

  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
 
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
 
end

[roborrr]

Following your tips I wrote a program for the model

y_calculate = A* (1 + (lambda * x).^m)+ C,

where it is necessary to calculate the best parameters

A,Lambda,m,C

or optimizing the objective functionЖ

mean( abs( y_calculate-y)/y) .

Here is the code that I wrote following your tips:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = 1;
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})
fval

function [fval,p]=objFun(theta, x,y)

  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);

  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
 
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
 
end

This code gives errors in my matlab. where am I making mistakes?

[Torsten]

You give an initial guess for only one parameter, but you try to access two parameters (lambda and m) in objFun.

[roborrr]

Dear Torsten, thanks for the useful hint.

[Matt J]

Yes, an initial guess is needed for all parameters that y_calculated has a nonlinear dependence on, in this case lambda and m.

[roborrr]

I corrected the code according to your tips and received:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})
fval
function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
 
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};

for the values ​​of the variables x,y specified in the code, the result is fine.

but I encountered values ​​of the variables x,y for which the code gives an error. for example, for the following values ​​of x,y:

 x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
 y = [21.58079135, 17.26809355, 15.3488609, 14.1074809, 13.81896005, 13.43865625, 13.49266905, 13.44, 13.44, 13.5349172, 13.4323965, 13.439328, 13.36635835, 13.365, 13.29, 13.29, 13.51, 13.365];

the code gives the following error:

Error using linprog (line 373)
A must be a real matrix.

Error in comment_3255504>minL1lin (line 147)
    [xz,~, varargout{3:nargout-1}]=linprog(f,A,b,Aeq,beq,LB,UB,args{:});

Error in comment_3255504>objFun (line 23)
 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));

Error in comment_3255504>@(theta)objFun(theta,x,y) (line 8)
error_function = @(theta) objFun(theta, x,y);

Error in fminsearch (line 336)
    fxr = funfcn(x,varargin{:});

how to make the code work correctly for all values ​​of the x,y variables?

[Torsten]

Use

Q=[ 1./(1 + (lambda^2 * x).^m), d]./y;

instead of

Q=[ 1./(1 + (lambda * x).^m), d]./y;

Of course, theta(2) must be squared to get the correct parameter "lambda" after the command

[A, theta,C]=deal(AthetaC{:})

[roborrr]

Use Q=[ 1./(1 + (lambda^2 * x).^m), d]./y;

great, this trick prevents complex numbers from appearing when negative lambda values ​​are used. I modified your idea a bit and instead of 'lambda^2' I used 'abs(lambda)' and got the lambda value directly, without the need to square it.

[Torsten]

I modified your idea a bit and instead of 'lambda^2' I used 'abs(lambda)' and got the lambda value directly, without the need to square it.

Usually, this is not a good idea since abs() is not differentiable, and most solvers for optimization need differentiability with respect to the parameters. But you are in luck: fminsearch does not belong to this class of solvers.

[roborrr]

Usually, this is not a good idea since abs() is not differentiable, and most solvers for optimization need differentiability with respect to the parameters. But you are in luck: fminsearch does not belong to this class of solvers.

Thank you, I understand you. I want the code to be universal and not depend on luck. So I will remake the code exactly as you advised.

And one more question: is it possible to impose restrictions on the parameters, for example, for the parameter C impose such a restriction:

C>0 and Abs(C-min(y))/min(y)<0.1

[Torsten]

In the call to minL1lin,

minL1lin(C,d,A,b,Aeq,beq,lb,ub,x0,options)

define the lower bound (lb) for C as 0.9*min(y) and the upper bound (ub) for C as 1.1*min(y).

[roborrr]

In the call to minL1lin,

minL1lin(C,d,A,b,Aeq,beq,lb,ub,x0,options)

define the lower bound (lb) for C as 0.9*min(y) and the upper bound (ub) for C as 1.1*min(y).

I already call 'minL1lin' once in the code

 [AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));

Do I need to call 'minL1lin' again? Where in the code should I call it minL1lin(C,d,A,b,Aeq,beq,lb,ub,x0,options)? Should I define lb and ub before or after the call?

[Matt J]

No you do not need to call minL1lin again. You must modify the call you have,

lb=[-inf,0.9,min(y)];
ub=[+inf, 1.1*min(y)];
[AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

[Torsten]

function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda^2 * x).^m), d]./y;
  %Q=[ 1./(1 + (lambda * x).^m), d]./y;
 
  %Since you also want C>0, the below setting assumes min(y)>0.
  lb = [-Inf,0.9*min(y)];
  ub = [Inf,1.1*min(y)];   
 [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));
 %[AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

[roborrr]

No you do not need to call minL1lin again. You must modify the call you have,

lb=[-inf,0.9,min(y)];

ub=[+inf, 1.1*min(y)];

[AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

Thank you for your answer. Please make changes to this design so that it also imposes restrictions on the rest of the parameters:

A>0, lambda0<lambda<lambda1, m0<m<m1

where the lambda0,lambda1,m0,m1 numbers will be preset in the code.

[Torsten]

Please make changes to this design so that it also imposes restrictions on the rest of the parameters:

A>0, lambda0<lambda<lambda1, m0<m<m1

where the lambda0,lambda1,m0,m1 numbers will be preset in the code.

For A>0, set

lb=[0,0.9*min(y)];
ub=[+inf,1.1*min(y)];

For the other two, switch from "fminsearch" to an optimizer that allows to set bounds on the parameters (e.g. "fmincon"). This way, the artificial solution to replace lambda by lambda^2 can also be avoided.

[roborrr]

For A>0, set

lb=[0,0.9*min(y)];
ub=[+inf,1.1*min(y)];

but, lb, ub we used to impose restrictions on parameter C, and where are the restrictions of parameter A? and can't restrictions of all parameters(A,lamda,m) be written into the following construction:

[AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

as Matt J did for parameter C?

[Torsten]

but, lb, ub we used to impose restrictions on parameter C, and where are the restrictions of parameter A?

The restrictions on A are the 0 and the +inf in

lb=[0,0.9*min(y)];
ub=[+inf,1.1*min(y)];

A and C are solved for by minL1lin - so you have to set bounds for them in the call to minL1lin.

lambda and m are solved for by fminsearch - so you have to set bounds for them in the call to fminsearch. But since fminsearch does not have the option to impose bounds on the parameters, you have to use a different optimizer here:

theta = fminsearch(error_function, initial_guess);

I suggested "fmincon":

theta = fmincon(error_function, initial_guess, [],[],[],[],[lambda0,m0],[lambda1,m1]);

You should try to understand the code you are using.

[roborrr]

A and C are solved for by minL1lin - so you have to set bounds for them in the call to minL1lin.

restrictions for parameter C (lb and ub) we wrote instead of five and six empty square brackets:

[AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

and where to write restrictions for parameter A? I don't know the syntax, how to write it.

[Torsten]

The code requires "fminsearchbnd" from the file exchange:

https://uk.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd-fminsearchcon

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = inf;
m0 = -inf;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
%theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})

A = 131.5017

theta = 1×2

0.2036 1.9539

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C = 73.9207

fval

fval = 0.0015

function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
  %Bounds on A and C
  A0 = 0;
  A1 = inf;
  C0 = 0.9*min(y);
  C1 = 1.1*min(y);
  lb=[A0,C0];
  ub=[A1,C1]; 
  %Call the optimizer
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));
 %[AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

[DGM]

@roborrr posted in a flag:

!!!

[DGM]

@roborrr: Flags are used to tell moderators that an item needs attention. I don't know what your flag is supposed to mean. If you want to tell @Torsten something, just post a comment.

[roborrr]

TThanks to Torsen for perfecting the code, it was trully amazing. Also thanks to Matt J who wrote the initial version of the code, that was also excellent.

The main beauty of this code is that it eliminates the occurrence of local minima during optimization, that could arise when using built-in Matlab functions.

In this code, the dependent variable 'y' depends only on one variable 'x' (one-dimensional problem). Is it possible to write code for a two-dimensional problem, when the dependent variable depends nonlinearly on two independent variables at once, using the scheme you used for the one-dimensional problem, which excludes the occurrence of local minima? Here is an example of such a problem:

The experiments were carried out on the following unique values of the variable x_un length of 5 and y_un length of 10 :

x_un = [0.4, 0.5, 0.6, 0.7, 0.8];
y_un = [283.15, 288.15, 290.65, 293.15, 295.15, 298.15, 300.65, 303.15, 313.15, 323.15];

based on x_un and y_un, new variables x and y length of 50 (5*10) are created so that the pairs (x(i),y(i)) i=1,50 exhaust all combinations of pairs (x_un, y_un):

x = [0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, ...
     0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, ...
     0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.5, 0.6, 0.7, 0.8, ...
     0.4, 0.5, 0.6, 0.7, 0.8];

y = [283.15, 283.15, 283.15, 283.15, 283.15, 288.15, 288.15, 288.15, 288.15, 288.15, ...
     290.65, 290.65, 290.65, 290.65, 290.65, 293.15, 293.15, 293.15, 293.15, 293.15, ...
     295.15, 295.15, 295.15, 295.15, 295.15, 298.15, 298.15, 298.15, 298.15, 298.15, ...
     300.65, 300.65, 300.65, 300.65, 300.65, 303.15, 303.15, 303.15, 303.15, 303.15, ...
     313.15, 313.15, 313.15, 313.15, 313.15, 323.15, 323.15, 323.15, 323.15, 323.15];

variable z is also created, the length of which is also 50, in such a way that experimental data are entered into z(i), which correspond to (x(i),y(i)) i=1,50:

z = [12.61323, 15.96843, 22.20326, 30.15335, 60.71139, 10.97453, 14.54062, 18.84638, 24.82916, ...
     53.02121, 9.57412, 13.34543, 17.50987, 22.78405, 46.60942, 8.74127, 11.59837, 14.62923, ...
     21.16607, 42.31916, 7.93542, 10.73998, 13.71756, 19.05432, 39.08479, 7.33475, 10.12368, ...
     12.67264, 18.14767, 35.87183, 6.52910, 9.07151, 11.60863, 16.87576, 32.95969, 5.43683, ...
     8.19517, 10.54588, 15.03453, 29.99945, 4.25370, 6.86437, 9.20690, 12.29640, 24.78083, ...
     2.99082, 5.05884, 6.68976, 9.54918, 18.76493];

approximation model is considered:

z_calculate = A * exp(E/ (8.314 * y)) * (1 + k * x^n);

and based on this model, an objective function is created, which is equal to the average value of the relative deviations between z and z_calculate:

error_function = mean(abs(z - z_calculate)) /z)

Our goal is to select parameters A,E,k,n in such a way that 'error_function' takes the value of the global minimum.

if possible, set lower and upper limits for parameters:

lb_A = 10^(-3); ub_A = 100;

lb_E = 2000; ub_E = 50000;

lb_k = 0.1; ub_k = 10;

lb_n = 1; ub_n = 3;

[Matt J]

Just modify the error function to take additional inputs,

error_function = @(theta) objFun(theta, x,y,z);

and construct Q this time using x,y, and z.

[roborrr]

and construct Q this time using x,y, and z.

in the linear part of the model z_calculate = A * exp(E/ (8.314 * y)) * (1 + k * x^n) there is only the coefficient A, and the free term C is missing. the nonlinear part of the model depends on two variables at once. how to construct Q? and how to replace those commands where C is present? I don't know the syntax, how to do this.

[Matt J]

The role of minL1lin is to solve for the linear parameters, for a fixed and given set of values of the nonlinear parameters. So, if A is your only linear parameter, then Q should contain coefficients only for A. Bounds on the nonlinear parameters should be given to fminsearchbnd as before.

function [fval,p]=objFun(theta, x,y,z)

  x=x(:);y=y(:);z=z(:); theta=theta(:)';
  E=theta(1); k=theta(2); n=theta(3);
  
  M=numel(z);
  
  d=ones(M,1);
  Q=  exp(E./(8.314 .* y)) .* (1 + k .* x.^n) ./z;

  lb_A=___;
  ub_A=___; 

  %Solve the linear part of the model
  [A,fval]=minL1lin(Q,d,[],[],[],[],lb_A,ub_A,[],optimoptions('linprog','Display','off'));

  fval=fval/M;

  p={A, theta};
  
  
end

[roborrr]

Thanks Matt J for the answer.

How do you optimize a model that doesn't have a linear part? For example, in the model

y = A.* (1 + (lambda * x).^m)+ C,

where A and C are vectors whose values ​​are already known and only the parameters of the nonlinear part - lambda and m - need to be optimized.

If we introduce the following designation:

y1 =( y- C)./A-1

then the original model is reduced to the model

y1= (lambda * x).^m), which doesn't have a linear part. What should we do in this case?

This model is a one-dimensional projection of a three-dimensional model onto the x-axis, so the values ​​in the vector x are repeated with different corresponding values ​​of y.

[Torsten]

How do you optimize a model that doesn't have a linear part?

In this case, you delete the lines of code in "objFun" in which you solve the linear part of the model and return e.g. sum((y1-(lambda*x).^m).^2) to "fminsearchbnd".

[roborrr]

Thanks Torsten for the answer, but I'm again confused by the code syntax. Here is the code you created for the model, which contains linear coefficients:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = inf;
m0 = -inf;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
%theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})
fval
function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
  %Bounds on A and C
  A0 = 0;
  A1 = inf;
  C0 = 0.9*min(y);
  C1 = 1.1*min(y);
  lb=[A0,C0];
  ub=[A1,C1]; 
  %Call the optimizer
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));
 %[AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

as I understand, you suggest to replace "minL1lin" with "fminsearchbnd" in the "objFun" subroutine. but "fminsearchbnd" is already used in the main part of the code to determine the boundaries of the "lambda" and "m" parameters. I don't know the syntax of how to combine these two uses of "fminsearchbnd", how to output the optimized values of "lambda" and "m" and the error function "fval", for a model that does not have a linear part.

If it is not too much trouble, Please remake the code you sent which is written for the model y = A.* (1 + (lambda * x).^m)+ C, for the model y1= (lambda * x).^m)

[Torsten]

Why do you still solve for A and C if you say their values are known ?

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = inf;
m0 = -inf;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:);
  lambda=theta(1); m=theta(2);
  % Set A and B to appropriate values !
  A = 1;    
  C = 1;
  fval = sum(abs((y - (A.* (1 + (lambda * x).^m)+ C))./y));
end

[roborrr]

You guessed correctly: I wanted only lambda and m to be optimized, with A and C known, And you sent the code i wanted.

but I tested your code as follows:

first I solved the optimization problem for the model y=A.*(1+(lambda*x).^m)+C using minL1lin, assuming first that A and C are unknown and should be optimized jointly lambda and m:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];

initial_guess = [1,1];

% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 =0;
lambda1=inf;


m0 =0;
m1 =inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
%theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})
fval
function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
  %Bounds on A and C
  A0 = 0;
  A1 = inf;
  C0 =  0.2*min(y);
  C1 = 1.1*min(y);
  lb=[A0,C0];
  ub=[A1,C1]; 
  %Call the optimizer
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));
 %[AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
 
 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

and I got the following results:

 A = 131.5161945547369  
 theta = 0.2037    1.9538
C = 73.9206126555175
fval = 0.0015

after that the optimized values A= 131.5161945547369 and C= 73.9206126555175 (obtained using minL1lin) I used as previously known coefficients and optimized only lambda and m parameters using fminsearchbnd:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = inf;
m0 = 0;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:);
  lambda=theta(1); m=theta(2);
  % Set A and B to appropriate values !
  A = 131.5161945547369;   
  C = 73.9206126555175;
  fval = sum(abs((y - (A.* (1 + (lambda * x).^m)+ C))./y));
  
  fval= fval/numel(y);
end

and got a much worse result for Fval:

fval = 3.3576

why the result of the first script - fval = 0.0015 and the result of the second script - fval = 3.3576 are not equal? where is the error in the second script?

[Torsten]

I think the error is in the first script. If you take

y=A.*(1+(lambda*x).^m)+C

as your model equation, you should take

Q=[ 1 + (lambda * x).^m, d]./y;

instead of

Q=[ 1./(1 + (lambda * x).^m), d]./y;

shouldn't you ?

Or you change the objective of your second script to

function fval=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:);
  lambda=theta(1); m=theta(2);
  % Set A and B to appropriate values !
  A = 131.5161945547369;   
  C = 73.9206126555175;
  fval = sum(abs((y - (A./(1 + (lambda * x).^m)+ C))./y));
  
  fval= fval/numel(y);
end

In this case, your model function is changed to

y=A./(1+(lambda*x).^m)+C

[roborrr]

thank you dear Torshen for correcting the error in my second script - instead of model function y=A.*(1+(lambda*x).^m)+C, there should indeed be function y=A./(1+(lambda*x).^m)+C:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 =0;
lambda1=inf;


m0 =0;
m1 =inf;
  A0 = 0;
  A1 = inf;
  C0 =  0.2*min(y);
  C1 = 1.1*min(y);
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:);
  lambda=theta(1); m=theta(2);
 
  % Set A and B to appropriate values !
   A = 131.5161945547369;   
   C = 73.9206126555175;
  fval = sum(abs((y - (A./(1 + (lambda * x).^m)+ C))./y));
  
  fval= fval/numel(y);
end

This construct does give the best result - Fval = 0.015, when the best coefficients A and C are known in advance. But I tried to use this construct to optimize all coefficients A, C, lambda and m and got the worst result - Fval = 0.0240. Here is the script I wrote:

clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter  theta
initial_guess = [1,1,1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 =0;
lambda1=inf;


m0 =0;
m1 =inf;
  A0 = 0;
  A1 = inf;
  C0 =  0.2*min(y);
  C1 = 1.1*min(y);
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0,A0,C0],[lambda1,m1,A1,C1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:);
  lambda=theta(1); m=theta(2);
  A=theta(3); 
  C=theta(4);
  % the best values ​​of A and C that should be obtained in this script
%    A = 131.5161945547369;   
%    C = 73.9206126555175;
  fval = sum(abs((y - (A./(1 + (lambda * x).^m)+ C))./y));
  
  fval= fval/numel(y);
end

why does the value - Fval = 0.015, obtained from the first script not match the value obtained from the second script - Fval = 0.0240? where am I making a mistake in the second script?

how universal is this construct? is it possible to use it for models of any complexity? or is the class of models limited?

[Torsten]

You tried two solution strategies: Solving for all 4 parameters simultaneously (first strategy) and solving for the 2 linear and the 2 nonlinear coefficients separately with result coupling (second strategy). The second strategy turned out to be better and more robust - especially because no initial guesses for A and C are needed (your values A = 1 and C = 1 used for the first strategy most probably are too far off).

That's life in numerical analysis: usually, you have to go the way of trial and error for each problem you want to solve, and generalizations are not possible.

[roborrr]

I use construct:

  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

matlab gives a message:

Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.012493 

What is the default maximum iteration count in aa and how to increase it?

and also: in some cases when using this construct matlab gives the following errors:

         Error using linprog (line 346)
LINPROG encountered an internal error that caused the solution to lose feasibility. We are sorry for the inconvenience.

Please contact technical support for assistance with your problem, quoting the code "-2_4".

Error in minL1lin (line 116)
    [xz,~, varargout{3:nargout-1}]=linprog(f,A,b,Aeq,beq,LB,UB,args{:});

Error in A_A_optimize_problem_A_B_C_lambda_m/objFun (line 42)
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

Error in A_A_optimize_problem_A_B_C_lambda_m>@(theta)objFun(theta,x,y) (line 9)
error_function = @(theta) objFun(theta, x,y);

Error in fminsearchbnd>intrafun (line 259)
fval = feval(params.fun,reshape(xtrans,params.xsize),params.args{:});

Error in fminsearch (line 324)
    x(:) = xr; fxr = funfcn(x,varargin{:});

Error in fminsearchbnd (line 229)
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

Error in A_A_optimize_problem_A_B_C_lambda_m (line 19)
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);

[Torsten]

Please include the full code to reproduce the error.

[Matt J]

Error in A_A_optimize_problem_A_B_C_lambda_m/objFun (line 42)
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

We need to examine Q,d, lb,ub when this happens.

[roborrr]

MJ: "We need to examine Q,d, lb,ub when this happens."

Below is the script:

clear
% Data
% the first series of x , y values, at which the script works successfully and gives the result fval=0.0015:

% x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
% y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];


%the second series of x , y values, at which the script also works successfully, but gives the message:

% x = [14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.6740, 66, ...
%      73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
% y = [8.80982634, 9.359883, 8.83235329, 9.05116997, 8.640265855, 8.87955602, 8.90198004, 8.351331895, ...
%      8.435, 8.23, 8.22833633, 8.15836833, 8.0604534, 8.15348853, 7.671930455, 8.346546765, 8.15];

% the third series of x, y values ​​at which the script fails and returns an error


     x = [14.6784; 22.0176; 29.3436; 36.6828; 44.0088; 51.3348; 58.6740; 66; ...
         73.3392; 80.6652; 88.0044; 95.3304; 102.6696; 109.9956; 117.3348; 124.6608; 132];
    y = [6.53; 6.714628295; 5.75798471; 5.75746671; 5.75884823; 6.17125225; 5.75928009; 6.4; ...
         6.0475162; 6.28375061; 5.75971201; 5.981722515; 6.17125225; 6.34975669; 6.325585625; ...
         6.304601135; 6.285626125];

% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = 3;
m0 = 0;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
%theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})
fval
function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
  %Bounds on A and C
  A0 = 0;
  A1 = inf;
  C0 = 0.9*min(y);
  C1 = 1.1*min(y);
  lb=[A0,C0];
  ub=[A1,C1]; 
  %Call the optimizer
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

In this script there are three series of values ​​of variables x, y. One of them is active and the rest are commented out.

when the first series is active, then the script works successfully and gives the result fval=0.0015.

when the second series is active, then the script also works successfully, but gives the message:

Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.015660 

when the third series is active, then the script fails and returns an error:

Error using linprog (line 346)
LINPROG encountered an internal error that caused the solution to lose feasibility. We are sorry for the inconvenience.

Please contact technical support for assistance with your problem, quoting the code "-2_4".

Error in minL1lin (line 116)
    [xz,~, varargout{3:nargout-1}]=linprog(f,A,b,Aeq,beq,LB,UB,args{:});

Error in A_A_example_opt>objFun (line 54)
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

Error in A_A_example_opt>@(theta)objFun(theta,x,y)

Error in fminsearchbnd>intrafun (line 259)
fval = feval(params.fun,reshape(xtrans,params.xsize),params.args{:});

Error in fminsearch (line 324)
    x(:) = xr; fxr = funfcn(x,varargin{:});

Error in fminsearchbnd (line 229)
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

Error in A_A_example_opt (line 36)
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);

[Torsten]

I don't get errors in all the three cases. You should update to the most recent MATLAB release.

clear
% Data
% the first series of x , y values, at which the script works successfully and gives the result fval=0.0015:

% x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
% y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];


%the second series of x , y values, at which the script also works successfully, but gives the message:

% x = [14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.6740, 66, ...
%      73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
% y = [8.80982634, 9.359883, 8.83235329, 9.05116997, 8.640265855, 8.87955602, 8.90198004, 8.351331895, ...
%      8.435, 8.23, 8.22833633, 8.15836833, 8.0604534, 8.15348853, 7.671930455, 8.346546765, 8.15];

% the third series of x, y values at which the script fails and returns an error


     x = [14.6784; 22.0176; 29.3436; 36.6828; 44.0088; 51.3348; 58.6740; 66; ...
         73.3392; 80.6652; 88.0044; 95.3304; 102.6696; 109.9956; 117.3348; 124.6608; 132];
    y = [6.53; 6.714628295; 5.75798471; 5.75746671; 5.75884823; 6.17125225; 5.75928009; 6.4; ...
         6.0475162; 6.28375061; 5.75971201; 5.981722515; 6.17125225; 6.34975669; 6.325585625; ...
         6.304601135; 6.285626125];

% Initial guess for parameter  theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m 
lambda0 = 0;
lambda1 = 3;
m0 = 0;
m1 = inf;
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
%theta = fminsearch(error_function, initial_guess);
[fval,AthetaC]=error_function(theta);
[A, theta,C]=deal(AthetaC{:})

A = 0.3587

theta = 1×2

0.0418 461.2531

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C = 6.1713

fval

fval = 0.0343

function [fval,p]=objFun(theta, x,y)
  x=x(:);y=y(:); theta=theta(:)';
  lambda=theta(1); m=theta(2);
  d=ones(size(y));
  Q=[ 1./(1 + (lambda * x).^m), d]./y;
  %Bounds on A and C
  A0 = 0;
  A1 = inf;
  C0 = 0.9*min(y);
  C1 = 1.1*min(y);
  lb=[A0,C0];
  ub=[A1,C1]; 
  %Call the optimizer
  [AC,fval]=minL1lin(Q,d,[],[],[],[],lb,ub,[],optimoptions('linprog','Display','off'));

 fval=fval/numel(y);
 
 A=AC(1); C=AC(2);
 
 p={A, theta, C};
end

function [x,varargout]=minL1lin(C,d,varargin)
% The submission minL1lin finds the  minimum L1-norm solution of the linear equations C*x=d, 
% optionally under linear constraints. It is similar to the Optimization Toolbox's lsqlin except that it minimizes with 
% respect to the L1 norm by reformulating the problem as a linear program. The input/output syntax,
%  
%    [x,resnorm,residual,exitflag,output,lambda] = 
%           minL1lin(C,d,A,b,Aeq,beq,lb,ub,x0,options)
%  
%  is the same as for lsqlin() and all the usual defaults for A,b,..,x0 from
%  lsqlin apply here. However, the "options" are those of linprog() which is used 
%  internally. And, of course, the minimization is instead over norm(C*x-d,1)
%  rather than norm(C*x-d,2).
%
%
%
%EXAMPLES:
%
% We first construct some data for the examples
% 
% C=rand(3);
% 
% Xtrue=(1:3)'; 
%  noise=rand(3,1)*.1;
%  
% d=C*Xtrue+noise;
%  dTrue=C*Xtrue;
% 
%     %EXAMPLE 1:
%     %
%     %  This is an example of an unconstrained problem, in which we check the optimality
%     %  of the solution by simulation. First, we find the minimum L1 solution
%     %  of C*X-d
% 
%             [X,resnormX] = minL1lin(C,d,[],[],[],[],[],[],[],optimset('Display','none'));
% 
%     %  Now, we randomly generate many neighboring X and test the value of the objective
%     %  function in this neighborhood,
% 
%             Xneighbors=bsxfun(@plus,X,randn(3,1e6)*1e-6); %explore neighborhood of X
% 
%             resnormNeighbors=sum(abs(  bsxfun(@minus,  C*Xneighbors, d) )); %L1 error of neighbors
% 
%      % A value of 1 indicates all neighbors have a sub-optimal objective
%      
%             isOptimal=~any(resnormNeighbors<resnormX);
%             
%      % Typically it returns 1, indicating that all the neighbors are less
%      % optimal.
%      
%      isOptimal=~any(resnormNeighbors<resnormX)
% 
% 
%     %%EXAMPLE 2: 
%     %
%     %  In this continuation of EXAMPLE 1, we add linear constraints, and test the
%     %  improvement that this gives in the L1 error relative to the
%     %  noise-free true data, dTrue. For variety's sake, we now use the dual simplex
%     %  algorithm for this example,
% 
%                 Aeq=[1,1,1];
%                 beq=6;
%                 lb=[1;1;1];
%                 ub=[3;3;3];
% 
%                  options=optimoptions(@linprog,'Algorithm','dual-simplex','Display','none');
% 
% 
%             Xcon = minL1lin(C,d,[],[],Aeq,beq,lb,ub,[],options);
% 
%             error_Xunc = norm(C*X - dTrue,1),
%             error_Xcon = norm(C*Xcon - dTrue,1),
%             
%             
%       %Typically, one will observe error_Xcon < error_Xunc 
% 
% 

if length(varargin)<8
   varargin(end+1:8)={[]}; 
end

[A,b,Aeq,beq,lb,ub,x0,options]=deal(varargin{:});

[m,n]=size(C);

    
    f=[zeros(1,n), ones(1,m)];
    
   
    Ax=[C,-speye(m);-C,-speye(m)];
    bx=[d;-d];

  [~,nx]=size(Ax);   
 
    LB(1:nx)=-inf; 
       LB(n+1:end)=0;
    UB(1:nx)=inf;

    LB(1:length(lb))=lb;
    UB(1:length(ub))=ub;
    
    if ~isempty(A),
       A(end,nx)=0;
    end
    
    if ~isempty(Aeq),
       Aeq(end,nx)=0;
    end
    
    A=[A;Ax]; b=[b;bx];
    
   try
    [xz,~, varargout{3:nargout-1}]=linprog(f,A,b,Aeq,beq,LB,UB,x0,options);
   catch
    args={}; if ~isempty(options), args{1}=options; end
    [xz,~, varargout{3:nargout-1}]=linprog(f,A,b,Aeq,beq,LB,UB,args{:});
   end
    
    if ~isempty(xz)
     x=xz(1:n);
    else
      x=xz;
      varargout(1:2)={[],[]}; 
      return
    end
    
    
    if nargout>=2 %truncate resnorm if requested
       residual=C*x-d;
       resnorm=norm(residual,1);
       varargout(1:2)={resnorm,residual}; 
    end
    
end   

function [x,fval,exitflag,output] = fminsearchbnd(fun,x0,LB,UB,options,varargin)
% FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
% usage: x=FMINSEARCHBND(fun,x0)
% usage: x=FMINSEARCHBND(fun,x0,LB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
% usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
% 
% arguments:
%  fun, x0, options - see the help for FMINSEARCH
%
%  LB - lower bound vector or array, must be the same size as x0
%
%       If no lower bounds exist for one of the variables, then
%       supply -inf for that variable.
%
%       If no lower bounds at all, then LB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
%  UB - upper bound vector or array, must be the same size as x0
%
%       If no upper bounds exist for one of the variables, then
%       supply +inf for that variable.
%
%       If no upper bounds at all, then UB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
% Notes:
%
%  If options is supplied, then TolX will apply to the transformed
%  variables. All other FMINSEARCH parameters should be unaffected.
%
%  Variables which are constrained by both a lower and an upper
%  bound will use a sin transformation. Those constrained by
%  only a lower or an upper bound will use a quadratic
%  transformation, and unconstrained variables will be left alone.
%
%  Variables may be fixed by setting their respective bounds equal.
%  In this case, the problem will be reduced in size for FMINSEARCH.
%
%  The bounds are inclusive inequalities, which admit the
%  boundary values themselves, but will not permit ANY function
%  evaluations outside the bounds. These constraints are strictly
%  followed.
%
%  If your problem has an EXCLUSIVE (strict) constraint which will
%  not admit evaluation at the bound itself, then you must provide
%  a slightly offset bound. An example of this is a function which
%  contains the log of one of its parameters. If you constrain the
%  variable to have a lower bound of zero, then FMINSEARCHBND may
%  try to evaluate the function exactly at zero.
%
%
% Example usage:
% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%
% fminsearch(rosen,[3 3])     % unconstrained
% ans =
%    1.0000    1.0000
%
% fminsearchbnd(rosen,[3 3],[2 2],[])     % constrained
% ans =
%    2.0000    4.0000
%
% See test_main.m for other examples of use.
%
%
% See also: fminsearch, fminspleas
%
%
% Author: John D'Errico
% E-mail: [e-mail]
% Release: 4
% Release date: 7/23/06

% size checks
xsize = size(x0);
x0 = x0(:);
n=length(x0);

if (nargin<3) || isempty(LB)
  LB = repmat(-inf,n,1);
else
  LB = LB(:);
end
if (nargin<4) || isempty(UB)
  UB = repmat(inf,n,1);
else
  UB = UB(:);
end

if (n~=length(LB)) || (n~=length(UB))
  error 'x0 is incompatible in size with either LB or UB.'
end

% set default options if necessary
if (nargin<5) || isempty(options)
  options = optimset('fminsearch');
end

% stuff into a struct to pass around
params.args = varargin;
params.LB = LB;
params.UB = UB;
params.fun = fun;
params.n = n;
% note that the number of parameters may actually vary if 
% a user has chosen to fix one or more parameters
params.xsize = xsize;
params.OutputFcn = [];

% 0 --> unconstrained variable
% 1 --> lower bound only
% 2 --> upper bound only
% 3 --> dual finite bounds
% 4 --> fixed variable
params.BoundClass = zeros(n,1);
for i=1:n
  k = isfinite(LB(i)) + 2*isfinite(UB(i));
  params.BoundClass(i) = k;
  if (k==3) && (LB(i)==UB(i))
    params.BoundClass(i) = 4;
  end
end

% transform starting values into their unconstrained
% surrogates. Check for infeasible starting guesses.
x0u = x0;
k=1;
for i = 1:n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      if x0(i)<=LB(i)
        % infeasible starting value. Use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(x0(i) - LB(i));
      end
      
      % increment k
      k=k+1;
    case 2
      % upper bound only
      if x0(i)>=UB(i)
        % infeasible starting value. use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(UB(i) - x0(i));
      end
      
      % increment k
      k=k+1;
    case 3
      % lower and upper bounds
      if x0(i)<=LB(i)
        % infeasible starting value
        x0u(k) = -pi/2;
      elseif x0(i)>=UB(i)
        % infeasible starting value
        x0u(k) = pi/2;
      else
        x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
        % shift by 2*pi to avoid problems at zero in fminsearch
        % otherwise, the initial simplex is vanishingly small
        x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
      end
      
      % increment k
      k=k+1;
    case 0
      % unconstrained variable. x0u(i) is set.
      x0u(k) = x0(i);
      
      % increment k
      k=k+1;
    case 4
      % fixed variable. drop it before fminsearch sees it.
      % k is not incremented for this variable.
  end
  
end
% if any of the unknowns were fixed, then we need to shorten
% x0u now.
if k<=n
  x0u(k:n) = [];
end

% were all the variables fixed?
if isempty(x0u)
  % All variables were fixed. quit immediately, setting the
  % appropriate parameters, then return.
  
  % undo the variable transformations into the original space
  x = xtransform(x0u,params);
  
  % final reshape
  x = reshape(x,xsize);
  
  % stuff fval with the final value
  fval = feval(params.fun,x,params.args{:});
  
  % fminsearchbnd was not called
  exitflag = 0;
  
  output.iterations = 0;
  output.funcCount = 1;
  output.algorithm = 'fminsearch';
  output.message = 'All variables were held fixed by the applied bounds';
  
  % return with no call at all to fminsearch
  return
end

% Check for an outputfcn. If there is any, then substitute my
% own wrapper function.
if ~isempty(options.OutputFcn)
  params.OutputFcn = options.OutputFcn;
  options.OutputFcn = @outfun_wrapper;
end

% now we can call fminsearch, but with our own
% intra-objective function.
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

% undo the variable transformations into the original space
x = xtransform(xu,params);

% final reshape to make sure the result has the proper shape
x = reshape(x,xsize);

% Use a nested function as the OutputFcn wrapper
  function stop = outfun_wrapper(x,varargin);
    % we need to transform x first
    xtrans = xtransform(x,params);
    
    % then call the user supplied OutputFcn
    stop = params.OutputFcn(xtrans,varargin{1:(end-1)});
    
  end

end % mainline end

% ======================================
% ========= begin subfunctions =========
% ======================================
function fval = intrafun(x,params)
% transform variables, then call original function

% transform
xtrans = xtransform(x,params);

% and call fun
fval = feval(params.fun,reshape(xtrans,params.xsize),params.args{:});

end % sub function intrafun end

% ======================================
function xtrans = xtransform(x,params)
% converts unconstrained variables into their original domains

xtrans = zeros(params.xsize);
% k allows some variables to be fixed, thus dropped from the
% optimization.
k=1;
for i = 1:params.n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      xtrans(i) = params.LB(i) + x(k).^2;
      
      k=k+1;
    case 2
      % upper bound only
      xtrans(i) = params.UB(i) - x(k).^2;
      
      k=k+1;
    case 3
      % lower and upper bounds
      xtrans(i) = (sin(x(k))+1)/2;
      xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
      % just in case of any floating point problems
      xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));
      
      k=k+1;
    case 4
      % fixed variable, bounds are equal, set it at either bound
      xtrans(i) = params.LB(i);
    case 0
      % unconstrained variable.
      xtrans(i) = x(k);
      
      k=k+1;
  end
end

end % sub function xtransform end

[roborrr]

Torsten: I don't get errors in all the three cases. You should update to the most recent MATLAB release.

Thanks for the advice. I will install the latest version of matlab.