Bilinear terms
Adapted from: (Floudas et al., 1999; Section 3.2) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 3.2).
using DynamicPolynomials
@polyvar x[1:8]
p = sum(x[1:3])
using SumOfSquares
K = @set 0.0025 * (x[4] + x[6]) <= 1 &&
0.0025 * (-x[4] + x[5] + x[7]) <= 1 &&
0.01 * (-x[5] + x[8]) <= 1 &&
100x[1] - x[1] * x[6] + 8333.33252x[4] <= 250000/3 &&
x[2] * x[4] - x[2] * x[7] - 1250x[4] + 1250x[5] <= 0 &&
x[3] * x[5] - x[3] * x[8] - 2500x[5] + 1250000 <= 0 &&
100 <= x[1] && x[1] <= 10000 &&
1000 <= x[2] && x[2] <= 10000 &&
1000 <= x[3] && x[3] <= 10000 &&
10 <= x[4] && x[4] <= 1000 &&
10 <= x[5] && x[5] <= 1000 &&
10 <= x[6] && x[6] <= 1000 &&
10 <= x[7] && x[7] <= 1000 &&
10 <= x[8] && x[8] <= 1000Basic semialgebraic Set defined by no equality
22 inequalities
1.0 - 0.0025*x[6] - 0.0025*x[4] ≥ 0
1.0 - 0.0025*x[7] - 0.0025*x[5] + 0.0025*x[4] ≥ 0
1.0 - 0.01*x[8] + 0.01*x[5] ≥ 0
83333.33333333333 - 8333.33252*x[4] - 100.0*x[1] + x[1]*x[6] ≥ 0
-1250.0*x[5] + 1250.0*x[4] + x[2]*x[7] - x[2]*x[4] ≥ 0
-1.25e6 + 2500.0*x[5] + x[3]*x[8] - x[3]*x[5] ≥ 0
-100.0 + x[1] ≥ 0
10000.0 - x[1] ≥ 0
-1000.0 + x[2] ≥ 0
10000.0 - x[2] ≥ 0
-1000.0 + x[3] ≥ 0
10000.0 - x[3] ≥ 0
-10.0 + x[4] ≥ 0
1000.0 - x[4] ≥ 0
-10.0 + x[5] ≥ 0
1000.0 - x[5] ≥ 0
-10.0 + x[6] ≥ 0
1000.0 - x[6] ≥ 0
-10.0 + x[7] ≥ 0
1000.0 - x[7] ≥ 0
-10.0 + x[8] ≥ 0
1000.0 - x[8] ≥ 0
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of 2100
model2 = solve(2)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 21
constraints = 29
nnz(P) = 0
nnz(A) = 64
cones (total) = 2
: Zero = 1, numel = 9
: Nonnegative = 1, numel = 20
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 -2.3512e+03 -2.3512e+03 1.93e-16 2.36e-03 2.84e-02 1.00e+00 9.80e+01 ------
1 -2.1101e+03 -2.1100e+03 9.17e-06 4.66e-05 5.21e-04 3.70e-02 1.86e+00 9.81e-01
2 -2.1001e+03 -2.1001e+03 9.66e-08 4.80e-07 5.34e-06 3.84e-04 1.93e-02 9.90e-01
3 -2.1000e+03 -2.1000e+03 9.66e-10 4.80e-09 5.34e-08 3.84e-06 1.93e-04 9.90e-01
4 -2.1000e+03 -2.1000e+03 9.66e-12 4.80e-11 5.34e-10 3.84e-08 1.93e-06 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 583μs
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : 2.10000e+03
Dual objective value : 2.10000e+03
* Work counters
Solve time (sec) : 5.83163e-04
Barrier iterations : 4This page was generated using Literate.jl.