MATLAB Answers

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

45 ビュー (過去 30 日間)
rotton
rotton 2018 年 8 月 15 日
編集済み: rotton 2018 年 8 月 15 日
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?

回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by