|
| 1 | +# # Goldstein-price function |
| 2 | + |
| 3 | +#md # [](@__BINDER_ROOT_URL__/generated/Polynomial Optimization/goldstein_price.ipynb) |
| 4 | +#md # [](@__NBVIEWER_ROOT_URL__/generated/Polynomial Optimization/goldstein_price.ipynb) |
| 5 | +# **Contributed by**: Benoît Legat |
| 6 | + |
| 7 | +# In this example, we consider the minimization of the [Goldstein-price function](https://en.wikipedia.org/wiki/Test_functions_for_optimization). |
| 8 | + |
| 9 | +using Test #src |
| 10 | +using SumOfSquares |
| 11 | +using DynamicPolynomials |
| 12 | + |
| 13 | +# Create *symbolic* variables (not JuMP *decision* variables) |
| 14 | + |
| 15 | +@polyvar x[1:2] |
| 16 | + |
| 17 | +# To use Sum-of-Squares Programming, we first need to pick an SDP solver, |
| 18 | +# see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices. |
| 19 | + |
| 20 | +import Clarabel |
| 21 | +using Dualization |
| 22 | +model = SOSModel(dual_optimizer(Clarabel.Optimizer)) |
| 23 | + |
| 24 | +# Create a JuMP decision variable for the lower bound |
| 25 | + |
| 26 | +@variable(model, γ) |
| 27 | + |
| 28 | +# `f(x)` is the Goldstein-Price function |
| 29 | + |
| 30 | +f1 = x[1] + x[2] + 1 |
| 31 | +f2 = 19 - 14*x[1] + 3*x[1]^2 - 14*x[2] + 6*x[1]*x[2] + 3*x[2]^2 |
| 32 | +f3 = 2*x[1] - 3*x[2] |
| 33 | +f4 = 18 - 32*x[1] + 12*x[1]^2 + 48*x[2] - 36*x[1]*x[2] + 27*x[2]^2 |
| 34 | +f = (1 + f1^2*f2) * (30 + f3^2*f4) |
| 35 | + |
| 36 | +# Constraints `f(x) - γ` to be a sum of squares |
| 37 | + |
| 38 | +con_ref = @constraint(model, f >= γ) |
| 39 | +@objective(model, Max, γ) |
| 40 | +optimize!(model) |
| 41 | + |
| 42 | +# The lower bound found is 3 |
| 43 | + |
| 44 | +@test objective_value(model) ≈ 3 rtol=1e-3 #src |
| 45 | +solution_summary(model) |
| 46 | + |
| 47 | +# The moment matrix is as follows, we can already see the global minimizer |
| 48 | +# `[0, -1]` from the entries `(2, 1)` and `(3, 1)`. |
| 49 | +# This heuristic way to obtain solutions to the polynomial optimization problem |
| 50 | +# is suggested in [L09, (6.15)]. |
| 51 | +# |
| 52 | +# [L09] Laurent, Monique. |
| 53 | +# *Sums of squares, moment matrices and optimization over polynomials.* |
| 54 | +# Emerging applications of algebraic geometry (2009): 157-270. |
| 55 | + |
| 56 | +ν = moment_matrix(con_ref) |
| 57 | + |
| 58 | +# Many entries of the matrix actually have the same moment. |
| 59 | +# We can obtain the following list of these moments without duplicates |
| 60 | +# (ignoring when difference of entries representing the same moments is below `1e-5`) |
| 61 | + |
| 62 | +μ = measure(ν, atol = 1e-5) |
| 63 | + |
| 64 | +# The truncated moment matrix can then be obtained as follows |
| 65 | + |
| 66 | +ν_truncated = moment_matrix(μ, monomials(x, 0:3)) |
| 67 | + |
| 68 | +# Let's check if the flatness property is satisfied. |
| 69 | +# The rank of `ν_truncated` seems to be 1: |
| 70 | + |
| 71 | +using LinearAlgebra |
| 72 | +LinearAlgebra.svdvals(Matrix(ν_truncated.Q)) |
| 73 | +LinearAlgebra.rank(Matrix(ν_truncated.Q), rtol = 1e-3) |
| 74 | +@test LinearAlgebra.rank(Matrix(ν_truncated.Q), rtol = 1e-3) == 1 #src |
| 75 | +svdvals(Matrix(ν_truncated.Q)) |
| 76 | + |
| 77 | +# The rank of `ν` is clearly higher than 1, closer to 3: |
| 78 | + |
| 79 | +@test 3 <= LinearAlgebra.rank(Matrix(ν.Q), rtol = 1e-3) <= 4 #src |
| 80 | +svdvals(Matrix(ν.Q)) |
| 81 | + |
| 82 | +# Even if the flatness property is not satisfied, we can |
| 83 | +# still try extracting the minimizer with a low rank decomposition of rank 3. |
| 84 | +# We find the optimal solution again doing so: |
| 85 | + |
| 86 | +atoms = atomic_measure(ν, FixedRank(3)) #src |
| 87 | +@test length(atoms.atoms) == 1 #src |
| 88 | +@test atoms.atoms[1].center ≈ [0, -1] rtol=1e-3 #src |
| 89 | +atomic_measure(ν, FixedRank(3)) |
0 commit comments