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

5 次查看（过去 30 天）
winkmal
winkmal 2018-8-15
移动：Matt J 2022-10-29
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?
  3 个评论
显示 1更早的评论
winkmal
winkmal 2022-10-27
移动：Matt J 2022-10-29
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.

请先登录，再进行评论。

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

回答（0 个）

类别

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

标签

产品


版本

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by