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 = hinfpartition(Gtrue, WS, WU, WT)
+
+# Check if the system is feasible for synythesis
+flag = hinfassumptions(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 = hinfpartition(G, WS, WU, WT)
+
+# Check if the problem is feasible
+flag = hinfassumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, C, gamma = hinfsynthesize(P)
+
+# Extract the transfer functions defining some signals of interest
+Pcl, S, CS, T = hinfsignals(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 = bilineard2c(ss(Gd))
+
+# Form the P in the LFT Fl(P,C) as a partitioned state-space object
+Pc = hinfpartition(Gc, WS, WU, WT)
+
+# Check if the problem is feasible
+flag = hinfassumptions(Pc)
+
+# Synthesize the H-infinity optimal controller
+flag, Cc, gamma = hinfsynthesize(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 = hinfsignals(
+  bilinearc2d(Pc, ts),
+  bilinearc2d(Gc, ts),
+  bilinearc2d(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_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 = hinfpartition(G, WS, WU, WT)
+
+# Check that the assumptions are satisfied
+flag = hinfassumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, C, gamma = hinfsynthesize(P)
+
+# Extract the transfer functions defining some signals of interest
+Pcl, S, CS, T = hinfsignals(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

example/hinf_example_MIT_discrete.jl

@@ -0,0 +1,60 @@
+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], however, this example is in discrete time
+
+[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
+  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.005
+Gd = hInf_bilinear_s2z(ss(tf([200], [0.025,1.0025,10.1,1])),ts)
+
+# 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 = []
+
+# Create continuous time approximation of the process
+Gc = bilineard2c(ss(Gd))
+
+# Form augmented P dynamics in state-space
+Pc = hinfpartition(Gc, WS, WU, WT)
+
+# Check that the assumptions are satisfied
+flag = hinfassumptions(Pc)
+
+# Synthesize the H-infinity optimal controller
+flag, Cc, gamma = hinfsynthesize(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 = hinfsignals(
+  bilinearc2d(Pc, ts),
+  bilinearc2d(Gc, ts),
+  bilinearc2d(Cc, ts)
+)
+
+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_tank.jl

@@ -0,0 +1,99 @@
+using ControlSystems
+using Plots
+using LinearAlgebra
+"""
+This is a simple SISO example with integrator dynamics corresponding to the
+quad tank process in the lab.
+
+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 proces parameters
+k1, k2, kc, g = 3.33, 3.35, 0.5, 981
+A1, A3, A2, A4 = 28, 28, 32, 32
+a1, a3, a2, a4= 0.071, 0.071, 0.057, 0.057
+h01, h02, h03, h04 = 12.4, 12.7, 1.8, 1.4
+T1, T2 = (A1/a1)*sqrt(2*h01/g), (A2/a2)*sqrt(2*h02/g)
+T3, T4 = (A3/a3)*sqrt(2*h03/g), (A4/a4)*sqrt(2*h04/g)
+c1, c2 = (T1*k1*kc/A1), (T2*k2*kc/A2)
+gamma1, gamma2 = 0.7, 0.6
+
+# Define the process dynamics
+A = [-1/T1     0 A3/(A1*T3)          0;
+     0     -1/T2          0 A4/(A2*T4);
+     0         0      -1/T3          0;
+     0         0          0      -1/T4];
+B = [gamma1*k1/A1     0;
+     0                gamma2*k2/A2;
+     0                (1-gamma2)*k2/A3;
+     (1-gamma1)*k1/A4 0              ];
+C = [kc 0 0 0;
+     0 kc 0 0];
+D = zeros(2,2)
+G = ss(A,B,C,D);
+
+# Sensitivity weight function
+WSelement = 100*tf([1,1],[1000,1])
+WS = [WSelement 0; 0 WSelement]
+iWSelement = 1/WSelement
+iWS = [iWSelement 0; 0 iWSelement]
+
+# Output sensitivity weight function
+WUelement = 5*tf([1,1],[0.1,1])
+WUelement = ss(0.01)
+WU = [WUelement 0; 0 WUelement]
+iWUelement = 1/WUelement
+iWU = [iWUelement 0; 0 iWUelement]
+
+# Complementary sensitivity weight function
+WTelement = tf([10,0.1],[1,1])
+WT  = [WTelement 0; 0 WTelement]
+iWTelement = 1/WTelement
+iWT = [iWTelement 0; 0 iWTelement]
+
+# Form augmented P dynamics in state-space
+P = hinfpartition(G, WS, WU, WT)
+
+# Check that the assumptions are satisfied
+flag = hinfassumptions(P)
+
+# Synthesize the H-infinity optimal controller
+flag, C, gamma = hinfsynthesize(P)
+
+Pcl, S, CS, T = hinfsignals(P, G, C)
+
+## Plot the specifications
+# TODO figure out why I get segmentation errors when using ss instead of tf for
+# the weighting functions, makes no sense at all
+if MakePlots
+  specificationplot([S, CS, T], [WSelement, 0.01, WTelement], gamma)
+  if SavePlots
+    savefig("example_tank_specifications.pdf")
+  end
+end
+
+## Plot the closed loop gain from w to z
+# TODO figure out why the legends don't seem to work in this case
+if MakePlots
+  specificationplot(Pcl, gamma; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
+  if SavePlots
+    savefig("example_tank_clgain.pdf")
+  end
+end
+
+if MakePlots
+  times = [i for i in range(0, stop=300, length=10000)]
+  stepy, stept, stepx = step(T, times)
+  pStep1=plot(stept, stepy[:,1,1], color = :blue, w=2, label="\$u_1\\rightarrow y_1\$")
+  pStep2=plot(stept, stepy[:,1,2], color = :blue, w=2, label="\$u_1\\rightarrow y_2\$", ylims = (-0.5,1.1))
+  pStep3=plot(stept, stepy[:,2,1], color = :blue, w=2, label="\$u_2\\rightarrow y_1\$", ylims = (-0.5,1.1))
+  pStep4=plot(stept, stepy[:,2,2], color = :blue, w=2, label="\$u_2\\rightarrow y_2\$")
+  l = @layout [ a b c d ]
+  plt=plot(pStep1, pStep2, pStep3, pStep4, layout=l, size=(1000,250))
+  if SavePlots
+    savefig("example_tank_stepresponse.pdf")
+  end
+end