Using the Positivstellensatz, you know that the set {(x, y) | f1(x, y) <= 0, f2(x, y) <= 0, x = a}
is empty if and only if there exists f
in the preorder generated by -f1, -f2
and h
in the ideal generated by x - a
such that f + h = -1
.
You can formulate it as a Sum-of-Squares program with
s_1 * f1 + s_2 * f2 + s_3 * f1 * f2 + p * (x - a) = -1
where s_i
are SOS polynomials and p
is any polynomial.
To solve it numerically, you will need to fix the maximal degree of the expression but you can build a hierarchy of SOS programs with higher and higher maximal degree.
An issue with this is that if the coefficients of p
are JuMP decision variables and a
is also a JuMP decision variable, then p * (x - a)
will be bilinear in the JuMP decision variables. Since you have SDP constraints (to constraint the s_i
to be SOS) and bilinear constraint, it becomes tricky to solve. You will need to mix an SDP solver with a meta-solver like GitHub - lanl-ansi/Alpine.jl: A Julia/JuMP-based Global Optimization Solver for Non-convex Programs to handle bilinearity.
However, if you fix a
, e.g. you are commuting g(a)
and use a separate method that find the root of g
given an oracle for g(a)
then you don’t have any bilinearity issues.
Here is how it would look like. First define your polynomials:
using DynamicPolynomials
@polyvar x y
f1 = -(1 - x)* y + 1/2* (1 - x - x^2 - x^3 - 2* (1 - x)* y)
f2 = (1 - 2*x)*(1 + x + x^2) - 2*(-1 + x)^2*(1 + x)*y
You want s * p
to have degree d
, hence s
to have degree d - deg(p)
. As s
is SOS, it has even degree. Taking all this into account, you get something like this to create s
:
function s(model, p, d)
deg = div(d - maxdegree(p, y) + 1, 2)
monos = monomials([y], 0:deg)
return @variable(model, variable_type=SOSPoly(monos))
end
Now you can combine all this to create your oracle for g
.
function g(a, d, solver = CSDP.Optimizer)
f_1 = subs(f1, x => a)
f_2 = subs(f2, x => a)
model = Model(solver)
set_silent(model)
p = sum(s(model, f, d) * f for f in [f_1, f_2, f_1 * f_2])
@constraint(model, p == -1)
optimize!(model)
status = termination_status(model)
if status == MOI.OPTIMAL
return false
elseif status == MOI.INFEASIBLE
return true
else
error("Termination status: $status.")
end
end
To test it, suppose you want to evaluate g
at the given values of x
: x_values = range(0, stop=1, length=10)
.
With degree d
, you get g.(x_values, 2)
which gives true
for all values except for x = 1
.
Climbing up the hierarchy, with d = 6
: g.(x_values, 6)
, you get the same except that you also get false
for an additional values before x = 1
.
With d = 10
, you get false
for the last three values of x
.
Accepting the status MOI.ALMOST_OPTIMAL
as false
, with d = 20
, you get false
for the last 4 values of g
.
EDIT: I just realized that in the implementation, I forgot we need to take -f1
and -f2
insted of f1
and f2
as mentioned in the beginning.