MATLAB: Issues about Using Function ‘fmincon’ solve Optimization Problem

fminconoptimizationspeed slow

Dear all, I am using 'fmincon' to solve the following optimization problem:
I need to figure out the optimal and correponding values which minimize the objective function.
My code is as below in which I use the symbolic mathematic tool to create the constraint conditions and their derivative for 'fmincon' and turn to 'fmincon'. In this code, I want to solve out the optimal respectively, e.g the optimal given .
clear
%% Parameters
%% Using symbolic mathematic tools to create Constraint conditions
syms Z_H psi_HH psi_HF psi_FF psi_FH N_H N_F L_H_M L_F_M w_H t_IH t_EH t_IH_M t_EH_M real
z = [Z_H;psi_HH;psi_HF;psi_FF;psi_FH;N_H;N_F;L_H_M;L_F_M;w_H;t_IH;t_EH;t_IH_M;t_EH_M];
% COND1-10 are the constraint conditions
t_HF     = t_EH*t_IF;
t_FH     = t_EF*t_IH;
t_HF_M   = t_EH_M*t_IF_M;
t_FH_M   = t_EF_M*t_IH_M;
Phi_H_1 = T_H*w_H^(-theta);
Phi_F_1 = T_F*w_F^(-theta);
Phi_H_2 = T_H*w_H^(-theta)+T_F*(w_F*d_M*t_FH_M)^(-theta);
Phi_F_2 = T_F*w_F^(-theta)+T_H*(w_H*d_M*t_HF_M)^(-theta);
phi_H   = (Phi_H_2/Phi_H_1)^((sigma-1)*(1-alpha)/theta);
phi_F   = (Phi_F_2/Phi_F_1)^((sigma-1)*(1-alpha)/theta);
psi_HH_S = psi_HH * (f_S/f_D)^(1/(sigma-1))*(phi_H-1)^(1/(1-sigma));
psi_FF_S = psi_FF * (f_S/f_D)^(1/(sigma-1))*(phi_F-1)^(1/(1-sigma));
m_H_S = (psi_HH/psi_HH_S)^(k);
m_HF  = (psi_HH/psi_HF)^(k);
m_F_S = (psi_FF/psi_FF_S)^(k);
m_FH  = (psi_FF/psi_FH)^(k);
N_HF  = N_H*m_HF;
N_FH  = N_F*m_FH;
N_H_S = N_H*m_H_S;
N_F_S = N_F*m_F_S;
beta_HH = Phi_H_1/Phi_H_2;
beta_FH = 1- beta_HH;
beta_FF = Phi_F_1/Phi_F_2;
beta_HF = 1- beta_FF;
X_HH = sigma*lambda*N_H*w_H*(f_D+m_H_S*f_S);  
X_FF = sigma*lambda*N_F*w_F*(f_D+m_F_S*f_S);
X_HF = sigma*lambda*N_HF*w_H*f_X;
X_FH = sigma*lambda*N_FH*w_F*f_X;
X_HH_M = lambda*(1-alpha)*(sigma-1)*N_H*w_H*(f_D+beta_HH*m_HF*f_X+((beta_HH*phi_H-1)/(phi_H-1))*m_H_S*f_S);
X_FF_M = lambda*(1-alpha)*(sigma-1)*N_F*w_F*(f_D+beta_FF*m_FH*f_X+((beta_FF*phi_F-1)/(phi_F-1))*m_F_S*f_S);
X_FH_M = beta_FH*lambda*(1-alpha)*(sigma-1)*w_H*(((N_H_S*f_S)/(1-phi_H^(-1)))+N_HF*f_X);
X_HF_M = beta_HF*lambda*(1-alpha)*(sigma-1)*w_F*(((N_F_S*f_S)/(1-phi_F^(-1)))+N_FH*f_X);
%E_H and P_H
E_H = X_HH + t_FH*X_FH;
P_HH = (lambda*C*N_H )^(1/(1-sigma))*(w_H^(alpha)*Gamma^(1-alpha)*Phi_H_1^((alpha-1)/theta))*(psi_HH^(sigma-1)+(phi_H-1)*m_H_S*psi_HH_S^(sigma-1))^(1/(1-sigma));
P_FH = (lambda*C*N_FH)^(1/(1-sigma))*(d*t_FH*w_F^(alpha)*Gamma^(1-alpha)*Phi_F_2^((alpha-1)/theta)*psi_FH^(-1));
P_H  = (P_HH^(1-sigma)+P_FH^(1-sigma))^(1/(1-sigma));
% Constraint conditions
COND1 = (psi_HH/psi_FH)^(sigma-1) - (t_FH^(-sigma)/d^(sigma-1))*(f_D/f_X)*(w_H/w_F)^(alpha*(sigma-1)+1)*(Phi_F_2/Phi_H_1)^((1-alpha)*(sigma-1)/theta);
COND2 = (psi_FF/psi_HF)^(sigma-1) - (t_HF^(-sigma)/d^(sigma-1))*(f_D/f_X)*(w_F/w_H)^(alpha*(sigma-1)+1)*(Phi_H_2/Phi_F_1)^((1-alpha)*(sigma-1)/theta);
COND3 = (lambda-1)*psi_min^(k)*(f_D*psi_HH^(-k)+f_S*psi_HH_S^(-k)+f_X*psi_HF^(-k))-f_E;
COND4 = (lambda-1)*psi_min^(k)*(f_D*psi_FF^(-k)+f_S*psi_FF_S^(-k)+f_X*psi_FH^(-k))-f_E;
COND5 =  w_H*(L_H-L_H_M) - (1/sigma+alpha*(sigma-1)/sigma)*(X_HH+X_HF);
COND6 =  w_F*(L_F-L_F_M) - (1/sigma+alpha*(sigma-1)/sigma)*(X_FF+X_FH);
COND7 =  w_H*L_H_M - (X_HH_M + X_HF_M/t_HF_M); 
COND8 =  w_F*L_F_M - (X_FF_M + X_FH_M/t_FH_M);
COND9 = t_EH*X_HF + X_HF_M/t_IF_M - t_EF*X_FH - X_FH_M/t_IH_M;
COND10 = E_H/P_H;
CSTRT1 = [COND1; COND2;COND3; COND4;COND5; COND6;COND7; COND8;COND9];
CSTRT2 = Z_H - COND10;
ceq  = [CSTRT1; CSTRT2];
gradceq = jacobian(ceq,z).';
constraint1 = matlabFunction([],ceq,[],gradceq,'File','derivative_cons1','vars',{z});
%Initial value and boundary of search area
x_0 = [9.1347e+03,2.0408,    4.3421,    2.0408,    4.3421,    6.7266,    6.7266,  251.6883,  251.6883,    1.0000];
length_b = 2;
UB = [repmat(length_b.*x_0,4,1),(diag(ones(1,4)).*(length_b-1))+ones(4,4)];
LB = [repmat((-length_b).*x_0,4,1),(diag(ones(1,4)).*(-length_b-1))+ones(4,4)];
x_0(1,11:14) = 1;
x_0 = x_0.';
%% Optimset for 'fmincon'
tol = 1.0E-13;
options = optimset( ...
    'Display',         'off', ...
    'GradObj',         'on', ...
    'GradConstr',      'on', ...
    'DerivativeCheck', 'off', ...
    'FinDiffType',     'central', ...
    'TolFun',          tol, ...
    'TolX',            tol, ...
    'TolCon',          tol, ...
    'algorithm',       'active-set', ...
    'MaxFunEvals',     inf, ...
    'MaxIter',         5000);
%% Using 'fmincon' to solve the optimization problem 
result = zeros(length(x_0),4);
% Corresponding to the first 3 situations
 for j = 1:3
   optimal =  fmincon(@(x)sample_objective(x),x_0,...
[],[],[],[],LB(j,:), UB(j,:),@(x)constraint1(x),options);
   result(:,j) = optimal;
end
% Corresponding to the last situations 
result(:,4) =  fmincon(@(x)sample_funcTariffOptimal(x),x_0,...
[],[],[],[],LB(4,:), UB(4,:),@(x)constraint1(x),options);
with 
function [obj,grad] = sample_objective(x)
obj = -x(1,1);  
     if nargout>1
       grad=zeros(size(x)); 
        grad(1,1)=-1;
     end
end
It takes 5s, 14s and 6s to solve out the optimal solution for the first three situations. However, I run the code for solving the optimal given for one day and the result have not been obtained. The role of and are quite similar in the model (although not totally symmetric).
What's the possible reason behind the issue about the last situation? How can I increase the speed of the code? Any code for me to see how it goes when the code is running?
Your comments are welcome.

Best Answer

tol = 1.0E-13;
This is an incredibly tight tolerance - possibly at the limit of double float precision, depending on the order of magnitude of of your functions. It's conceviable that it would be hard to reach.
Also, I recommend using
>> dbstop if naninf
to see if your constraints are generating non-finite values.
Related Question