Open baggepinnen opened 1 year ago
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
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))
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?
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 ?
I'm looking for functionality akin to
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?