diff --git a/.gitignore b/.gitignore index b067edd..ba39cc5 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1 @@ -/Manifest.toml +Manifest.toml diff --git a/README.md b/README.md index e1dad2a..4a73036 100644 --- a/README.md +++ b/README.md @@ -35,11 +35,12 @@ zs = zero_skip_next(sys, K, max_misses) To run the Bounded Runs Iteration algorithm for this automaton, first create the initial set, then run the algorithm for a given per-iteration run -length `n` and number of iterations `t`: +length `n` to produce the reachable set of length `H`. The number of iterations +is the smallest integer `t` such that `n`×`t`≥`H`. ```julia augbounds = augment(automaton, bounds) -all_bounds = bounded_runs_iter(automaton, augbounds, n, t) +all_bounds = bounded_runs_iter(automaton, augbounds, n, H) ``` The deviation at each time step can then be computed thusly: diff --git a/docs/src/index.md b/docs/src/index.md index 4aa10a1..b66f2f7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -31,11 +31,12 @@ zs = zero_skip_next(sys, K, max_misses) To run the Bounded Runs Iteration algorithm for this automaton, first create the initial set, then run the algorithm for a given per-iteration run -length `n` and number of iterations `t`: +length `n` to produce the reachable set of length `H`. The number of iterations +is the smallest integer `t` such that `n`×`t`≥`H`. ```julia augbounds = augment(automaton, bounds) -all_bounds = bounded_runs_iter(automaton, augbounds, n, t) +all_bounds = bounded_runs_iter(automaton, augbounds, n, H) ``` The deviation at each time step can then be computed thusly: diff --git a/docs/src/schedule_synthesis.md b/docs/src/schedule_synthesis.md index 0fc9983..9fd95e1 100644 --- a/docs/src/schedule_synthesis.md +++ b/docs/src/schedule_synthesis.md @@ -1,14 +1,30 @@ # Schedule Synthesis -The main functionality of this section is divided into two functions +The main functionality of this section is divided into three functions 1. [`synthesize_constraints`](@ref): given the dynamics and safety margin of **one** control task, find all weakly-hard constraints (up to a maximum window size `maxwindow`) - with which the system behaves safely, and -2. [`schedule_xghtc`](@ref): given the lists of safe weakly-hard constraints for **all** + with which the system behaves safely, +2. [`estimate_constraints`](@ref): similar to [`synthesize_constraints`](@ref), but uses + the statistical method [`estimate_deviation`](@ref) to find the deviation upper bound of + each weakly-hard constraints. +3. [`schedule_xghtc`](@ref): given the lists of safe weakly-hard constraints for **all** control tasks, synthesizes a schedule so that they all behave safely (if such schedules exist). ```@docs synthesize_constraints schedule_xghtc +estimate_constraints +``` + +## Utility Functions + +A few utility functions are exported as well. These support the implementation +of the above functions, and are potentially useful to their callers as well. + +```@docs +schedule_to_sequence +verify_schedule +devub +devest ``` diff --git a/notebooks/control_timing_safety.jl b/notebooks/control_timing_safety.jl index 2449670..0672e0b 100644 --- a/notebooks/control_timing_safety.jl +++ b/notebooks/control_timing_safety.jl @@ -370,11 +370,10 @@ This section illustrates the Bounded Runs Iteration algorithm, used to efficient # ╔═╡ efc6aa00-cfdb-45fb-809b-f99c663690f6 begin - t_4 = div(time_4, n_4) all_bounds = let automaton = strat_map[sim_strategy_4](sysd, K, MissRow(max_miss_4)) augbounds = augment(automaton, bounds) - bounded_runs_iter(automaton, augbounds, n_4, t_4) + bounded_runs_iter(automaton, augbounds, n_4, time_4) end md""" diff --git a/src/ControlTimingSafety.jl b/src/ControlTimingSafety.jl index 6f0bb99..c63eab0 100644 --- a/src/ControlTimingSafety.jl +++ b/src/ControlTimingSafety.jl @@ -19,6 +19,8 @@ export maximum_deviation_random, estimate_deviation include("probablesafety.jl") export schedule_xghtc, synthesize_constraints +export estimate_constraints, verify_schedule, schedule_to_sequence +export devub, devest include("schedule_synthesis.jl") end diff --git a/src/safety.jl b/src/safety.jl index e3e624d..a3fea76 100644 --- a/src/safety.jl +++ b/src/safety.jl @@ -135,17 +135,19 @@ function bounded_runs(a::Automaton, z_0::AbstractVecOrMat, n::Integer) end """ - bounded_runs_iter(a, z_0, n, t) + bounded_runs_iter(a, z_0, n, H) -Iterate [`bounded_runs`](@ref)`(a, z_0, n)` for `t` iterations, returning the reachable -set at each of the `n`×`t+1` time steps. +Iterate [`bounded_runs`](@ref)`(a, z_0, n)` multiple times, returning the reachable set at each +of the `H` time stpes. The number of iterations is the smallest integer `t` such that +`n`×`t`≥`H`. See also [`deviation`](@ref), which can be called with the result of this function to find the deviation from a nominal trajectory. """ -function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, t::Integer; safety_margin::Float64=Inf) +function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, H::Integer; safety_margin::Real=Inf) + t = ceil(Int64, H / n) # Dimensions: time, augmented state, min/max - all_bounds = Array{Float64}(undef, n*(t+1)+1, a.nz, 2) + all_bounds = Array{Float64}(undef, n*t+1, a.nz, 2) if isa(z_0, AbstractVector) all_bounds[1,:,:] = [z_0 z_0] else @@ -176,7 +178,7 @@ function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, t::I # Dimensions: initial location, final location, time, augmented state, min/max new_bounds = Array{Float64}(undef, nlocations(a), nlocations(a), n+1, a.nz, 2) - for i in 1:t + for i in 2:t # Simulate each box from previous iteration Threads.@threads for i in a.L new_bounds[i,:,:,:,:] = bounded_runs(A[i], bounds[i,end,:,:], n) @@ -186,16 +188,16 @@ function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, t::I merge_bounds!(view(bounds, i, :, :, :), view(new_bounds, :, i, :, :, :)) end # Save the bounds - merge_bounds!(view(all_bounds, n*i+1:n*(i+1)+1, :, :), bounds) + merge_bounds!(view(all_bounds, n*(i-1)+1:n*i+1, :, :), bounds) if isfinite(safety_margin) - d = deviation(a, z_0, all_bounds[n*i+2:n*(i+1)+1,:,:], nominal_trajectory=nom[:,:,n*i+2:n*(i+1)+1]) + d = deviation(a, z_0, all_bounds[n*(i-1)+2:n*i+1,:,:], nominal_trajectory=nom[:,:,n*(i-1)+2:n*i+1]) if maximum(d) > safety_margin - return all_bounds[1:n*(i+1)+1,:,:] + return all_bounds[1:n*i+1,:,:] end end end - all_bounds + all_bounds[1:H+1,:,:] end """ diff --git a/src/schedule_synthesis.jl b/src/schedule_synthesis.jl index ce5522f..91f5537 100644 --- a/src/schedule_synthesis.jl +++ b/src/schedule_synthesis.jl @@ -2,7 +2,7 @@ using RealTimeScheduling using DataStructures """ - schedule_xghtc(constraints; slotsize=1, H=100, work_conserving=false) + schedule_xghtc(constraints, H; slotsize=1, work_conserving=false) Generate a schedule for a set of weakly hard constraints. The schedule returned has the type Matrix{Int64}, where the first dimension iterates through tasks, and the second @@ -20,7 +20,8 @@ Constraints." ASPDAC 2023. DOI: [10.1145/3566097.3567848](https://doi.org/10.1145/3566097.3567848) """ -function schedule_xghtc(constraints::Vector{<:MeetAny}; slotsize::Int64=1, H::Int64=100, work_conserving::Bool=false) +function schedule_xghtc(constraints::Vector{<:MeetAny}, H::Integer; + slotsize::Integer=1, work_conserving::Bool=false) # Check if "utilization" is greater than available slot size utilization = sum(c -> c.meet/c.window, constraints) if utilization > slotsize @@ -72,41 +73,181 @@ function schedule_xghtc(constraints::Vector{<:MeetAny}; slotsize::Int64=1, H::In end """ - synthesize_constraints(sysd, K, z_0, d_max, maxwindow, n, t) + schedule_xghtc(allconstraints, H; slotsize=1, work_conserving=false) + +Given a list of safe constraints for each system, synthesize a schedule using one safe +constraint from each list. The first successful synthesis result is returned. +""" +function schedule_xghtc(allconstraints::Vector{<:Vector{<:MeetAny}}, H::Integer; + slotsize::Integer=1, work_conserving::Bool=false) + for constraints in Iterators.product(allconstraints...) + sch = schedule_xghtc(collect(constraints), H, + slotsize=slotsize, work_conserving=work_conserving) + if length(sch) > 0 + println(constraints) + return sch + end + end + return zeros(Int64, 0, 0) +end + +""" + verify_schedule(sysd, K, z_0, σ) + +Given a discrete state-space model, feedback gain, the initial state, and a sequence of +deadline hits and misses, returns the maximum deviaiton between the resulting trajectory +and the nominal trajectory for that system where all deadlines are met. +""" +function verify_schedule(sysd::AbstractStateSpace{<:Discrete}, + K::AbstractMatrix{<:Real}, z_0::AbstractVecOrMat, σ::AbstractVector{<:Integer}) + a = hold_kill(sysd, K) + + # Convert the actions σ to: 1=hit, 2=miss (instead of 0=miss) + σ = [i == 0 ? 2 : 1 for i in σ] + z = evol(a, z_0, σ) + + # Convert trajectory to reachable sets by repeating the state twice for min/max + reachable = cat(z, z, dims=3) + maximum(deviation(a, z_0, reachable)) +end + +""" + schedule_to_sequence(schedule, task, H) + +Given a schedule in the matrix form, the position of the task in the schedule and the desired +length of the sequence, returns a sequence of 0s and 1s of length H representing the deadline +hits and misses for that task under the given schedule. +""" +function schedule_to_sequence(schedule::Matrix{<:Integer}, task::Integer, H::Integer) + σ = schedule[task, :] + [repeat(σ, H ÷ length(σ)); σ[1:H % length(σ)]] +end + +""" + synthesize_constraints(sysd, K, z_0, d_max, maxwindow, n, H) Find all `MeetAny` weakly hard constraints with window size at most `maxwindow` that guarantees the deviation upper bound is at most `d_max`. The system is specified by -[`Automaton`](@ref) `a` and initial state is `z_0`. `n` and `t` are as in +[`Automaton`](@ref) `a` and initial state is `z_0`. `n` and `H` are as in [`bounded_runs_iter`](@ref). """ function synthesize_constraints(sysd::AbstractStateSpace{<:Discrete}, - K::AbstractMatrix{Float64}, z_0::AbstractVecOrMat, d_max::Float64, - maxwindow::Int64, n::Int64, t::Int64) - - safe_constraints = MeetAny[] - - # Do not need to go through all O(maxwindow^2) constraints, - # see paper for optimization argument - meet = 1 - for window in 2:maxwindow - while meet < window - constraint = MeetAny(meet, window) - a = hold_kill(sysd, K, constraint) - # Check if the deviation bound is within the safety margin - reachable = bounded_runs_iter(a, z_0, n, t) - m = maximum(deviation(a, z_0, reachable)) - if m <= d_max - # All constraints with (m, window) where m >= meet are valid - for i in meet:window-1 - push!(safe_constraints, MeetAny(i, window)) + K::AbstractMatrix{<:Real}, z_0::AbstractVecOrMat, d_max::Real, + maxwindow::Integer, n::Integer, H::Integer; fullresults=false, + nominal::Union{Matrix{<:Real}, Nothing}=nothing) + + devs = fill(Inf, maxwindow, maxwindow) + safe_constraints = [MeetAny(1, 1)] + + if fullresults + for window = 1:maxwindow, meet=1:window + devs[window, meet] = devub(meet, window, sysd, K, z_0, d_max, n, H) + if devs[window, meet] <= d_max && meet < window + push!(safe_constraints, MeetAny(meet, window)) + end + end + else + # Do not need to go through all O(maxwindow^2) constraints, + # see paper for optimization argument + meet = 1 + for window in 1:maxwindow + while meet <= window + devs[window, meet] = devub(meet, window, sysd, K, z_0, d_max, n, H, nominal) + if devs[window, meet] <= d_max + # All constraints with (m, window) where m >= meet are valid + for i in meet:window-1 + push!(safe_constraints, MeetAny(i, window)) + end + break end - break + meet += 1 end - meet += 1 end end - safe_constraints + safe_constraints, devs +end + +""" + devub(meet, window, sysd, K, z_0, d_max, n, H) + +Compute the deviation upper bound for a given system under the weakly hard constraint +(meet, window), over a specified time horizon. Uses [`bounded_runs_iter`](@ref). +""" +function devub(meet::Integer, window::Integer, sysd::AbstractStateSpace{<:Discrete}, + K::AbstractMatrix{<:Real}, z_0::AbstractVecOrMat, d_max::Real, n::Integer, + H::Integer, nominal::Union{Matrix{<:Real}, Nothing}=nothing) + @boundscheck nominal === nothing || size(nominal, 1) == H+1 || throw(ArgumentError("nominal and H mismatch")) + if meet == window && nominal === nothing + return 0. + end + constraint = MeetAny(meet, window) + a = hold_kill(sysd, K, constraint) + reachable = bounded_runs_iter(a, z_0, n, H, safety_margin=d_max) + if size(reachable, 1) != H+1 + @info "Size mismatch" + return d_max + end + # @info "Data" meet window size(reachable, 1) argmax(deviation(a, z_0, reachable)) + maximum(deviation(a, z_0, reachable, nominal_trajectory=nominal)) +end + +""" + estimate_constraints(sysd, K, z_0, d_max, maxwindow, c, B, H) + +Find all `MeetAny` weakly hard constraints with window size at most `maxwindow` that +guarantees the deviation upper bound is at most `d_max`. The system is specified by +[`Automaton`](@ref) `a` and initial state is `z_0`. This function uses SHT for +estimating the deviation upper bound as in [`estimate_deviation`](@ref). +""" +function estimate_constraints(sysd::AbstractStateSpace{<:Discrete}, + K::AbstractMatrix{<:Real}, z_0::AbstractVecOrMat, d_max::Real, + maxwindow::Integer, c::Real, B::Real, H::Integer; fullresults=false) + + devs = fill(Inf, maxwindow, maxwindow) + safe_constraints = [MeetAny(1, 1)] + + if fullresults + for window = 1:maxwindow, meet=1:window + devs[window, meet] = devest(meet, window, sysd, K, z_0, c, B, H) + if devs[window, meet] <= d_max && meet < window + push!(safe_constraints, MeetAny(meet, window)) + end + end + else + meet = 1 + for window in 1:maxwindow + while meet <= window + devs[window, meet] = devest(meet, window, sysd, K, z_0, c, B, H) + if devs[window, meet] <= d_max + for i in meet:window-1 + push!(safe_constraints, MeetAny(i, window)) + end + break + end + meet += 1 + end + end + end + + safe_constraints, devs +end + +""" + devest(meet, window, sysd, K, z_0, c, B, H) + +Estimate the deviation upper bound for a given system under the weakly hard constraint +(meet, window), over a specified time horizon. Uses [`estimate_deviation`](@ref). +""" +function devest(meet::Integer, window::Integer, sysd::AbstractStateSpace{<:Discrete}, + K::AbstractMatrix{<:Real}, z_0::AbstractVecOrMat, c::Real, B::Real, H::Integer) + if meet == window + return 0. + end + constraint = MeetAny(meet, window) + a = hold_kill(sysd, K, constraint) + sampler = SamplerUniformMeetAny(constraint, H) + estimate_deviation(a, sampler, z_0, c, B) end # Helper functions @@ -120,7 +261,7 @@ julia> ControlTimingSafety._undigit([1, 0, 0]) 4 ``` """ -function _undigit(d::Vector{Int64}; base=2) +function _undigit(d::Vector{<:Integer}; base=2) s = 0 mult = 1 for val in reverse(d) @@ -136,7 +277,7 @@ end Digits **b**ase **2**, **r**everse A shortcut for `digits(x, base=2, pad=pad) |> reverse` """ -_digits_b2r(x::Int64; pad::Int64=0) = digits(x, base=2, pad=pad) |> reverse +_digits_b2r(x::Integer; pad::Integer=0) = digits(x, base=2, pad=pad) |> reverse """ _state_separation(l, B[, indigits=false]) @@ -159,7 +300,7 @@ The explanation of the example is as follows [[1], [1, 0]] -> [1, 2] ``` """ -function _state_separation(l::Int64, B::Vector{Int64}; indigits=false) +function _state_separation(l::Integer, B::Vector{<:Integer}; indigits=false) @boundscheck l < 2^sum(B) || throw(ArgumentError("l has more bits than the sum of B")) bits = digits(l, base=2, pad=sum(B)) |> reverse @@ -241,7 +382,7 @@ end """ Build a scheduler automaton from a given list of controller automata """ -function _SynthesizedAutomaton(controllers::Vector{_ConstraintAutomaton}; slotsize::Int64=1, work_conserving::Bool=false) +function _SynthesizedAutomaton(controllers::Vector{_ConstraintAutomaton}; slotsize::Integer=1, work_conserving::Bool=false) # Converting tuple to array with collect() N = length(controllers) B = map(a -> ceil(Int64, log2(a.L)), controllers) |> collect @@ -277,7 +418,7 @@ function _SynthesizedAutomaton(controllers::Vector{_ConstraintAutomaton}; slotsi _SynthesizedAutomaton(N, B, L, Σ, T, L-1, Q_f) end -function _path_to_schedule(path::Union{LinkedList{Int64}, Vector{Int64}}, AS::_SynthesizedAutomaton) +function _path_to_schedule(path::Union{LinkedList{<:Integer}, Vector{<:Integer}}, AS::_SynthesizedAutomaton) # Convert path to Vector path = collect(path) diff --git a/test/safety_margin_issue.jl b/test/safety_margin_issue.jl new file mode 100644 index 0000000..81759bf --- /dev/null +++ b/test/safety_margin_issue.jl @@ -0,0 +1,35 @@ +using ControlTimingSafety +using ControlSystemsBase + +sys = let + r_1 = 100000 + r_2 = 500000 + r_3 = 200000 + c_1 = 0.000002 + c_2 = 0.000010 + A = [-1/c_1 * (1/r_1 + 1/r_2) 1/(r_2*c_1) + 1/(r_2*c_2) -1/c_2 * (1/r_2 + 1/r_3)] + B = [1/(r_1*c_1) + 1/(r_3*c_2)] + C = [1 0; 0 1] + D = 0 + + ss(A, B, C, D) +end +K = [0.16007906919762957 0.21425701514482584 0.022142131526773] + +p = 0.023 + +sysd = c2d(sys, p) + +z0 = [100., 100., 0.] + +d_max = 1000. +maxwindow = 6 +n = 10 +H = 100 + +# display(synthesize_constraints(sysd, K, z0, d_max, maxwindow, n, H, fullresults=true)[2]) +for i in 1:10 + println(devub(1, 2, sysd, K, z0, d_max, n, H)) +end diff --git a/test/schedule_synthesis_tests.jl b/test/schedule_synthesis_tests.jl index 783c481..5e5f13d 100644 --- a/test/schedule_synthesis_tests.jl +++ b/test/schedule_synthesis_tests.jl @@ -36,6 +36,7 @@ using RealTimeScheduling end @test synthesize_constraints(sysd, K, repeat([10.0], 3), 5.5, 6, 10, 10) == [ + RealTimeScheduling.MeetAny(1, 1), RealTimeScheduling.MeetAny(1, 2), RealTimeScheduling.MeetAny(2, 3), RealTimeScheduling.MeetAny(3, 4), @@ -49,7 +50,7 @@ using RealTimeScheduling RealTimeScheduling.MeetAny(1, 4), RealTimeScheduling.MeetAny(1, 4), RealTimeScheduling.MeetAny(1, 2) - ], slotsize=2) == [ + ], 100, slotsize=2) == [ 0 1 0 0 1 0 0 1 0 0 1 0 1 0 1 1 0 1 1 0 1 1 0 1 1 0 0 0 1 0 0 0 1 0 0 0