WIP: Robust toolbox dev

.gitignore

@@ -1,2 +1,3 @@
 docs/build/
 docs/site/
+.DS_Store

docs/make.jl

@@ -1,4 +1,4 @@
-using Documenter, ControlSystems, Plots, LinearAlgebra, DSP
+using Documenter, ControlSystems, Plots
 import GR # Bug with world age in Plots.jl, see https://github.com/JuliaPlots/Plots.jl/issues/1047
 
 include("src/makeplots.jl")
@@ -23,6 +23,6 @@ end
 deploydocs(
     repo = "github.com/JuliaControl/ControlSystems.jl.git",
     latest = "master",
-    julia = "1.0",
+    julia = "0.6",
     deps = myDeps
 )

docs/src/makeplots.jl

@@ -16,7 +16,7 @@ Q = Matrix{Float64}(I,2,2)
 R = Matrix{Float64}(I,1,1)
 L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used
 
-u(x,t)  = -L*x .+ 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5
+u(x,t)  = -L*x + 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5
 t=0:h:5
 x0 = [1,0]
 y, t, x, uout = lsim(sys,u,t,x0=x0)
@@ -71,7 +71,7 @@ Plots.savefig(plotsDir*"/ppgofplot.svg")
 
 P1(s) = exp(-sqrt(s))
 f1 = stabregionPID(P1,exp10.(range(-5, stop=1, length=1000))); Plots.savefig(plotsDir*"/stab1.svg")
-P2 = s -> 100*(s+6).^2 ./(s.*(s+1).^2 .*(s+50).^2)
+P2 = s -> 100*(s+6).^2./(s.*(s+1).^2.*(s+50).^2)
 f2 = stabregionPID(P2,exp10.(range(-5, stop=2, length=1000))); Plots.savefig(plotsDir*"/stab2.svg")
 P3 = tf(1,[1,1])^4
 f3 = stabregionPID(P3,exp10.(range(-5, stop=0, length=1000))); Plots.savefig(plotsDir*"/stab3.svg")
@@ -84,15 +84,15 @@ f3 = stabregionPID(P3,exp10.(range(-5, stop=0, length=1000))); Plots.savefig(plo
 P = tf([1.],[1., 1])
 ζ = 0.5 # Desired damping
 ws = exp10.(range(-1, stop=2, length=8)) # A vector of closed-loop bandwidths
-kp = 2*ζ*ws.-1 # Simple pole placement with PI given the closed-loop bandwidth, the poles are placed in a butterworth pattern
+kp = 2*ζ*ws-1 # Simple pole placement with PI given the closed-loop bandwidth, the poles are placed in a butterworth pattern
 ki = ws.^2
 pidplots(P,:nyquist,;kps=kp,kis=ki, ω= exp10.(range(-2, stop=2, length=500)))
 Plots.savefig(plotsDir*"/pidplotsnyquist1.svg")
 pidplots(P,:gof,;kps=kp,kis=ki, ω= exp10.(range(-2, stop=2, length=500)))
 Plots.savefig(plotsDir*"/pidplotgof1.svg")
 
 kp = range(-1, stop=1, length=8) # Now try a different strategy, where we have specified a gain crossover frequency of 0.1 rad/s
-ki = sqrt.(1 .-kp.^2)./10
+ki = sqrt.(1-kp.^2)./10
 pidplots(P,:nyquist,;kps=kp,kis=ki)
 Plots.savefig(plotsDir*"/pidplotsnyquist2.svg")
 pidplots(P,:gof,;kps=kp,kis=ki)

example/hinf_example_DC.jl

@@ -0,0 +1,64 @@
+using Plots
+using ControlSystems
+"""
+This is a simple SISO example with a pole in the origin, corresponding to the
+DC servos used in the Lund laboratories. It serves to exeplify how the syntheis
+can be done for simple SISO systems, and also demonstrates how we chan verify
+if the problem is feasible to solve using the ARE method.
+
+The example can be set to visualize and save plots using the variables
+  MakePlots - true/false (true if plots are to be generated, false for testing)
+  SavePlots - true/false (true if plots are to be saved, false for testing)
+"""
+MakePlots, SavePlots = false, false
+
+# Define the process
+Gtrue   = tf([11.2], [1, 0.12,0])
+
+# Sensitivity weight function
+M, wB, A = 1.5, 20.0, 1e-8
+WS = tf([1/M, wB],[1, wB*A])
+
+# Output sensitivity weight function
+WU = ss(1)
+
+# Complementary sensitivity weight function
+WT = []
+
+# Form the P in the LFT F_l(P,C) as a partitioned state-space object
+P = hInf_partition(Gtrue, WS, WU, WT)
+
+# Check if the system is feasible for synythesis
+flag = hInf_assumptions(P)
+
+# Since it is not, modify the plant desciption
+epsilon = 1e-5
+G = tf([11.2], [1, 0.12]) * tf([1], [1, epsilon])
+
+# Form the P in the LFT Fl(P,C) as a partitioned state-space object
+P = hInf_partition(G, WS, WU, WT)
+
+# Check if the problem is feasible
+flag = hInf_assumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, C, gamma = hInf_synthesize(P)
+
+# Extract the transfer functions defining some signals of interest
+Pcl, S, CS, T = hInf_signals(P, G, C)
+
+## Plot the specifications
+if MakePlots
+  specificationplot([S, CS, T], [WS, WU, WT], gamma)
+  if SavePlots
+    savefig("example_DC_specifications.pdf")
+  end
+end
+
+## Plot the closed loop gain from w to z
+if MakePlots
+  specificationplot(Pcl, gamma; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
+  if SavePlots
+    savefig("example_DC_clgain.pdf")
+  end
+end

example/hinf_example_DC_discrete.jl

@@ -0,0 +1,60 @@
+using Plots
+using ControlSystems
+"""
+This is a simple SISO example with a pole in the origin, corresponding to the
+DC servos used in the Lund laboratories. It serves to exeplify how the syntheis
+can be done for simple SISO systems, and also demonstrates how we chan verify
+if the problem is feasible to solve using the ARE method.
+
+The example can be set to visualize and save plots using the variables
+  MakePlots - true/false (true if plots are to be generated, false for testing)
+  SavePlots - true/false (true if plots are to be saved, false for testing)
+"""
+MakePlots = true
+
+# Define the process
+ts = 0.01
+epsilon = 1e-5
+Gd = ss(c2d(tf([11.2], [1, 0.12]) * tf([1], [1, epsilon]), ts))
+
+# Sensitivity weight function
+M, wB, A = 1.5, 20.0, 1e-8
+WS = tf([1/M, wB],[1, wB*A])
+
+# Output sensitivity weight function
+WU = ss(1)
+
+# Complementary sensitivity weight function
+WT = []
+
+# Create continuous time approximation of the process
+Gc = hInf_bilinear_z2s(ss(Gd))
+
+# Form the P in the LFT Fl(P,C) as a partitioned state-space object
+Pc = hInf_partition(Gc, WS, WU, WT)
+
+# Check if the problem is feasible
+flag = hInf_assumptions(Pc)
+
+# Synthesize the H-infinity optimal controller
+flag, Cc, gamma = hInf_synthesize(Pc)
+
+# Extract the transfer functions defining some signals of interest, but do so
+# using discrete equivalent of the continuous time objects Pc, Cc and Gc
+PclD, SD, CSD, TD = hInf_signals(
+  hInf_bilinear_s2z(Pc, ts),
+  hInf_bilinear_s2z(Gc, ts),
+  hInf_bilinear_s2z(Cc, ts)
+)
+
+# This solution is a bit hacky and should be revised
+Pcl = ss(PclD.A, PclD.B, PclD.C, PclD.D, ts)
+S   = ss(SD.A, SD.B, SD.C, SD.D, ts)
+CS  = ss(CSD.A, CSD.B, CSD.C, CSD.D, ts)
+T   = ss(TD.A, TD.B, TD.C, TD.D, ts)
+
+# Visualize results
+if MakePlots
+  specificationplot([S, CS, T], [WS, WU, WT], gamma)
+  specificationplot(Pcl, gamma; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
+end

example/hinf_example_MIMO.jl

@@ -0,0 +1,109 @@
+using ControlSystems
+using Plots
+
+"""
+This is very dense and horrible example used for debugging. CUrrently doesn't
+work due to the absence of enforcing a loopshift in h´the solver - I need to
+figure out how to doe this prolperly.
+"""
+
+# Define the process dynamics
+nx = 30
+nu = 8
+ny = 8
+A = 10*rand(nx,nx)
+A -= (abs(maximum(real(eigvals(A))))+0.1) * Matrix{Float64}(I, size(A,1), size(A,2))
+B = rand(nx,nu)
+C = rand(ny,nx)
+D = rand(ny,nu)
+
+G = hInf_bilinear_z2s(hInf_bilinear_s2z(ss(A,B,C,D), 0.005))
+
+# Sensitivity weight function
+WSe = tf([0.1,1],[1000,1])
+#WS = [WSe 0 0 ; 0 WSe 0 ; 0 0 WSe]
+WS = [WSe 0 0 0 ; 0 WSe 0 0 ; 0 0 WSe 0 ; 0 0 0 WSe]
+WS = [WS 0*WS; WS*0 WS]
+
+# Output sensitivity weight function
+WUe = 0.0001*tf([1,1],[0.1,1])
+WU = [WUe 0 0 0 ; 0 WUe 0 0 ; 0 0 WUe 0 ; 0 0 0 WUe]
+WU = [WU 0*WU; WU*0 WU]
+
+#WU = [WUe 0 0 ; 0 WUe 0 ; 0 0 WUe ]
+# Complementary sensitivity weight function
+WTe = tf([10,1],[1,1])
+WT = [WTe 0 0 0 ; 0 WTe 0 0 ; 0 0 WTe 0 ; 0 0 0 WTe]
+WT = [WT 0*WT; WT*0 WT]
+
+# Form augmented P dynamics in state-space
+P = hInf_partition(G, WS, WU, WT)
+
+# Check that the assumptions are satisfied
+flag = hInf_assumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, K, gamma = hInf_synthesize(P; maxIter=20)
+
+Pcl, S, KS, T = hInf_signals(P, G, K)
+
+MakePlots = true
+SavePlots = false
+
+
+fmin = -8
+fmax =  8
+# Plot the singular values  of the system
+f = [10^i for i in range(fmin, stop=fmax, length=10001)]
+
+valPcl  = sigma(Pcl, f)[1];
+valS    = sigma(S, f)[1];
+valKS   = sigma(KS, f)[1];
+valT    = sigma(T, f)[1];
+
+# Visualize the close loop gain
+pGain = plot(f, valPcl[:,1], xscale = :log10, yscale = :log10, color = :black, w=2, label="\$\\sigma(P_{w\\rightarrow z}(j\\omega))\$")
+plot!(f, gamma*ones(size(f)), xscale = :log10, yscale = :log10, color = :black, w=3, style=:dot, label="\$\\gamma\$")
+if size(valPcl,2) > 1
+  for ii = 2:size(valPcl,2)
+    plot!(f, valPcl[:,ii], xscale = :log10, yscale = :log10, color = :black, w=1)
+  end
+end
+
+# Plot the sensitivity functions of interest
+pSigma = plot(f, valS[:,1], xscale = :log10, yscale = :log10, color = :red, w=3, label="\$\\sigma(S(j\\omega))\$")
+if size(valS,2) > 1
+  for ii = 2:size(valS,2)
+    plot!(f, valS[:,ii], xscale = :log10, yscale = :log10, color = :red, w=1, label="")
+  end
+end
+plot!(f, valKS[:,1], xscale = :log10, yscale = :log10, color = :green, w=3, label="\$\\sigma(C(j\\omega)S(j\\omega)))\$")
+if size(valKS,2) > 1
+  for ii = 2:size(valKS,2)
+    plot!(f, valKS[:,ii], xscale = :log10, yscale = :log10, color = :green, w=1, label="")
+  end
+end
+plot!(f, valT[:,1], xscale = :log10, yscale = :log10, color = :blue, w=3, label="\$\\sigma(T(j\\omega)))\$")
+if size(valT,2) > 1
+  for ii = 2:size(valT,2)
+    plot!(f, valT[:,ii], xscale = :log10, yscale = :log10, color = :blue, w=1, label="")
+  end
+end
+
+# Visualize the weighting functions
+if isa(WS, LTISystem) || isa(WS, Number)
+  valiWS = sigma(gamma/WS[1,1], f)[1]
+  plot!(f, valiWS, xscale = :log10, yscale = :log10, color = :red, w=2, style=:dot, label="\$\\gamma\\sigma(W_S(j\\omega)^{-1})\$")
+end
+if isa(WU, LTISystem) || isa(WU, Number)
+  valiWU = sigma(gamma/WU[1,1], f)[1]
+  plot!(f, valiWU, xscale = :log10, yscale = :log10, color = :green, w=2, style=:dot, label="\$\\gamma\\sigma(W_U(j\\omega)^{-1})\$")
+end
+if isa(WT, LTISystem) || isa(WT, Number)
+  valiWT = sigma(gamma/WT[1,1], f)[1]
+  plot!(f, valiWT, xscale = :log10, yscale = :log10, color = :blue, w=2, style=:dot, label="\$\\gamma\\sigma(W_T(j\\omega)^{-1})\$")
+end
+
+l = @layout [ a ; b ]
+plt=plot(pGain, pSigma, layout=l, size=(600,600))
+gui()

example/hinf_example_MIT.jl

@@ -0,0 +1,58 @@
+using ControlSystems
+using Plots
+
+"""
+This is a simple SISO example which was used for debugging the implementation,
+as this exact example was use in the lecture notes of the "Principles of Optimal
+Control" cours of the MIT OpenCourseWare [1], where the choice of weighting
+functions and dynamics gave rise to an H-infinity optimal cotnroller with a
+γ of approximately 1.36, where, in our case, we get a controller at 0.93
+
+[1] https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-323-principles-of-optimal-control-spring-2008/lecture-notes/lec15.pdf
+
+The example can be set to visualize and save plots using the two variables
+  ShowPlots - true/false (true if plots are to be generated, false for testing)
+  filename  - Set to string if files are to be saved, otherwise set a empty list
+"""
+MakePlots, SavePlots = false, false
+
+# Define the process
+G   = tf([200], [0.025,1.0025,10.1,1])
+
+# Sensitivity weight function
+M, wB, A = 1.5, 10, 1e-4
+WS = tf([1/M, wB],[1, wB*A])
+
+# Output sensitivity weight function
+WU = ss(0.1)
+
+# Complementary sensitivity weight function
+WT = []
+
+# Form augmented P dynamics in state-space
+P = hInf_partition(G, WS, WU, WT)
+
+# Check that the assumptions are satisfied
+flag = hInf_assumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, C, gamma = hInf_synthesize(P)
+
+# Extract the transfer functions defining some signals of interest
+Pcl, S, CS, T = hInf_signals(P, G, C)
+
+## Plot the specifications
+if MakePlots
+  specificationplot([S, CS, T], [ss(WS), WU, WT], gamma)
+  if SavePlots
+    savefig("example_MIT_specifications.pdf")
+  end
+end
+
+## Plot the closed loop gain from w to z
+if MakePlots
+  specificationplot(Pcl, gamma; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
+  if SavePlots
+    savefig("example_MIT_clgain.pdf")
+  end
+end