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

26 views (last 30 days)
winkmal on 15 Aug 2018
Moved: Matt J on 29 Oct 2022
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
% 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 on 29 Oct 2022
Moved: Matt J on 29 Oct 2022
thank you so much

Categories

Find more on Linear Programming and Mixed-Integer Linear Programming in Help Center and File Exchange

R2017b

Community Treasure Hunt

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

Start Hunting!