diff --git a/examples/Controlled Invariant Set.ipynb b/examples/Controlled Invariant Set.ipynb index 017ed38..09489f3 100644 --- a/examples/Controlled Invariant Set.ipynb +++ b/examples/Controlled Invariant Set.ipynb @@ -39,7 +39,7 @@ "0 & 0\n", "\\end{bmatrix}.\n", "$$\n", - "As shown in [LTJ18], a set $S$ is controlled invariant if\n", + "As shown in [LTJ18], a set $S$ is controlled invariant iff\n", "$$\n", "\\begin{bmatrix}\n", "1 & \\Delta t\n", @@ -59,7 +59,9 @@ { "cell_type": "code", "execution_count": 2, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "using Polyhedra\n", @@ -77,7 +79,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "using SCS\n", @@ -88,7 +92,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "using CSDP # CSDP is less appropriate than SCS and Mosek because it does not natively support\n", @@ -99,7 +105,9 @@ { "cell_type": "code", "execution_count": 3, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "using MosekTools\n", @@ -1942,7 +1950,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.1.0" + "version": "1.1.1" } }, "nbformat": 4, diff --git a/src/Sets/ellipsoids.jl b/src/Sets/ellipsoids.jl index 0c619c7..51eb942 100644 --- a/src/Sets/ellipsoids.jl +++ b/src/Sets/ellipsoids.jl @@ -92,7 +92,7 @@ function _HPH(D, d, δ) d D] end -_HPH(ell::EllipsoidAtOrigin) = _HPH(ell.Q, zeros(size(q.Q, 1)), -1.0) +_HPH(ell::EllipsoidAtOrigin) = _HPH(ell.Q, zeros(size(ell.Q, 1)), -1.0) """ struct ShiftedEllipsoid{T} diff --git a/src/Sets/transformations.jl b/src/Sets/transformations.jl index db7b4b7..44abaea 100644 --- a/src/Sets/transformations.jl +++ b/src/Sets/transformations.jl @@ -7,6 +7,16 @@ perspective_variable(li::LinearImage) = perspective_variable(li.set) space_variables(li::LinearImage) = nothing dimension(li::LinearImage) = size(li.A, 1) +# A^{-1} * S + b +struct AffineImage{S, T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractSet{T} + set::S + A::MT + b::VT +end +perspective_variable(ai::AffineImage) = perspective_variable(ai.set) +space_variables(ai::AffineImage) = nothing +dimension(ai::AffineImage) = size(ai.A, 1) + # A^{-1} * S struct LinearPreImage{S, T, MT <: AbstractMatrix{T}} <: AbstractSet{T} set::S diff --git a/src/affine.jl b/src/affine.jl index 71394ab..e8a5559 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -8,6 +8,10 @@ struct AffineExpression{T, F <: AbstractScalarFunction} <: AbstractScalarFunctio constant::T end +function Base.copy(aff::AffineExpression) + return AffineExpression(copy(aff.terms), aff.constant) +end + function Base.:(+)(f::F, g::F) where F <: AbstractScalarFunction return AffineExpression([AffineTerm(1.0, f), AffineTerm(1.0, g)], 0.0) end diff --git a/src/inclusion.jl b/src/inclusion.jl index 1b2814d..4be6096 100644 --- a/src/inclusion.jl +++ b/src/inclusion.jl @@ -50,6 +50,10 @@ function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage, <:Sets.LinearImage}) return set_space(space, DualSpace) end +function set_space(space::Space, ::InclusionConstraint{<:Sets.AffineImage, + <:Sets.LinearImage}) + return set_space(space, DualSpace) +end # AS ⊆ S <=> S ⊆ A^{-1}S so PrimalSpace work # AS ⊆ S <=> A^{-T}S∘ ⊆ S∘ so DualSpace work function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage, @@ -140,14 +144,23 @@ function JuMP.build_constraint(_error::Function, end # See [LTJ18] +function _apply_maps(_error, subset, supset; kws...) + dim = Sets.dimension(subset) + @polyvar x[1:dim] + JuMP.build_constraint(_error, apply_map(subset, x), + PowerSet(apply_map(supset, x)); kws...) +end function JuMP.build_constraint(_error::Function, subset::Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}}, sup_powerset::PowerSet{<:Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}}}; kws...) - dim = Sets.dimension(subset) - @polyvar x[1:dim] - JuMP.build_constraint(_error, apply_map(subset, x), - PowerSet(apply_map(sup_powerset.set, x)); kws...) + _apply_maps(_error, subset, sup_powerset.set; kws...) +end +function JuMP.build_constraint(_error::Function, + subset::Sets.AffineImage{<:Sets.PerspectiveDual}, + sup_powerset::PowerSet{<:Sets.LinearImage{<:Sets.PerspectiveDual}}; + kws...) + _apply_maps(_error, subset, sup_powerset.set; kws...) end function JuMP.build_constraint(_error::Function, subset::Sets.AbstractSet, diff --git a/src/map.jl b/src/map.jl index 97dae27..6fbad34 100644 --- a/src/map.jl +++ b/src/map.jl @@ -2,6 +2,11 @@ need_variablify(lm::Sets.LinearImage) = need_variablify(lm.set) function variablify(lm::Sets.LinearImage) return Sets.LinearImage(variablify(lm.set), lm.A) end +need_variablify(lm::Sets.AffineImage) = need_variablify(lm.set) +function variablify(lm::Sets.AffineImage) + return Sets.AffineImage(variablify(lm.set), lm.A, lm.b) +end + """ apply_map(li::Sets.LinearImage{<:Sets.PolarOf{<:Sets.EllipsoidAtOrigin}}, new_vars) @@ -64,7 +69,28 @@ function apply_map(li::Sets.LinearImage{<:Sets.PerspectiveDualOf{<:Sets.Househol return Sets.perspective_dual(dual) end +# [x'] = [A b] [x] +# [z'] = [0 1] [z] +# The transpose is +# [A' 0] +# [b' 1] +# So we substitute x for A' x' and then z for b' x' + z' +function apply_map(li::Sets.AffineImage{<:Sets.PerspectiveDualOf{<:Sets.Householder{T}}}, + new_vars) where T + z = Sets.perspective_variable(li.set) + old_vars = Sets.space_variables(li.set) + q = subs(Sets.perspective_gauge0(li.set.set), old_vars => li.A' * new_vars) + # We do this substitution after in case `new_vars` and `old_vars` have an intersection + p = subs(q, z => z + dot(li.b, new_vars)) + dual = Sets.Householder(Sets.UnknownSet{T}(), p, li.A * li.set.set.h + li.b, + z, new_vars) + return Sets.perspective_dual(dual) +end + # FIXME, for Sets.AbstractSet, we should apply it directly function Base.:(*)(A::AbstractMatrix, set::Union{SetVariableRef, Sets.AbstractSet}) return Sets.LinearImage(set, A) end +function Base.:(+)(li::Sets.LinearImage, b::AbstractVector) + return Sets.AffineImage(li.set, li.A, b) +end diff --git a/src/variables.jl b/src/variables.jl index fe8f8ee..8cdab18 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -1,3 +1,5 @@ +export InteriorPoint, CenterPoint + abstract type AbstractVariable <: JuMP.AbstractVariable end abstract type HintPoint end