Minimization of a function using fmincon with no constraints vs. using fminbnd

I am trying to minimize an equation F(r), using both 'fmincon' and 'fminbnd'.

By using 'fminbnd', I got certain r values optimizing F(r) (r=r_opt) and had no special problem.

Trying to solve the function with constaints, I used 'fmincon'. I added a constraints function [c,ceq] = heightconst(r), with nonlinear inequality constraints c(r)<=0, but it appeared to return solutions that do not converge well with the expected ones.

Therefore, in order to isolate to problem, I used 'fmincon', but this time using no constraints at all, relying on the following syntax:

x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)

with A = [], b = [], Aeq =[], beq = [], and appropriate values of lb, ub.

I expected to get similiar values of r, such the ones I got while using non-constrained minimization ('fminbnd'), but I failed to. Instead, I got the very same values which were produced using the syntax:

x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

i.e. the same values of the minimization with the non-linear inequality.

Is there any explanation for this?

Attached is the script I used:

function single_bnd_con_comp

%Parameters for the system:

Rgc=100;
kappa=20;
a =  linspace(0.0000001,4*pi*Rgc^2,200);
kappa_t = 50;
r0=0;
theta_0 = 0.5*pi;
lambda = 5;

%__________________________________________________________________________________________________________
%options for minimization functions:

opts = optimset('MaxIter',2e10,'MaxFunEvals',2e10,'Display','off');
options = optimoptions(@fmincon,'MaxIter',2e10,'MaxFunEvals', 2e10,'StepTolerance',1e-20,'Display','off');
%__________________________________________________________________________________________________________
%Preallocation of arrays:
r_opt_mat=zeros(length(kappa_t),length(a));
height_mat=zeros(length(kappa_t),length(a));
zmax_mat=zeros(length(kappa_t),length(a));
const_r_opt_mat=zeros(length(kappa_t),length(a));
const_height_mat=zeros(length(kappa_t),length(a));
%__________________________________________________________________________________________________________

    function [c,ceq] = heightconst(r)% constraint function.
        z_max_func=(Rgc+sqrt(Rgc^2-r^2));%
        z_func=sqrt(A/pi-r^2);
        c = z_func-z_max_func;%non linear inequality.
        ceq = [];
    end
%__________________________________________________________________________________________________________
    % nested_f - A function used to minimize F(r).
    %Returns normzlied variables which minimize F(r).

    function [norm_r_opt,z,const_norm_r_opt,const_z,z_max_const]=nested_f(lambda, A , k_t , theta_0)
        
        %The function to be minimized:
        F=@(r,lambda, A , k_t , theta_0) 2.*pi.*(lambda/kappa).*r -8*pi*(pi*r^2 /A - 1)+(k_t/kappa)*pi*r*(acos(2*pi*r.^2 /(A)-1)-theta_0).^2;
        func=@(r)F(r,lambda, A , k_t , theta_0);% function handle (in single variable).
%__________________________________________________________________________________________________________
        
        r_disc=sqrt(A/pi);
        
        %Setting lower and upper boundaries:
        
        lb=r0;%an initial point satisfying all the constraints.
        if r_disc<Rgc
            ub=r_disc;
        else
            ub=Rgc;
        end
        
        %Minimization - without any constraints: 
        r_opt = fminbnd(func,lb,ub,opts);
        z = sqrt(A/pi-r_opt^2);
%__________________________________________________________________________________________________________
        %Minimization with constraints:

        %There are no linear constraints so A,b, Aeq and beq are null arrays:
        G = [];
        b = [];
        Geq = [];
        beq = [];
        nonlcon=@heightconst;
        
        const_r_opt = fmincon(func,r0,G,b,Geq,beq,lb,ub,nonlcon,options);
        const_z = sqrt(A/pi-const_r_opt^2);
        z_max_const=Rgc+sqrt(Rgc^2-const_r_opt^2);

        %normalized variables:
        const_norm_r_opt=const_r_opt/r_disc;
        norm_r_opt=r_opt/r_disc;
    end
%__________________________________________________________________________________________________________


% Finding solutions for F(r) for each A:

for i=[1:1:length(a)]% an column index for the area.
    k_t = kappa_t;
    A =a(i);
    
    %Calling the minimization function:
    [norm_r_opt,z,const_norm_r_opt,const_z,z_max_const]=nested_f(lambda, A , k_t , theta_0);

    %Assigning outputs of the non-constrained function:
    height_mat(i)=z;
    r_opt_mat(i)=norm_r_opt;%Calculate the free energy, using indiced of each vectors as described.

    %Assigning outputs of the constrained function:
  const_height_mat(i)=const_z;
    zmax_mat(i)=z_max_const;
    const_r_opt_mat(i)=const_norm_r_opt;%Calculate the free energy, using indiced of each vectors as described.
end
%__________________________________________________________________________________________________________

norm_a=a./(pi*Rgc^2);

%Plot z_opt:
figure(1)
plot(norm_a,zmax_mat,'+',norm_a,height_mat,norm_a,const_height_mat);
xlabel('A (\pi*r_{gc}^2 [nm^{2}])','FontWeight','bold');
ylabel('z [nm]','FontWeight','bold');
legend('z_{max}','z_{unconst.}','z_{const}');

title({...
    ['\fontsize{15}  z'];...
    ['\fontsize{12}\color{black} \kappa_{m}=\color{red}' num2str(kappa) '\color{black}  [kT] ,',...
    '  \color{black}\lambda=\color{red}' num2str(lambda),'\color{black}  [kT/nm] ,',...
    '  \color{black}\theta_{0}=\color{red}' num2str(theta_0/pi),'\pi',...
    ]})

%Plot r_opt:
figure(2)
plot(norm_a,const_r_opt_mat,norm_a,r_opt_mat);
xlabel('A (\pi*r_{gc}^2 [nm^{2}])','FontWeight','bold');
ylabel('r_{opt} [nm]','FontWeight','bold');
legend('r_{opt}^{const.}','r_{opt}^{unconst.}');

title({...
    ['\fontsize{15}  r_{opt}'];...
    ['\fontsize{12}\color{black} \kappa_{m}=\color{red}' num2str(kappa) '\color{black}  [kT] ,',...
    '  \color{black}\lambda=\color{red}' num2str(lambda),'\color{black}  [kT/nm] ,',...
    '  \color{black}\theta_{0}=\color{red}' num2str(theta_0/pi),'\pi',...
    ]})

end

Answer 1

[Swatantra Mahato]

Hi,

Kindly refer to the limitations of the fmincon function below:

https://www.mathworks.com/help/optim/ug/fmincon.html#busp5fq-7

In particular, the first order derivatives of both the objective as well as constraint functions must be continuous.

It may be helpful to look into the Symbolic Math Toolbox to calculate the derivative of a function.

The points of discontinuity of a function f of the symbolic variable x can be calculated as

feval(symengine,'discont',f,x)

Additionally, I observed that varying r0 and lb gives different results for fmincon with and without constraints.

For example,

nonlcon=@heightconst; %with constraints
%nonlcon=[]; %without constraints
const_r_opt = fmincon(func,1e-7,G,b,Geq,beq,-Inf,ub,nonlcon,options);

gives different plots for const_r_opt_mat based on the value of nonlcon.

Hope this helps