From 177f97560f6a98d41e553c1f25bc9322980748e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 30 Sep 2026 17:57:21 +0200 Subject: [PATCH 1/6] [ADD] templates: setprog template behind an extension SetProgTemplate uses the sets of SetProg as node functions: ellipsoids, polynomial sets, and what the other templates do not have, piecewise semi-ellipsoids and piecewise polytopes. Stability only, like the sum-of-squares template. Needs SetProg with support for constant supersets, blegat/SetProg.jl branch bl/constant_superset. --- CLAUDE.md | 5 +- Project.toml | 3 + docs/src/manual/support.md | 8 +- docs/src/manual/templates.md | 8 +- docs/src/reference/templates.md | 10 + ext/PathCompleteCertificatesSetProgExt.jl | 236 ++++++++++++++++++++++ src/PathCompleteCertificates.jl | 3 + src/templates/set_prog.jl | 76 +++++++ test/Project.toml | 4 +- test/runtests.jl | 1 + test/templates/set_prog.jl | 205 +++++++++++++++++++ 11 files changed, 552 insertions(+), 7 deletions(-) create mode 100644 ext/PathCompleteCertificatesSetProgExt.jl create mode 100644 src/templates/set_prog.jl create mode 100644 test/templates/set_prog.jl diff --git a/CLAUDE.md b/CLAUDE.md index aed643b..14d0414 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -189,7 +189,8 @@ src/ │ ├── quadratic.jl │ ├── polyhedral.jl symmetric 2n-face, fixed facets │ ├── conic_polyhedral.jl free facets, plus the partition that linearises them -│ └── sum_of_squares.jl the type only -- its methods are in ext/ +│ ├── sum_of_squares.jl the type only -- its methods are in ext/ +│ └── set_prog.jl the type only -- its methods are in ext/ ├── problems/ AXIS 2 — what the edge inequality says │ ├── abstract.jl AbstractProblem, add_edge_constraint!, shared statuses, and │ │ what refine reads: best_certificate, objective, optimality_gap @@ -216,7 +217,7 @@ src/ | Path | What it is | | :--- | :--- | -| `ext/` | Optional interop, one extension per weak dependency. Today: the plotting recipe for `Trajectory` behind `RecipesBase`, and the sum-of-squares template behind `SumOfSquares`. **A template may be split this way** — type in `src/templates/`, primitives in `ext/` — when its stack is too heavy to load for an LP | +| `ext/` | Optional interop, one extension per weak dependency. Today: the plotting recipe for `Trajectory` behind `RecipesBase`, the sum-of-squares template behind `SumOfSquares`, and the SetProg template behind `SetProg`. **A template may be split this way** — type in `src/templates/`, primitives in `ext/` — when its stack is too heavy to load for an LP | | `test/` | Mirrors `src/` **including its subdirectories**. Entry point `test/runtests.jl`; every file is standalone-runnable | | `docs/` | The manual, the examples and these developer docs | | `docs/src/examples/` | Runnable scripts, **executed by the docs build** via Literate. Run one with `--project=docs`, which is where `Clarabel` and `Plots` live — never a dependency of the package or of the environment CI instantiates. There is no top-level `examples/`: it was folded in here so an example cannot go stale unnoticed | diff --git a/Project.toml b/Project.toml index bc5f891..0196fa1 100644 --- a/Project.toml +++ b/Project.toml @@ -18,14 +18,17 @@ MathematicalSystems = "0.14" Random = "1" MultivariatePolynomials = "0.5" RecipesBase = "1" +SetProg = "0.4" SumOfSquares = "0.7, 0.8" julia = "1.10" [weakdeps] MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +SetProg = "39881422-4b75-5582-a5c7-0feb14562a65" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" [extensions] PathCompleteCertificatesRecipesBaseExt = "RecipesBase" +PathCompleteCertificatesSetProgExt = "SetProg" PathCompleteCertificatesSumOfSquaresExt = ["SumOfSquares", "MultivariatePolynomials"] diff --git a/docs/src/manual/support.md b/docs/src/manual/support.md index c659da3..5fbbaad 100644 --- a/docs/src/manual/support.md +++ b/docs/src/manual/support.md @@ -15,11 +15,12 @@ three problems are still quadratic-only. | [`PolyhedralTemplate`](@ref) | ✓ | ✗ | ✗ | | [`ConicPolyhedralTemplate`](@ref) | ✓ | ✗ | ✗ | | [`SumOfSquaresTemplate`](@ref) | ✓ ² | ✗ | ✗ | +| [`SetProgTemplate`](@ref) | ✓ ² | ✗ | ✗ | ¹ Entrywise nonnegative mode matrices only. -² Behind a package extension: `using SumOfSquares` first, or its primitives are -not defined. +² Behind a package extension: `using SumOfSquares` first, or `using SetProg` +for [`SetProgTemplate`](@ref), or its primitives are not defined. A ✗ throws an `ArgumentError` before any model is built. Optimal control additionally requires a **complete** graph, not merely a path-complete one. @@ -44,8 +45,9 @@ tight: | [`PolyhedralTemplate`](@ref) | ✓ | map; relation on the duals when the node matrices coincide | | [`ConicPolyhedralTemplate`](@ref) | ✓ | map | | [`SumOfSquaresTemplate`](@ref) | ✗ ³ | map, linear program | +| [`SetProgTemplate`](@ref) | ✗ ³ | map | -³ No `domination_slack` yet, so the loop cannot read where that template is +³ No `domination_slack` yet, so the loop cannot read where these templates are tight. The second column is what the template's closures entitle diff --git a/docs/src/manual/templates.md b/docs/src/manual/templates.md index 2c9a556..b4706e8 100644 --- a/docs/src/manual/templates.md +++ b/docs/src/manual/templates.md @@ -13,7 +13,7 @@ member. data, [`PolyhedralTemplate`](@ref) holds one matrix per node, and a type has nowhere to put it. -## The six templates +## The seven templates | Template | ``V_s(x)`` | Degree | Solved for | Source | | :-- | :-- | :-: | :-- | :-- | @@ -23,6 +23,7 @@ member. | [`PolyhedralTemplate`](@ref) | ``\max_k \lvert (G_s x)_k \rvert / w_{s,k}`` | 1 | weights ``w_s`` | [athanasopoulos2019polyhedral](@cite) | | [`ConicPolyhedralTemplate`](@ref) | ``\max_i \lvert p_{s,i}^\top x \rvert`` | 1 | facets ``p_{s,i}`` | unpublished ¹ | | [`SumOfSquaresTemplate`](@ref) | ``z(x)^\top Q_s z(x)`` | ``2d`` | ``Q_s \succeq 0`` | [parrilo2008approximation](@cite) ² | +| [`SetProgTemplate`](@ref) | the gauge of a SetProg set | its degree | the set | [SetProg](https://github.com/blegat/SetProg.jl) ³ | ¹ The free-facet construction is not yet published; it is under submission. @@ -30,6 +31,11 @@ member. Raising `d` tightens the bound monotonically, and `d = 1` is exactly [`QuadraticTemplate`](@ref), a useful thing to check a change against. +³ Behind a package extension too. An ellipsoid is [`QuadraticTemplate`](@ref) +and a polynomial set [`SumOfSquaresTemplate`](@ref); what it adds are the +piecewise semi-ellipsoids, quadratic on each cone of a partition, which improve +on the quadratic template as the partition is refined. + The degree is [`rate_exponent`](@ref): ``V(cx) = c^d V(x)``. It is accounted for when rates are computed, so bounds from different templates are comparable. diff --git a/docs/src/reference/templates.md b/docs/src/reference/templates.md index 54130dc..061527e 100644 --- a/docs/src/reference/templates.md +++ b/docs/src/reference/templates.md @@ -58,3 +58,13 @@ required before this template can be solved with. Modules = [PathCompleteCertificates] Pages = ["templates/sum_of_squares.jl"] ``` + +## SetProg sets + +The primitives are defined in a package extension, so `using SetProg` is +required before this template can be solved with. + +```@autodocs +Modules = [PathCompleteCertificates] +Pages = ["templates/set_prog.jl"] +``` diff --git a/ext/PathCompleteCertificatesSetProgExt.jl b/ext/PathCompleteCertificatesSetProgExt.jl new file mode 100644 index 0000000..664b268 --- /dev/null +++ b/ext/PathCompleteCertificatesSetProgExt.jl @@ -0,0 +1,236 @@ +# The SetProg template, as a package extension. +# +# SetProg pulls SumOfSquares, Polyhedra and the polynomial stack, so it stays +# weak for the same reason as the sum-of-squares template. What it adds over +# that one is the sets SetProg already knows how to constrain, the piecewise +# semi-ellipsoids and the piecewise polytopes in particular. +# +# A node function is the gauge-like function of a set, so every primitive is a +# set inclusion. All of them are modeled in SetProg's primal space, where a set +# is represented by that function: the normalization `S ⊆ unit` needs it (a +# constant polynomial or piecewise superset has no dual representation) and +# the edge inclusions accept either space. +module PathCompleteCertificatesSetProgExt + +import JuMP +import LinearAlgebra +import SetProg +import PathCompleteCertificates as PCC + +const Sets = SetProg.Sets +const MP = SetProg.MP + +function _variable(template::PCC.SetProgTemplate, node::Integer) + template.shared && return only(template.variables) + + node <= length(template.variables) || throw( + ArgumentError( + "the template carries $(length(template.variables)) set variables " * + "but the graph has at least $node nodes", + ), + ) + + return template.variables[node] +end + +_degree(::SetProg.Ellipsoid) = 2 +_degree(variable::SetProg.PolySet) = variable.degree +_degree(::SetProg.Polytope) = 1 + +function PCC.rate_exponent(template::PCC.SetProgTemplate) + degrees = unique(_degree.(template.variables)) + + # One rate for the whole graph needs one degree: the edge condition scales + # `V_src` by gamma^d and `V_dst` would need a different power otherwise. + length(degrees) == 1 || + throw(ArgumentError("the set variables of every node must have the same degree")) + + return only(degrees) +end + +# The variable with its dimension filled in, so that the normalization, which +# is not given the dimension, can build the unit set. +function _with_dimension(variable::SetProg.Ellipsoid, dimension) + variable.point === nothing && variable.superset === nothing || + throw(ArgumentError("`point` and `superset` are not supported by SetProgTemplate")) + variable.symmetric || + throw(ArgumentError("SetProgTemplate needs `Ellipsoid(symmetric = true)`")) + _check_dimension(variable.piecewise, dimension) + + return SetProg.Ellipsoid(; + symmetric = true, + dimension = dimension, + piecewise = variable.piecewise, + ) +end + +function _with_dimension(variable::SetProg.Polytope, dimension) + variable.symmetric && variable.piecewise !== nothing || throw( + ArgumentError( + "SetProgTemplate needs `Polytope(symmetric = true, piecewise = P)`: " * + "without a partition the function is linear, not a gauge", + ), + ) + _check_dimension(variable.piecewise, dimension) + + return SetProg.Polytope(; + symmetric = true, + dimension = dimension, + piecewise = variable.piecewise, + ) +end + +function _with_dimension(variable::SetProg.PolySet, dimension) + variable.point === nothing && variable.superset === nothing || + throw(ArgumentError("`point` and `superset` are not supported by SetProgTemplate")) + variable.symmetric || + throw(ArgumentError("SetProgTemplate needs `PolySet(symmetric = true)`")) + + # The unit set is a polynomial, and it has to be written in the variables of + # the set variable for the two to be compared. + variable.variables === nothing && throw( + ArgumentError( + "SetProgTemplate needs the polynomial variables, " * + "use `PolySet(variables = x, ...)`", + ), + ) + length(variable.variables) == dimension || throw( + ArgumentError( + "the set variable carries $(length(variable.variables)) polynomial " * + "variables but the system has dimension $dimension", + ), + ) + + return SetProg.PolySet(; + symmetric = true, + degree = variable.degree, + dimension = dimension, + convex = variable.convex, + variables = variable.variables, + basis = variable.basis, + ) +end + +_check_dimension(::Nothing, _) = nothing + +function _check_dimension(partition, dimension) + SetProg.Polyhedra.fulldim(partition) == dimension || throw( + ArgumentError( + "the partition lives in R^$(SetProg.Polyhedra.fulldim(partition)) " * + "but the system is $dimension-dimensional", + ), + ) + + return nothing +end + +function PCC.add_function_variables!( + model::JuMP.Model, + template::PCC.SetProgTemplate, + dimension::Integer, + node::Integer, +) + @assert dimension > 0 + variable = _with_dimension(_variable(template, node), dimension) + + return JuMP.@variable(model, variable_type = variable, base_name = "S_$(node)") +end + +_identity(dimension) = + LinearAlgebra.Symmetric(Matrix(1.0LinearAlgebra.I, dimension, dimension)) + +# The 1-sublevel set of `‖x‖₂^d`, or of the gauge of the partition for a +# polytope, in the representation of the set variable. +function _unit(variable::SetProg.Ellipsoid) + ball = Sets.Ellipsoid(_identity(variable.dimension)) + variable.piecewise === nothing && return ball + + partition = variable.piecewise + + return Sets.Piecewise( + [ball for _ in SetProg.Polyhedra.halfspaces(partition)], + partition, + ) +end + +function _unit(variable::SetProg.Polytope) + partition = variable.piecewise + facets = [Sets.PolarPoint(h.a / h.β) for h in SetProg.Polyhedra.halfspaces(partition)] + + return Sets.Piecewise(facets, partition) +end + +function _unit(variable::SetProg.PolySet) + half_degree = div(variable.degree, 2) + monomials = MP.monomials(variable.variables, half_degree) + + # (x₁² + ⋯ + xₙ²)ᵈ = Σ_{|α| = d} d! / α! x^(2α), a diagonal Gram matrix. + coefficients = [ + factorial(half_degree) / prod(factorial, MP.exponents(monomial)) for + monomial in monomials + ] + gram = SetProg.GramMatrix(Matrix(LinearAlgebra.Diagonal(coefficients)), monomials) + + return Sets.PolySet(variable.degree, gram) +end + +function PCC.add_normalization!( + model::JuMP.Model, + ::PCC.SetProgTemplate, + S::SetProg.SetVariableRef, +) + # `S ⊆ unit` is `V ≥ ‖x‖₂^d`. Any bounded set containing the origin in its + # interior can be scaled into it, so this excludes `V ≡ 0` and nothing else. + JuMP.@constraint(model, S ⊆ _unit(S.set)) + + return nothing +end + +function PCC.add_domination!( + model::JuMP.Model, + template::PCC.SetProgTemplate, + S_src::SetProg.SetVariableRef, + S_dst::SetProg.SetVariableRef, + map::AbstractMatrix; + scale = 1, + margin = 0, +) + margin == 0 || + throw(ArgumentError("SetProgTemplate does not support a positive margin")) + + # The scale is folded into the map, `V_dst(A x) ≤ γᵈ V_src(x)` being + # `(A / γ) S_src ⊆ S_dst`, and a fractional root of a JuMP expression is + # not expressible. + scale isa Real || + throw(ArgumentError("SetProgTemplate needs a constant scale, got $scale")) + scale > 0 || throw(ArgumentError("the scale must be positive, got $scale")) + + gamma = scale^(1 / PCC.rate_exponent(template)) + + JuMP.@constraint(model, (map ./ gamma) * S_src ⊆ S_dst) + + return nothing +end + +PCC.solution_value(::PCC.SetProgTemplate, S::SetProg.SetVariableRef) = + PCC.SetProgFunction(JuMP.value(S)) + +(V::PCC.SetProgFunction)(x::AbstractVector{<:Real}) = _evaluate(V.set, x) + +_evaluate(set::Sets.Ellipsoid, x) = LinearAlgebra.dot(x, set.Q * x) + +_evaluate(set::Union{Sets.PolySet, Sets.ConvexPolySet}, x) = + MP.polynomial(set.p)(Sets.space_variables(set) => collect(x)) + +_evaluate(set::Sets.PolarPoint, x) = LinearAlgebra.dot(set.a, x) + +function _evaluate(set::Sets.Piecewise, x) + # The function is continuous, so a point on the boundary between two pieces + # can be evaluated in either of them. + piece = findfirst(piece -> x in piece, set.pieces) + piece === nothing && throw(ArgumentError("$x is in none of the pieces")) + + return _evaluate(set.sets[piece], x) +end + +end diff --git a/src/PathCompleteCertificates.jl b/src/PathCompleteCertificates.jl index b5dca79..08141e3 100644 --- a/src/PathCompleteCertificates.jl +++ b/src/PathCompleteCertificates.jl @@ -58,6 +58,9 @@ include("templates/conic_polyhedral.jl") # The methods live in ext/PathCompleteCertificatesSumOfSquaresExt.jl: the type # is cheap, the polynomial stack behind it is not. include("templates/sum_of_squares.jl") +# Same split, behind `SetProg`: its sets reach templates the others do not have, +# piecewise semi-ellipsoids among them. +include("templates/set_prog.jl") # --- Axis 2: one file per problem, each composing those primitives ------------ include("problems/stability.jl") diff --git a/src/templates/set_prog.jl b/src/templates/set_prog.jl new file mode 100644 index 0000000..13dd68a --- /dev/null +++ b/src/templates/set_prog.jl @@ -0,0 +1,76 @@ +""" + SetProgTemplate(variable) + SetProgTemplate(variables::AbstractVector) + +The template of the sets of [SetProg](https://github.com/blegat/SetProg.jl): +node ``s`` carries a set ``S_s`` and its node function is the gauge-like +function ``V_s`` whose 1-sublevel set is ``S_s``. For homogeneous functions of +degree ``d``, the edge inequality ``V_{dst}(A x) \\le \\gamma^d V_{src}(x)`` is +the set inclusion ``(A / \\gamma) S_{src} \\subseteq S_{dst}``, which is what +this template asks SetProg for. + +`variable` is a SetProg set variable, one of + +| `variable` | ``V_s(x)`` | Degree | +| :-- | :-- | :-: | +| `Ellipsoid(symmetric = true)` | ``x^\\top Q_s x`` | 2 | +| `PolySet(symmetric = true, degree = 2d, variables = x)` | a polynomial ``p_s(x)`` | ``2d`` | +| `Ellipsoid(symmetric = true, piecewise = P)` | ``x^\\top Q_{s,i} x`` on the ``i``th cone | 2 | +| `Polytope(symmetric = true, piecewise = P)` | ``a_{s,i}^\\top x`` on the ``i``th cone | 1 | + +The piecewise ones are what this template adds to the others: a piecewise +semi-ellipsoid, or a polytope, whose function is quadratic, or linear, on each +cone over a facet of the polytope `P` (which must contain the origin in its +interior). Pass a vector with one variable per node to give each node its own +partition. + +```julia +using SetProg # before this is usable +template = SetProgTemplate(Ellipsoid(symmetric = true, piecewise = P)) +``` + +The node functions are normalized by ``V_s(x) \\ge \\lVert x \\rVert_2^d``, +where ``d`` is the degree, or by ``V_s(x) \\ge`` the gauge of `P` for a +polytope. The edge conditions are homogeneous, so this only fixes the scale. + +!!! note "The methods live in a package extension" + Everything this template does is behind `SetProg`, which pulls the + polynomial and polyhedral stacks. Without `using SetProg` the type exists + but its primitives are not defined. + +!!! warning "Use a partition with exact facets" + Rounding errors in the facet normals of `P`, such as the `1e-16` of a + `cos(π / 2)`, can make the solver report an almost optimal solution that + violates the edge inequalities. Round the normals, or use exact numbers. + +[`domination_slack`](@ref) is not implemented for this template, so neither +are [`edge_slacks`](@ref) and what builds on it. +""" +struct SetProgTemplate{V <: AbstractVector} <: AbstractTemplate + variables::V + shared::Bool + + SetProgTemplate(variables::AbstractVector) = new{typeof(variables)}(variables, false) + SetProgTemplate(variable) = new{Vector{typeof(variable)}}([variable], true) +end + +""" + SetProgFunction(set) + +A fitted [`SetProgTemplate`](@ref) node function: the solved SetProg set, +callable as ``V(x)``, the function whose 1-sublevel set is `set`. +""" +struct SetProgFunction{S} + set::S +end + +# The sets are constructed in their cone (PSD matrices, SOS polynomials), and +# the normalization keeps them positive, so there is nothing to add. +add_nonnegativity!(::JuMP.Model, ::SetProgTemplate, _) = nothing + +node_value( + ::SetProgTemplate, + ::AbstractProblem, + V::SetProgFunction, + x::AbstractVector{<:Real}, +) = V(x) diff --git a/test/Project.toml b/test/Project.toml index 8bcefdf..c9031c6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,15 +1,17 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" +DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathematicalSystems = "d14a8603-c872-5ed3-9ece-53e0e82e39da" PathCompleteCertificates = "b391c707-6ba7-482f-97f4-86d517cfbcc9" -DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SetProg = "39881422-4b75-5582-a5c7-0feb14562a65" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [sources] PathCompleteCertificates = {path = ".."} +SetProg = {rev = "bl/constant_superset", url = "https://github.com/blegat/SetProg.jl"} diff --git a/test/runtests.jl b/test/runtests.jl index cfc7e23..5f30f3b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,6 +28,7 @@ const TEST_FILES = [ ("./templates/polyhedral.jl", :slow), ("./templates/conic_polyhedral.jl", :slow), ("./templates/sum_of_squares.jl", :slow), + ("./templates/set_prog.jl", :slow), ("./problems/certify.jl", :slow), ("./problems/tight_edges.jl", :slow), ("./problems/stability.jl", :slow), diff --git a/test/templates/set_prog.jl b/test/templates/set_prog.jl new file mode 100644 index 0000000..f4ed4ea --- /dev/null +++ b/test/templates/set_prog.jl @@ -0,0 +1,205 @@ +module TestSetProg + +# The SetProg template lives in a package extension, so this file is the only +# one that loads SetProg. +# +# The anchors are the templates it overlaps with: an ellipsoid is +# `QuadraticTemplate`, a polynomial set is `SumOfSquaresTemplate`, and a +# polytope that is linear on the cones over the facets of a partition is +# `ConicPolyhedralTemplate` on the same cones. What it adds, the piecewise +# semi-ellipsoids, is checked against these and against the edge inequality. + +using Test +using LinearAlgebra +using HybridSystems +using SetProg +using SumOfSquares, DynamicPolynomials +import PathCompleteCertificates as PCC +import Clarabel + +const OPTIMIZER = Clarabel.Optimizer + +# A rotation and a shear: no quadratic attains the joint spectral radius. +const A = [0.69 * [0.0 1.0; -1.0 0.0], 0.69 * [1.0 1.0; 0.0 1.0]] +const PROBLEM = PCC.StabilityProblem(PCC.switched_system(A)) +# Every product of modes bounds the joint spectral radius from below. +const LOWER_BOUND = maximum(Ai -> maximum(abs, eigvals(Ai)), A) + +const MEMORYLESS = GraphAutomaton(1) +add_transition!(MEMORYLESS, 1, 1, 1) +add_transition!(MEMORYLESS, 1, 1, 2) + +# The Polyhedra library needs an LP solver to detect the linearity of the pieces +const LIB = Polyhedra.DefaultLibrary{Float64}( + JuMP.optimizer_with_attributes(Clarabel.Optimizer, JuMP.MOI.Silent() => true), +) + +# Partitions with exact facets, see the warning in the docstring. The cones over +# the facets of the octagon are the cones of `planar_conic_partition(1)`: its +# vertices are on the rays ±(1, 0), ±(1, 1), ±(0, 1) and ±(-1, 1). +const SQUARE = polyhedron( + hrep([ + 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), + ]), + LIB, +) +const OCTAGON = polyhedron( + hrep([ + HalfSpace(a, 6.0) for a in [ + [2.0, 1.0], + [1.0, 2.0], + [-1.0, 2.0], + [-2.0, 1.0], + [-2.0, -1.0], + [-1.0, -2.0], + [1.0, -2.0], + [2.0, -1.0], + ] + ]), + LIB, +) + +@polyvar x[1:2] + +rate(template, graph = MEMORYLESS) = + PCC.jsr_bound(template, graph, PROBLEM; optimizer = OPTIMIZER, rtol = 1e-4) + +# The edge inequality `V_dst(A x) ≤ γᵈ V_src(x)`, sampled on the circle. +function test_edges(certificate, template, graph) + V = PCC.functions(certificate) + scale = certificate.rate^PCC.rate_exponent(template) + + for edge in PCC.edges(graph), θ in range(0, 2π; length = 97) + point = [cos(θ), sin(θ)] + src = V[PCC.source(edge)] + dst = V[PCC.dest(edge)] + @test scale * src(point) - dst(A[PCC.label(graph, edge)] * point) >= -1e-5 + end +end + +@testset "the extension is what supplies the methods" begin + @test Base.get_extension(PCC, :PathCompleteCertificatesSetProgExt) !== nothing +end + +@testset "the degree is the one of the set variable" begin + @test PCC.rate_exponent(PCC.SetProgTemplate(Ellipsoid(symmetric = true))) == 2 + @test PCC.rate_exponent( + PCC.SetProgTemplate(PolySet(symmetric = true, degree = 6, variables = x)), + ) == 6 + @test PCC.rate_exponent( + PCC.SetProgTemplate(Polytope(symmetric = true, piecewise = SQUARE)), + ) == 1 + mixed = PCC.SetProgTemplate([ + Ellipsoid(symmetric = true), + PolySet(symmetric = true, degree = 4, variables = x), + ]) + @test_throws ArgumentError PCC.rate_exponent(mixed) +end + +@testset "an ellipsoid reproduces the quadratic template" begin + quadratic = rate(PCC.QuadraticTemplate()) + ellipsoid = rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true))) + + @test PCC.is_feasible(ellipsoid) + @test abs(ellipsoid.rate - quadratic.rate) <= 1e-3 + + V = only(PCC.functions(ellipsoid)) + @test V isa PCC.SetProgFunction + # V(cx) = c² V(x), and the normalization `V ≥ ‖x‖²` + @test V([2.0, 3.0]) ≈ 4 * V([1.0, 1.5]) rtol = 1e-6 + @test V([1.0, 1.0]) >= 2 - 1e-6 + test_edges(ellipsoid, PCC.SetProgTemplate(Ellipsoid(symmetric = true)), MEMORYLESS) +end + +@testset "a polynomial set reproduces the sum-of-squares template" begin + sos = rate(PCC.SumOfSquaresTemplate(2, x)) + template = PCC.SetProgTemplate(PolySet(symmetric = true, degree = 4, variables = x)) + polyset = rate(template) + + @test PCC.is_feasible(polyset) + @test abs(polyset.rate - sos.rate) <= 1e-3 + test_edges(polyset, template, MEMORYLESS) +end + +@testset "a piecewise polytope reproduces the conic polyhedral template" begin + conic = rate(PCC.ConicPolyhedralTemplate([PCC.planar_conic_partition(1)])) + template = PCC.SetProgTemplate(Polytope(symmetric = true, piecewise = OCTAGON)) + polytope = rate(template) + + @test PCC.is_feasible(polytope) + @test abs(polytope.rate - conic.rate) <= 1e-3 + test_edges(polytope, template, MEMORYLESS) +end + +@testset "piecewise semi-ellipsoids" begin + ellipsoid = rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true))).rate + + @testset "$name" for (name, partition) in [("square", SQUARE), ("octagon", OCTAGON)] + template = PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = partition)) + certificate = rate(template) + + @test PCC.is_feasible(certificate) + # An ellipsoid is a piecewise semi-ellipsoid whose pieces are equal + @test certificate.rate <= ellipsoid + 1e-4 + @test certificate.rate >= LOWER_BOUND - 1e-4 + + V = only(PCC.functions(certificate)) + @test V.set isa SetProg.Sets.Piecewise + test_edges(certificate, template, MEMORYLESS) + end + + @testset "more pieces, lower rate" begin + octagon = + rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = OCTAGON))) + @test octagon.rate < ellipsoid - 1e-2 + end + + @testset "one partition per node" begin + graph = PCC.de_bruijn(1, 2) + template = PCC.SetProgTemplate([ + Ellipsoid(symmetric = true, piecewise = SQUARE), + Ellipsoid(symmetric = true, piecewise = OCTAGON), + ]) + certificate = rate(template, graph) + + @test PCC.is_feasible(certificate) + @test certificate.rate >= LOWER_BOUND - 1e-4 + test_edges(certificate, template, graph) + end +end + +@testset "what the template does not support" begin + @test_throws ArgumentError rate(PCC.SetProgTemplate(Ellipsoid(dimension = 2))) + @test_throws ArgumentError rate( + PCC.SetProgTemplate(PolySet(symmetric = true, degree = 4)), + ) + @test_throws ArgumentError rate( + PCC.SetProgTemplate(Polytope(symmetric = true, dimension = 2)), + ) + @polyvar y[1:3] + @test_throws ArgumentError rate( + PCC.SetProgTemplate(PolySet(symmetric = true, degree = 4, variables = y)), + ) + @test_throws ArgumentError rate( + PCC.SetProgTemplate([Ellipsoid(symmetric = true)]), + PCC.de_bruijn(1, 2), + ) + + model = JuMP.Model() + template = PCC.SetProgTemplate(Ellipsoid(symmetric = true)) + S = PCC.add_function_variables!(model, template, 2, 1) + @test_throws ArgumentError PCC.add_domination!(model, template, S, S, A[1]; margin = 1) + @test_throws ArgumentError PCC.add_domination!( + model, + template, + S, + S, + A[1]; + scale = JuMP.@variable(model), + ) +end + +end From 79b6ac5a8d609d4f310f01ac94e5adc1c1366bfa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 8 Oct 2026 12:02:44 +0200 Subject: [PATCH 2/6] [ADD] templates: reindex the setprog template on a lift --- src/templates/set_prog.jl | 6 ++++++ test/templates/set_prog.jl | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/templates/set_prog.jl b/src/templates/set_prog.jl index 13dd68a..4461c4b 100644 --- a/src/templates/set_prog.jl +++ b/src/templates/set_prog.jl @@ -64,6 +64,12 @@ struct SetProgFunction{S} set::S end +# A shared variable follows any lift. One per node gives each new node the +# variable of its origin. +reindex(template::SetProgTemplate, origins::AbstractVector) = + template.shared ? template : + SetProgTemplate(_reindexed_data(template.variables, origins, "set variables")) + # The sets are constructed in their cone (PSD matrices, SOS polynomials), and # the normalization keeps them positive, so there is nothing to add. add_nonnegativity!(::JuMP.Model, ::SetProgTemplate, _) = nothing diff --git a/test/templates/set_prog.jl b/test/templates/set_prog.jl index f4ed4ea..21621be 100644 --- a/test/templates/set_prog.jl +++ b/test/templates/set_prog.jl @@ -171,6 +171,27 @@ end end end +@testset "the template follows a lift" begin + shared = PCC.SetProgTemplate(Ellipsoid(symmetric = true)) + @test PCC.reindex(shared, [1, 1, 1]) === shared + + square = Ellipsoid(symmetric = true, piecewise = SQUARE) + octagon = Ellipsoid(symmetric = true, piecewise = OCTAGON) + per_node = PCC.SetProgTemplate([square, octagon]) + + # Each new node gets the variable of its origin. + lifted = PCC.reindex(per_node, [2, 1, 2]) + @test lifted.variables[1] === octagon + @test lifted.variables[2] === square + @test lifted.variables[3] === octagon + @test !lifted.shared + + # Merging nodes is fine when they carry the same variable, not otherwise. + @test only(PCC.reindex(PCC.SetProgTemplate([square, square]), [[1, 2]]).variables) === + square + @test_throws ArgumentError PCC.reindex(per_node, [[1, 2]]) +end + @testset "what the template does not support" begin @test_throws ArgumentError rate(PCC.SetProgTemplate(Ellipsoid(dimension = 2))) @test_throws ArgumentError rate( From 7c04df7c895c4089fd2f03987cd92a172f6c8ab2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 8 Oct 2026 12:10:14 +0200 Subject: [PATCH 3/6] [IMP] meta: require setprog 0.4.4 --- Project.toml | 2 +- test/Project.toml | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0196fa1..8e52db2 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,7 @@ MathematicalSystems = "0.14" Random = "1" MultivariatePolynomials = "0.5" RecipesBase = "1" -SetProg = "0.4" +SetProg = "0.4.4" SumOfSquares = "0.7, 0.8" julia = "1.10" diff --git a/test/Project.toml b/test/Project.toml index c9031c6..c2856de 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -14,4 +14,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [sources] PathCompleteCertificates = {path = ".."} -SetProg = {rev = "bl/constant_superset", url = "https://github.com/blegat/SetProg.jl"} From dd84b592a16c9a499414715acc087362c842538b Mon Sep 17 00:00:00 2001 From: JulienCalbert Date: Mon, 12 Oct 2026 09:54:37 +0200 Subject: [PATCH 4/6] IMP test: setprog bounds shared, coarser, partition on HiGHS --- src/templates/catalogue/set_prog.jl | 4 +++ test/Project.toml | 1 + test/templates/set_prog.jl | 41 ++++++++++++++++++----------- 3 files changed, 30 insertions(+), 16 deletions(-) diff --git a/src/templates/catalogue/set_prog.jl b/src/templates/catalogue/set_prog.jl index f1774b0..8b381c3 100644 --- a/src/templates/catalogue/set_prog.jl +++ b/src/templates/catalogue/set_prog.jl @@ -43,6 +43,10 @@ polytope. The edge conditions are homogeneous, so this only fixes the scale. `cos(π / 2)`, can make the solver report an almost optimal solution that violates the edge inequalities. Round the normals, or use exact numbers. +Build `P` with a Polyhedra library whose solver is a simplex one, such as +HiGHS: SetProg runs many small linear programs on the partition at every +solve, and an interior-point solver makes each of them slow. + [`domination_slack`](@ref) is not implemented for this template, so neither are [`edge_slacks`](@ref) and what builds on it. """ diff --git a/test/Project.toml b/test/Project.toml index c2856de..edf0683 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/test/templates/set_prog.jl b/test/templates/set_prog.jl index 846062f..8a6c658 100644 --- a/test/templates/set_prog.jl +++ b/test/templates/set_prog.jl @@ -16,6 +16,7 @@ using SetProg using SumOfSquares, DynamicPolynomials import PathCompleteCertificates as PCC import Clarabel +import HiGHS const OPTIMIZER = Clarabel.Optimizer @@ -29,9 +30,11 @@ const MEMORYLESS = GraphAutomaton(1) add_transition!(MEMORYLESS, 1, 1, 1) add_transition!(MEMORYLESS, 1, 1, 2) -# The Polyhedra library needs an LP solver to detect the linearity of the pieces +# The Polyhedra library needs an LP solver to detect the linearity of the +# pieces, and SetProg runs many small programs on it at every solve: a simplex +# solver keeps them cheap, where Clarabel would triple the time of each solve. const LIB = Polyhedra.DefaultLibrary{Float64}( - JuMP.optimizer_with_attributes(Clarabel.Optimizer, JuMP.MOI.Silent() => true), + JuMP.optimizer_with_attributes(HiGHS.Optimizer, JuMP.MOI.Silent() => true), ) # Partitions with exact facets, see the warning in the docstring. The cones over @@ -64,8 +67,10 @@ const OCTAGON = polyhedron( @polyvar x[1:2] +# Every bound is bisected to `rtol`, so two of them agree to about twice that. rate(template, graph = MEMORYLESS) = - PCC.jsr_bound(template, graph, PROBLEM; optimizer = OPTIMIZER, rtol = 1e-4) + PCC.jsr_bound(template, graph, PROBLEM; optimizer = OPTIMIZER, rtol = 1e-3) +const AGREE = 2e-3 # The edge inequality `V_dst(A x) ≤ γᵈ V_src(x)`, sampled on the circle. function test_edges(certificate, template, graph) @@ -99,19 +104,21 @@ end @test_throws ArgumentError PCC.homogeneity_degree(mixed) end +# The one-node ellipsoid bound, shared by the tests below. +const ELLIPSOID = rate(PCC.SetProgTemplate(Ellipsoid(; symmetric = true))) + @testset "an ellipsoid reproduces the quadratic template" begin quadratic = rate(PCC.QuadraticTemplate()) - ellipsoid = rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true))) - @test PCC.is_feasible(ellipsoid) - @test abs(ellipsoid.rate - quadratic.rate) <= 1e-3 + @test PCC.is_feasible(ELLIPSOID) + @test abs(ELLIPSOID.rate - quadratic.rate) <= AGREE - V = only(PCC.functions(ellipsoid)) + V = only(PCC.functions(ELLIPSOID)) @test V isa PCC.SetProgFunction # V(cx) = c² V(x), and the normalization `V ≥ ‖x‖²` @test V([2.0, 3.0]) ≈ 4 * V([1.0, 1.5]) rtol = 1e-6 @test V([1.0, 1.0]) >= 2 - 1e-6 - test_edges(ellipsoid, PCC.SetProgTemplate(Ellipsoid(symmetric = true)), MEMORYLESS) + test_edges(ELLIPSOID, PCC.SetProgTemplate(Ellipsoid(symmetric = true)), MEMORYLESS) end @testset "a polynomial set reproduces the sum-of-squares template" begin @@ -120,7 +127,7 @@ end polyset = rate(template) @test PCC.is_feasible(polyset) - @test abs(polyset.rate - sos.rate) <= 1e-3 + @test abs(polyset.rate - sos.rate) <= AGREE test_edges(polyset, template, MEMORYLESS) end @@ -130,20 +137,24 @@ end polytope = rate(template) @test PCC.is_feasible(polytope) - @test abs(polytope.rate - conic.rate) <= 1e-3 + @test abs(polytope.rate - conic.rate) <= AGREE test_edges(polytope, template, MEMORYLESS) end @testset "piecewise semi-ellipsoids" begin - ellipsoid = rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true))).rate + piecewise = Dict( + name => rate( + PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = partition)), + ) for (name, partition) in [("square", SQUARE), ("octagon", OCTAGON)] + ) @testset "$name" for (name, partition) in [("square", SQUARE), ("octagon", OCTAGON)] template = PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = partition)) - certificate = rate(template) + certificate = piecewise[name] @test PCC.is_feasible(certificate) # An ellipsoid is a piecewise semi-ellipsoid whose pieces are equal - @test certificate.rate <= ellipsoid + 1e-4 + @test certificate.rate <= ELLIPSOID.rate + AGREE @test certificate.rate >= LOWER_BOUND - 1e-4 V = only(PCC.functions(certificate)) @@ -152,9 +163,7 @@ end end @testset "more pieces, lower rate" begin - octagon = - rate(PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = OCTAGON))) - @test octagon.rate < ellipsoid - 1e-2 + @test piecewise["octagon"].rate < ELLIPSOID.rate - 1e-2 end @testset "one partition per node" begin From 48108c94d639a2aec37d7183c1ffb4aad9bee2ce Mon Sep 17 00:00:00 2001 From: JulienCalbert Date: Mon, 12 Oct 2026 09:55:04 +0200 Subject: [PATCH 5/6] ADD templates: setprog closed under addition on one variable --- docs/src/manual/support.md | 2 +- docs/src/manual/templates.md | 2 +- src/templates/catalogue/set_prog.jl | 6 ++++++ test/templates/set_prog.jl | 16 ++++++++++++++++ 4 files changed, 24 insertions(+), 2 deletions(-) diff --git a/docs/src/manual/support.md b/docs/src/manual/support.md index fe48a04..5fa6adf 100644 --- a/docs/src/manual/support.md +++ b/docs/src/manual/support.md @@ -43,7 +43,7 @@ how the loop reads where a certificate is tight: | [`PolyhedralTemplate`](@ref) | ✓ | identity; max when the node matrices coincide | | [`ConicPolyhedralTemplate`](@ref) | ✓ | identity | | [`SumOfSquaresTemplate`](@ref) | ✗ ³ | identity, sum; composition if invertible | -| [`SetProgTemplate`](@ref) | ✗ ³ | identity | +| [`SetProgTemplate`](@ref) | ✗ ³ | identity; sum when the nodes share the set variable | ³ No `domination_slack` yet, so the loop cannot read where these templates are tight. diff --git a/docs/src/manual/templates.md b/docs/src/manual/templates.md index 2080887..36a3cd2 100644 --- a/docs/src/manual/templates.md +++ b/docs/src/manual/templates.md @@ -76,7 +76,7 @@ hold for it. A template declares its closures with [`is_closed_under`](@ref): | [`DualCopositiveTemplate`](@ref) | | ✓ | | | | [`PolyhedralTemplate`](@ref) | | ✓ if the node matrices coincide | | | | [`ConicPolyhedralTemplate`](@ref) | | | | | -| [`SetProgTemplate`](@ref) | | | | | +| [`SetProgTemplate`](@ref) | ✓ if the nodes share the set variable | | | | ¹ Not closed as a set of functions, yet the min lift is valid for it on positive linear systems ([debauche2024thesis](@cite), Thm. 7.43), the one diff --git a/src/templates/catalogue/set_prog.jl b/src/templates/catalogue/set_prog.jl index 8b381c3..c034cd3 100644 --- a/src/templates/catalogue/set_prog.jl +++ b/src/templates/catalogue/set_prog.jl @@ -77,3 +77,9 @@ reindex(template::SetProgTemplate, origins::AbstractVector) = # The sets are constructed in their cone (PSD matrices, SOS polynomials), and # the normalization keeps them positive, so there is nothing to add. add_nonnegativity!(::JuMP.Model, ::SetProgTemplate, _) = nothing + +# The sum of two functions of one family, on one partition, is in the family +# again: closed under addition exactly when every node carries the same set +# variable, as the fixed-facet polyhedral template is under maximum. +is_closed_under(template::SetProgTemplate, ::Addition, system) = + template.shared || allequal(template.variables) diff --git a/test/templates/set_prog.jl b/test/templates/set_prog.jl index 8a6c658..79ffa3e 100644 --- a/test/templates/set_prog.jl +++ b/test/templates/set_prog.jl @@ -201,6 +201,22 @@ end @test_throws ArgumentError PCC.reindex(per_node, [[1, 2]]) end +@testset "closed under addition when the nodes share the variable" begin + system = PROBLEM.system + shared = PCC.SetProgTemplate(Ellipsoid(symmetric = true)) + @test PCC.is_closed_under(shared, PCC.Addition(), system) + @test PCC.is_valid(PCC.SumLift(2), shared, system) + + octagon = Ellipsoid(symmetric = true, piecewise = OCTAGON) + same = PCC.SetProgTemplate([octagon, octagon]) + @test PCC.is_closed_under(same, PCC.Addition(), system) + + mixed = PCC.SetProgTemplate([Ellipsoid(symmetric = true, piecewise = SQUARE), octagon]) + @test !PCC.is_closed_under(mixed, PCC.Addition(), system) + @test !PCC.is_valid(PCC.SumLift(2), mixed, system) + @test !PCC.is_closed_under(shared, PCC.Maximum(), system) +end + @testset "what the template does not support" begin @test_throws ArgumentError rate(PCC.SetProgTemplate(Ellipsoid(dimension = 2))) @test_throws ArgumentError rate( From 7fafe82d349c9502d26a636d27e4fc673376e1b2 Mon Sep 17 00:00:00 2001 From: JulienCalbert Date: Mon, 12 Oct 2026 10:20:41 +0200 Subject: [PATCH 6/6] ADD docs: the setprog example, four sets on one node --- docs/Project.toml | 9 +- docs/make.jl | 4 +- docs/src/examples/set_prog.jl | 154 ++++++++++++++++++++++++++++++++++ docs/src/index.md | 4 +- docs/src/manual/templates.md | 3 +- 5 files changed, 167 insertions(+), 7 deletions(-) create mode 100644 docs/src/examples/set_prog.jl diff --git a/docs/Project.toml b/docs/Project.toml index 58bd82f..5a59971 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,24 +1,27 @@ # The examples are executed as part of the docs build, so this environment needs -# a solver, a plotting backend and the polynomial stack on top of the -# documentation tooling. None of them is a dependency of the package, nor of the -# environment CI instantiates for the tests. +# a solver, a plotting backend and the polynomial and set-programming stacks on +# top of the documentation tooling. None of them is a dependency of the package, +# nor of the environment CI instantiates for the tests. [deps] Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" PathCompleteCertificates = "b391c707-6ba7-482f-97f4-86d517cfbcc9" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +SetProg = "39881422-4b75-5582-a5c7-0feb14562a65" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" [compat] Documenter = "1" DocumenterCitations = "1" Literate = "2" +SetProg = "0.4.4" SumOfSquares = "0.7, 0.8" [sources] diff --git a/docs/make.jl b/docs/make.jl index 11d11dc..c733fa8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -21,6 +21,7 @@ const ORDER = [ "switched_systems", "stability", "templates_and_graphs", + "set_prog", "safety", "optimal_control", "comparing", @@ -146,8 +147,7 @@ makedocs(; style = :authoryear, ), ], - # The tool paper has six pages including references, so it cannot explain - # the package -- these docs have to. `checkdocs = :all` ignores the export + # `checkdocs = :all` ignores the export # list, so it covers every docstring in the module and *fails the build* on # one that no `@autodocs` block under `src/reference/` picks up. Adding a # file to `src/` therefore means adding it to a `Pages` list as well. diff --git a/docs/src/examples/set_prog.jl b/docs/src/examples/set_prog.jl new file mode 100644 index 0000000..909f81d --- /dev/null +++ b/docs/src/examples/set_prog.jl @@ -0,0 +1,154 @@ +# # Piecewise semi-ellipsoids: a template from SetProg +# +# [SetProg](https://github.com/blegat/SetProg.jl) writes optimisation problems +# over sets: ellipsoids, sublevel sets of polynomials, polytopes, and sets that +# are one of these on each cone of a partition. [`SetProgTemplate`](@ref) +# makes such a set the node function of a node, the function whose unit +# sublevel set it is, and the stability inequality +# ``V_d(A_i x) \le \gamma^d V_s(x)`` becomes the inclusion +# ``(A_i / \gamma)\, S_s \subseteq S_d``, which SetProg imposes. Three of its +# sets are templates this package already has, and the example checks that; +# the fourth, the piecewise semi-ellipsoid, is new. + +import PathCompleteCertificates as PCC +import Clarabel, HiGHS, JuMP +using SetProg +using DynamicPolynomials +using LinearAlgebra +using Plots + +const OPTIMIZER = Clarabel.Optimizer + +# The rotation and shear of the other examples, on the one-node graph, so that +# the template alone decides the rate. + +A = [0.69 * [0.0 1.0; -1.0 0.0], 0.69 * [1.0 1.0; 0.0 1.0]] +problem = PCC.StabilityProblem(PCC.switched_system(A)) +graph = PCC.common_graph(problem.system) + +bound(template) = + PCC.jsr_bound(template, graph, problem; optimizer = OPTIMIZER, rtol = 1e-3) + +# A partition is a polytope with the origin inside: its pieces are the cones +# over its facets. The polytope is built with Polyhedra on a linear programming +# solver, which SetProg works on the pieces with; a simplex solver such as +# HiGHS keeps that cheap. + +LIBRARY = Polyhedra.DefaultLibrary{Float64}( + JuMP.optimizer_with_attributes(HiGHS.Optimizer, JuMP.MOI.Silent() => true), +) + +polygon(normals) = polyhedron(hrep([HalfSpace(a, norm(a)) for a in normals]), LIBRARY) + +octagon = polygon([[2, 1], [1, 2], [-1, 2], [-2, 1], [-2, -1], [-1, -2], [1, -2], [2, -1]]) + +# ## Three sets that are templates already +# +# An ellipsoid is a quadratic form, the sublevel set of a polynomial is a sum +# of squares, and a polytope that is linear on each cone of a partition is the +# free-facet polyhedral template on the same cones; the cones over the facets +# of this octagon are those of [`planar_conic_partition`](@ref)`(1)`. Each +# pair certifies the same rate. + +@polyvar x[1:2] + +pairs = [ + ( + "quadratic", + PCC.QuadraticTemplate(), + PCC.SetProgTemplate(Ellipsoid(; symmetric = true)), + ), + ( + "sum of squares, degree 4", + PCC.SumOfSquaresTemplate(2, x), + PCC.SetProgTemplate(PolySet(; symmetric = true, degree = 4, variables = x)), + ), + ( + "free facets, 4 cones", + PCC.ConicPolyhedralTemplate([PCC.planar_conic_partition(1)]), + PCC.SetProgTemplate(Polytope(; symmetric = true, piecewise = octagon)), + ), +] + +[ + (name, round(bound(ours).rate; digits = 4), round(bound(theirs).rate; digits = 4)) for + (name, ours, theirs) in pairs +] + +# ## The piecewise semi-ellipsoid +# +# `Ellipsoid(symmetric = true, piecewise = P)` is a function quadratic on each +# cone of `P` and continuous across them, a template none of the others gives. +# Refining the partition can only help, since a function on the coarser one is +# a function on the finer one too. + +square = polygon([[1, 0], [0, 1], [-1, 0], [0, -1]]) +sixteen = polygon([ + [1, 0], + [2, 1], + [1, 1], + [1, 2], + [0, 1], + [-1, 2], + [-1, 1], + [-2, 1], + [-1, 0], + [-2, -1], + [-1, -1], + [-1, -2], + [0, -1], + [1, -2], + [1, -1], + [2, -1], +]) + +partitions = [("square", square), ("octagon", octagon), ("16-gon", sixteen)] + +piecewise = [ + (name, bound(PCC.SetProgTemplate(Ellipsoid(; symmetric = true, piecewise = P)))) for + (name, P) in partitions +] + +[ + (name, nhalfspaces(P), round(certificate.rate; digits = 4)) for + ((name, P), (_, certificate)) in zip(partitions, piecewise) +] + +# On the four cones of the square the piecewise function certifies the rate of +# the plain ellipsoid; the octagon's eight cones improve on it, and sixteen +# cones of degree two beat the polynomial of degree four above. + +# ## Drawing it +# +# The sublevel set ``\{x : V(x) \le 1\}`` of each function, normalised as in +# [Templates and graphs: the two axes](@ref). The piecewise one is a union of +# arcs of ellipses, meeting along the rays of the partition. + +function sublevel_boundary(certificate; samples = 721) + degree = PCC.homogeneity_degree(PCC.template(certificate)) + angles = range(0, 2pi; length = samples) + radii = [certificate([cos(a), sin(a)])^(-1 / degree) for a in angles] + radii ./= maximum(radii) + + return radii .* cos.(angles), radii .* sin.(angles) +end + +figure = plot(; + title = "one node, three partitions", + aspect_ratio = :equal, + legend = :outerbottom, + framestyle = :origin, +) + +for (name, certificate) in piecewise + xs, ys = sublevel_boundary(certificate) + plot!( + figure, + xs, + ys; + label = "$name, rate $(round(certificate.rate; digits = 3))", + lw = 2, + ) +end + +figure diff --git a/docs/src/index.md b/docs/src/index.md index c350e25..7ece0c2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -64,7 +64,9 @@ In reading order: The two other problems have an example each: [Safety: a barrier certificate](@ref) and -[Optimal control: a feedback per node](@ref). +[Optimal control: a feedback per node](@ref). A template from another +package has one too: +[Piecewise semi-ellipsoids: a template from SetProg](@ref). ## Scope diff --git a/docs/src/manual/templates.md b/docs/src/manual/templates.md index 36a3cd2..c0d3888 100644 --- a/docs/src/manual/templates.md +++ b/docs/src/manual/templates.md @@ -38,7 +38,8 @@ the bound monotonically, and ``d = 1`` is exactly [`QuadraticTemplate`](@ref). quadratic template, a polynomial set the sum-of-squares one and a piecewise polytope the free-facet one on the same cones; what it adds are the piecewise semi-ellipsoids, quadratic on each cone of a partition, which improve on the -quadratic template as the partition is refined. +quadratic template as the partition is refined; +[Piecewise semi-ellipsoids: a template from SetProg](@ref) shows all four. The **homogeneity degree**, [`homogeneity_degree`](@ref), is the ``d`` in ``V(cx) = c^d V(x)``. Stability scales its edge inequality by ``\gamma^d``, so