From 1c356561de6cccbfecf90a83d45ecebecab5ea9d Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 15 Mar 2023 17:44:27 +0100 Subject: [PATCH] add simple model reduction for uncertain models --- docs/make.jl | 1 + docs/src/uncertainty.md | 39 ++++++++++++++++ src/RobustAndOptimalControl.jl | 4 ++ src/mcm_nugap.jl | 84 ++++++++++++++++++++++++++++++++++ test/runtests.jl | 5 ++ test/test_mcm_nugap.jl | 32 +++++++++++++ 6 files changed, 165 insertions(+) create mode 100644 src/mcm_nugap.jl create mode 100644 test/test_mcm_nugap.jl diff --git a/docs/make.jl b/docs/make.jl index 7ea7456a..7e347e34 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,6 +10,7 @@ makedocs( sitename = "RobustAndOptimalControl Documentation", doctest = true, modules = [RobustAndOptimalControl, ControlSystemsBase], + pagesonly = false, pages = [ "Home" => "index.md", "Uncertainty modeling" => "uncertainty.md", diff --git a/docs/src/uncertainty.md b/docs/src/uncertainty.md index f18d94ec..4c9f9b08 100644 --- a/docs/src/uncertainty.md +++ b/docs/src/uncertainty.md @@ -4,6 +4,11 @@ We provide two general means of modeling uncertainty, the traditional $M\Delta$ Both the $M\Delta$ framework and parametric-uncertainty approaches are illustrated below. +```@contents +Pages = ["uncertainty.md"] +Depth = 3 +``` + ## Uncertainty API - [`δc`](@ref) Creates an uncertain complex parameter. - [`δr`](@ref) Creates an uncertain real parameter. @@ -183,6 +188,40 @@ sigmaplot!(So, w, c=2, lab="So") The sensitivity at the plant output is enormous. A low sensitivity with the nominal system does not guarantee robustness! +## Model-order reduction for uncertain models +The ``\nu``-gap metric is a measure of distance between models when they are used in a feedback loop. This metric has the nice property that a controller designed for a process ``P`` that achieves a normalized coprime factor margin ([`ncfmargin`](@ref)) of ``m``, will stabilize all models that are within a ``\nu``-gap distance of ``m`` from ``P``. This can be used to reduce the number of uncertain realizations for a model represented with `Particles` like above. Say that we have a plant model ``P`` +```@example MCM_NUGAP +using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements +ω = with_nominal(0.9 .. 1.1, 1) +ζ = with_nominal(0.5 ± 0.01, 0.5) +P = tf([ω^2], [1, 2*ζ*ω, ω^2]) |> ss +``` +represented by 2000 samples. We can compute the ``\nu``-gap metric between each realization in ``P`` and the nominal value (encoded using `with_nominal` above): +```@example MCM_NUGAP +gap = nugap(P) +``` +The worst-case gap is: +```@example MCM_NUGAP +pmaximum(gap) # p for "particle" maximum +``` +That means that if we design a controller for the nominal ``P`` without any uncertainty, and make sure that it achieves an [`ncfmargin`](@ref) of at least this value, it will stabilize all realizations in ``P``. + +We can also reduce the number of realizations in ``P`` by discarding those that are close in the ``\nu``-gap sense to the nominal value: +```@example MCM_NUGAP +Pr = nu_reduction(P, 0.1) +``` +here, all realizations that were within a ``\nu``-gap distance of 0.1 from the nominal value were discarded. [`nu_reduction`](@ref) usually reduces the number of realizations substantially. + +We can reduce the number of realizations even further using [`nu_reduction_recursive`](@ref): +```@example MCM_NUGAP +Gr = nu_reduction_recursive(P, 0.1) +``` +The algorithm used in [`nu_reduction_recursive`](@ref) has a worst-case complexity of ``O(N^2)`` where ``N`` is the number of realizations (particles) in ``P``, but this complexity is only problematic for small gaps and large number of realizations, say, ``\nu < 0.02`` and ``N > 50``. + +!!! warn "Stochastic interpretation" + If `P` has a stochastic interpretation, i.e., the coefficients come from some distribution, this interpretation will be lost after reduction, mean values and standard deviations will not be preserved. The reduced system should instead be interpreted as preserving worst-case uncertainty. + +The reduced model is useful for plotting and manual loop-shaping etc. ## Using the $M\Delta$ framework The examples above never bothered with things like the "structured singular value", $\mu$ or linear-fractional transforms. We do, however, provide some elementary support for this modeling framework. diff --git a/src/RobustAndOptimalControl.jl b/src/RobustAndOptimalControl.jl index 8a3b9604..3c74b810 100644 --- a/src/RobustAndOptimalControl.jl +++ b/src/RobustAndOptimalControl.jl @@ -79,4 +79,8 @@ export diskmargin, Diskmargin, Disk, sim_diskmargin, loop_diskmargin, structured include("diskmargin.jl") include("mimo_diskmargin.jl") + +export nu_reduction, nu_reduction_recursive +include("mcm_nugap.jl") + end diff --git a/src/mcm_nugap.jl b/src/mcm_nugap.jl new file mode 100644 index 00000000..344acdf5 --- /dev/null +++ b/src/mcm_nugap.jl @@ -0,0 +1,84 @@ +""" + nugap(G; map = map) + +Compute the νgap between the nominal system `Gₙ` represented by the first particle index in `G`, and all other systems in `G`. +Returns a `Particles` object with the νgap for each system in `G`. + +See `with_nominal` to endow uncertain values with a nominal value, and `nominal` to extract the nominal value. + +The value returned by this function, `νᵧ` is useful for robust synthesis, by designing a controller for the nominal system `Gₙ`, that achieves an [`ncfmargin`](@ref) of at least `νᵧ` is guaranteed to stabilize all realizations within `G`. + +To speed up computation for large systems, a threaded or distributed `map` function can be supplied, e.g., `ThreadTools.tmap` or `Distributed.pmap`. +""" +function RobustAndOptimalControl.nugap(G, i = 1; map = map) + ControlSystemsBase.numeric_type(G) <: MonteCarloMeasurements.AbstractParticles || + throw(ArgumentError("When calling nugap with a single argument, G must be an uncertain system with particles where the nominal value is the first particle.")) + Gn = RobustAndOptimalControl.sys_from_particles(G, i) + N = nparticles(G.A) + gaps = map(1:N) do i + Gi = RobustAndOptimalControl.sys_from_particles(G, i) + nugap(Gn, Gi)[1] + end + Particles(gaps) +end + + + +## + +""" + nu_reduction(G, g=0.1; gap = nugap(G)) + +Reduce the number of particles in an uncertain system `G` by removing all particles that are within the νgap `g` of the nominal system `Gₙ`. + +Note: If `G` has a stochastic interpretation, i.e., the coefficients come from some distribution, this interpretation will be lost after reduction, mean values and standard deviations will not be preserved. The reduced system should instead be interpreted as preserving worst-case uncertainty. + +If the `gap = nugap(G)` has already been precomputed, it can be supplied as an argument to avoid potentially costly recomputaiton. +""" +function nu_reduction(G, g=0.1; gap = nugap(G)) + nind = argmin(gap.particles) + keepinds = [nind; findall(gap.particles .> g)] # keep also nominal + syss = [RobustAndOptimalControl.sys_from_particles(G, i) for i in keepinds] + RobustAndOptimalControl.ss2particles(syss) +end + + +""" + nu_reduction_recursive(G, g = 0.1; gap = nugap(G), keepinds = Set{Int}(1), verbose = false) + +Find a νgap cover of balls of radius `g` (in the νgap metric) that contain all realizations in `G`. + +If the `gap = nugap(G)` has already been precomputed, it can be supplied as an argument to avoid potentially costly recomputaiton. If a manually computed `gap` is supplied, you must also supply `keepinds=Set{Int}(index)` where `index` is the index of the nominal system in `G` used to compute `gap`. + +The returned cover `Gr` is of the same type as `G`, but with a smaller number of particles. A controller designed for `Gr` that achieves a [`ncfmargin`](@ref) of at least `g` for all realizations in `Gr` will stabilize all realizations in the original `G`. The extreme case cover where `Gr = Gnominal` is a single realization only can be computed by calling `g = nugap(G, i)` where `i` is the index of the nominal system in `G`. + +# Arguments: +- `G`: An uncertain model in the form of a `StateSpace{TE, Particles}` (a multi-model). +- `g`: The radius of the balls in the νgap cover. +- `gap`: An optional precomputed gap +- `verbose`: Print progress +""" +function nu_reduction_recursive(G, g=0.1; gap = nugap(G), discardinds=Set{Int}(), keepinds=Set{Int}(1), verbose = false) + nparticles(gap) >= 10 && g <= 0.01 && @warn "This algorithm has a complexity of O(N²) for small `g`" maxlog=1 + discardinds = union(discardinds, Set(findall(gap.particles .< g))) + verbose && println("Covered $(length(discardinds)) particles") + maxind, _ = argmax(enumerate(gap.particles)) do (i,g) + if i ∈ discardinds + zero(g) + else + g + end + end + # @show maxind + if maxind ∈ discardinds + N = nparticles(G.A) + keepers = sort(unique([collect(keepinds); [i for i in 1:N if i ∉ discardinds]])) + syss = [RobustAndOptimalControl.sys_from_particles(G, i) for i in keepers] + return RobustAndOptimalControl.ss2particles(syss) + else + push!(keepinds, maxind) + end + gap = nugap(G, maxind) + nu_reduction_recursive(G, g; gap, discardinds, keepinds, verbose) +end + diff --git a/test/runtests.jl b/test/runtests.jl index 37523849..2ff24916 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -121,4 +121,9 @@ using Test include("../examples/uncertain.jl") end + @testset "mcm_nugap" begin + @info "Testing mcm_nugap" + include("test_mcm_nugap.jl") + end + end diff --git a/test/test_mcm_nugap.jl b/test/test_mcm_nugap.jl new file mode 100644 index 00000000..0f00d3eb --- /dev/null +++ b/test/test_mcm_nugap.jl @@ -0,0 +1,32 @@ +using ControlSystemsBase +using RobustAndOptimalControl +using MonteCarloMeasurements + +ω = with_nominal(0.9 .. 1.1, 1) +ζ = with_nominal(0.5 ± 0.01, 0.5) +G = tf([ω^2], [1, 2*ζ*ω, ω^2]) |> ss + +@test nominal(G) == ss(tf([1], [1, 2*0.5*1, 1])) + +gap = nugap(G) +@test gap isa typeof(ω) + +# plot(gap) + +# unsafe_comparisons(true) +# w = exp10.(LinRange(-2, 2, 200)) +# bodeplot(G) + +Gr = nu_reduction(G, 0.05) +@test 10 < nparticles(Gr.A) < nparticles(G.A) + +Gr = nu_reduction_recursive(G, 0.05; verbose=true) +@test 1 < nparticles(Gr.A) < nparticles(G.A) + +if isinteractive() + unsafe_comparisons(true) + bodeplot(G); bodeplot!(Gr, c=2) + # bodeplot(Gr, c=2) +end + +