|
| 1 | +# # Nonconvex QP |
| 2 | + |
| 3 | +#md # [](@__BINDER_ROOT_URL__/generated/Polynomial Optimization/qp.ipynb) |
| 4 | +#md # [](@__NBVIEWER_ROOT_URL__/generated/Polynomial Optimization/qp.ipynb) |
| 5 | +# **Adapted from**: Section 2.2 of [F99], Table 5.1 of [Las09] |
| 6 | +# |
| 7 | +# [F99] Floudas, Christodoulos A. et al. |
| 8 | +# *Handbook of Test Problems in Local and Global Optimization.* |
| 9 | +# Nonconvex Optimization and Its Applications (NOIA, volume 33). |
| 10 | +# |
| 11 | +# [Las09] Lasserre, J. B. |
| 12 | +# *Moments, positive polynomials and their applications* |
| 13 | +# World Scientific, **2009**. |
| 14 | + |
| 15 | +# ## Introduction |
| 16 | + |
| 17 | +# Consider the nonconvex Quadratic Program (QP) [F99, Section 2.2] |
| 18 | +# that minimizes the *concave* function $c^\top x - x^\top Qx / 2$ |
| 19 | +# over the polyhedron obtained by intersecting the hypercube $[0, 1]^5$ |
| 20 | +# with the halfspace $10x_1 + 12x_2 + 11x_3 + 7x_4 + 4x_5 \le 40$. |
| 21 | + |
| 22 | +using Test #src |
| 23 | + |
| 24 | +using LinearAlgebra |
| 25 | +c = [42, 44, 45, 47, 47.5] |
| 26 | +Q = 100I |
| 27 | + |
| 28 | +using DynamicPolynomials |
| 29 | +@polyvar x[1:5] |
| 30 | +p = c'x - x' * Q * x / 2 |
| 31 | +using SumOfSquares |
| 32 | +K = @set x[1] >= 0 && x[1] <= 1 && |
| 33 | + x[2] >= 0 && x[2] <= 1 && |
| 34 | + x[3] >= 0 && x[3] <= 1 && |
| 35 | + x[4] >= 0 && x[4] <= 1 && |
| 36 | + x[5] >= 0 && x[5] <= 1 && |
| 37 | + 10x[1] + 12x[2] + 11x[3] + 7x[4] + 4x[5] <= 40 |
| 38 | + |
| 39 | +# We will now see how to find the optimal solution using Sum of Squares Programming. |
| 40 | +# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices. |
| 41 | + |
| 42 | +import Clarabel |
| 43 | +solver = Clarabel.Optimizer |
| 44 | + |
| 45 | +# A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. |
| 46 | +# The following function searches for the largest lower bound and finds zero using the `d`th level of the hierarchy`. |
| 47 | + |
| 48 | +function solve(d) |
| 49 | + model = SOSModel(solver) |
| 50 | + @variable(model, α) |
| 51 | + @objective(model, Max, α) |
| 52 | + @constraint(model, c, p >= α, domain = K, maxdegree = d) |
| 53 | + optimize!(model) |
| 54 | + println(solution_summary(model)) |
| 55 | + return model |
| 56 | +end |
| 57 | + |
| 58 | +# The first level of the hierarchy does not give any lower bound: |
| 59 | + |
| 60 | +model2 = solve(2) |
| 61 | +nothing # hide |
| 62 | +@test termination_status(model2) == MOI.INFEASIBLE #src |
| 63 | + |
| 64 | +# Indeed, as the constraints have degree 1 and their multipliers are SOS |
| 65 | +# so they have an even degree, with `maxdegree` 2 we can only use degree 0 |
| 66 | +# multipliers hence constants. The terms of maximal degree in resulting |
| 67 | +# sum will therefore only be in `-x' * Q * x/2` hence it is not SOS whatever |
| 68 | +# is the value of the multipliers. Let's try with `maxdegree` 3 so that the |
| 69 | +# multipliers can be quadratic. |
| 70 | +# This second level is now feasible and gives a lower bound of `-22`. |
| 71 | + |
| 72 | +model3 = solve(3) |
| 73 | +nothing # hide |
| 74 | +@test objective_value(model3) ≈ -22 rtol=1e-4 #src |
| 75 | +@test termination_status(model3) == MOI.OPTIMAL #src |
0 commit comments