Option to treat equality constraints in domain differently

[tweisser]

By default SumOfSquares handles equality constraints in the domain by working over the quotient ring R[x]/I where I is the ideal generated by the equality constraints.
This can lead to unexpected behaviour when dealing with noisy data, e.g., when the equality constraints are the result of some preprocessing step which involved floating point computations. A common source for such rounding errors could emmerge from normalizing data.

using DynamicPolynomials
using SemialgebraicSets

@polyvar vr[1:3]
@polyvar vi[1:3]
@polyvar pg[1:2]
@polyvar qg[1:2]
@polyvar pd[1]

p = [
vr[2]^2 + vi[2]^2 - 1.098778957631452,
 vr[3]^2 - 1.0,
 0.08879023307436182*vr[2]^2 - 0.04439511653718091*vr[2]*vr[3] - 0.04439511653718091*vr[2]*vr[1] - 1.3318534961154274*vr[2]*vi[1] + 1.3318534961154274*vr[3]*vi[2] + 1.3318534961154274*vr[1]*vi[2] + 0.08879023307436182*vi[2]^2 - 0.04439511653718091*vi[2]*vi[1] - 0.4117054263835753 ,
 1.9637069922308548*vr[2]^2 - 1.3318534961154274*vr[2]*vr[3] - 1.3318534961154274*vr[2]*vr[1] + 0.04439511653718091*vr[2]*vi[1] - 0.04439511653718091*vr[3]*vi[2] - 0.04439511653718091*vr[1]*vi[2] + 1.9637069922308548*vi[2]^2 - 1.3318534961154274*vi[2]*vi[1] - qg[1] + 0.5,
 -0.04439511653718091*vr[2]*vr[3] + 0.08879023307436182*vr[3]^2 - 0.04439511653718091*vr[3]*vr[1] - 1.3318534961154274*vr[3]*vi[2] - 1.3318534961154274*vr[3]*vi[1] - pg[2] + 1.0,
 -1.3318534961154274*vr[2]*vr[3] + 1.9637069922308548*vr[3]^2 - 1.3318534961154274*vr[3]*vr[1] + 0.04439511653718091*vr[3]*vi[2] + 0.04439511653718091*vr[3]*vi[1] - qg[2] + 0.5,
 -0.04439511653718091*vr[2]*vr[1] + 1.3318534961154274*vr[2]*vi[1] - 0.04439511653718091*vr[3]*vr[1] + 1.3318534961154274*vr[3]*vi[1] + 0.08879023307436182*vr[1]^2 - 1.3318534961154274*vr[1]*vi[2] - 0.04439511653718091*vi[2]*vi[1] + 0.08879023307436182*vi[1]^2 + pd[1],
 -1.3318534961154274*vr[2]*vr[1] - 0.04439511653718091*vr[2]*vi[1] - 1.3318534961154274*vr[3]*vr[1] - 0.04439511653718091*vr[3]*vi[1] + 1.9637069922308548*vr[1]^2 + 0.04439511653718091*vr[1]*vi[2] - 1.3318534961154274*vi[2]*vi[1] + 1.9637069922308548*vi[1]^2 + pd[1] - 0.5]

I = ideal(p)

Symbolically the ideal is empty:

SemialgebraicSets.computegröbnerbasis!(I)
println(I.p)
# Polynomial{true,Float64}[1.0]

However, if we allow for some numerical error the following point can be seen as a solution:

sub = Dict(
vr[1] =>  0.955906,
vr[2] =>  1.04674,
vr[3] => 1.0,
vi[1] => -0.395553,
vi[2] => -0.0558258,
vi[3] => 1.88079e-37,
pg[1] => 1.41171,
pg[2] =>1.60105,
qg[1] =>-0.111999,
qg[2] =>-0.223562,
pd[1] =>1.00000
)

subs.(p, sub...)

If we minimize pd[1] over the ideal intersected with (pd[1]-0.5)*(1.5-pd[1]) >=0 we could expect that the mimum is between 0.5 and 1. However, the problem is infeasible.

using SumOfSquares
using MosekTools

m = SOSModel(Mosek.Optimizer)
@variable m t
@objective m Max t
@constraint(m,
    pd[1]-t >=0,
    domain = basicsemialgebraicset(algebraicset(p), [(pd[1]-0.5)*(1.5-pd[1])]),
    maxdegree = 4
            )

optimize!(m)
termination_status(m)
# DUAL_INFEASIBLE::TerminationStatusCode = 3

This problem could be solved by introducing an option to encode equality constraints differently. A common way is to introduce polynomial multipliers for the equality (such as SOSPoly multipliers are introduced for inequalities).

vars = sort!(collect(keys(sub)); rev = true)
m = SOSModel(Mosek.Optimizer)
@variable m t
@objective m Max t
lhs = let 
    lhs = pd[1] - t
    for poly in p 
        pi = @variable m variable_type =  Poly(monomials(vars, 0:4-maxdegree(poly)))
    lhs -= pi*poly
    end
    lhs 
end
@constraint(m,
    lhs >=0,
    domain = @set((pd[1]-0.5)*(1.5-pd[1])>=0),
    maxdegree = 4
            )

optimize!(m)
termination_status(m)
# OPTIMAL::TerminationStatusCode = 1

Now the problem solves with optimal and actually returns an expected solutions.