Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated XFZ18under
). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.
In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using XFZ18under
. This is a Julia implementation that relies on the JuMP ecosystem (JuMP
, PolyJuMP
, SumOfSquares
, MathProgInterface
, MathOptInterfaceMosek
) and the JuliaAlgebra
ecosystem (MultivariatePolynomials
, DynamicPolynomials
). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.
References:
Quick links to documentation:
Consider the following van-der-Pol system:
$$ \dot{x}_1 = x_2 \\ \dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2 $$using MultivariatePolynomials,
JuMP,
PolyJuMP,
SumOfSquares,
DynamicPolynomials,
MathOptInterfaceMosek,
MathematicalSystems,
MosekTools
const ∂ = differentiate
differentiate (generic function with 18 methods)
# symbolic variables
@polyvar x₁ x₂ t
# time duration (scaled, see dynamics below)
T = 1.0
# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂]
# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25
# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 25 - x₁^2 - x₂^2
# degree of the relaxation
k = 12
# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)
# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(Mosek.Optimizer))
# add unknown Φ to the model
@variable(model, Φ, Poly(XT))
# jacobian
∂t = α -> ∂(α, t)
∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
LΦ = ∂t(Φ) + ∂xf(Φ)
# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)
# scalar variable
@variable(model, ϵ)
dom1 = @set t*(T-t) >= 0 && g >= 0
dom2 = @set g >= 0
@constraint(model, ϵ >= 0.)
@constraint(model, LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)
@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)
@objective(model, Min, ϵ)
optimize!(model)
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 1543 Cones : 0 Scalar variables : 456 Matrix variables : 10 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 1543 Cones : 0 Scalar variables : 456 Matrix variables : 10 Integer variables : 0 Optimizer - threads : 8 Optimizer - solved problem : the primal Optimizer - Constraints : 1542 Optimizer - Cones : 1 Optimizer - Scalar variables : 457 conic : 456 Optimizer - Semi-definite variables: 10 scalarized : 30074 Factor - setup time : 0.16 dense det. time : 0.00 Factor - ML order time : 0.06 GP order time : 0.00 Factor - nonzeros before factor : 4.80e+05 after factor : 7.71e+05 Factor - dense dim. : 2 flops : 9.25e+08 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.0e+00 8.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.18 1 9.1e-01 3.6e+00 6.2e-01 -7.85e-06 3.887272364e-02 -7.277377899e-04 4.6e-01 0.41 2 2.7e-01 1.1e+00 4.6e-01 9.08e-01 7.784730830e-01 4.303286309e-03 1.3e-01 0.65 3 6.2e-02 2.5e-01 4.5e-01 1.70e+00 2.542790179e-01 2.166139301e-02 3.1e-02 0.87 4 1.7e-02 6.7e-02 2.7e-01 1.87e+00 1.039863645e-01 6.446008794e-02 8.4e-03 1.10 5 9.8e-03 3.9e-02 2.1e-01 1.38e+00 1.294203141e-01 1.090778203e-01 4.9e-03 1.33 6 5.5e-03 2.2e-02 1.6e-01 1.15e+00 5.083213108e-02 4.002533190e-02 2.8e-03 1.55 7 1.9e-03 7.7e-03 1.0e-01 1.25e+00 2.537306998e-02 2.192283756e-02 9.7e-04 1.78 8 7.9e-04 3.2e-03 6.5e-02 1.18e+00 1.234485829e-02 1.103229702e-02 4.0e-04 2.00 9 1.7e-04 6.8e-04 3.2e-02 1.12e+00 3.993774249e-03 3.722721925e-03 8.5e-05 2.23 10 3.9e-05 1.6e-04 1.5e-02 1.04e+00 1.425668588e-03 1.364297146e-03 2.0e-05 2.45 11 1.1e-05 4.5e-05 7.1e-03 8.67e-01 8.002658418e-04 7.817030218e-04 5.7e-06 2.67 12 4.4e-06 1.8e-05 3.5e-03 6.29e-01 7.247429610e-04 7.167569170e-04 2.2e-06 2.89 13 1.1e-06 4.6e-06 9.2e-04 2.91e-01 1.130477537e-03 1.129671794e-03 5.7e-07 3.12 14 5.1e-07 2.0e-06 3.4e-04 -2.08e-01 1.970391763e-03 1.973234438e-03 2.5e-07 3.37 15 1.8e-07 6.0e-07 7.3e-05 -3.88e-01 4.513283862e-03 4.525581099e-03 7.3e-08 3.60 16 1.0e-07 1.7e-07 1.5e-05 -4.09e-01 9.367423708e-03 9.394061798e-03 2.0e-08 3.84 17 4.0e-08 6.3e-08 4.5e-06 -4.44e-01 1.682736656e-02 1.687283866e-02 7.6e-09 4.18 18 3.1e-08 3.2e-08 2.2e-06 -2.46e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.61 19 3.1e-08 3.2e-08 2.2e-06 -3.04e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.98 Optimizer terminated. Time: 5.25 Relaxation order : k = 12 STALL JuMP.termination_status(model) = SLOW_PROGRESS JuMP.primal_status(model) = FEASIBLE_POINT JuMP.dual_status(model) = FEASIBLE_POINT STALL JuMP.objective_bound(model) = 0.0 STALL JuMP.objective_value(model) = 0.02305835379792861
Package | k | Constraints | Scalar variables | Matrix variables | Time(s) |
---|---|---|---|---|---|
JuMP v0.18.5 | 2 | 245 | 175 | 8 | < 1 |
JuMP v0.19.0 | 2 | 83 | 13 | 8 | < 1 |
YALMIP | 2 | 152 | 63 | 10 | < 1 |
JuMP v0.18.5 | 3 | 893 | 715 | 10 | ~ 1 |
JuMP v0.19.0 | 3 | 199 | 21 | 10 | < 1 |
YALMIP | 3 | 254 | 121 | 10 | ~ 1 |
JuMP v0.18.5 | 4 | 893 | 730 | 10 | ~1 |
JuMP v0.19.0 | 4 | 199 | 36 | 10 | < 1 |
YALMIP | 4 | 394 | 206 | 10 | 1.18 |
JuMP v0.18.5 | 5 | 2639 | 2309 | 10 | 1.17 |
JuMP v0.19.0 | 5 | 387 | 57 | 10 | < 1 |
YALMIP | 5 | 578 | 323 | 10 | 0.11 |
JuMP v0.18.5 | 6 | 2639 | 2337 | 10 | 0.90 |
JuMP v0.19.0 | 6 | 387 | 85 | 10 | < 1 |
YALMIP | 6 | 812 | 477 | 10 | 1.10 |
JuMP v0.18.5 | 7 | 6725 | 6183 | 10 | 5.62 |
JuMP v0.19.0 | 7 | 663 | 121 | 10 | < 1 |
YALMIP | 7 | 1102 | 673 | 10 | 1.52 |
JuMP v0.18.5 | 8 | 6725 | 6228 | 10 | 6.06 |
JuMP v0.19.0 | 8 | 663 | 166 | 10 | < 1 |
YALMIP | 8 | 1454 | 916 | 10 | 1.10 |
JuMP v0.18.5 | 9 | 15269 | 14447 | 10 | 41.2 |
JuMP v0.19.0 | 9 | 1043 | 221 | 10 | 1.70 |
YALMIP | 9 | 1874 | 1211 | 10 | 1.58 |
JuMP v0.18.5 | 10 | 15269 | 14513 | 10 | 38.1 |
JuMP v0.19.0 | 10 | 1043 | 287 | 10 | 1.67 |
YALMIP | 10 | 2368 | 1563 | 10 | 2.73 |
JuMP v0.18.5 | 11 | 31617 | 30439 | 10 | 244 |
JuMP v0.19.0 | 11 | 1543 | 365 | 10 | 4.88 |
YALMIP | 11 | 2942 | 1977 | 10 | 2.30 |
JuMP v0.18.5 | 12 | 31617 | 30530 | 10 | 231 |
JuMP v0.19.0 | 12 | 1543 | 456 | 10 | 5.02 |
YALMIP | 12 | 3602 | 2458 | 10 | 6.57 |
Related: https://github.com/JuliaOpt/MosekTools.jl/issues/11.
MOI.get(model, MOI.SolveTime())
STALL
5.254102945327759
# Recovering the solution:
ϵopt = JuMP.objective_value(model)
# Punder <= 0
Punder = subs(JuMP.value(model[:Φ]), t => T)
STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL
# Pover <= 0
Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)
STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL STALL
See notebook Constraints.ipynb
.
# no necesito convertir, si utilizo el IntervalContractors#fix_zero_power
using ModelingToolkit
# TODO: temporary definition here, needs convert(Variable, z^2 + t^2)
vars = ModelingToolkit.@variables x1 x2
function tomtexpr(poly::Polynomial)
global vars
expr = 0.0
for (i, term) in enumerate(poly.x)
if iszero(exponents(term))
expr += poly.a[i]
else
expr += poly.a[i] * prod([vars[j]^value for (j, value) in enumerate(exponents(term)) if value != 0])
end
end
return expr
end
Punder_mt = tomtexpr(Punder)
include("constraints.jl")
using Plots
S = Separator(-Inf..0.0, Punder_mt)
p = pave(S, IntervalBox(-2..2, 2), 0.1)
S = Separator(-Inf..0.0, Punder_mt)
p = pave(S, IntervalBox(-3..3, 2), 0.1)
plot(p.inner, lw=0, aspectratio=1, label="inner")
plot!(p.boundary, lw=0, label="boundary", legend=:topleft)
Punder(2.0, 0.5)
S = Separator(-Inf..0.0, Punder_mt)
p = pave(S, IntervalBox(-2..2, 2), 0.1)
Pover_mt = tomtexpr(Pover)
Sover = Separator(-Inf..0.0, Pover_mt)
pover = pave(Sover, IntervalBox(-2..2, 2), 0.1)
plot(p.inner, lw=0, aspectratio=1, label="under", alpha=.4)
plot!(pover.inner, lw=0, label="over", legend=:topleft, alpha=.4)
using LinearAlgebra, StaticPolynomials
# filter out the smallest coefficients?
filter_coeffs = false
if filter_coeffs
max_coeff = maximum(abs.(Punder.a))
# plots are *very* sensitive to small values
# make tol bigger than 0 to filter out some coefficients
tol = 0.0 # max_coeff / 1e19
kept_coeffs = map(c -> abs(c) > tol, Punder.a)
Punder = dot(Punder.a[kept_coeffs], Punder.x[kept_coeffs])
Pover = dot(Pover.a[kept_coeffs], Pover.x[kept_coeffs]);
end
# we now convert to a static polynomial for faster evaluation
Punder_st = StaticPolynomials.Polynomial(Punder)
Pover_st = StaticPolynomials.Polynomial(Pover);
_Punder(x, y) = Punder_st([x, y])
_Pover(x, y) = Pover_st([x, y])
@btime Punder(1.0, 1.0)
@btime Punder_st([1.0, 1.0])
2e-6 / 96e-9
Punder_st
using ImplicitEquations, Plots
gr()
G = plot()
plot!(G, _Punder ⩵ 0., xlims=(-4, 4), ylims=(-4, 4), color="red")
plot!(G, _Pover ⩵ 0., xlims=(-4, 4), ylims=(-4, 4), color="blue")
G
G = plot()
Gu = plot(_Punder ≪ 0., xlims=(-2, 2), ylims=(-2, 2))
Go = plot(_Pover ≪ 0., xlims=(-2, 2), ylims=(-2, 2))
plot(Gu, Go)
Punder(2.0, 1.5)
using IntervalConstraintProgramming, ValidatedNumerics
@function Funder(x₁, x₂) = _Punder(x₁, x₂)
Bx₁x₂ = IntervalBox(-10..10, -10..10)
Cunder = IntervalConstraintProgramming.@constraint Funder(x₁, x₂) <= 0.0 # no
paving = pave(Cunder, Bx₁x₂)
# symbolic variables
@polyvar x₁ x₂ t
# time duration (scaled, see dynamics below)
T = 1.0
# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂]
# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25
# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 25 - x₁^2 - x₂^2
# degree of the relaxation
k = 12
# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)
# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(Mosek.Optimizer))
# add unknown Φ to the model
@variable(model, Φ, Poly(XT))
# jacobian
∂t = α -> ∂(α, t)
∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
LΦ = ∂t(Φ) + ∂xf(Φ)
# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)
# scalar variable
@variable(model, ϵ)
dom1 = @set t*(T-t) >= 0 && g >= 0
dom2 = @set g >= 0
@constraint(model, ϵ >= 0.)
@constraint(model, LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)
@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)
# add a safety constarint
#@constraint(model, x₂ <= 100.)
@objective(model, Min, ϵ)
optimize!(model)
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 1543 Cones : 0 Scalar variables : 456 Matrix variables : 10 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 1543 Cones : 0 Scalar variables : 456 Matrix variables : 10 Integer variables : 0 Optimizer - threads : 8 Optimizer - solved problem : the primal Optimizer - Constraints : 1542 Optimizer - Cones : 1 Optimizer - Scalar variables : 457 conic : 456 Optimizer - Semi-definite variables: 10 scalarized : 30074 Factor - setup time : 0.15 dense det. time : 0.00 Factor - ML order time : 0.06 GP order time : 0.00 Factor - nonzeros before factor : 4.80e+05 after factor : 7.71e+05 Factor - dense dim. : 2 flops : 9.25e+08 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.0e+00 8.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.17 1 9.1e-01 3.6e+00 6.2e-01 -7.85e-06 3.887272364e-02 -7.277377899e-04 4.6e-01 0.39 2 2.7e-01 1.1e+00 4.6e-01 9.08e-01 7.784730830e-01 4.303286309e-03 1.3e-01 0.61 3 6.2e-02 2.5e-01 4.5e-01 1.70e+00 2.542790179e-01 2.166139301e-02 3.1e-02 0.81 4 1.7e-02 6.7e-02 2.7e-01 1.87e+00 1.039863645e-01 6.446008794e-02 8.4e-03 1.08 5 9.8e-03 3.9e-02 2.1e-01 1.38e+00 1.294203141e-01 1.090778203e-01 4.9e-03 1.30 6 5.5e-03 2.2e-02 1.6e-01 1.15e+00 5.083213108e-02 4.002533190e-02 2.8e-03 1.54 7 1.9e-03 7.7e-03 1.0e-01 1.25e+00 2.537306998e-02 2.192283756e-02 9.7e-04 1.78 8 7.9e-04 3.2e-03 6.5e-02 1.18e+00 1.234485829e-02 1.103229702e-02 4.0e-04 2.01 9 1.7e-04 6.8e-04 3.2e-02 1.12e+00 3.993774249e-03 3.722721925e-03 8.5e-05 2.22 10 3.9e-05 1.6e-04 1.5e-02 1.04e+00 1.425668588e-03 1.364297146e-03 2.0e-05 2.48 11 1.1e-05 4.5e-05 7.1e-03 8.67e-01 8.002658418e-04 7.817030218e-04 5.7e-06 2.69 12 4.4e-06 1.8e-05 3.5e-03 6.29e-01 7.247429610e-04 7.167569170e-04 2.2e-06 2.94 13 1.1e-06 4.6e-06 9.2e-04 2.91e-01 1.130477537e-03 1.129671794e-03 5.7e-07 3.15 14 5.1e-07 2.0e-06 3.4e-04 -2.08e-01 1.970391763e-03 1.973234438e-03 2.5e-07 3.38 15 1.8e-07 6.0e-07 7.3e-05 -3.88e-01 4.513283862e-03 4.525581099e-03 7.3e-08 3.58 16 1.0e-07 1.7e-07 1.5e-05 -4.09e-01 9.367423708e-03 9.394061798e-03 2.0e-08 3.79 17 4.0e-08 6.3e-08 4.5e-06 -4.44e-01 1.682736656e-02 1.687283866e-02 7.6e-09 4.08 18 3.1e-08 3.2e-08 2.2e-06 -2.46e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.48 19 3.1e-08 3.2e-08 2.2e-06 -3.04e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.83 Optimizer terminated. Time: 5.01 Relaxation order : k = 12 STALL JuMP.termination_status(model) = SLOW_PROGRESS JuMP.primal_status(model) = FEASIBLE_POINT JuMP.dual_status(model) = FEASIBLE_POINT STALL JuMP.objective_bound(model) = 0.0 STALL JuMP.objective_value(model) = 0.02305835379792861
# symbolic variables
@polyvar x₁ x₂ t
# time duration (scaled, see dynamics below)
T = 1.0
# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂]
# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25
# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 25 - x₁^2 - x₂^2
# degree of the relaxation
k = 12
# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)
# create a SOS JuMP model to solve with Mosek
model_verif = SOSModel(with_optimizer(Mosek.Optimizer))
# add unknown Φ to the model
@variable(model, Φ, Poly(XT))
# jacobian
∂t = α -> ∂(α, t)
∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
LΦ = ∂t(Φ) + ∂xf(Φ)
# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)
# scalar variable
@variable(model, ϵ)
dom1 = @set t*(T-t) >= 0 && g >= 0
dom2 = @set g >= 0
@constraint(model, ϵ >= 0.)
@constraint(model, LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)
@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)
# add a safety constarint
#@constraint(model, x₂ <= 100.)
@objective(model, Min, ϵ)
Punder
model_verif = Model(with_optimizer(Mosek.Optimizer))
@variable(model_verif, t)
@variable(model_verif, x₁)
@variable(model_verif, x₂)
@constraint(model_verif, Punder <= 0)
@objective(model_verif, Max, x₂)
optimize!(model_verif)
PolyJuMP is just a JuMP extension for modelling Polynomial Optimization: it does not implement any reformulation. To use automatic sums of squares (SOS) reformulations, install the SumOfSquares Julia package and try `using SumOfSquares` and `setpolymodule!(SumOfSquares)` or use `SOSModel` instead of `Model`. Stacktrace: [1] error(::String) at ./error.jl:33 [2] get_default(::Nothing) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:20 [3] getdefault(::PolyJuMP.Data, ::Type{PolyJuMP.NonNegPoly}) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:12 [4] getdefault(::PolyJuMP.Data, ::PolyJuMP.NonNegPoly) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:11 [5] getdefault(::Model, ::PolyJuMP.NonNegPoly) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:7 [6] add_constraint(::Model, ::PolyJuMP.Constraint{getfield(JuMP, Symbol("#_error#55")){Symbol},Polynomial{true,Float64},PolyJuMP.NonNegPoly,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::String) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/constraint.jl:166 [7] top-level scope at /Users/forets/.julia/packages/JuMP/jnmGG/src/macros.jl:621 [8] top-level scope at In[29]:6