# Sum of squares optimization

I have the following optimization problem that I want to solve by the sum of squares approach. It has a linear objective function and I am looking for the best set of (a, b, and c) that satisfies a specific polynomial.

min (1/3ax_max^3+1/2bx_max^2+cx_max+sin(x_max)- 1/3ax_min^3+1/2bx_min^2+cx_min+sin(x_min)); (x_min=-pi/2, x_max=pi/2)

subject to:
x<=pi/2;
x>=-pi/2;
I also have the following constraint that enforces (ax^2+bx+c) to lay above the Taylor expansion of cosine function which is non-negative if we enforce x to vary within [-pi/2, pi/2].
ax^2+bx+c > = 1-x^2/2! +x^4/4! -x^6/6! +x^8/8! -x^10/10! + epsilon;
(epsilon: a very small number)
I run the following code to solve the problem.

using MultivariatePolynomials
using JuMP
using PolyJuMP
using SumOfSquares
using DynamicPolynomials
using Mosek

@polyvar x
m = SOSModel(solver = MosekSolver())
@variable m a
@variable m b
@variable m c
P= ax^2+bx+c-1+x^2/factorial(2) -x^4/factorial(4) +x^6/factorial(6) -x^8/factorial(8) +x^10/factorial(10) – epsilon;
S = @set x >= -pi/2 && x <= pi/2
x_max=pi/2;
x_min=-pi/2;
@objective m Min (1/3ax_max^3+1/2bx_max^2+cx_max+sin(x_max)- 1/3ax_min^3+1/2bx_min^2+cx_min+sin(x_min));
@constraint( m, p >= 0, domain = S)

I got following error when I executed the code:
Can you please let me know how I can debug this. Thanks in advance!

You need to do `using SemialgebraicSets`

As a tip, you can use triple backticks (```) to enclose code to format it.

This:

```
println(“example”)
```

becomes this:

``````println("example")
``````

It makes your code easier to read 2 Likes

I run the following piece if code to obtain an understimator for cosine function. I have replace cosine function with its corresponding Taylor series.

``````using MultivariatePolynomials
using JuMP
using PolyJuMP
using SumOfSquares
using DynamicPolynomials
using Mosek
using SemialgebraicSets

@polyvar x
m = SOSModel(solver = MosekSolver())
@variable m a
@variable m b
@variable m c

epsilon = 1.387895246507398e-07
x_min = -3*pi/2
x_max = -pi/2

p = -1 + (x+pi)^2/factorial(2) - (x+pi)^4/factorial(4) + (x+pi)^6/factorial(6) - (x+pi)^8/factorial(8) + (x+pi)^10/factorial(10) - (x+pi)^12/factorial(12) + (x+pi)^14/factorial(14) - (x+pi)^16/factorial(16) + (x+pi)^18/factorial(18) - (x+pi)^20/factorial(20) + (x+pi)^20/factorial(20) - epsilon - a*x^2 - b*x - c
#p =  1 - x^2/factorial(2) + x^4/factorial(4) - x^6/factorial(6) + x^8/factorial(8) - x^10/factorial(10) + x^12/factorial(12) - x^14/factorial(14) + x^16/factorial(16) - x^18/factorial(18) + x^20/factorial(20) - epsilon -a*x^2 - b*x - c
s = @set x >= x_min && x <= x_max

@objective m Min   (sin(x_max) - 1/3*a*x_max^3 - 1/2*b*x_max^2 - c*x_max - sin(x_min) + 1/3*a*x_min^3 + 1/2*b*x_min^2 + c*x_min )
@constraint( m, p >= 0, domain = s)

solve(m)

println("a =", getvalue(a))
println("b =", getvalue(b))
println("c =", getvalue(c))
``````

After running the code I got following “WARNING: Not solved to optimality, status: Stall”
and the result do not provide a reasonable understimator for the cosine function within [-3pi/2, -pi/2] span. I will appreciate any help.

@Rasoul You might try to change the polynomial basis to have better numerical behavior. Note that this is currently only available on SumOfSquares master branch.

Thank you so much!

1 Like