Why is nonlinear inequality constraint not respected by fmincon?

3 views (last 30 days)
winkmal on 9 Apr 2019
Commented: winkmal on 10 Apr 2019
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?
winkmal on 10 Apr 2019
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?

Categories

Find more on Get Started with Optimization Toolbox in Help Center and File Exchange

R2018a

Community Treasure Hunt

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

Start Hunting!