From 7b673bcd379f74369c6500b771ae8c64743b520c Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Wed, 29 Sep 2021 12:36:09 -0300 Subject: [PATCH] #133 - Add transmission line test matrix (#181) * refactor tests * add tline modelc add tline test * Update Project.toml --- Project.toml | 1 + test/arithmetic.jl | 43 +++++++++++ test/constructors.jl | 24 ++++++ test/exponential.jl | 54 ++++++++++++++ test/models.jl | 21 ++++++ test/runtests.jl | 170 ++----------------------------------------- test/setops.jl | 51 +++++++++++++ 7 files changed, 200 insertions(+), 164 deletions(-) create mode 100644 test/arithmetic.jl create mode 100644 test/constructors.jl create mode 100644 test/exponential.jl create mode 100644 test/models.jl create mode 100644 test/setops.jl diff --git a/Project.toml b/Project.toml index 973bff7..c460fac 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19" diff --git a/test/arithmetic.jl b/test/arithmetic.jl new file mode 100644 index 0000000..1d25960 --- /dev/null +++ b/test/arithmetic.jl @@ -0,0 +1,43 @@ +@testset "Interval arithmetic" begin + a = -1.5 ± 0.5 + b = -1..1 + @test a + b == -3.0..0.0 + @test a * b == -2.0..2.0 + @test a * b + a == -4.0..1.0 + @test a * (b + 1) == -4.0..0.0 +end + +@testset "Interval matrix arithmetic" begin + a = 1.0..1.3; b = 2.0..3.5; c = -0.5 ± 0.5; d = 0.0 ± 0.1 + a₊ = 2.0..2.6; b₊ = 4.0..7.0; c₊ = -2.0..0.0; d₊ = -0.2..0.2 + a₋ = -0.3..0.3; b₋ = -1.5..1.5; c₋ = -1.0..1.0; d₋ = -0.2..0.2 + A = IntervalMatrix([a b; c d]) + + B = A + A + @test B isa IntervalMatrix && B == IntervalMatrix([a₊ b₊; c₊ d₊]) + + B = A - A + @test B isa IntervalMatrix && B == IntervalMatrix([a₋ b₋; c₋ d₋]) + + B = A * A + @test B isa IntervalMatrix + + # multiply interval and interval matrix + x = 0.0..2.0 + @test x * A == A * x == + IntervalMatrix([0.0..2.6 0.0..7.0; -2.0..0.0 -0.2..0.2]) + # multiply scalar and interval matrix + x = 1.0 + for A2 in [x * A, A * x, A / x] + @test A2 == A && typeof(A2) == typeof(A) + end + + # arithmetic closure using interval matrices and non-interval matrices + Ainf = inf(A) + @test A + Ainf isa IntervalMatrix + @test Ainf + A isa IntervalMatrix + @test A - Ainf isa IntervalMatrix + @test Ainf - A isa IntervalMatrix + @test A * Ainf isa IntervalMatrix + @test Ainf * A isa IntervalMatrix +end diff --git a/test/constructors.jl b/test/constructors.jl new file mode 100644 index 0000000..c74b657 --- /dev/null +++ b/test/constructors.jl @@ -0,0 +1,24 @@ +@testset "Interval matrix construction" begin + m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9]) + m2 = IntervalMatrix{Float64}(undef, 2, 2) + @test m2 isa IntervalMatrix{Float64} && size(m2) == (2, 2) + m3 = similar(m1) + @test m3 isa IntervalMatrix{Float64} && size(m3) == size(m1) + m = [1.0 2.0; 3.0 4.0] + mint = IntervalMatrix([Interval(1) Interval(2); Interval(3) Interval(4)]) + @test IntervalMatrix(m) == mint + + A = [1 2; 3 4] + B = [1 2; 4 5] + + @test IntervalMatrix(A, B) == IntervalMatrix([1..1 2..2; 3..4 4..5]) + + @test A ± B == IntervalMatrix([0..2 0..4;-1..7 -1..9]) +end + +@testset "Interval matrix split" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + c, s = split(m) + @test c ≈ [-0.1 -4.; 4. -0.1] + @test s ≈ [1. 0.1; 0.1 1.] +end diff --git a/test/exponential.jl b/test/exponential.jl new file mode 100644 index 0000000..c0ba3ef --- /dev/null +++ b/test/exponential.jl @@ -0,0 +1,54 @@ +@testset "Interval matrix exponential" begin + + @test quadratic_expansion(-3..3, 1.0, 2.0) == Interval(-0.125, 21) + + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + + for i in 0:4 + _truncated_exponential_series(m, 1.0, i) + end + + overapp1 = exp_overapproximation(m, 1.0, 4) + overapp2 = horner(m, 10) + overapp3 = scale_and_square(m, 5, 1.0, 4) + underapp = exp_underapproximation(m, 1.0, 4) + + @test underapp isa IntervalMatrix + for overapp in [overapp1, overapp2, overapp3] + @test overapp isa IntervalMatrix + end +end + +@testset "Interval matrix correction terms" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + f = correction_hull(m, 1e-3, 5) + f2 = input_correction(m, 1e-3, 5) + f = correction_hull(mid(m), 1e-3, 5) +end + +@testset "Interval matrix square" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + a = m * m + b = square(m) + @test all(inf(a) .<= inf(b)) && all(sup(b) .>= sup(a)) +end + +@testset "Interval matrix power" begin + m = IntervalMatrix([2.0..2.0 2.0..3.0; 0.0..0.0 -1.0..1.0]) + pow = IntervalMatrixPower(m) + increment!(pow) # power of two + increment!(pow, algorithm="multiply") + increment!(pow) # power of two + increment!(pow, algorithm="power") + increment!(pow, algorithm="decompose_binary") + increment!(pow, algorithm="intersect") +end + +@testset "Transmission line model" begin + A = transmission_line() |> Matrix + expA = exp(A) + @test opnorm(expA, Inf)< 1e-6 + Aint = IntervalMatrix(interval.(A)) + expA_int = exp_overapproximation(Aint, 1.0, 10) + @test opnorm(expA_int, Inf) > 1e122 +end diff --git a/test/models.jl b/test/models.jl new file mode 100644 index 0000000..17ad62c --- /dev/null +++ b/test/models.jl @@ -0,0 +1,21 @@ +#= +State matrix `A` from the transmission line model [1]; we refer to that paper +for the interpretation of the physical parameters in this matrix. Here `A` is a +sparse, structured, square matrix of order 2η, where η ≥ 2 is a parameter +(optional, default: 2). For any choice of (positive) parameters, this matrix +is stable in the sense of Hurwitz (i.e. its eigenvalues have strictly negative +real-part). A parameter `scale` (optional, default: 1e-9) to rescale the physical +quantities is applied to `A`. + +[1] M. Althoff, B. H. Krogh, and O. Stursberg. Analyzing Reachability of +Linear Dynamic Systems with Parametric Uncertainties. Modeling, Design, +and Simulation of Systems with Uncertainties. 2011. +=# +function transmission_line(; η=2, R₀=1.00, Rd₀=10.0, L₀=1e-10, C₀=1e-13 * 4.00, scale=1e-9) + A₁₁ = zeros(η, η) + A₁₂ = Bidiagonal(fill(-1/C₀, η), fill(1/C₀, η-1), :U) + A₂₁ = Bidiagonal(fill(1/L₀, η), fill(-1/L₀, η-1), :L) + A₂₂ = Diagonal(vcat(-Rd₀/L₀, fill(-R₀/L₀, η-1))) + A = [A₁₁ A₁₂; A₂₁ A₂₂] + return A .* scale +end diff --git a/test/runtests.jl b/test/runtests.jl index 27d3bd2..78a497f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,171 +1,13 @@ -using IntervalMatrices, Test, LinearAlgebra +using IntervalMatrices, Test, LinearAlgebra, SparseArrays using IntervalMatrices: _truncated_exponential_series, horner, scale_and_square, correction_hull, input_correction +include("models.jl") -@testset "Interval arithmetic" begin - a = -1.5 ± 0.5 - b = -1..1 - @test a + b == -3.0..0.0 - @test a * b == -2.0..2.0 - @test a * b + a == -4.0..1.0 - @test a * (b + 1) == -4.0..0.0 -end - -@testset "Interval matrix construction" begin - m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9]) - m2 = IntervalMatrix{Float64}(undef, 2, 2) - @test m2 isa IntervalMatrix{Float64} && size(m2) == (2, 2) - m3 = similar(m1) - @test m3 isa IntervalMatrix{Float64} && size(m3) == size(m1) - m = [1.0 2.0; 3.0 4.0] - mint = IntervalMatrix([Interval(1) Interval(2); Interval(3) Interval(4)]) - @test IntervalMatrix(m) == mint - - A = [1 2; 3 4] - B = [1 2; 4 5] - - @test IntervalMatrix(A, B) == IntervalMatrix([1..1 2..2; 3..4 4..5]) - - @test A ± B == IntervalMatrix([0..2 0..4;-1..7 -1..9]) -end - -@testset "Interval matrix arithmetic" begin - a = 1.0..1.3; b = 2.0..3.5; c = -0.5 ± 0.5; d = 0.0 ± 0.1 - a₊ = 2.0..2.6; b₊ = 4.0..7.0; c₊ = -2.0..0.0; d₊ = -0.2..0.2 - a₋ = -0.3..0.3; b₋ = -1.5..1.5; c₋ = -1.0..1.0; d₋ = -0.2..0.2 - A = IntervalMatrix([a b; c d]) - - B = A + A - @test B isa IntervalMatrix && B == IntervalMatrix([a₊ b₊; c₊ d₊]) - - B = A - A - @test B isa IntervalMatrix && B == IntervalMatrix([a₋ b₋; c₋ d₋]) - - B = A * A - @test B isa IntervalMatrix - - # multiply interval and interval matrix - x = 0.0..2.0 - @test x * A == A * x == - IntervalMatrix([0.0..2.6 0.0..7.0; -2.0..0.0 -0.2..0.2]) - # multiply scalar and interval matrix - x = 1.0 - for A2 in [x * A, A * x, A / x] - @test A2 == A && typeof(A2) == typeof(A) - end - - # arithmetic closure using interval matrices and non-interval matrices - Ainf = inf(A) - @test A + Ainf isa IntervalMatrix - @test Ainf + A isa IntervalMatrix - @test A - Ainf isa IntervalMatrix - @test Ainf - A isa IntervalMatrix - @test A * Ainf isa IntervalMatrix - @test Ainf * A isa IntervalMatrix -end - -@testset "Interval matrix methods" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9]) - @test opnorm(m) == opnorm(m, Inf) ≈ 5.2 - @test opnorm(m, 1) ≈ 5.3 - l = inf(m) - r = sup(m) - c = mid(m) - d = diam(m) - rad = radius(m) - m2 = copy(m) - @test m2 isa IntervalMatrix && m.mat == m2.mat - @test l == inf.(m) && r == sup.(m) && c == mid.(m) - @test d ≈ r - l - @test rad ≈ d/2 - sm = scale(m, 2.0) - @test sm == 2.0 .* m - @test sm ≠ m - scale!(m, 2.0) # in-place - @test sm == m - m3 = IntervalMatrix([-2.0..2.0 -2.0..0.0; 0.0..2.0 -1.0..1.0]) - m4 = IntervalMatrix([-1.0..1.0 -1.0..1.0; -1.0..1.0 -2.0..2.0]) - @test m3 ∩ m4 == IntervalMatrix([-1.0..1.0 -1.0..0.0; 0.0..1.0 -1.0..1.0]) - @test m3 ∪ m4 == IntervalMatrix([-2.0..2.0 -2.0..1.0; -1.0..2.0 -2.0..2.0]) - @test diam_norm(m3) ≈ 6.0 # default diameter p-norm is Inf - @test diam_norm(m3, 1) ≈ 6.0 -end - -@testset "Interval matrix exponential" begin - - @test quadratic_expansion(-3..3, 1.0, 2.0) == Interval(-0.125, 21) - - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - - for i in 0:4 - _truncated_exponential_series(m, 1.0, i) - end - - overapp1 = exp_overapproximation(m, 1.0, 4) - overapp2 = horner(m, 10) - overapp3 = scale_and_square(m, 5, 1.0, 4) - underapp = exp_underapproximation(m, 1.0, 4) - - @test underapp isa IntervalMatrix - for overapp in [overapp1, overapp2, overapp3] - @test overapp isa IntervalMatrix - end -end - -@testset "Interval matrix split" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - c, s = split(m) - @test c ≈ [-0.1 -4.; 4. -0.1] - @test s ≈ [1. 0.1; 0.1 1.] -end - -@testset "Interval matrix membership" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - a1 = [0. -4.; 4. 0.] - a2 = [0. -3.; 4. 0.] - @test a1 ∈ m - @test a2 ∉ m -end - -@testset "Interval matrix inclusion" begin - m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - m2 = IntervalMatrix([-1.2..1.0 -4.1.. -3.9; 3.9..4.2 -1.2..0.9]) - @test m1 ⊆ m2 && !(m2 ⊆ m1) - -end - -@testset "Interval matrix rand" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - for i in 1:3 - @test rand(m) ∈ m - end -end - -@testset "Interval matrix correction terms" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - f = correction_hull(m, 1e-3, 5) - f2 = input_correction(m, 1e-3, 5) - f = correction_hull(mid(m), 1e-3, 5) -end - -@testset "Interval matrix square" begin - m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) - a = m * m - b = square(m) - @test all(inf(a) .<= inf(b)) && all(sup(b) .>= sup(a)) -end - -@testset "Interval matrix power" begin - m = IntervalMatrix([2.0..2.0 2.0..3.0; 0.0..0.0 -1.0..1.0]) - pow = IntervalMatrixPower(m) - increment!(pow) # power of two - increment!(pow, algorithm="multiply") - increment!(pow) # power of two - increment!(pow, algorithm="power") - increment!(pow, algorithm="decompose_binary") - increment!(pow, algorithm="intersect") -end +include("constructors.jl") +include("arithmetic.jl") +include("setops.jl") +include("exponential.jl") diff --git a/test/setops.jl b/test/setops.jl new file mode 100644 index 0000000..a289fa8 --- /dev/null +++ b/test/setops.jl @@ -0,0 +1,51 @@ +@testset "Interval matrix methods" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9]) + @test opnorm(m) == opnorm(m, Inf) ≈ 5.2 + @test opnorm(m, 1) ≈ 5.3 + + l = inf(m) + r = sup(m) + c = mid(m) + d = diam(m) + rad = radius(m) + m2 = copy(m) + @test m2 isa IntervalMatrix && m.mat == m2.mat + @test l == inf.(m) && r == sup.(m) && c == mid.(m) + @test d ≈ r - l + @test rad ≈ d/2 + + sm = scale(m, 2.0) + @test sm == 2.0 .* m + @test sm ≠ m + scale!(m, 2.0) # in-place + @test sm == m + + m3 = IntervalMatrix([-2.0..2.0 -2.0..0.0; 0.0..2.0 -1.0..1.0]) + m4 = IntervalMatrix([-1.0..1.0 -1.0..1.0; -1.0..1.0 -2.0..2.0]) + @test m3 ∩ m4 == IntervalMatrix([-1.0..1.0 -1.0..0.0; 0.0..1.0 -1.0..1.0]) + @test m3 ∪ m4 == IntervalMatrix([-2.0..2.0 -2.0..1.0; -1.0..2.0 -2.0..2.0]) + @test diam_norm(m3) ≈ 6.0 # default diameter p-norm is Inf + @test diam_norm(m3, 1) ≈ 6.0 +end + +@testset "Interval matrix membership" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + a1 = [0. -4.; 4. 0.] + a2 = [0. -3.; 4. 0.] + @test a1 ∈ m + @test a2 ∉ m +end + +@testset "Interval matrix inclusion" begin + m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + m2 = IntervalMatrix([-1.2..1.0 -4.1.. -3.9; 3.9..4.2 -1.2..0.9]) + @test m1 ⊆ m2 && !(m2 ⊆ m1) + +end + +@testset "Interval matrix rand" begin + m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9]) + for i in 1:3 + @test rand(m) ∈ m + end +end