Why do MINLP solvers have such a hard time with this simple problem?

I am trying to solve a global optimization problem, belonging to the category of MINLP (mixed-integer nonlinear programming).

The objective is to maximize a revenue, which depends on a switch (12 binary states). So far, so easy. But there is a constraint of the availability of a substrate (let's call it "volume"), which must never fall below zero, and should not exceed the storage limit. Thus, the volume condition can either be imposed as a nonlinear constraint function or introduced as a penalty term to the objective function, which is the better option as far as I have tested.

Consider this script (using either ga or MINLP solvers from OPTI toolbox ):

% Script to construct an optimization problem as MINLP
% 12 slots (i.e., 12 hours of a switch) which can be either 1 or 0,
% constrained by available volume
outputTime = (0:1:11)';
totalTime = 0:1:12;
initialVolume = 0.33;
% Discrete input array (initial guess, to be optimized)
inputVector = [0 0 1 1 1 0 0 0 1 1 0 0]';
revenuePerHour = [27; 30; 33; 38; 40; 37; 27; 29; 39; 43; 34; 25];
% bestHours = find(revenuePerHour > median(revenuePerHour));

%%Simulate model with initial guess
[totalRevenue, V] = objectiveFunction(inputVector, initialVolume, revenuePerHour);

%%Build optimization problem
numOfVars = length(outputTime);   %numberOfVariables
objFun = @(inputVector) objectiveFunction(inputVector, initialVolume, revenuePerHour);

% Inequality Constraint: Sum of hours between 4 and 8 (commented out with NOMAD)
A = ones(2, numOfVars);
b = [4, 8];
% Bounds (lb <= x <= ub), i.e., decision variables can only range between 0 and 1
lowerBound = zeros(1, numOfVars);
upperBound = ones(1, numOfVars);
% Binary Constraints
decisionVariableType = repmat('B', 1, numOfVars);
% Solve Problem: Optimal value is -230.02
optimizationRoutine = 'nomad';
% Distinguish several cases
switch optimizationRoutine
case 'ga'
    optionen = optimoptions('ga', 'Display', 'iter', 'FunctionTolerance', 1e-2);
    tic; [inputVectorOpt, fval] = ga(objFun, numOfVars, [], [], [], [],...
    lowerBound, upperBound, [], 1:numOfVars, optionen); toc
    inputVectorOpt = inputVectorOpt';
otherwise
    % Build OPTI Object
    optionen = optiset('solver', optimizationRoutine,...    % nomad, bonmin
        'display', 'iter', 'tolrfun', 1e-2, 'maxiter', 1e4,...
        'maxfeval', 4100, 'maxtime', 1800);
    Opt = opti('fun', objFun,... %'ineq', A, b,...
        'bounds', lowerBound, upperBound, 'xtype', decisionVariableType,...
        'x0', inputVector, 'options', optionen);
    tic; [inputVectorOpt, fval] = solve(Opt); toc
end

%%Rerun simulation with optimized values
[totalRevenue, V] = objectiveFunction(inputVectorOpt, initialVolume, revenuePerHour);

%%Plot results in Matlab
figure();
hold on; grid on;
stairs(outputTime, inputVector, 'DisplayName', 'Orignal schedule');
stairs(outputTime, revenuePerHour/100, 'DisplayName', 'Revenue per hour');
stairs(outputTime, inputVectorOpt, '-.', 'DisplayName', 'Adapted schedule');
plot(totalTime, [initialVolume; V], '-^', 'MarkerSize', 4, 'DisplayName', 'Volume 1');
legend('Location', 'Best');
xlim([0 12.5]);
ylim([-0.05 1.1]);
xlabel('Time (hours)');
set(gca, 'FontSize', 14);
hold off;

function [objFunValue, V, totalRevenueUnconstrained, penaltyTerm_low,...
penaltyTerm_up] = objectiveFunction(inputVector, initialValue_V, revenuePerHour)
%objectiveFunction (Non)linear objective function for (OPTI) optimization problem
%   Calculate revenue +/- penalty terms
    productionRate = 5.5/24;    % Constant value instead of vector (dynamic production)
  consumptionRate = 2*productionRate;
    upperThreshold_V = 1;

      % Check if inputVector is a column vector; if not, convert it to column vector
      % (necessary for Global Optimization Toolbox)
      if (size(inputVector, 1) == 1) && (size(inputVector, 2) > 1)
          inputVector = inputVector';
      end
      % Calculate volume change (rate) per hour
      dV = (productionRate - inputVector*consumptionRate);
      V = zeros(size(inputVector));
      V(1) = initialValue_V + dV(1);
      for hI = 2:numel(inputVector)
          V(hI) = V(hI - 1) + dV(hI);
      end
      % Total revenue + penalty term, with negative sign since its called by a minimizer
      totalRevenueUnconstrained = sum(inputVector.*revenuePerHour);
      % Negative and excessive volume is penalized (instead of given as constraint)
      % since practically infeasible, high weight
      penaltyTerm_low = 200*sum(abs(V) - V);
      % since practically feasible, low weight
      penaltyTerm_up = 28*sum(abs(V - upperThreshold_V) + (V - upperThreshold_V));
      objFunValue = -totalRevenueUnconstrained + penaltyTerm_low + penaltyTerm_up;
  end

The optimal solution would consist in the second 1s block to be prolonged by 1, and the first to be shifted forward, giving the optimal solution

inputVectorOpt = [0 0 0 1 1 1 0 0 1 1 1 0]';

with the objective function value -230.02.

Unfortunately, NOMAD needs very long to find the solution, usually taking 4096 function evaluations, which is equivalent to simulate all possible cases (2^12 = 4096). BONMIN does not even give a solution, stating:

Cbc0006I The LP relaxation is infeasible or too expensive

ga is the only one which finds the optimal solution in very short time.

Which leads me to the question: Why do two out of three solvers have such a hard time with this relatively simple problem? Am I missing something here? Do you have any suggestions?

[Xuehan Zhang]

I also face the same problem, could you let me know how you solved it? thank you so much

[winkmal]

Hey, not really, and it's been a while, as you can see.

I was, however, able to replace the for-loop by cumsum, which might make the calculation of the objective faster.

Also, and more importantly, I found a way to replace the nonlinear constraint by a matrix multiplication. Thus, it could be handled as a MILP, rather than MINLP, and solved easily with intlinprog.

The Matlab code is not published, but a very similar Julia code can be found here. It's actually a translation of my former Matlab code. Hope this helps.

[Xuehan Zhang]

thank you so much