Discrete-time Invariant Set

In this example we compute the maximal (resp. minimal) invariant set contained in (resp. containing) the square with vertices $(\pm 1, \pm 1)$ for the system

\[\begin{aligned} x_{k+1} & = -y_k\\ y_{k+1} & = x_k. \end{aligned}\]

The system is $x_{k+1} = A x_k$ where

\[A = \begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix}.\]

A set $S$ is invariant if $A S \subseteq S$.

using Polyhedra
h = HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1)
□ = polyhedron(h)

Polyhedron DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}}:
4-element iterator of HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([1.0, 0.0], 1.0)
 HalfSpace([-1.0, 0.0], 1.0)
 HalfSpace([0.0, 1.0], 1.0)
 HalfSpace([0.0, -1.0], 1.0)

We need to pick an SDP solver, see here for a list of available ones.

using SetProg
using SumOfSquares
import Hypatia
sdp_solver = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)

A = [0.0 -1.0
     1.0  0.0]

2×2 Matrix{Float64}:
 0.0  -1.0
 1.0   0.0

Ellipsoidal template

We first look for the largest ellipsoid contained in the square that is invariant under A.

model = Model(sdp_solver)
@variable(model, S, Ellipsoid(symmetric=true))
@constraint(model, S ⊆ □)
@constraint(model, A * S ⊆ S)
@objective(model, Max, nth_root(volume(S)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
maximal = value(S)

SetProg.Sets.Polar{Float64, SetProg.Sets.Ellipsoid{Float64}}(SetProg.Sets.Ellipsoid{Float64}([0.9999999634184679 0.0; 0.0 0.9999999634184679]))

We then look for the smallest invariant ellipsoid containing the square.

model = Model(sdp_solver)
@variable(model, S, Ellipsoid(symmetric=true))
@constraint(model, □ ⊆ S)
@constraint(model, A * S ⊆ S)
@objective(model, Min, nth_root(volume(S)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
minimal = value(S)

SetProg.Sets.Ellipsoid{Float64}([0.4999999856232981 2.2146125131708216e-20; 2.2146125131708216e-20 0.4999999856232981])

We plot the two ellipsoids together with the square.

using Plots
plot(ratio=:equal)
plot!(minimal)
plot!(□)
plot!(maximal)

We can also look for convex sublevel sets of a quartic polynomial. The basis keyword selects the Gram basis exposed by SumOfSquares. Here we use Chebyshev polynomials instead of the default monomials. SetProg keeps the defining polynomial homogeneous of degree four in either basis. Convexity and homogeneity ensure nonnegativity, so psd = false omits the redundant SOS constraint on the defining polynomial.

model = Model(sdp_solver)
@variable(
    model,
    S,
    PolySet(degree = 4, convex = true, symmetric = true, psd = false, basis = Chebyshev)
)
@constraint(model, S ⊆ □)
@constraint(model, A * S ⊆ S)
@objective(model, Max, nth_root(volume(S)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
maximal_convex = value(S)

SetProg.Sets.Polar{Float64, SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}}(SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}(4, GramMatrix with row/column basis:
 SubBasis{ChebyshevFirstKind}([1, x[2]^2, x[1]*x[2], x[1]^2])
And entries in a 4×4 MultivariateMoments.SymMatrix{Float64}:
 1.230557414768613      0.6152787073843071     …  0.6152787073843063
 0.6152787073843071     0.24999999517293853       0.3652787122113685
 6.98337310895222e-16  -4.290391410676423e-14     4.3144962744459745e-14
 0.6152787073843063     0.3652787122113685        0.24999999517293836, [11.999999768301054 -5.206426592995398e-13 -5.144002638630726e-13 8.485280921489254; -5.206426592995398e-13 5.656853759428914 2.8284265973685843 5.174335295366242e-13; -5.144002638630726e-13 2.8284265973685843 5.656853759428917 5.196024730456617e-13; 8.485280921489254 5.174335295366242e-13 5.196024730456617e-13 11.999999768301048]))

And the minimal containing the square:

model = Model(sdp_solver)
@variable(
    model,
    S,
    PolySet(degree = 4, convex = true, symmetric = true, psd = false, basis = Chebyshev)
)
@constraint(model, □ ⊆ S)
@constraint(model, A * S ⊆ S)
@objective(model, Min, nth_root(volume(S)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
minimal_convex = value(S)

SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.ChebyshevFirstKind, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}(4, GramMatrix with row/column basis:
 SubBasis{ChebyshevFirstKind}([1, x[2]^2, x[1]*x[2], x[1]^2])
And entries in a 4×4 MultivariateMoments.SymMatrix{Float64}:
 0.2627240620677149   0.13136203103385757     …   0.13136203103385746
 0.13136203103385757  0.07544417330775251         0.05591785772610506
 1.7575534687432e-16  2.3481484153169438e-11     -2.3481452384745062e-11
 0.13136203103385746  0.05591785772610506         0.07544417330775252, [3.6213203187721232 2.817764860884378e-10 2.8177832643899257e-10 1.4999999984432044; 2.817764860884378e-10 0.7928932190351994 0.08578643962719408 -2.8177786607934655e-10; 2.8177832643899257e-10 0.08578643962719408 0.7928932190351989 -2.817770748731098e-10; 1.4999999984432044 -2.8177786607934655e-10 -2.817770748731098e-10 3.6213203187721223])

We plot the two quartic sublevel sets together with the square.

plot(ratio = :equal)
plot!(minimal_convex)
plot!(□)
plot!(maximal_convex)

This page was generated using Literate.jl.

Discrete-time Invariant Set

Ellipsoidal template

Quartic template