add simple model reduction for uncertain models

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",

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.

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

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
+

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

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
+
+