Use IntervalArithmetic package for intervals

.github/workflows/CI.yml

@@ -19,15 +19,15 @@ jobs:
       matrix:
         version:
           - '1.6'
-          - '1.8'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/documentation.yml

@@ -13,10 +13,10 @@ jobs:
       contents: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.6'
+          version: '1.10'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -7,6 +7,8 @@ version = "0.1.0"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RealTimeScheduling = "f1524149-62f5-490d-8249-9bfac76a7111"
@@ -15,6 +17,7 @@ RealTimeScheduling = "f1524149-62f5-490d-8249-9bfac76a7111"
 ControlSystemsBase = "1"
 DataStructures = "0.18"
 Distances = "0.10"
+IntervalArithmetic = "0.20"
 RealTimeScheduling = "0.3"
 julia = "^1.6"
 

README.md

@@ -54,8 +54,8 @@ notebook, `notebooks/control_timing_safety.jl`.
 ## Citing
 
 If you use this code in your research, we ask that you consider citing the
-relevant papers.  For the `bounded_runs_iter` function for overapproximating
-the maximum deviation, please cite:
+relevant papers.  For the `uncertain_linear_system` and `bounded_runs_iter`
+functions for overapproximating the maximum deviation, please cite:
 
 > Clara Hobbs, Bineet Ghosh, Shengjie Xu, Parasara Sridhar Duggirala, and
 > Samarjit Chakraborty.

docs/make.jl

@@ -7,6 +7,7 @@ makedocs(
     sitename = "ControlTimingSafety",
     format = Documenter.HTML(),
     modules = [ControlTimingSafety],
+    checkdocs = :exports,
     pages = [
         "index.md",
         "automata.md",

docs/src/safety.md

@@ -1,13 +1,14 @@
 # Checking Safety
 
-The main safety-checking function of ControlTimingSafety is
-[`bounded_runs_iter`](@ref), which implements the Bounded Runs Iteration
-algorithm from "Safety Analysis of Embedded Controllers under Implementation
-Platform Timing Uncertainties."  This implementation permits an initial point,
-or an axis-aligned box initial set.  The helper function
+ControlTimingSafety implements two methods from "Safety Analysis of Embedded
+Controllers under Implementation Platform Timing Uncertainties": the Uncertain
+Linear Systems method in [`uncertain_linear_system`](@ref), and the Bounded
+Runs Iteration method in [`bounded_runs_iter`](@ref).  The implementation
+permits an axis-aligned box initial set.  The helper function
 [`bounded_runs`](@ref), which computes a single iteration, is also exported.
 
 ```@docs
+uncertain_linear_system
 bounded_runs_iter
 bounded_runs
 ```

notebooks/control_timing_safety.jl

@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.19.14
+# v0.19.22
 
 #> [frontmatter]
 #> title = "ControlTimingSafety.jl Demo"
@@ -29,12 +29,12 @@ begin
     Pkg.instantiate()
 
 	using Revise
-	using ControlTimingSafety
+	using ControlTimingSafety, IntervalArithmetic
 	using RealTimeScheduling
     using Plots, PlutoUI, LaTeXStrings
 	using Random, Distributions
 	using Distances
-	using ControlSystems, LinearAlgebra
+	using ControlSystemsBase, LinearAlgebra
 	using Printf
 end
 
@@ -150,7 +150,7 @@ let
 	for i = size(hit_prob, 1):-1:1
 	    seq = (rand(sim_time) .>= hit_prob[i]) .+ 1
 		z = evol(automaton, z_0, seq)
-	    plot!(z[:,1], z[:,2], label=L"P(\mathrm{hit}) = %$(hit_prob[i])", linecolor=line_colors[i])
+	    plot!(z[1,:], z[2,:], label=L"P(\mathrm{hit}) = %$(hit_prob[i])", linecolor=line_colors[i])
 	end
 	plot!()
 end
@@ -208,7 +208,7 @@ trj_2 = let
 	# Compute random trajectories
 	sp = SamplerUniformMissRow(MissRow(maxmiss_2), H_2)
 	seq = ones(Int64, H_2)
-	trj = Array{Float64}(undef, n_samples_2, H_2+1, size(T.Φ[1], 1))
+	trj = Array{Float64}(undef, n_samples_2, size(T.Φ[1], 1), H_2+1)
 	max_possible = (H_2 / (maxmiss_2 + 1)) * maxmiss_2 + H_2 % (maxmiss_2 + 1)
 	for i = axes(trj, 1)
 		accepted = false
@@ -232,7 +232,7 @@ We next compute the distance between each point of each random trajectory and th
 dist_2 = let	
 	dist = Array{Float64}(undef, n_samples_2, H_2+1)
 	for i in axes(trj_2, 1)
-		dist[i,:] = colwise(Euclidean(), trj_2[i,:,1:2]', nominal_2[:,1:2]')
+		dist[i,:] = colwise(Euclidean(), trj_2[i,1:2,:], nominal_2[1:2,:])
 	end
 	dist
 end
@@ -256,16 +256,16 @@ let
 	pipe_radius = sort(dev_2, dims=1)[end-1] - 0.01
 	circ_x = cos.(θ) * pipe_radius
 	circ_y = sin.(θ) * pipe_radius
-	for i in axes(nominal_2, 1)
-		plot!(circ_x .+ nominal_2[i,1], circ_y .+ nominal_2[i,2], repeat([i,], size(θ,1)), label=(i == 1) ? "Safety pipe" : "", seriestype=[:shape,], c=:lightblue, linecolor=:lightblue)
+	for i in axes(nominal_2, 2)
+		plot!(circ_x .+ nominal_2[1,i], circ_y .+ nominal_2[2,i], repeat([i,], size(θ,1)), label=(i == 1) ? "Safety pipe" : "", seriestype=[:shape,], c=:lightblue, linecolor=:lightblue)
 	end
 
 	# Plot random trajectories
 	label_good = false
 	# Plot the good trajectories first
 	for i = axes(trj_2, 1)
 		if dev_2[i] < pipe_radius
-			plot!(trj_2[i,:,1], trj_2[i,:,2], 1:H_2+1, label=(!label_good) ? "Random" : "", color=:green, opacity=10/n_samples_2)
+			plot!(trj_2[i,1,:], trj_2[i,2,:], 1:H_2+1, label=(!label_good) ? "Random" : "", color=:green, opacity=10/n_samples_2)
 			label_good = true
 		end
 	end
@@ -277,14 +277,14 @@ let
 			t_viol = argmax(dist_2[i,:], dims=1)
 			#label_bad || plot!(cos.(θ) * pipe_radius .+ nominal_2[t_viol,1], sin.(θ) * pipe_radius .+ nominal_2[t_viol,2], repeat([t_viol,], size(θ,1)), label="Radius violated", style=:dot, linecolor=:gray)
 			crosses = [(d >= pipe_radius) ? 1 : 0 for d in dist_2[i,:]]
-			plot!(trj_2[i,:,1], trj_2[i,:,2], 1:H_2+1, label=(!label_bad) ? "Violation" : "", color=:red, markershape=:x, markeralpha=crosses)
+			plot!(trj_2[i,1,:], trj_2[i,2,:], 1:H_2+1, label=(!label_bad) ? "Violation" : "", color=:red, markershape=:x, markeralpha=crosses)
 			nom_crosses .+= crosses
 			label_bad = true
 		end
 	end
 
 	# Finally, plot the nominal trajectory on top of everything else
-	plot!(nominal_2[:,1], nominal_2[:,2], 1:H_2+1, label="Nominal", color=:black, linewidth=3, marker=:x, markeralpha=nom_crosses)
+	plot!(nominal_2[1,:], nominal_2[2,:], 1:H_2+1, label="Nominal", color=:black, linewidth=3, marker=:x, markeralpha=nom_crosses)
 end
   ╠═╡ =#
 
@@ -314,7 +314,7 @@ This section illustrates the Bounded Runs algorithm, used to explore all possibl
 """
 
 # ╔═╡ 0843dd0e-0091-402f-9f28-fb112232f140
-bounds = [2; 2;; 2.5; 2.5]
+bounds = IntervalBox(2..2.5, 2..2.5)
 
 # ╔═╡ ef358d22-a530-4881-a282-b850d3655427
 begin
@@ -338,7 +338,7 @@ let
 	plot(title="Reachable set for $(n_bounded_runs) time steps", xlabel=L"x_1", ylabel=L"x_2", legend=:bottomright)
 	merged = merge_bounds(points_bounded_runs)
 	for k = 0:n_bounded_runs
-		corners = corners_from_bounds(merged[begin+k,:,:], cycle=true, dims=1:2)
+		corners = corners_from_bounds(merged[begin+k], cycle=true, dims=1:2)
 		plot!(corners[1,:], corners[2,:], label=L"x[%$(k)]")
 	end
 
@@ -347,7 +347,7 @@ let
 		corners = augment(hsn, corners_from_bounds(bounds, dims=[1,2]))
 		for c in eachcol(corners)
 	    	x = evol(hsn, c, ones(Int64, n_bounded_runs))
-			plot!(x[:,1], x[:,2], label="Sample", linecolor=:blue, marker=:circle)
+			plot!(x[1,:], x[2,:], label="Sample", linecolor=:blue, marker=:circle)
 		end
 	end
 	plot!()
@@ -390,7 +390,7 @@ end
 let
 	plot(title="Reachable set for $(time_4) time steps, $(n_4) per tree", xlabel=L"x_1", ylabel=L"x_2", legend=:bottomright)
 	for t in 1:time_4+1
-		corners = corners_from_bounds(all_bounds[t,:,:], cycle=true, dims=1:2)
+		corners = corners_from_bounds(all_bounds[t], cycle=true, dims=1:2)
 		plot!(corners[1,:], corners[2,:], label=L"x[%$(t-1)]")
 	end
 
@@ -399,7 +399,7 @@ let
 		corners = augment(hsn, corners_from_bounds(bounds, dims=[1,2]))
 		for c in eachcol(corners)
 	    	x = evol(hsn, c, ones(Int64, n_4*t_4))
-			plot!(x[:,1], x[:,2], label="Sample", linecolor=:blue, marker=:circle)
+			plot!(x[1,:], x[2,:], label="Sample", linecolor=:blue, marker=:circle)
 		end
 	end
 	plot!(legend=false)
@@ -413,7 +413,7 @@ Using these computed reachable sets, we finally compute the maximum deviation th
 # ╔═╡ 7afec811-9ccc-4e93-9cc8-a3c6b6ee5c49
 let
 	automaton = strat_map[sim_strategy_4](sysd, K, MissRow(max_miss_4))
-	d = deviation(automaton, augment(automaton, bounds), all_bounds, dims=1:2)[1:time_4+1]
+	d = deviation(automaton, augment(automaton, bounds), all_bounds)[1:time_4+1]
 	i = argmax(d)
 	v = maximum(d)
 

src/ControlTimingSafety.jl

@@ -4,13 +4,15 @@ using LinearAlgebra
 using Random
 using ControlSystemsBase
 using Distances
+using IntervalArithmetic
 
 export Automaton, nlocations, nactions
 export hold_skip_next, zero_skip_next, hold_kill, zero_kill
 export strat_map, strat_names
 export evol_final, evol, evol_final!, evol!, augment
 include("automata.jl")
 
+export uncertain_linear_system
 export corners_from_bounds, merge_bounds, merge_bounds!
 export bounded_runs, bounded_runs_iter, deviation
 include("safety.jl")

src/automata.jl

@@ -397,19 +397,47 @@ and `l` is the final location in the automaton.
 """
 function evol_final(a::Automaton, z_0::AbstractVector{Float64}, input::AbstractVector{Int64})
     @boundscheck length(z_0) == a.nz || throw(DimensionMismatch("z_0 must have length a.nz"))
-    z = zeros(length(input) + 1, a.nz)
-    z[1,:] = z_0
+    z = zeros(a.nz, length(input) + 1)
+    z[:,1] = z_0
+    evol_final!(a, z, input)
+end
+
+function evol_final(a::Automaton, z_0::IntervalBox, input::AbstractVector{Int64})
+    @boundscheck length(z_0) == a.nz || throw(DimensionMismatch("z_0 must have length a.nz"))
+    z = Vector{IntervalBox{length(z_0), eltype(inf.(z_0))}}(undef, length(input) + 1)
+    z[1] = z_0
     evol_final!(a, z, input)
 end
 
 """
     evol_final!(a, z, input)
 
-Same as [`evol_final`](@ref), but writes to the input matrix `z`, whose first row is `z_0`.
+Same as [`evol_final`](@ref), but writes to the input matrix `z`, whose first column is `z_0`.
+`z` may also be a vector of `IntervalBox`es, whose first element is `z_0`.
 """
 function evol_final!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVector{Int64})
-    @boundscheck size(z,1) == length(input)+1 || throw(DimensionMismatch("z must have size (length(input)+1, a.nz)"))
-    @boundscheck size(z,2) == a.nz || throw(DimensionMismatch("z must have size (length(input)+1, a.nz)"))
+    @boundscheck size(z,1) == a.nz || throw(DimensionMismatch("z must have size (a.nz, length(input)+1)"))
+    @boundscheck size(z,2) == length(input)+1 || throw(DimensionMismatch("z must have size (a.nz, length(input)+1)"))
+    l = a.l_int
+    # For each time step and input action
+    for (t, in) in enumerate(input)
+        # Get the dynamics matrix
+        μ = a.μ[l, in]
+        # If we hit a missing transition, return the states that we reached,
+        # and a missing final location to signal the problem to the caller.
+        if ismissing(μ)
+            return z[:,1:t], missing
+        end
+        # Apply the dynamics
+        mul!(view(z,:,t+1), a.Φ[μ], view(z,:,t))
+        # Transition to the new location
+        l = a.T[l, in]
+    end
+    z, l
+end
+
+function evol_final!(a::Automaton, z::AbstractVector{IntervalBox{S, T}}, input::AbstractVector{Int64}) where {S, T}
+    @boundscheck length(z) == length(input)+1 || throw(DimensionMismatch("z must have length length(input)+1"))
     l = a.l_int
     # For each time step and input action
     for (t, in) in enumerate(input)
@@ -418,10 +446,10 @@ function evol_final!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVe
         # If we hit a missing transition, return the states that we reached,
         # and a missing final location to signal the problem to the caller.
         if ismissing(μ)
-            return z[1:t,:], missing
+            return z[1:t], missing
         end
         # Apply the dynamics
-        mul!(view(z,t+1,:), a.Φ[μ], view(z,t,:))
+        z[t+1] = a.Φ[μ] * z[t]
         # Transition to the new location
         l = a.T[l, in]
     end
@@ -439,14 +467,14 @@ Returns `z`, a matrix of states over time.
 See also [`evol_final`](@ref), which additionally returns the final location in the
 automaton.
 """
-evol(a::Automaton, z_0::AbstractVector{Float64}, input::AbstractVector{Int64}) = evol_final(a, z_0, input)[1]
+evol(a::Automaton, z_0::Union{AbstractVector{Float64}, IntervalBox{S, T}}, input::AbstractVector{Int64}) where {S, T} = evol_final(a, z_0, input)[1]
 
 """
     evol!(a, z, input)
 
-Same as [`evol`](@ref), but writes to the input matrix `z`, whose first row is `z_0`.
+Same as [`evol`](@ref), but writes to the input matrix `z`, whose first column is `z_0`.
 """
-evol!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVector{Int64}) = evol_final!(a, z, input)[1]
+evol!(a::Automaton, z::Union{AbstractMatrix{Float64}, AbstractVector{IntervalBox{S, T}}}, input::AbstractVector{Int64}) where {S, T} = evol_final!(a, z, input)[1]
 
 """
     augment(a::Automaton, x)
@@ -456,3 +484,4 @@ augmented state space of `a`.
 """
 augment(a::Automaton, x::AbstractMatrix) = [x; zeros(a.nz - size(x, 1), size(x, 2))]
 augment(a::Automaton, x::AbstractVector) = [x; zeros(a.nz - size(x, 1))]
+augment(a::Automaton, x::IntervalBox) = IntervalBox([x...; zeros(a.nz - length(x))])