Why is nonlinear inequality constraint not respected by fmincon?

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?

[Torsten]

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]

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?