Merge pull request #55 from JuliaControl/CS1

bump compat ControlSystems v1

Project.toml

@@ -1,7 +1,7 @@
 name = "RobustAndOptimalControl"
 uuid = "21fd56a4-db03-40ee-82ee-a87907bee541"
 authors = ["Fredrik Bagge Carlson", "Marcus Greiff"]
-version = "0.4.11"
+version = "0.4.12"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -25,7 +25,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ChainRulesCore = "1"
 ComponentArrays = "0.9, 0.10, 0.11, 0.12"
-ControlSystems = "0.12.16"
+ControlSystems = "1.0.1"
 DescriptorSystems = "1.2"
 Distributions = "0.25"
 GenericSchur = "0.5.2"

docs/src/uncertainty.md

@@ -314,7 +314,7 @@ Ps.Δ # Ps.delta also works
 We can evaluate the frequency response of $M$ and calculate the structured singular value $\mu$
 
 ```@example satellite
-M = freqresp(lft(Ps, -K).M, w).parent # -K to get negative feedback
+M = freqresp(lft(Ps, -K).M, w) # -K to get negative feedback
 μ = structured_singular_value(M)
 plot(w, μ, xscale=:log10)
 ```

src/ExtendedStateSpace.jl

@@ -305,6 +305,11 @@ function Base.:*(s1::ExtendedStateSpace, s2::ExtendedStateSpace)
 end
 
 function Base.:*(s1::ExtendedStateSpace, s2::Number)
+    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(s1)
+    ss(A, s2*B1, B2, C1, C2, s2*D11, D12, s2*D21, D22, s1.timeevol)
+end
+
+function Base.:*(s2::Number, s1::ExtendedStateSpace)
     A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(s1)
     ss(A, B1, B2, s2*C1, C2, s2*D11, s2*D12, D21, D22, s1.timeevol)
 end

src/diskmargin.jl

@@ -423,16 +423,12 @@ passivityplot
         s.ny == s.nu || throw(ArgumentError("passivity_index only defined for square systems"))
         Is = ss(I(s.ny), s.timeevol)
         G = (Is-s)feedback(Is, s)
-        sv = sigma(G, w)[1]
-        for j in 1:size(sv, 2)
-            for i in 1:size(sv, 3)
-                @series begin
-                    xscale --> :log10
-                    yscale --> :log10
-                    label --> "System $si"
-                    to1series(ws, sv)
-                end
-            end
+        sv = sigma(G, w)[1] |> permutedims
+        @series begin
+            xscale --> :log10
+            yscale --> :log10
+            label --> "System $si"
+            to1series(ws, sv)
         end
     end
     @series begin

src/mimo_diskmargin.jl

@@ -107,7 +107,7 @@ using Optim
 #     X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
 #     M = starprod(X,L)
 #     # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
-#     M0 = freqresp(M, w).parent
+#     M0 = freqresp(M, w)
 #     M0, D
 # end
 
@@ -123,7 +123,7 @@ using Optim
 #     X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
 #     M = starprod(X,L)
 #     # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
-#     M0 = freqresp(M, w).parent
+#     M0 = freqresp(M, w)
 #     M0, D
 # end
 
@@ -237,12 +237,12 @@ robstab(M0::UncertainSS, args...; kwargs...) = 1/norm(structured_singular_value(
 function structured_singular_value(M0::LTISystem, w::AbstractVector; kwargs...)
     if M0 isa UncertainSS
         if length(M0.Δ) == 1 && M0.Δ[] isa δDyn
-            return sigma(M0.M, w)[1][:, 1]
+            return sigma(M0.M, w)[1][1,:]
         end
         all(d isa δ{<:Complex} for d in M0.Δ) || error("Structured singular value only supported for diagonal complex perturbations")
         M0 = M0.M # not it's own method due to method ambiguity misery
     end
-    M = freqresp(M0, w).parent
+    M = freqresp(M0, w)
     μ = structured_singular_value(M; kwargs...)
 end
 
@@ -264,7 +264,7 @@ function sim_diskmargin(L::LTISystem,σ::Real,w::AbstractVector)
     n = L.ny
     X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
     M = starprod(X,L)
-    M0 = freqresp(M, w).parent
+    M0 = freqresp(M, w)
     mu = structured_singular_value(M0)
     imu = inv.(structured_singular_value(M0))
     simultaneous = [Diskmargin(imu, σ; ω0 = w, L) for (imu, w) in zip(imu,w)]

src/plotting.jl

@@ -69,14 +69,14 @@ specificationplot
             if ControlSystems._PlotScale == "dB"
                 singval = 20 * log10.(singval)
             end
-            for i = 1:min(size(singval, 2), nsigma)
+            for i = 1:min(size(singval, 1), nsigma)
                 @series begin
                     xscale --> :log10
                     yscale --> ControlSystems._PlotScaleFunc
                     linestyle --> :solid
                     color --> colors[mod(index - 1, 3)+1]
                     label --> (i == 1 ? s_labels[mod(index - 1, 3)+1] : "")
-                    w ./ (hz ? 2pi : 1), singval[:, i]
+                    w ./ (hz ? 2pi : 1), singval[i, :]
                 end
             end
         end
@@ -98,7 +98,7 @@ specificationplot
             if ControlSystems._PlotScale == "dB"
                 singval = 20 * log10.(singval)
             end
-            for i = 1:size(singval, 2)
+            for i = 1:size(singval, 1)
                 weightlabel = (i == 1 ? w_labels[mod(index - 1, 3)+1] : "")
                 @series begin
                     xscale --> :log10
@@ -107,7 +107,7 @@ specificationplot
                     color --> colors[mod(index - 1, 3)+1]
                     linewidth --> 2
                     label --> (i == 1 ? w_labels[mod(index - 1, 3)+1] : "")
-                    w ./ (hz ? 2pi : 1), singval[:, i]
+                    w ./ (hz ? 2pi : 1), singval[i, :]
                 end
             end
         end
@@ -158,7 +158,7 @@ mvnyquistplot
     framestyle --> :zerolines
     θ = range(0, stop=2π, length=100)
     S, C = sin.(θ), cos.(θ)
-    L = freqresp(sys, w).parent
+    L = freqresp(sys, w)
     dets = map(axes(L, 3)) do i
         det(I + L[:,:,i])
     end
@@ -230,7 +230,7 @@ muplot
     xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
     yguide --> "μ Singular Values $(ControlSystems._PlotScaleStr)"
     @views for (si, s) in enumerate(systems)
-        M = freqresp(s, w).parent
+        M = freqresp(s, w)
         mu, D = structured_singular_value(M, scalings=true, tol=1e-6, dynamic=true)
         Di = Diagonal(D)
         sv = map(axes(M, 3)) do i

src/uncertainty_interface.jl

@@ -67,9 +67,9 @@ end
 
 Base.:(+)(d2::δ, n::Number) = +(n, d2)
 
-function Base.:(*)(n::Number, d2::δ)
-    δ(*(n, d2.val), abs(n)*d2.radius, d2.name)
-end
+Base.:(*)(n::Number, d2::δ) = δ(*(n, d2.val), abs(n)*d2.radius, d2.name)
+
+Base.:(*)(d2::δ, n::Number) = δ(*(d2.val, n), abs(n)*d2.radius, d2.name)
 
 function Base.:(-)(d2::δ)
     δ(-d2.val, d2.radius, d2.name)
@@ -263,6 +263,11 @@ function Base.:*(s1::UncertainSS, n::Number)  # reverse case is handled by switc
     UncertainSS(sys, s1.Δ)
 end
 
+function Base.:*(n::Number, s1::UncertainSS)  # reverse case is handled by switching
+    sys = invert_mappings(n*invert_mappings(s1.sys))
+    UncertainSS(sys, s1.Δ)
+end
+
 
 function Base.:+(s1::UncertainSS, s2::UncertainSS) 
     sys1 = s1.sys

src/weights.jl

@@ -118,7 +118,7 @@ For each frequency in `w`, fit a circle in the complex plane that contains all m
 If `realtive = true`, circles encompassing `|(P - Pn)/Pn|` will be returned (multiplicative/relative uncertainty).
 If `realtive = false`, circles encompassing `|P - Pn|` will be returned (additive uncertainty).
 
-If `nominal = :mean`, the mean of `P` will be used as nominal model. If `nominal = :first`, the first particle will be used. See `MonteCarloMeasurements.with_nominal` to set the nominal value in the first particle. 
+If `nominal = :mean`, the mean of `P` will be used as nominal model. If `nominal = :first`, the first particle will be used. See `MonteCarloMeasurements.with_nominal` to set the nominal value in the first particle. If `nominal = :center`, the middle point `(pmaximum(ri)+pminimum(ri))/2` will be used.
 
 See also [`nyquistcircles`](@ref) to plot circles (only if relative=false).
 """
@@ -127,16 +127,32 @@ function fit_complex_perturbations(P, w; relative=true, nominal=:mean)
     r = freqresp(P, w)
     centers = Vector{Complex{eltype(w)}}(undef, length(w))
     radii = Vector{eltype(w)}(undef, length(w))
-    for i in axes(r,1)
-        ri = r[i,1,1] # NOTE: assumes SISO
-        c = nominal === :mean ? pmean(ri) : MonteCarloMeasurements.nominal(ri)
-        rad = abs(relative ? pmaximum((ri - c)/c) : pmaximum(ri - c))
+    for i in axes(r,3)
+        ri = r[1,1,i] # assumes SISO
+        c = if nominal === :mean
+            pmean(ri)
+        elseif nominal === :median
+            complex(pmedian(real(ri)), pmedian(imag(ri)))
+        elseif nominal === :center
+            center(ri)
+        elseif nominal ∈ (:nominal, :first)
+            MonteCarloMeasurements.nominal(ri)
+        else
+            error("Unknown nominal option $nominal")
+        end
+        rad = abs(relative ? pmaximum(abs((ri - c)/c)) : pmaximum(abs(ri - c)))
         centers[i] = c
         radii[i] = rad
     end
     centers, radii
 end
 
+function center(ri)
+    mir, mar = pextrema(real(ri))
+    mic, mac = pextrema(imag(ri))
+    complex((mar+mir)/2, (mac+mic)/2)
+end
+
 
 """
     nyquistcircles(w, centers, radii)

test/test_diskmargin.jl

@@ -191,15 +191,15 @@ mu = 1e5*[2.961061403927230
 0.000003489834242
 0.000003488286208]
 
-mu2 = structured_singular_value(permutedims(freqresp(L3, w), (2,3,1)))
+mu2 = structured_singular_value(freqresp(L3, w))
 @test mu2 ≈ mu[2:2:end]
 
 
 ##
 w = exp10.(LinRange(-2, 2, 300))
 K = ss(I(3), L3.timeevol)
 L = feedback(L3, K)
-mu = structured_singular_value(permutedims(freqresp(L, w), (2,3,1)))
+mu = structured_singular_value(freqresp(L, w))
 plot(mu)
 @test mu[1] ≈ 1 rtol=1e-3
 @test maximum(mu) ≈ 2.503148529400597 rtol=1e-3

test/test_lqg.jl

@@ -35,7 +35,7 @@ C = observer_controller(G)
 @test C.B == kalman(G)
 @test all(iszero, C.D)
 RobustAndOptimalControl.gangoffourplot(sys, C)
-RobustAndOptimalControl.gangofsevenplot(sys, C, tf(1, [1, 1]))
+RobustAndOptimalControl.gangofsevenplot(sys, C, tf(1, [1, 1]) .* I(2))
 
 Ce = extended_controller(G)
 @test Ce.ny == C.ny

test/test_uncertainty.jl

@@ -226,7 +226,7 @@ G = lft(P, -K)
 hn = norm(G, Inf)
 
 w = exp10.(LinRange(-5, 1, 100))
-M = freqresp(G.M, w).parent
+M = freqresp(G.M, w)
 # mu = mussv_DG(M)
 # maximum(mu)
 # # maximum(structured_singular_value(M))