HiGHS solver >100x slower than Gurobi with bilinear big-m relaxation

I’m not sure where to put this as it’s not specifically an issue, and I’m not sure who would improve it.
For my research I need to solve a bilinear model and I’d like to use an open-source solver ideally.

Rough example of my constraints

My constraint is like xy==0. I relax this using big_M to roughly

@variable(model, 0 <= x[1:n])
@variable(model, y[1:n] <= 0)
@variable(model, 0 <= u[1:n] <= 1, binary=true)
@constraint(model, big_m_x, x .≤ 10 .* (1 .- u))
@constraint(model, big_m_y,  -10 .* u .≤ y)

I’ve formulated a MWE example (see below), in this scenario Gurobi solves the problem in 0.4s while HiGHS takes 40-60 seconds which is >100 times slower. As my actual problem is much bigger than the example below HiGHs become prohibitively slow.

I appreciate Gurobi is often quicker than HiGHS, but not usually by this much, it seems like my HiGHs implementation is missing a trick?

MWE

Note I haven’t gone through this MWE thoroughly, it just serves to demonstrate the timing scenario and is roughly similar to my actual problem

using JuMP, HiGHS, Gurobi

function create_model(n, nl, solver)
    if solver==:H
        model = Model(HiGHS.Optimizer)
    elseif solver==:G
        model = Model(Gurobi.Optimizer)
    end
    if nl
        @variable(model, 0 <= x[1:n] <= 10)
        @variable(model, -10 <= y[1:n] <= 0)
        @constraint(model, [i=1:n], x[i]*y[i] ==0)

    else
        @variable(model, 0 <= x[1:n])
        @variable(model, y[1:n] <= 0)
        @variable(model, 0 <= u[1:n] <= 1, binary=true)
        @constraint(model, big_m_x, x .≤ 10 .* (1 .- u))
        @constraint(model, big_m_y,  -10 .* u .≤ y)

    end

    @objective(model, Min, sum(x[i] + -2*y[i] for i in 1:n))
    @constraint(model, [i=1:n], x[i] - y[i] >= 1)
    @constraint(model, [i=1:n], x[i] - y[i] <= 15)

    return model
end

n = 100_000
model = create_model(n,false,:G)

optimize!(model)

# Gurobi = 0.4
# HiGHS = 58.6
solve_time(model)
objective_value(model)
sum(value.(model[:x]))
sum(value.(model[:y]))

with

(jl_MPgoHK) pkg> st
Status `/private/var/folders/63/lgx7112n3b58qj9r923y2v680000gp/T/jl_MPgoHK/Project.toml`
  [2e9cd046] Gurobi v1.2.3
  [87dc4568] HiGHS v1.9.0
  [4076af6c] JuMP v1.22.1

This is probably expected behavior.

You should expect that HiGHS is 10-100× slower than Gurobi.

It will be missing a trick, but because Gurobi is closed source, we don’t know which one(s).

I see, usually I’ve found it to be more like 10x, so found this up to 150x difference noteworthy. But as you say, we cannot find what is causing this disparity.

Unfortunately, I do not have a good enough understanding of HiGHS to try investigating beyond changing parameters, which has made little difference.