Forbid transformation of max(x,y) into MILP

Question

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*}

Answer 1

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(2*n,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, A*x <= b], c'*x)

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

Answer 2

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