parf
October 9, 2022, 4:14pm
1
I am trying to solve the following optimization problem:
using DynamicPolynomials
using SumOfSquares
using SDPAFamily
χ=1.0
@polyvar x1 x2 x3 x4
f = [2*x1*x4 - 2*x2*x3 - x1,
-2*x1*x3 - 2*x2*x4 - x2,
2*x1*x2,
x2*x2 - x1*x1]
model = SOSModel(() -> SDPAFamily.Optimizer{Float64}( # JuMP only supports Float64
variant = :sdpa_gmp, # use the arbitrary precision variant
params = SDPAFamily.STABLE_BUT_SLOW,
presolve=true,
#params = ( epsilonStar = 1e-10, # constraint tolerance
# epsilonDash = 1e-10, # normalized duality gap tolerance
# precision = 200) # arithmetric precision used in sdpa_gmp
))
X = monomials([x1, x2, x3, x4], 1:4);
@variable(model, V, Poly(X));
LHS = differentiate(V, x1)*f[1]+differentiate(V, x2)*f[2] + differentiate(V, x3)*f[3] +
+ differentiate(V, x4)*f[4] + x3^2 + x4^2 + differentiate(differentiate(V, x3),x3)/(4*χ) +
+ differentiate(differentiate(V, x4),x4)/(4*χ);
@variable(model, L);
@objective(model, Max, L);
@constraint(model, LHS >= L);
JuMP.optimize!(model)
The problem is solved if X contains only low powers of monomials (from 1 to 5), but with an increase in this parameter, the solver ceases to produce even the solution that it found before. I tried a lot of different solvers, but without success.
Maybe someone can suggest what can be done to increase the number of monomials in X?
1 Like
odow
October 9, 2022, 8:35pm
2
what can be done to increase the number of monomials in X ?
What powers do you want to get up to?
I tried a lot of different solvers, but without success.
Which ones? Have you tried Mosek?
parf
October 9, 2022, 9:19pm
3
What powers do you want to get up to?
The bigger, the better. I would expect that as the power of X increases, the value of L will also increase. From an independent stochastic simulation, I know the “theoretical limit” of L and I want to understand how close I can get using Sum of Squares programming.
Which ones? Have you tried Mosek?
Yes, I tried Mosek, CDCS, Clarabel, COSMO, CSDP, Hypatia, ProxSDP, SCS, SeDuMi, SDPA-GMP. However, I always used the Float64 data type because JUMP does not support BigFloat. Can I solve this problem with julia but using the BigFloat type and how? Perhaps then I will be able to increase the degree of X.
odow
October 10, 2022, 11:08pm
4
@blegat might have some suggestions for numerical improvements. He is the SumOfSquares.jl expert.
blegat
October 11, 2022, 1:38pm
5
What do you mean by “cease to produce any solution” ? Does it take too much time or does it finish but with an inconclusive status ?
parf
October 11, 2022, 4:43pm
6
It finishes with an inconclusive status. For example, the output for monomials of order 8 (Mosek solver):
Problem
Name :
Objective sense : maximize
Type : CONIC (conic optimization problem)
Constraints : 688
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 0
Scalar variables : 495
Matrix variables : 1
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 82
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 17
Presolve terminated. Time: 0.00
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name :
Objective sense : maximize
Type : CONIC (conic optimization problem)
Constraints : 688
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 0
Scalar variables : 495
Matrix variables : 1
Integer variables : 0
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 462
Optimizer - Cones : 1
Optimizer - Scalar variables : 319 conic : 319
Optimizer - Semi-definite variables: 1 scalarized : 2145
Factor - setup time : 0.02 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 1.07e+05 after factor : 1.07e+05
Factor - dense dim. : 0 flops : 4.35e+07
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.02
1 2.8e-01 2.8e-01 1.3e-01 8.41e-01 5.658786273e-02 1.292251824e-01 2.8e-01 0.02
2 1.2e-01 1.2e-01 3.1e-02 8.56e-01 2.971659068e-01 3.385846790e-01 1.2e-01 0.03
3 2.2e-02 2.2e-02 1.8e-03 1.24e+00 4.004710586e-01 4.118891854e-01 2.2e-02 0.03
4 3.7e-03 3.7e-03 1.6e-04 1.13e+00 3.933816935e-01 3.947614139e-01 3.7e-03 0.03
5 1.0e-03 1.0e-03 3.4e-05 6.84e-01 3.845216114e-01 3.846743932e-01 1.0e-03 0.03
6 2.8e-04 2.8e-04 8.1e-06 2.92e-01 3.631527616e-01 3.628386926e-01 2.8e-04 0.05
7 9.4e-05 9.4e-05 2.5e-06 1.94e-01 3.486957620e-01 3.482766828e-01 9.4e-05 0.05
8 2.6e-05 2.6e-05 6.6e-07 7.92e-02 3.312750693e-01 3.307510053e-01 2.6e-05 0.05
9 7.5e-06 7.5e-06 1.7e-07 2.59e-01 3.203259811e-01 3.198739211e-01 7.5e-06 0.06
10 2.3e-06 2.3e-06 5.1e-08 3.58e-02 3.085929591e-01 3.081301685e-01 2.3e-06 0.06
11 7.2e-07 7.2e-07 1.4e-08 1.88e-01 3.004379714e-01 3.000531242e-01 7.2e-07 0.06
12 2.0e-07 2.0e-07 3.9e-09 7.59e-02 2.916395607e-01 2.912761518e-01 2.0e-07 0.08
13 5.8e-08 5.8e-08 1.0e-09 2.22e-01 2.855252086e-01 2.852327366e-01 5.8e-08 0.08
14 2.5e-08 1.6e-08 2.6e-10 8.99e-02 2.793028964e-01 2.790338667e-01 1.6e-08 0.08
15 8.6e-09 4.2e-09 5.8e-11 3.29e-01 2.749540167e-01 2.747579127e-01 4.2e-09 0.09
16 3.4e-09 1.6e-09 2.5e-11 -1.04e-02 2.703414885e-01 2.700904299e-01 1.6e-09 0.11
17 1.3e-09 5.8e-10 6.6e-12 5.23e-01 2.689030845e-01 2.687708780e-01 5.8e-10 0.11
18 1.2e-09 4.9e-10 5.5e-12 3.23e-01 2.684525884e-01 2.683226175e-01 4.8e-10 0.11
19 1.1e-09 3.3e-10 3.7e-12 2.58e-01 2.673859138e-01 2.672588002e-01 3.2e-10 0.12
20 1.1e-09 3.2e-10 3.6e-12 4.80e-01 2.673412672e-01 2.672160643e-01 3.1e-10 0.12
21 1.1e-09 3.2e-10 3.6e-12 4.80e-01 2.673412672e-01 2.672160643e-01 3.1e-10 0.14
22 1.1e-09 3.2e-10 3.6e-12 4.80e-01 2.673412672e-01 2.672160643e-01 3.1e-10 0.16
23 1.1e-09 3.2e-10 3.6e-12 1.00e+00 2.673411582e-01 2.672159702e-01 3.1e-10 0.17
24 1.1e-09 3.2e-10 3.6e-12 1.00e+00 2.673411037e-01 2.672159232e-01 3.1e-10 0.17
Optimizer terminated. Time: 0.20
JuMP.primal_status(model)
UNKNOWN_RESULT_STATUS::ResultStatusCode = 8
blegat
October 12, 2022, 11:28am
7
It looks quite close to feasible, maybe you can relax the feasibility tolerance : MSK_DPAR_INTPNT_CO_TOL_PFEASConstraints · SumOfSquares .
parf
October 13, 2022, 7:31am
8
Thank you @blegat , I will try.
parf
October 14, 2022, 11:28am
9
Hi, @blegat . Actually, your suggestions helped me a lot. I rescaled the initial system of equations and reduced the basis in the space of which the optimization takes place (based on some physical reasoning). Now the problem becomes:
using DynamicPolynomials
using SumOfSquares
using MosekTools
χ=1.0
@polyvar x1 x2 x3 x4
f = [2*x1*x4 - 2*x2*x3 - x1*χ^(1/3),
-2*x1*x3 - 2*x2*x4 - x2*χ^(1/3),
2*x1*x2,
x2*x2 - x1*x1]
model = SOSModel(Mosek.Optimizer)
set_optimizer_attribute(model, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-3)
deg = 20
@polyvar n1 n2 J
arr = Polynomial{true, Float64}[]
for i in 0:ceil(deg/2)
for j in 0:ceil(deg/2)
for k in 0:ceil(deg/3)
if 1<=2*i+2*j+3*k<=deg
p = (n1^Int(i))*(n2^Int(j))*(J^Int(k))
push!(arr, p)
end
end
end
end
arr2 = subs(arr, n1=>x1*x1+x2*x2, n2=>x3*x3+x4*x4, J=>-4*x1*x2*x3+2*x1*x1*x4-2*x2*x2*x4)
basis = FixedPolynomialBasis(arr2)
@variable(model, V, Poly(basis))
LHS = differentiate(V, x1)*f[1]+differentiate(V, x2)*f[2] + differentiate(V, x3)*f[3] +
+ differentiate(V, x4)*f[4] + x3^2 + x4^2 + differentiate(differentiate(V, x3),x3)/(4) +
+ differentiate(differentiate(V, x4),x4)/(4);
@variable(model, L);
@objective(model, Max, L);
@constraint(model, LHS >= L);
JuMP.optimize!(model);
println(JuMP.primal_status(model))
println(objective_value(model)/χ^(2/3))
The output for deg=20 is following:
Problem
Name :
Objective sense : maximize
Type : CONIC (conic optimization problem)
Constraints : 11015
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 0
Scalar variables : 192
Matrix variables : 1
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 2
Eliminator terminated.
Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 152
Presolve terminated. Time: 0.00
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name :
Objective sense : maximize
Type : CONIC (conic optimization problem)
Constraints : 11015
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 0
Scalar variables : 192
Matrix variables : 1
Integer variables : 0
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 10565
Optimizer - Cones : 1
Optimizer - Scalar variables : 191 conic : 191
Optimizer - Semi-definite variables: 1 scalarized : 490545
Factor - setup time : 12.44 dense det. time : 4.50
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 5.57e+07 after factor : 5.57e+07
Factor - dense dim. : 0 flops : 8.77e+11
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 13.11
1 4.8e-01 4.8e-01 3.3e-01 9.44e-01 4.550285069e-02 6.429942074e-02 4.8e-01 30.56
2 1.3e-01 1.3e-01 9.7e-02 5.14e-01 4.848911013e-01 1.575909731e-01 1.3e-01 49.72
3 4.0e-02 4.0e-02 5.1e-02 -6.23e-01 1.500290024e+00 1.720124574e-01 4.0e-02 70.61
4 1.4e-02 1.4e-02 2.6e-02 -6.48e-01 3.393740475e+00 2.310617171e-01 1.4e-02 90.11
5 4.7e-03 4.7e-03 1.2e-02 -6.54e-01 6.762067251e+00 2.970061982e-01 4.7e-03 108.58
6 1.4e-03 1.4e-03 5.2e-03 -6.23e-01 1.345494511e+01 3.677029505e-01 1.4e-03 128.56
7 5.3e-04 5.3e-04 2.4e-03 -4.72e-01 2.170525715e+01 4.316049177e-01 5.3e-04 147.66
8 3.2e-04 3.2e-04 1.3e-03 -1.14e-01 1.631440577e+01 4.815503283e-01 3.0e-04 166.56
9 5.2e-05 5.2e-05 2.4e-04 -1.70e-01 2.130817677e+01 5.856556728e-01 5.4e-05 187.61
10 2.1e-05 2.1e-05 8.5e-05 1.96e-01 1.664418491e+01 6.194935873e-01 2.1e-05 206.39
11 3.5e-06 3.5e-06 1.0e-05 3.19e-01 9.105502953e+00 6.693353513e-01 3.8e-06 261.75
12 7.9e-07 7.9e-07 1.3e-06 6.24e-01 3.315519459e+00 6.989735069e-01 8.5e-07 300.72
13 3.2e-07 1.4e-07 9.2e-08 8.62e-01 1.136843993e+00 7.058032411e-01 1.4e-07 322.11
14 7.1e-07 1.1e-07 6.1e-08 1.01e+00 1.036732868e+00 7.054957416e-01 1.1e-07 357.44
15 9.9e-07 9.6e-08 5.3e-08 1.00e+00 1.005906292e+00 7.055728351e-01 1.0e-07 377.91
16 9.8e-07 4.3e-08 1.6e-08 9.76e-01 8.441681925e-01 7.062274679e-01 4.5e-08 397.52
17 3.1e-07 1.2e-08 2.6e-09 9.83e-01 7.497089483e-01 7.073477438e-01 1.3e-08 417.08
18 5.7e-07 5.5e-09 7.7e-10 9.98e-01 7.267646158e-01 7.074757932e-01 6.0e-09 434.61
19 4.8e-07 1.4e-09 9.8e-11 9.97e-01 7.125862404e-01 7.076853996e-01 1.5e-09 452.86
20 1.1e-06 1.3e-09 8.4e-11 1.00e+00 7.121236355e-01 7.076887856e-01 1.4e-09 489.89
21 9.3e-07 6.5e-10 1.8e-11 1.00e+00 7.092900816e-01 7.076978265e-01 5.0e-10 509.56
22 4.5e-07 2.9e-10 2.6e-12 1.00e+00 7.081567134e-01 7.077133249e-01 1.4e-10 529.94
23 7.8e-07 1.9e-09 3.2e-13 1.00e+00 7.078313046e-01 7.077185736e-01 3.4e-11 548.22
24 7.8e-07 1.9e-09 3.2e-13 1.00e+00 7.078313046e-01 7.077185736e-01 3.4e-11 588.75
Optimizer terminated. Time: 612.20
UNKNOWN_RESULT_STATUS
0.7078313045648063
and as I understand the result is again quite close to the feseability. Maybe you can tell me what other parameters of the solver should be paid attention to in order to try to improve its performance? I want to play around with the parameters a bit, but there are so many of them in the documentation that I can’t pick out the key ones.
You can see the parameters grouped by topic here:https://docs.mosek.com/10.0/capi/param-groups.html#doc-param-groups
The likely ones relevant to you all start with MSK_DPAR_INTPNT_CO_TOL_ under ‘Conic interior-point method’.
2 Likes
blegat
October 27, 2022, 1:45pm
11
I just added an example on how to change the basis used for the Gram matrix in Mminimization of a univariate polynomial · SumOfSquares
1 Like