Numerical issues

[blegat]

function test(ztol=Base.rtoldefault(Float64))
    X = [0.5459627556242905, 1.7288950489507429, 0.7167681447476535]
    Y = [-54.06002080721971, 173.77393714162503, 71.48154370522498]
    n = 3
    @polyvar W[1:n] α β
    I = @set(sum([W[i] * X[i] * (β * X[i] + α - Y[i]) for i = 1:n]) == 0, library = Buchberger(ztol))
    for i in 1:n
        addequality!(I, W[i] - W[i]^2)
    end
    SemialgebraicSets.computegröbnerbasis!(I.I)
    I
end

The Groebner basis found does not seem to be correct, it gives incorrect results with
jump-dev/SumOfSquares.jl#269
as it does no match the lagrangian version

using SumOfSquares
using DynamicPolynomials
using MosekTools
using LinearAlgebra
using SCS
using JuMP

function test(solver, d, obj, I, lagrangian_monos)
    model = SOSModel(solver)
    @variable(model, t)
    if lagrangian_monos isa Float64
		display(I)
        c = @constraint(model, t >= obj, domain = I, maxdegree = d)
		@show moi_set(constraint_object(c)).certificate
		@show moi_set(constraint_object(c)).domain
		@show I.I.gröbnerbasis
		display(I)
    else
        p = equalities(I)
        @variable(model, multipliers[eachindex(p)], Poly(lagrangian_monos))
        @constraint(model, t >= obj + multipliers ⋅ p)
    end
    @objective(model, Min, t)
    optimize!(model)
    solution_summary(model)
end

function sol(solver, d, ztol)
    X = [0.5459627556242905, 1.7288950489507429, 0.7167681447476535]

    Y = [-54.06002080721971, 173.77393714162503, 71.48154370522498]

    n = 3

    @polyvar W[1:n] α β

	I = @set(sum([W[i] * X[i] * (β * X[i] + α - Y[i]) for i = 1:n]) == 0, library = Buchberger(ztol))
    #I = @set sum([W[i] * X[i] * (β * X[i] + α - Y[i]) for i = 1:n]) == 0
	for i in 1:n
		addequality!(I, W[i] - W[i]^2)
	end

    obj = sum(W)

	if ztol === nothing
		monos = monomials([α, β], 0:d)
	else
		monos = ztol
	end

	test(solver, d, obj, I, monos)
end

Probably need take into account Section 10.1.3 of Stetter's "Numerical Polynomial Algebra".