6

The function $\max(x,y)$ can be linearized by making use of additional binary variables. I suppose global optimisers are implemented to perform this transformation automatically.

Is there a global optimiser that allows this function to be treated as a non-linear function (I mean not transforming it into a MILP equivalent)?

Does e.g. SCIP, Octeract, BARON, ... have an option to forbid the transformation into an MILP?

The model looks like

\begin{align*} \max&\quad\text{linear function}\\\text{s.t.}&\quad z_i = \max(x_i,y_i),\quad i > 1 \\&\quad\text{linear block of constraints}\\&\quad x_i \ge 0, \hspace{0.1 cm} y_i \ge0, \hspace{0.1 cm} z_i \ge 0 \end{align*}

RobPratt
  • 32,006
  • 1
  • 44
  • 84
Clement
  • 2,252
  • 7
  • 13

2 Answers2

6

In YALMIP you can use arbitrary black-box operators to circumvent modelling

n = 5;
c = randn(3*n,1);
A = randn(10*n,3*n);
b = rand(10*n,1);

% MILP model
x = sdpvar(2n,1);
Domain= [0 <= x <= 1]; % Big-M should have some help
for i = 1:n
    x = [x;max([x(i),x(i+n)])];
end
optimize([Domain, Ax <= b], c'*x)


% Black-box model (sdpfun is renamed to blackbox in next version...)
x = sdpvar(2n,1);
Domain= [0 <= x <= 1];
for i = 1:n
    x = [x;sdpfun(x([i i+n]),@max)];
end
% ...with local nonlinear solver
optimize([Domain, Ax <= b], c'*x)


% ...with global solver
%    which appears to have issues as it doesn't understand the thin envelopes

%    (I will look into that)
%    However, using a blackbox inside the global solver is shaky and should not be 
%    done and I would expect it to fail. A dedicated operator should be added
optimize([Domain, A*x <= b], c'*x, sdpsettings('solver','bmibnb'))

The important question is why though? A global solver will likely have to work hard to find the logic which is nicely described by the binary structure.

EDIT: Appears the global solver thinks I am crazy sending a multivariate black-box function so it does not do any bound propagation etc. Changing to a univariate works better, and that can be done by $\max(x,y) = \frac{x+y+|x-y|}{2}$. Absolute values are MILP-represented so once again black-box needed. Still though, this is not something you want to do unless you have some very particular reason. In particular at the moment with a sampling-based blackbox to which you can have absolutely no trust (next version will allow you to append envelope generators to a black box which improves the situation)

x = sdpvar(2*n,1);
Domain = [0 <= x <= 1];
for i = 1:n        
    z = (x(i)+x(i+n)+sdpfun(x(i)-x(i+n),@abs))/2;
    x = [x;z];
end
optimize([Domain, A*x <= b], c'*x,sdpsettings('solver','bmibnb'))
Johan Löfberg
  • 1,702
  • 5
  • 9
  • Hi Johan! I made the exprerience with a model of my, that y exp(x), where y binary and x continuous, is performing better when it is left it as it is, in comparison to when it is transformed into an MINLP (exp(x) = z and reformulation of y z) . Of course chances are high, that this happened by chance. I would, however, be happy if I experience the same with the problem in question (even if it happens by chance = because the whole model is as it is). – Clement May 25 '21 at 13:18
  • I implemented stronger black-box support to test this now, and it performed rather well. Much faster than a MILP solution using Mosek, and slightly slower than a MILP solution via Gurobi (for $n=25$). If you are solving a MINLP involving $ye^{x}$, I don't see why you would want to introduce an intermediate variable $z$ and then introduce a nonconvex equality to the model. I am not surprised that $ye^{x}$ performs better. – Johan Löfberg May 25 '21 at 13:44
  • I am not sure I understand you. $ye^x$ is nonconvex itself. Why should it be obvious that $z = e^x$ is harder to cope with than $y e^x$? – Clement May 25 '21 at 16:23
  • You might have $ye^x$ entering your model in a convex way for fixed $y$ meaning all nodes when $y$ is fixed will be convex, while if you add an equality involving $e^x$ your nodes will always be nonconvex as the model then involves a nonlinear equality. – Johan Löfberg May 25 '21 at 17:00
  • There you are of course, right. – Clement May 25 '21 at 18:27
4

The default Octeract Engine behaviour is to handle this automatically without reformulating to an MILP.

Nikos Kazazakis
  • 12,121
  • 17
  • 59