Skip to content

Commit

Permalink
Merge pull request #689 from JuliaReach/schillic/schur
Browse files Browse the repository at this point in the history
Fix crashes in Schur transformation
schillic authored Oct 16, 2019
2 parents 464fdcb + af48d59 commit d0eae9f
Showing 2 changed files with 62 additions and 12 deletions.
69 changes: 58 additions & 11 deletions src/Utils/transformations.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,33 @@
export transform
using LinearAlgebra

export transform,
backtransform

"""
backtransform(Rsets::ReachSolution, options::Options)
Undo a coordinate transformation.
### Input
- `Rsets` -- flowpipe
- `option` -- problem options containing an `:transformation_matrix` entry
### Output
A new flowpipe where each reach set has been transformed.
### Notes
The transformation is implemented with a lazy `LinearMap`.
"""
function backtransform(Rsets, options::Options)
transformation_matrix = options[:transformation_matrix]
if transformation_matrix == nothing
return Rsets
end
return project(Rsets, transformation_matrix)
end

"""
transform(problem::InitialValueProblem, options::Options)
@@ -25,8 +54,8 @@ function transform(problem::InitialValueProblem, options::Options)
if method == ""
nothing # no-op
elseif method == "schur"
options[:transformation_matrix] = Tm
problem = schur_transform(problem)
problem, T_inverse = schur_transform(problem)
options[:transformation_matrix] = T_inverse
else
error("the transformation method $method is undefined")
end
@@ -40,7 +69,7 @@ Applies a Schur transformation to a discrete or continuous initial-value problem
### Input
- `S` -- discrete or continuous initial-value problem
- `problem` -- discrete or continuous initial-value problem
### Output
@@ -58,19 +87,26 @@ function schur_transform(problem::InitialValueProblem{PT, ST}
n = size(problem.s.A, 1)

# full (dense) matrix is required by schur
A_new, T_new = schur(Matrix(problem.s.A))
# result S is a struct such that
# - S.Schur is in Schur form and
# - A == S.vectors * S.Schur * S.vectors'
S = schur(Matrix(problem.s.A))

# recall that for Schur matrices, inv(T) == T'
Z_inverse = copy(transpose(T_new))
Z_inverse = copy(transpose(S.vectors))

# obtain transformed system matrix
A_new = S.Schur

# apply transformation to the initial states
X0_new = Z_inverse * problem.x0

# apply transformation to the inputs
U_new = Z_inverse * problem.s.U
B_new = Z_inverse * problem.s.B
U_new = problem.s.U

# obtain the transformation matrix for reverting the transformation again
T_inverse = F[:vectors]
# matrix for reverting the transformation again
T_inverse = S.vectors

# apply transformation to the invariant
if hasmethod(stateset, Tuple{typeof(problem.s)})
@@ -79,6 +115,17 @@ function schur_transform(problem::InitialValueProblem{PT, ST}
invariant_new = Universe(n)
end

system_new = (A_new, I(n), invariant_new, U_new)
return InitialValueProblem(PT(system_new), X0_new)
system_new = _wrap_system(PT, A_new, B_new, invariant_new, U_new)
problem_new = InitialValueProblem(system_new, X0_new)
return problem_new, T_inverse
end

function _wrap_system(PT::Type{<:ConstrainedLinearControlDiscreteSystem},
A, B, invariant, U)
return ConstrainedLinearControlDiscreteSystem(A, B, invariant, U)
end

function _wrap_system(PT::Type{<:ConstrainedLinearControlContinuousSystem},
A, B, invariant, U)
return ConstrainedLinearControlContinuousSystem(A, B, invariant, U)
end
5 changes: 4 additions & 1 deletion src/solve.jl
Original file line number Diff line number Diff line change
@@ -71,7 +71,10 @@ function solve!(problem::InitialValueProblem,
problem, options = transform(problem, options)

# run the continuous-post operator
post(op, problem, options)
res = post(op, problem, options)

# undo a coordinate transformation
res = backtransform(res, options)
end

"""

0 comments on commit d0eae9f

Please sign in to comment.