From 5f72622e0e2193cacae6f39265ce33ce26036942 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 1 Oct 2019 09:21:07 +0200 Subject: [PATCH 1/2] fix Schur transformation --- src/Utils/transformations.jl | 40 +++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/src/Utils/transformations.jl b/src/Utils/transformations.jl index 8f97860a..04e4a779 100644 --- a/src/Utils/transformations.jl +++ b/src/Utils/transformations.jl @@ -1,3 +1,5 @@ +using LinearAlgebra + export transform """ @@ -25,8 +27,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 +42,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 +60,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 +88,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 From af48d590bea3ef7e97c68842f291ca976f851801 Mon Sep 17 00:00:00 2001 From: schillic Date: Sun, 13 Oct 2019 15:50:46 +0200 Subject: [PATCH 2/2] add code to undo a coordinate transformation --- src/Utils/transformations.jl | 29 ++++++++++++++++++++++++++++- src/solve.jl | 5 ++++- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/Utils/transformations.jl b/src/Utils/transformations.jl index 04e4a779..24044360 100644 --- a/src/Utils/transformations.jl +++ b/src/Utils/transformations.jl @@ -1,6 +1,33 @@ using LinearAlgebra -export transform +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) diff --git a/src/solve.jl b/src/solve.jl index c1bcd4e1..1c1d045d 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -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 """