the 'fmincon' optimisation doesn't stop at minimum with all possible algorithms

I've got an constrained optimisation problem, thus I adopted the function ‘fmincon’.

The curve of the f(x) is like the figure below (This is just 1 specific case as the coefficients of f(x) would vary from case to case but the form is the same) , and we can see the minimum point is around the labelled one at x=7.1e-12.

zoom in →

However, the optimisation would stop at x=5e-12 (with initial x=1e-11) because the ‘step tolerance’ is satisfied with the default algorithm 'interior-point’. After that, no matter how I decrease the step tolerance (shown in the figure below), the function just yields the same results (x=5e-12).

Then I tried algorithm ’sqp’, the step length became too small at the 3rd row, and the function just cannot get out from ‘fval=1.017079’ after that (don't forget the f_min=0.09793 approximately at x=7.1e-12), and it seems the optimisation would just never stop (shown below).

Finally, I tried 'active-set’ algorithm, it also gives me x=5e-12. Could anyone please tell me how to push the optimisation closer to the optimal point (x=7.1e-12)? Peronally the initial point is OK. Did I do something wrong in the optimisation?

%% 29.10.2020 optimisation codes
clear
syms x;
size=1e-11;
C_Gmm1=[441400000,298900000000.000];
C_Gmm2=[4140000,34650000000.0000];
fun_2Gmm=@(t) C_Gmm1(1).*t.*exp(-C_Gmm1(2).*t) + C_Gmm2(1).*t.*exp(-C_Gmm2(2).*t);% double gamma function

term_1=@(x) 2.1e-6 + 2.433093685001734e+06.*x;
term_2=@(x) integral(fun_2Gmm,0,x)./size;
f_x=@(x) 2.*term_1(x)./(4.*term_1(x)+term_2(x));% function to be optimised w.r.t. x
%% plot the function w.r.t. the x
fig_1=figure();
fplot(@(x) f_x(x),[0 5.*size],'b');% see where the minimum is before optimisation
%% fmincon optimisation
x0_con=size;

A_fmin=[];
b=[];Aeq=[];beq=[];
lb=0; % x should be bigger than 0
ub=[];nonlcon=[];
options_con = optimoptions('fmincon','Algorithm','interior-point','Display','iter-detailed','OptimalityTolerance', 1e-26,'ConstraintTolerance', 1e-7, 'StepTolerance', 1e-34, 'FunctionTolerance',1e-15,'MaxFunctionEvaluations', 1e10, 'MaxIterations', 1e11);

[x_optmzd_con,Fval_con,exitflag_con,output_con] =fmincon(f_x,x0_con,A_fmin,b,Aeq,beq,lb,ub,nonlcon,options_con);

[Walter Roberson]

Can you attach your code and data for testing purposes?

[Shuangfeng Jiang]

Okay, I've just attached my codes, and it can be run diretly.

Answer 1

[Matt J]

Change the units of your unknowns to a less extreme order of magnitude. Also, fmincon is way overkill for a 1D problem. Use fminbnd instead,

size=1e-11;
C_Gmm1=[441400000,298900000000.000];
C_Gmm2=[4140000,34650000000.0000];
fun_2Gmm=@(t) C_Gmm1(1).*t.*exp(-C_Gmm1(2).*t) + C_Gmm2(1).*t.*exp(-C_Gmm2(2).*t);% double gamma function

term_1=@(x) 2.1e-6 + 2.433093685001734e+06.*x;
term_2=@(x) integral(fun_2Gmm,0,x)./size;
f_x=@(x) 2.*term_1(x)./(4.*term_1(x)+term_2(x));% function to be optimised w.r.t. x


s=1e11; %unit-changing scale factor
[x_optmzd_con,Fval_con,exitflag_con,output_con] =fminbnd( @(x)f_x(x/s), 0,5);

x_optmzd_con=x_optmzd_con/s

x_optmzd_con = 7.0782e-12

Answer 2

[Walter Roberson]

Change the FiniteDifferenceStepSize

%% 29.10.2020 optimisation codes
close all
syms x;
size=1e-11;
C_Gmm1=[441400000,298900000000.000];
C_Gmm2=[4140000,34650000000.0000];
fun_2Gmm=@(t) C_Gmm1(1).*t.*exp(-C_Gmm1(2).*t) + C_Gmm2(1).*t.*exp(-C_Gmm2(2).*t);% double gamma function

term_1=@(x) 2.1e-6 + 2.433093685001734e+06.*x;
%fplot is for some reason invoking the first time with [-5 0 5]
%and then it complains about the nan that shows up for negative values
term_2=@(x) arrayfun(@(X)integral(fun_2Gmm,0,X),max(0,x))./size;

f_x=@(x) 2.*term_1(x)./(4.*term_1(x)+term_2(x));% function to be optimised w.r.t. x
%% plot the function w.r.t. the x
fig_1=figure();
fplot(@(x) f_x(x),[0 5.*size],'b');% see where the minimum is before optimisation
%% fmincon optimisation
x0_con=size;

A_fmin=[];
b=[];Aeq=[];beq=[];
lb=0; % x should be bigger than 0
ub=5e-11;nonlcon=[];
options_con = optimoptions('fmincon', ...
    'Algorithm', 'interior-point',    ...
    'Display', 'iter-detailed',       ...
    'OptimalityTolerance', 1e-16,     ...
    'ConstraintTolerance', 1e-8,      ...
    'StepTolerance', 1e-20,           ...
    'FunctionTolerance',1e-17,        ...
    'MaxFunctionEvaluations', 1e10,   ...
    'MaxIterations', 1e11,            ...
    'FiniteDifferenceStepSize', 1e-13); ... 'Diagnostics', 'off');

[x_optmzd_con,Fval_con,exitflag_con,output_con] =fmincon(f_x,x0_con,A_fmin,b,Aeq,beq,lb,ub,nonlcon,options_con);
disp(x_optmzd_con)
disp(Fval_con)

[Shuangfeng Jiang]

Thanks for your reply, and it works after dereasing the FiniteDifferenceStepSize!

Cold you pease explain the differences between 'FiniteDifferenceStepSize' and 'StepTolerance'? Does the optimisation strictly follow the 'FiniteDifferenceStepSize' we've set at each iteration?

[Walter Roberson]

FiniteDifferenceStepSize only applies during the Finite Differences phase that is used to estimate the gradient.

Unfortunately for a few years now, the key implementing routines are compiled in, which perhaps is more efficient but means that I cannot dig into how the variables fit together.

[Shuangfeng Jiang]

Well, do you mean there're different phases in an optimisation process where specfic options are called? Could you please tell me where to find the corresponding help page for the explanation of the different phases?