Allow relaxation of ReLUQuadratic

src/predictors/ReLU.jl

@@ -300,7 +300,16 @@ ReLUQuadratic()
   └ moai_ReLU[2]*moai_z[2] = 0
 ```
 """
-struct ReLUQuadratic <: AbstractPredictor end
+struct ReLUQuadratic <: AbstractPredictor
+    relax_equality::Bool
+    relaxation_parameter::Float64
+    function ReLUQuadratic(;
+        relax_equality::Bool = false,
+        relaxation_parameter::Float64 = 0.0,
+    )
+        return new(relax_equality, relaxation_parameter)
+    end
+end
 
 function add_predictor(
     model::JuMP.AbstractModel,
@@ -314,6 +323,11 @@ function add_predictor(
     z = JuMP.@variable(model, [1:m], base_name = "moai_z")
     _set_bounds_if_finite.(Ref(cons), z, 0, max.(0, -first.(bounds)))
     append!(cons, JuMP.@constraint(model, x .== y - z))
-    append!(cons, JuMP.@constraint(model, y .* z .== 0))
+    if predictor.relax_equality
+        ϵ = predictor.relaxation_parameter
+        append!(cons, JuMP.@constraint(model, y .* z .<= ϵ))
+    else
+        append!(cons, JuMP.@constraint(model, y .* z .== 0))
+    end
     return y, Formulation(predictor, Any[y; z], cons)
 end

test/test_predictors.jl

@@ -237,6 +237,31 @@ function test_ReLU_Quadratic()
     return
 end
 
+function test_ReLU_Quadratic_relaxed()
+    model = Model(Ipopt.Optimizer)
+    set_silent(model)
+    @variable(model, x[1:2])
+    f = MathOptAI.ReLUQuadratic(;
+        relax_equality = true,
+        relaxation_parameter = 1e-4,
+    )
+    y, formulation = MathOptAI.add_predictor(model, f, x)
+    # Maximize sum of all variables to exercise the ReLU relaxation
+    @objective(model, Max, sum(formulation.variables))
+    @test length(y) == 2
+    @test num_variables(model) == 6
+    @test num_constraints(model, AffExpr, MOI.EqualTo{Float64}) == 2
+    @test num_constraints(model, QuadExpr, MOI.LessThan{Float64}) == 2
+    fix.(x, [-1, 2])
+    optimize!(model)
+    @assert is_solved_and_feasible(model)
+    # We do not satisfy equality to a tight tolerance
+    @test !isapprox(value.(y), [0.0, 2.0]; atol = 1e-6)
+    # But we satisfy equality to a loose tolerance
+    @test isapprox(value.(y), [0.0, 2.0]; atol = 1e-2)
+    return
+end
+
 function test_Sigmoid()
     model = Model(Ipopt.Optimizer)
     set_silent(model)