Nonconvex quadratically constrained quadratic programs
Adapted from: (Hesse, 1973), (Floudas et al., 1999; Section 3.4), (Laurent, 2008; Example 6.22) and (Lasserre, 2009; Table 5.1)
We consider the nonconvex Quadratically Constrained Quadratic Programs (QCQP) introduced in [H73]. Consider now the polynomial optimization problem (Laurent, 2008; Example 6.22) of maximizing the convex quadratic function (hence nonconvex since convex programs should either maximize concave functions or minimize convex functions) $25(x_1 - 2)^2 + (x_2 - 2)^2 + (x_3 - 1)^2 + (x_4 - 4)^2 + (x_5 - 1)^2 + (x_6 - 4)^2$ over the basic semialgebraic set defined by the nonconvex quadratic inequalities $(x_3 - 3)^2 + x_4 \ge 4$, $(x_5 - 3)^2 + x_6 \ge 4$, and linear inequalities $x_1 - 3x_2 \le 2$, $-x_1 + x_2 \le 2$, $2 \le x_1 + x_2 \le 6$, $0 \le x_1, x_2$, $1 \le x_3 \le 5$, $0 \le x_4 \le 6$, $1 \le x_5 \le 5$, $0 \le x_6 \le 10$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$,
using DynamicPolynomials
@polyvar x[1:6]
centers = [2, 2, 1, 4, 1, 4]
weights = [25, 1, 1, 1, 1, 1]
p = -weights' * (x .- centers).^2
using SumOfSquares
K = @set x[1] >= 0 && x[2] >= 0 &&
x[3] >= 1 && x[3] <= 5 &&
x[4] >= 0 && x[4] <= 6 &&
x[5] >= 1 && x[5] <= 5 &&
x[6] >= 0 && x[6] <= 10 &&
(x[3] - 3)^2 + x[4] >= 4 &&
(x[5] - 3)^2 + x[6] >= 4 &&
x[1] - 3x[2] <= 2 &&
-x[1] + x[2] <= 2 &&
x[1] + x[2] <= 6 &&
x[1] + x[2] >= 2Basic semialgebraic Set defined by no equality
16 inequalities
x[1] ≥ 0
x[2] ≥ 0
-1 + x[3] ≥ 0
5 - x[3] ≥ 0
x[4] ≥ 0
6 - x[4] ≥ 0
-1 + x[5] ≥ 0
5 - x[5] ≥ 0
x[6] ≥ 0
10 - x[6] ≥ 0
5 + x[4] - 6*x[3] + x[3]^2 ≥ 0
5 + x[6] - 6*x[5] + x[5]^2 ≥ 0
2 + 3*x[2] - x[1] ≥ 0
2 - x[2] + x[1] ≥ 0
6 - x[2] - x[1] ≥ 0
-2 + x[2] + x[1] ≥ 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 cannot find any lower bound.
model2 = solve(2)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 45
constraints = 72
nnz(P) = 0
nnz(A) = 109
cones (total) = 3
: Zero = 1, numel = 28
: Nonnegative = 1, numel = 16
: PSDTriangle = 1, numel = 28
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 7.6836e+01 7.6836e+01 3.70e-16 3.89e-01 1.31e-01 1.00e+00 1.04e+02 ------
1 2.3712e+04 2.9567e+04 2.47e-01 2.41e-01 3.59e-02 5.92e+03 6.57e+01 9.27e-01
2 -2.2453e+02 9.6600e+02 5.30e+00 2.61e-01 5.17e-02 1.19e+03 5.27e+01 5.22e-01
3 1.9334e+03 7.3247e+03 2.79e+00 6.41e-02 7.49e-03 5.39e+03 1.13e+01 8.20e-01
4 1.7396e+04 3.1882e+04 8.33e-01 6.95e-03 7.47e-04 1.45e+04 2.37e+00 9.55e-01
5 1.2768e+05 2.1916e+05 7.16e-01 9.16e-04 9.60e-05 9.15e+04 3.33e-01 8.92e-01
6 5.0101e+06 8.4176e+06 6.80e-01 2.18e-05 2.27e-06 3.41e+06 8.11e-03 9.88e-01
7 4.3548e+08 7.3167e+08 6.80e-01 2.54e-07 2.63e-08 2.96e+08 9.44e-05 9.88e-01
8 2.2618e+10 3.8076e+10 6.83e-01 5.07e-09 5.26e-10 1.55e+10 1.88e-06 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = primal infeasible
solve time = 1.78ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : INFEASIBLE
Message from the solver:
"PRIMAL_INFEASIBLE"
* Candidate solution (result #1)
Primal status : INFEASIBLE_POINT
Dual status : INFEASIBILITY_CERTIFICATE
Objective value : -2.26177e+10
Dual objective value : -3.80764e+10
* Work counters
Solve time (sec) : 1.77532e-03
Barrier iterations : 8The second level of the hierarchy finds the lower bound of -310.
model3 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 423
constraints = 506
nnz(P) = 0
nnz(A) = 1243
cones (total) = 17
: Zero = 1, numel = 84
: Nonnegative = 1, numel = 2
: PSDTriangle = 15, numel = (28,28,28,28,...,28)
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 1.3550e+02 1.3550e+02 2.10e-16 1.24e-01 2.31e-01 1.00e+00 1.28e+01 ------
1 1.3966e+02 1.4303e+02 2.41e-02 6.58e-02 1.20e-01 3.94e+00 4.70e+00 7.91e-01
2 1.6240e+02 1.7566e+02 8.17e-02 4.50e-02 5.52e-02 1.38e+01 1.71e+00 8.78e-01
3 1.5720e+02 1.6295e+02 3.65e-02 1.79e-02 1.95e-02 5.95e+00 6.72e-01 6.85e-01
4 1.4675e+02 1.5215e+02 3.68e-02 1.25e-02 1.24e-02 5.53e+00 4.15e-01 6.37e-01
5 1.8037e+02 1.8247e+02 1.17e-02 3.28e-03 2.76e-03 2.14e+00 9.73e-02 7.90e-01
6 2.1477e+02 2.1622e+02 6.73e-03 1.98e-03 1.53e-03 1.47e+00 5.73e-02 5.20e-01
7 2.5401e+02 2.5440e+02 1.55e-03 4.18e-04 2.80e-04 3.99e-01 1.15e-02 8.97e-01
8 2.7007e+02 2.7034e+02 9.95e-04 1.66e-04 8.76e-05 2.71e-01 3.47e-03 8.17e-01
9 2.8259e+02 2.8277e+02 6.14e-04 7.86e-05 3.41e-05 1.75e-01 1.31e-03 8.24e-01
10 2.8775e+02 2.8799e+02 8.37e-04 7.38e-05 1.95e-05 2.42e-01 8.16e-04 6.61e-01
11 3.0102e+02 3.0109e+02 2.34e-04 2.35e-05 5.75e-06 7.09e-02 2.48e-04 9.90e-01
12 3.0953e+02 3.0953e+02 1.29e-05 1.25e-06 2.67e-07 4.01e-03 1.22e-05 9.54e-01
13 3.0999e+02 3.0999e+02 2.13e-07 2.07e-08 4.41e-09 6.64e-05 2.02e-07 9.83e-01
14 3.1000e+02 3.1000e+02 5.77e-09 5.60e-10 1.19e-10 1.80e-06 5.46e-09 9.73e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 38.6ms
* 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 : -3.10000e+02
Dual objective value : -3.10000e+02
* Work counters
Solve time (sec) : 3.86480e-02
Barrier iterations : 14This page was generated using Literate.jl.