Integrate polynomial over specified domain

[baggepinnen]

I'm looking for functionality akin to

integrate(poly, x, lb, ub)

i.e., to express the integral
$$\int_{l}^{u} p(x) dx$$ where $p(x)$ is a polynomial. Similar to YALMIP int

Scanning through the docs, I can't find this available. Is this supported by some of the polynomial packages that are supported by SumOfSquares.jl?

[baggepinnen]

Here's something that appears to work to some extent. I'm completely unfamiliar with this ecosystem of types, so I might not use the correct functions and types everywhere.

function integrate(p::Polynomial{C, T}, x::PolyVar{C}, l, u) where {C, T}
    # grlex order preserved
    i = something(findfirst(isequal(x), DynamicPolynomials._vars(p)), 0)
    S = typeof(zero(T) * 0)

    Z = copy.(p.x.Z)
    a = Vector{S}(undef, length(Z))
    for j in 1:length(Z)
        Z[j][i] += 1
        a[j] = p.a[j] / Z[j][i]
    end
    intpoly = Polynomial(a, MonomialVector(DynamicPolynomials._vars(p), Z))
    subs(intpoly, x=>u) - subs(intpoly, x=>l)
end

My remaining problem is that I'd like to use the integral of a polynomial as a cost function for optimization (maximize the area under a Lyapunov function) and I just can't figure out how to express this yet

[blegat]

You can find here code to integrate a multivariate polynomial over a hyperrectangle or arbitrary polytope: https://github.com/blegat/SetProg.jl/blob/master/src/L1_heuristic.jl
The idea is as follows: suppose p is a polynomial for which the coefficients are JuMP expressions and you have a function integrate that can integrate a monomial.
Then you can do sum(coefficient(term) * integrate(monomial(term)) for term in terms(p))

[baggepinnen]

Thanks for the link and the clever trick :)

The problem I had with having the integral of the polynomial as the objective could be worked around by this trivial trick

# @objective(model, Max, ∫J) # What I would like to do
@variable(model, t)
@objective(model, Max, t)
@constraint(model, ∫J == t)

it seems JuMP needs to be taught how to have a polynomial as objective?

[blegat]

With SumOfSquares, the decision variables are the JuMP variables that you create with @variable, not the polynomial variable that you create with @polyvar. When you write ∫J == t, it means that you want t to be equal to ∫J for all values of the polynomial variables. Does ∫J still contain polynomial variables after integration ?