Issues about Using Function 'fmincon' solve Optimization Problem

1 次查看（过去 30 天）
sxj
sxj 2019-7-5
评论： sxj 2019-7-11
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.

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2019-7-5
编辑：Matt J 2019-7-5
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.
  7 个评论
显示 5更早的评论
sxj
sxj 2019-7-8
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
sxj 2019-7-11
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Optimization Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by