Issues about Using Function 'fmincon' solve Optimization Problem

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.

Answer 1

[Matt J]

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.

[sxj]

Thank you for the comement, Matt.

(1) I guess it won't be issue about the tolerance because it works for the other three situations, but I will try increasing the tol to see how it goes.

(2) You mean I should put the code ' dbstop if naninf' before the code about the last situation? I have restrict the search area of endogenous variables to be finite, so you mean it is the value of constraint achieve to be infinity?

The code is still running, and I will let it run until tomorrow morning to see what happens and fix it according to your suggestions.

[Matt J]

You mean I should put the code ' dbstop if naninf' before the code about the last situation?

You can run it at the command line

>> dbstop if naninf

and then run the optimization.

I have restrict the search area of endogenous variables to be finite, so you mean it is the value of constraint achieve to be infinity?

Yes. I'm thinking that fmincon may be struggling to find a region where your constraints are not NaN or Inf. That may be why the optimization is hanging, even though you have limited it to only 5000 iterations. Do you know how long a single iteration should take?

Another thing you could do is use the PlotFcn option of fmincon to plot the progress of the algorithm and see where it is getting stuck.

[sxj]

Hey, Matt, here are the progresses:

(1) I make the tolerance be smaller as 1.0e-8, but it still doesn't work, so I guess it is not the problem of tol.

(2) I tried the code 'dbstop if naninf' as you suggests, and the code stops in somewhere and the following information comes.

NaN/Inf breakpoint hit for fmincon.p on line 136.

136 'DiffMaxChange',Inf, ...

and MATLAB open the file of function 'fmincon' between line 131 and 141 as below

                defaultopt = struct( ...
    'Algorithm','interior-point', ...
    'AlwaysHonorConstraints','bounds', ...
    'DerivativeCheck','off', ...
    'Diagnostics','off', ...
    'DiffMaxChange',Inf, ...
    'DiffMinChange',0, ...
    'Display','final', ...
    'FinDiffRelStep', [], ...
    'FinDiffType','forward', ...
    'ProblemdefOptions', struct, ...

It seems to be the issue of 'DiffmaxChange' here. Does it means that there exists NaN or inf in somewhere as you suggest?

Besides, I also change the search area from [-2,2] to [-10,10] and also put 'PlotFcns', {@optimplotfval,@optimplotconstrviolation} and 'FunValCheck','on' in the optimset line. The following results and error messages come out:

Error using optimfcnchk/checkfun (line 244)

Supplied function '@(x)constraint1(x)' returned a complex value when evaluated;

FMINCON cannot continue.

Error in nlconst (line 755)

[nctmp,nceqtmp] = feval(confcn{3},x,varargin{:});

Error in fmincon (line 775)

nlconst(funfcn,X,l,u,full(A),B,full(Aeq),Beq,confcn,options,defaultopt, ...

Does it mean that it is because of the returning of complex value for constraints which make the fmincon stop continuing? (my objective function is simple as where is one of the endogenous variables, so it is less likely the complex comes from the objective) Any suggestions about fixing the complex value issue?

[Matt J]

Yes, fmincon requires both your constraint function and your objective function to return real, finite values at all points in the search space. The messages you are getting clearly show that your constraint function @(x)constraint1(x) is not returning real and finite values.

This is not entirely surprising, because I see lots of operations like z^(-k) in your constraint functions, with k=4.25. These operations will return complex values whenever z<0.

>> (-2)^(-4.25)

ans =

   0.0372 - 0.0372i

In any case, when dbstop stops your code, it should be a simple matter to go inside your constraint function and see what expression is generating the bad values.

[sxj]

Hey, Matt, thank you for the quick reply.

I also found other discussions mention that the complex issues in 'fmincon' usually come from the objective or constaint function which contain logs or non-integer powers of variables that take negative values. (as you mentioned above) I will restrict the search area to avoid the possible issue firstly and consider other possible solution like SQP algorithm as you mentioned in other questions.

I will let you know if there is any new progress.

[sxj]

Hey, Matt, it works well now after I restrict the search area for endogenous variables to be non-negative (the non-negative searching area still consistent with the reality of my model). No complex value issue comes out now. Thank you for the comments.

[sxj]

Hey, Matt, I found that if I change the search area by changing the value of 'b', the optimization result for the 4th situation changes, while the optimization result remain unchanged under different 'b'. I give details in my new question and here is the link for it. https://www.mathworks.com/matlabcentral/answers/471167-different-bounds-global-or-local-minimimum-in-fmincon

Your comments are appreciated.