Why is nonlinear inequality constraint not respected by fmincon?

4 次查看（过去 30 天）
winkmal
winkmal 2019-4-9
评论： winkmal 2019-4-10
The following script uses fmincon for minimization and works as expected:
% Feeding optimization over time, with algebraic model
% Parameters
substr_0 = 6.720;
substr_in = 4.800;
k_1 = 0.5;
f_1 = 43;
phi_1 = 0.5;
initialVolume = 47;
daysTotal = 7;
hoursTotal = daysTotal*24;
%% Consumer schedule per hour
onOff = [0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 ...
1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 ...
1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 ...
1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0]';
% Initial value: Feeding only during "Consumer ON" time
substr_in_time = substr_in/12.*onOff;
%% Setup optimization problem
parameters = [substr_0, k_1, f_1, phi_1, initialVolume];
% Objective function: Standard deviation of storage volume
objFun = @(feedingVector) calculateStorageRes(feedingVector, onOff, parameters);
nlcon = @(feedingVector) nonLinearConstraintFcn(feedingVector, onOff, parameters);
% Inequality Constraint: Limit sum of fed substrate
A = ones(1, hoursTotal);
b = 1.00*sum(daysTotal*substr_in);
% Bounds (lb <= x <= ub), i.e., decision variables can only range between 0 and 6.3
lowerBound = zeros(1, hoursTotal);
upperBound = 0.95*substr_in.*ones(1, hoursTotal);
% Possible solvers: FilterSD, Ipopt, fmincon (local)
optimOpts = optimoptions('fmincon', 'Display', 'iter', 'FunctionTolerance', 1e-2,...
'MaxIterations', 1e4*daysTotal, 'MaxFunctionEvaluations', 1e5*daysTotal);
% 'HonorBounds', true, 'Algorithm', 'sqp'
%% Solve optimization problem
tic; [feeding_Opt, fval] = fmincon(objFun, substr_in_time, A, b, [], [],...
lowerBound, upperBound, [], optimOpts); toc
%% Calculate model with optimal values
feeding_Opt(feeding_Opt < 1e-15) = 0; % Overwrite very small values with 0
[~, interm1, storageVolume] = substr_Model(feeding_Opt, onOff, parameters);
%% Plot results
figure();
plot(1:hoursTotal, interm1, '-o', 'LineWidth', 1.5, 'MarkerSize', 4, ...
'DisplayName', 'Interm production');
hold on; grid on
stairs(0:hoursTotal-1, onOff, 'LineWidth', 1.5, 'DisplayName', 'Consumption schedule');
plot(1:hoursTotal, storageVolume/10, '-^', 'LineWidth', 1.5, 'MarkerSize', 4,...
'DisplayName', 'Storage volume [1e-1]');
stairs(0:hoursTotal-1, feeding_Opt, 'LineWidth', 1.5, 'DisplayName', 'Feeding schedule');
xlim([0 hoursTotal]);
xtickVector = 0:12:hoursTotal;
xticks(xtickVector);
xticklabels(xtickVector/24);
xlabel('Time (days)');
legend('Location', 'Best');
set(gca, 'FontSize', 18);
hold off
function storageRes = calculateStorageRes(substr_feeding, onOff, parameters)
hoursTotal = length(substr_feeding);
[~, ~, storageVolume] = substr_Model(substr_feeding, onOff, parameters);
% Formula of standard deviation
storageRes = sqrt(1/(hoursTotal-1)*sum((storageVolume - mean(storageVolume)).^2));
end
function [storageConstraints, ceq] = nonLinearConstraintFcn(substr_feeding, onOff, parameters)
%nonLinearConstraintFcn Nonlinear constraint function for optimization problem
% --> Could this be linearized?
storageThreshold_min = 1.5;
storageThreshold_max = 105.5;
[~, ~, storageVolume] = substr_Model(substr_feeding, onOff, parameters);
minVolume = min(storageVolume);
maxVolume = max(storageVolume);
lowerStorageConstraint = storageThreshold_min - minVolume;
upperStorageConstraint = maxVolume - storageThreshold_max;
storageConstraints = [lowerStorageConstraint, upperStorageConstraint];
ceq = [];
end
function [substr, interm1, storageVolume] = substr_Model(substr_feeding, onOff, parameters)
hoursTotal = length(substr_feeding);
substr_0 = parameters(1);
k_1 = parameters(2);
f_1 = parameters(3);
phi_1 = parameters(4);
initialVolume = parameters(5);
% Simple for-loop with feeding as vector per hour
substr = zeros(hoursTotal, 1);
substr(1) = (substr_0 + substr_feeding(1))*exp(-k_1/24);
for hour = 2:hoursTotal
substr(hour) = (substr(hour-1) + substr_feeding(hour))*exp(-k_1/24);
end
% Calculate resulting interm1 production and storage volume
interm1 = phi_1*f_1*substr/24;
interm1ConsumptionRate = 17.700;
storageVolume = initialVolume + cumsum(interm1 - onOff*interm1ConsumptionRate);
end
But, as soon as I tell fmincon to use nlcon,
tic; [feeding_Opt, fval] = fmincon(objFun, substr_in_time, A, b, [], [],...
lowerBound, upperBound, nlcon, optimOpts); toc
results become weird, and though I get convergence, conditions are not met:
>>[max(storageVolume), min(storageVolume)]
ans =
108.2434 -0.3176
>>nlcon(feeding_Opt)
ans =
1.8176 2.7434
(Should both be below 0)
Am I doing something wrong? How could I reformulate my code for fmincon to find the desired minimum?
  2 个评论
Torsten
Torsten 2019-4-9
编辑：Torsten 2019-4-9
You shouldn't work with min(storageVolume) and max(storageVolume), but constrain
storageThreshold_min <= storageVolume(hour) <= storageThreshold_max
for 1 <= hour <= hoursTotal.
Note that min() and max() operations generate non-differentiable constraint functions, but fmincon requires differentiabiliy.
winkmal
winkmal 2019-4-10
So you mean, instead of...
minVolume = min(storageVolume);
maxVolume = max(storageVolume);
lowerStorageConstraint = storageThreshold_min - minVolume;
upperStorageConstraint = maxVolume - storageThreshold_max;
storageConstraints = [lowerStorageConstraint, upperStorageConstraint];
I should use
lowerStorageConstraint = storageThreshold_min - storageVolume;
upperStorageConstraint = storageVolume - storageThreshold_max;
storageConstraints = [lowerStorageConstraint, upperStorageConstraint];
I tried, but optimum is still not within constraints... If I increase storageThreshold_max from 105.5 to 115.5, then it works. Does this mean that the problem is simply infeasible with the smaller constraint?

请先登录，再进行评论。

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

回答（0 个）

类别

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

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by