From 582e8914eda530d269144eb8e048ed367ff26132 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 15 Jun 2024 23:43:30 +0200 Subject: [PATCH] revise discretization module --- src/Discretization/CorrectionHullModule.jl | 7 +- src/Discretization/DiscretizationModule.jl | 8 +- src/Discretization/Exponentiation.jl | 205 +++++++++--------- .../FirstOrderZonotopeModule.jl | 6 +- src/Discretization/ForwardBackwardModule.jl | 2 +- src/Discretization/ForwardModule.jl | 15 +- 6 files changed, 120 insertions(+), 123 deletions(-) diff --git a/src/Discretization/CorrectionHullModule.jl b/src/Discretization/CorrectionHullModule.jl index 85f237f77..5cf632712 100644 --- a/src/Discretization/CorrectionHullModule.jl +++ b/src/Discretization/CorrectionHullModule.jl @@ -4,9 +4,10 @@ module CorrectionHullModule using ..DiscretizationModule -using ..Exponentiation: _exp, _alias, IntervalExpAlg, interval_matrix +using ..Exponentiation: _exp, _alias, IntervalExpAlg using ..Overapproximate: _convert_or_overapproximate, _overapproximate -using IntervalMatrices: IntervalMatrix, correction_hull, input_correction, _exp_remainder +using IntervalMatrices: AbstractIntervalMatrix, IntervalMatrix, correction_hull, + input_correction, _exp_remainder using MathematicalSystems using LazySets using LazySets: LinearMap @@ -137,7 +138,7 @@ function discretize(ivp::IVP{<:CLCCS,<:LazySet}, δ, alg::CorrectionHull) Φ = _exp(A, δ, alg.exp) - A_interval = interval_matrix(A) + A_interval = A isa AbstractIntervalMatrix ? A : IntervalMatrix(A) origin_not_contained_in_U = zeros(dim(U)) ∉ Uz diff --git a/src/Discretization/DiscretizationModule.jl b/src/Discretization/DiscretizationModule.jl index 67d431074..5f720a8b3 100644 --- a/src/Discretization/DiscretizationModule.jl +++ b/src/Discretization/DiscretizationModule.jl @@ -22,7 +22,7 @@ next_set(inputs::AbstractInput, state::Int64) = collect(nextinput(inputs, state) abstract type AbstractApproximationModel end -function _default_approximation_model(ivp::IVP{<:AbstractContinuousSystem}) +function _default_approximation_model(::IVP{<:AbstractContinuousSystem}) return Forward() end @@ -39,17 +39,13 @@ isinterval(A::IntervalMatrix{N,IT}) where {N,IT<:IA.Interval{N}} = true isinterval(A::AbstractMatrix{IT}) where {IT<:IA.Interval} = true # options for a-posteriori transformation of a discretized set -# valid options are: -# AbstractDirections, Val{:lazy}, Val{:concrete}, Val{:vrep}, Val{:zono}, Val{:zonotope} -# _alias(setops) = setops # no-op - _alias(setops::AbstractDirections) = setops _alias(setops::Val{:lazy}) = setops _alias(setops::Val{:concrete}) = setops _alias(setops::Val{:vrep}) = setops _alias(setops::Val{:box}) = setops _alias(setops::Val{:zono}) = setops -_alias(setops::Val{:zonotope}) = Val(:zono) +_alias(::Val{:zonotope}) = Val(:zono) """ discretize(ivp::IVP, δ, alg::AbstractApproximationModel) diff --git a/src/Discretization/Exponentiation.jl b/src/Discretization/Exponentiation.jl index a76f3303c..59c5d2237 100644 --- a/src/Discretization/Exponentiation.jl +++ b/src/Discretization/Exponentiation.jl @@ -1,8 +1,8 @@ """ Interface to matrix exponential backends of different kinds. -Includes common integral computations arising in the discretization of linear differential equations -using matrix methods. For applications see e.g. [1] and references therein. +Includes common integral computations arising in the discretization of linear differential +equations using matrix methods. For applications see e.g. [1] and references therein. [1] Conservative Time Discretization: A Comparative Study. Marcelo Forets and Christian Schilling (2022). @@ -20,8 +20,7 @@ using IntervalArithmetic using IntervalMatrices using LazySets -export BaseExp, BaseExpAlg, IntervalExpAlg, LazyExpAlg, PadeExpAlg, elementwise_abs, Φ₂, Φ₁, Φ₁_u, - interval_matrix +export BaseExp, BaseExpAlg, IntervalExpAlg, LazyExpAlg, PadeExpAlg, elementwise_abs, Φ₂, Φ₁, Φ₁_u # ------------------------- # Exponentiation interface @@ -43,23 +42,22 @@ _alias(alg::Symbol) = _alias(Val(alg)) """ BaseExpAlg <: AbstractExpAlg -Matrix exponential using the scaling and squaring algorithm implemented -in Julia Base. +Matrix exponential using the scaling-and-squaring algorithm implemented in Julia +Base. ### Notes -If the array is static, the method implemented in `StaticArrays.jl` is applied. -This algorithm admits the alias `:base`. +The alias for this algorithm is `:base`. """ struct BaseExpAlg <: AbstractExpAlg end const BaseExp = BaseExpAlg() -_alias(alg::Val{:base}) = BaseExp +_alias(::Val{:base}) = BaseExp """ LazyExpAlg <: AbstractExpAlg -Lazy wrapper for the matrix exponential, operations defined in `LazySets.jl`. +Matrix exponential computed in a lazy fashion. ### Fields @@ -68,18 +66,20 @@ Lazy wrapper for the matrix exponential, operations defined in `LazySets.jl`. ### Notes -This algorithm admits the alias `:lazy` and also `:krylov`(using default options). +The aliases for this algorithm are `:lazy` and `:krylov`. + +The operations are defined in the package `LazySets.jl` (`SparseMatrixExp`). """ struct LazyExpAlg <: AbstractExpAlg - m::Int # size of the Krylov subspace - tol::Float64 # tolerance + m::Int + tol::Float64 end # default constructor LazyExpAlg() = LazyExpAlg(30, 1e-10) const LazyExp = LazyExpAlg() -_alias(alg::Union{Val{:lazy},Val{:krylov}}) = LazyExp +_alias(::Union{Val{:lazy},Val{:krylov}}) = LazyExp """ IntervalExpAlg <: AbstractExpAlg @@ -92,9 +92,10 @@ Matrix exponential using an interval enclosure of the Taylor series remainder. ### Notes -This method allows to overapproximate ``exp(Aδ)`` with an interval matrix. -It also accepts an interval matrix ``A`` given as input. -This method admits the alias `:interval` and `:taylor`. +The aliases for this algorithm are `:interval` and `:taylor`. + +This algorithm allows to overapproximate ``exp(Aδ)`` with an interval matrix. +It also accepts an interval matrix ``A`` as input. """ struct IntervalExpAlg <: AbstractExpAlg order::Int @@ -104,8 +105,7 @@ end IntervalExpAlg() = IntervalExpAlg(10) const IntervalExp = IntervalExpAlg() -_alias(alg::Val{:interval}) = IntervalExp -_alias(alg::Val{:taylor}) = IntervalExp +_alias(::Union{Val{:interval},Val{:taylor}}) = IntervalExp """ PadeExpAlg <: AbstractExpAlg @@ -114,55 +114,14 @@ Matrix exponential for sparse matrices using Pade approximants. ### Notes -Requires `Expokit.jl`. This algorithm admits the alias `:pade`. +The alias for this algorithm is `:pade`. + +This algorithm requires to load the package `Expokit.jl`. """ struct PadeExpAlg <: AbstractExpAlg end const PadeExp = PadeExpAlg() -_alias(alg::Val{:pade}) = PadeExp - -# general case: convert to Matrix -_exp(A::AbstractMatrix, ::BaseExpAlg) = exp(Matrix(A)) -_exp(A::Matrix, ::BaseExpAlg) = exp(A) - -# static arrays have their own exp method -_exp(A::StaticArray, ::BaseExpAlg) = exp(A) - -# exponential of an identity multiple, defined in MathematicalSystems.jl -_exp(A::IdentityMultiple, ::BaseExpAlg) = exp(A) - -# lazy wrapper of the matrix exponential -_exp(A::AbstractMatrix, ::LazyExpAlg) = SparseMatrixExp(A) - -# pade approximants (requires Expokit.jl) -_exp(::AbstractMatrix, ::PadeExpAlg) = throw(ArgumentError("algorithm requires a sparse matrix")) - -function _exp(A::SparseMatrixCSC, alg::PadeExpAlg) - require(@__MODULE__, :Expokit) - return _exp_pade(A, alg) -end - -function load_expokit_pade() - return quote - _exp_pade(A::SparseMatrixCSC, ::PadeExpAlg) = Expokit.padm(A) - end -end # quote / load_expokit_pade - -function _exp(A::AbstractMatrix, alg::IntervalExpAlg) - return exp_overapproximation(interval_matrix(A), one(eltype(A)), alg.order) -end - -function _exp(A::AbstractMatrix, δ, alg::IntervalExpAlg) - return exp_overapproximation(interval_matrix(A), δ, alg.order) -end - -# TODO Refactor to IntervalMatrices.jl? -function interval_matrix(A::AbstractIntervalMatrix) - return A -end -function interval_matrix(A::AbstractMatrix) - return IntervalMatrix(A) -end +_alias(::Val{:pade}) = PadeExp """ _exp(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp) @@ -174,8 +133,9 @@ Compute the matrix exponential ``e^{Aδ}``. - `A` -- matrix - `δ` -- step size - `alg` -- (optional, default: `BaseExp`) the algorithm used to take the matrix - exponential of `Aδ`, possible options are `BaseExp`, `LazyExp`, `PadeExp` and `IntervalExp` - see details in the *Algorithm* section below + exponential of `Aδ`, possible options are `BaseExp`, `LazyExp`, + `PadeExp` and `IntervalExp` (see details in the *Algorithm* section + below) ### Output @@ -183,18 +143,18 @@ A matrix or a lazy wrapper of the matrix exponential, depending on `alg`. ### Algorithm -- `BaseExp` -- (alias: `:base`) use Higham's scaling and squaring method implemented - in Julia standard library; see `?exp` for details; if `A` is a static array, +- `BaseExp` -- (alias: `:base`) use Higham's scaling-and-squaring method implemented + in Julia's standard library; see `?exp` for details; if `A` is a static array, uses the implementation in `StaticArrays.jl` -- `LazyExp` -- (alias: `:lazy`) return a lazy wrapper type around the matrix exponential +- `LazyExp` -- (alias: `:lazy`) return a lazy wrapper around the matrix exponential using the implementation `LazySets.SparseMatrixExp` -- `PadeExp` -- (alias: `pade`) apply Pade approximant method to compute the matrix +- `PadeExp` -- (alias: `pade`) apply the Padé approximant method to compute the matrix exponential of a sparse matrix (requires `Expokit.jl`) - `IntervalExp` -- (alias: `interval`, `taylor`) apply the Taylor series expansion of the matrix - exponential with an interval remainder; works if `A` is an interval matrix + exponential with an interval remainder; works if `A` is an interval matrix ### Notes @@ -203,32 +163,53 @@ evaluated with an external library such as `ExponentialUtilities.jl` or `Expokit.jl`. """ function _exp(A::AbstractMatrix, δ, alg::AbstractExpAlg=BaseExp) - n = checksquare(A) + checksquare(A) return _exp(A * δ, alg) end -@inline function _P_2n(A::AbstractMatrix{N}, δ, n) where {N} - return [Matrix(A * δ) Matrix(δ * I, n, n); - zeros(n, 2 * n)]::Matrix{N} -end +# general case: convert to Matrix +_exp(A::AbstractMatrix, ::BaseExpAlg) = exp(Matrix(A)) -@inline function _P_2n(A::SparseMatrixCSC{N,M}, δ, n) where {N,M} - return [sparse(A * δ) sparse(δ * I, n, n); - spzeros(n, 2 * n)]::SparseMatrixCSC{N,M} -end +# Base's `exp` +_exp(A::Matrix, ::BaseExpAlg) = exp(A) -@inline function _P_3n(A::AbstractMatrix{N}, δ, n) where {N} - return [Matrix(A * δ) Matrix(δ * I, n, n) zeros(n, n); - zeros(n, 2 * n) Matrix(δ * I, n, n); - zeros(n, 3 * n)]::Matrix{N} +# static arrays have their own exp method +_exp(A::StaticArray, ::BaseExpAlg) = exp(A) + +# exponential of an identity multiple, defined in MathematicalSystems.jl +_exp(A::IdentityMultiple, ::BaseExpAlg) = exp(A) + +# lazy wrapper of the matrix exponential +_exp(A::AbstractMatrix, ::LazyExpAlg) = SparseMatrixExp(A) + +# pade approximants (requires Expokit.jl) +_exp(::AbstractMatrix, ::PadeExpAlg) = throw(ArgumentError("algorithm requires a sparse matrix")) + +function _exp(A::SparseMatrixCSC, ::PadeExpAlg) + require(@__MODULE__, :Expokit) + return _exp_pade(A) end -@inline function _P_3n(A::SparseMatrixCSC{N,M}, δ, n) where {N,M} - return [sparse(A * δ) sparse(δ * I, n, n) spzeros(n, n); - spzeros(n, 2 * n) sparse(δ * I, n, n); - spzeros(n, 3 * n)]::SparseMatrixCSC{N,M} +function load_expokit_pade() + return quote + _exp_pade(A::SparseMatrixCSC) = Expokit.padm(A) + end +end # quote / load_expokit_pade + +function _exp(A::AbstractIntervalMatrix, δ, alg::IntervalExpAlg) + return exp_overapproximation(A, δ, alg.order) end +# convert to IntervalMatrix +_exp(A::AbstractMatrix, δ, alg::IntervalExpAlg) = _exp(IntervalMatrix(A), δ, alg) + +# add δ explicitly +_exp(A::AbstractMatrix, alg::IntervalExpAlg) = _exp(A, one(eltype(A)), alg) + +# ------------------------ +# Discretization functions +# ------------------------ + """ Φ₁(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp, [isinv]::Bool=false, [Φ]=nothing) @@ -282,13 +263,13 @@ It can be shown that where ``Φ(A, δ) = e^{Aδ}``. In particular, `Φ₁(A, δ) = P[1:n, (n+1):2*n]`. This method can be found in [[FRE11]](@ref). """ -function Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg=BaseExp, isinv=false, Φ=nothing) +function Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg=BaseExp, isinv::Bool=false, Φ=nothing) return Φ₁(A, δ, alg, Val(isinv), Φ) end # dispatch for discretization methods -Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, isinv::Val{:false}, Φ) = _Φ₁_blk(A, δ, alg) -Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, isinv::Val{:true}, Φ) = _Φ₁_inv(A, δ, alg, Φ) +Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, ::Val{:false}, Φ) = _Φ₁_blk(A, δ, alg) +Φ₁(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, ::Val{:true}, Φ) = _Φ₁_inv(A, δ, alg, Φ) # evaluate the series Φ₁(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+1}}{(i+1)!}A^i # without assuming invertibility of A, by taking the exponential of a 2n x 2n matrix @@ -299,6 +280,16 @@ function _Φ₁_blk(A, δ, alg) return _P₁_blk(Q, n, alg) end +@inline function _P_2n(A::AbstractMatrix{N}, δ, n) where {N} + return [Matrix(A * δ) Matrix(δ * I, n, n); + zeros(n, 2 * n)]::Matrix{N} +end + +@inline function _P_2n(A::SparseMatrixCSC{N,M}, δ, n) where {N,M} + return [sparse(A * δ) sparse(δ * I, n, n); + spzeros(n, 2 * n)]::SparseMatrixCSC{N,M} +end + @inline _P₁_blk(P, n::Int, ::AbstractExpAlg) = P[1:n, (n + 1):(2 * n)] @inline _P₁_blk(P, n::Int, ::LazyExpAlg) = sparse(get_columns(P, (n + 1):(2 * n))[1:n, :]) @@ -324,11 +315,11 @@ function _Φ₁_inv(A::IdentityMultiple, δ, alg, Φ=nothing) end # method to compute Φ₁ * u, where u is a vector -function Φ₁_u(A, δ, alg, isinv::Val{:true}, u::AbstractVector, Φ=nothing) +function Φ₁_u(A, δ, alg, ::Val{:true}, u::AbstractVector, Φ=nothing) return Φ₁_u(A, δ, alg, Φ, u) end -function Φ₁_u(A, δ, alg, isinv::Val{:false}, u::AbstractVector, Φ=nothing) +function Φ₁_u(A, δ, alg, ::Val{:false}, u::AbstractVector, Φ=nothing) M = _Φ₁_blk(A, δ, alg) return M * u end @@ -406,16 +397,13 @@ It can be shown that where ``Φ(A, δ) = e^{Aδ}``. In particular, `Φ₂ = P_{3n}[1:n, (2*n+1):3*n]`. This method can be found in [[FRE11]](@ref). """ -function Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg=BaseExp, isinv=false, Φ=nothing) +function Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg=BaseExp, isinv::Bool=false, Φ=nothing) return Φ₂(A, δ, alg, Val(isinv), Φ) end # dispatch for discretization methods -Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, isinv::Val{:false}, Φ) = _Φ₂_blk(A, δ, alg) -Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, isinv::Val{:true}, Φ) = _Φ₂_inv(A, δ, alg, Φ) - -@inline _P₂_blk(P, n, ::AbstractExpAlg) = P[1:n, (2 * n + 1):(3 * n)] -@inline _P₂_blk(P, n, ::LazyExpAlg) = sparse(get_columns(P, (2 * n + 1):(3 * n))[1:n, :]) +Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, ::Val{:false}, Φ) = _Φ₂_blk(A, δ, alg) +Φ₂(A::AbstractMatrix, δ::Real, alg::AbstractExpAlg, ::Val{:true}, Φ) = _Φ₂_inv(A, δ, alg, Φ) # evaluate the series Φ₂(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+2}}{(i+2)!}A^i # without assuming invertibility of A, by taking the exponential of a 3n x 3n matrix @@ -426,9 +414,23 @@ function _Φ₂_blk(A, δ, alg) return _P₂_blk(P, n, alg) end +@inline function _P_3n(A::AbstractMatrix{N}, δ, n) where {N} + return [Matrix(A * δ) Matrix(δ * I, n, n) zeros(n, n); + zeros(n, 2 * n) Matrix(δ * I, n, n); + zeros(n, 3 * n)]::Matrix{N} +end + +@inline function _P_3n(A::SparseMatrixCSC{N,M}, δ, n) where {N,M} + return [sparse(A * δ) sparse(δ * I, n, n) spzeros(n, n); + spzeros(n, 2 * n) sparse(δ * I, n, n); + spzeros(n, 3 * n)]::SparseMatrixCSC{N,M} +end + +@inline _P₂_blk(P, n, ::AbstractExpAlg) = P[1:n, (2 * n + 1):(3 * n)] +@inline _P₂_blk(P, n, ::LazyExpAlg) = sparse(get_columns(P, (2 * n + 1):(3 * n))[1:n, :]) + # compute the matrix Φ₂ = A^{-2} (exp(A*δ) - I - A*δ) assuming that A is invertible # and explicitly computing inv(A); this function optionally receives Φ = exp(Aδ) -# TODO don't pass algorithm since it is ignored if Φ is given # Φ2 = ReachabilityAnalysis._Φ₂_inv(abs.(A), δ, ReachabilityAnalysis.BaseExp, Φ) function _Φ₂_inv(A::AbstractMatrix, δ, alg, Φ=nothing) Aδ = A * δ @@ -463,11 +465,12 @@ function _Eplus(A::SparseMatrixCSC{N,D}, X0::AbstractHyperrectangle{N}, δt; m=m require(@__MODULE__, :ExponentialUtilities) Pv = _phiv(Aabs, v, 1, δt; m, tol) - return E⁺ = Hyperrectangle(zeros(n), Pv) + return Hyperrectangle(zeros(n), Pv) end # --------------- # Absolute values +# --------------- elementwise_abs(A::AbstractMatrix) = abs.(A) function elementwise_abs(A::SparseMatrixCSC) diff --git a/src/Discretization/FirstOrderZonotopeModule.jl b/src/Discretization/FirstOrderZonotopeModule.jl index b137e6a1d..f1a1d84de 100644 --- a/src/Discretization/FirstOrderZonotopeModule.jl +++ b/src/Discretization/FirstOrderZonotopeModule.jl @@ -4,7 +4,7 @@ module FirstOrderZonotopeModule using ..DiscretizationModule -using ..Exponentiation: _exp, _alias, BaseExp +using ..Exponentiation: _exp, BaseExp using LinearAlgebra using MathematicalSystems using LazySets @@ -123,10 +123,6 @@ end # ----------- function _discretize_zonotope(Φ, X0, ::FirstOrderZonotope, δ, norm_A) - c = center(X0) - G = genmat(X0) - n, p = size(G) - # compute X(δ) Xδ = linear_map(Φ, X0) diff --git a/src/Discretization/ForwardBackwardModule.jl b/src/Discretization/ForwardBackwardModule.jl index f5d49df5e..ba3d1a32a 100644 --- a/src/Discretization/ForwardBackwardModule.jl +++ b/src/Discretization/ForwardBackwardModule.jl @@ -155,7 +155,7 @@ function _get_e(d, E) end # homogeneous case -function ω(λ, d, Φᵀ, X0, U::ZeroSet, Eψ, δ, e⁺, e⁻) +function ω(λ, d, Φᵀ, X0, ::ZeroSet, Eψ, δ, e⁺, e⁻) aux1h = (1 - λ) * ρ(d, X0) + λ * ρ(Φᵀ * d, X0) aux2 = sum(min(λ * e⁺[i], (1 - λ) * e⁻[i]) * abs(d[i]) for i in eachindex(d)) return aux1h + aux2 diff --git a/src/Discretization/ForwardModule.jl b/src/Discretization/ForwardModule.jl index cff6aff4d..aeccebf09 100644 --- a/src/Discretization/ForwardModule.jl +++ b/src/Discretization/ForwardModule.jl @@ -94,7 +94,7 @@ function discretize(ivp::IVP{<:CLCS,<:LazySet}, δ, alg::Forward) return InitialValueProblem(Sdis, Ω0) end -function discretize(ivp::IVP{<:CLCS,Interval{N}}, δ, alg::Forward) where {N} +function discretize(ivp::IVP{<:CLCS,Interval{N}}, δ, ::Forward) where {N} A = state_matrix(ivp) @assert size(A, 1) == 1 X0 = initial_state(ivp) @@ -102,19 +102,19 @@ function discretize(ivp::IVP{<:CLCS,Interval{N}}, δ, alg::Forward) where {N} a = A[1, 1] aδ = a * δ Φ = exp(aδ) - A_abs = abs(a) # use inverse method @assert !iszero(a) "the given matrix should be invertible" + # A_abs = abs(a) # a_sqr = a * a #P2A_abs = (1/a_sqr) * (Φ - one(N) - aδ) - #Einit = (P2A_abs * a_sqr) * RA._symmetric_interval_hull(X0).dat + #E⁺ = (P2A_abs * a_sqr) * RA._symmetric_interval_hull(X0).dat #P2A_abs = (1/a_sqr) * (Φ - one(N) - aδ) - Einit = (Φ - one(N) - aδ) * convert(Interval, symmetric_interval_hull(X0)).dat + E⁺ = (Φ - one(N) - aδ) * convert(Interval, symmetric_interval_hull(X0)).dat - Ω0 = Interval(hull(X0.dat, Φ * X0.dat + Einit)) + Ω0 = Interval(hull(X0.dat, Φ * X0.dat + E⁺)) X = stateset(ivp) # the system constructor creates a matrix Sdis = ConstrainedLinearDiscreteSystem(Φ, X) @@ -135,7 +135,8 @@ function discretize(ivp::IVP{<:CLCCS,<:LazySet}, δ, alg::Forward) Φcache = sum(A) == abs(sum(A)) ? Φ : nothing P2A_abs = Φ₂(A_abs, δ, alg.exp, alg.inv, Φcache) - Einit = sih(P2A_abs * sih((A * A) * X0, alg.sih), alg.sih) + # TODO outsource to Exponentiation module and merge with _Eplus + E⁺ = sih(P2A_abs * sih((A * A) * X0, alg.sih), alg.sih) U = next_set(inputset(ivp), 1) Eψ0 = sih(P2A_abs * sih(A * U, alg.sih), alg.sih) @@ -143,7 +144,7 @@ function discretize(ivp::IVP{<:CLCCS,<:LazySet}, δ, alg::Forward) Ud = δ * U ⊕ Eψ0 In = IdentityMultiple(one(eltype(A)), size(A, 1)) - Ω0 = ConvexHull(X0, Φ * X0 ⊕ Ud ⊕ Einit) + Ω0 = ConvexHull(X0, Φ * X0 ⊕ Ud ⊕ E⁺) Ω0 = _apply_setops(Ω0, alg.setops) X = stateset(ivp) Sdis = ConstrainedLinearControlDiscreteSystem(Φ, In, X, Ud)