using BallAndBeam, ControlSystems, JLD2, LabProcesses, Plots # @load "workspace.jld" # Run this command to restore a saved workspace # All biases between -0.05 and 0.05, perform experiments to find out yours! P = LabProcesses.Beam(bias = 0.0) # Replace for BeamSimulator() to run simulations h = sampletime(P) settling_time = 1 nbr_of_periods = 3 # Below we define some frequency vectors (using logarithmic spacing) # and run three experiments. You may modify the freqency vectors # any way you want and add/remove experiments as needed. w1_100 = logspace(log10(1),log10(300),8) G1 = fra(P, w1_100, amplitude=1.0, nbr_of_periods=nbr_of_periods, settling_time=settling_time) w1_200 = logspace(log10(50),log10(100),10) G2 = fra(P, w1_200, amplitude=1.0, nbr_of_periods=nbr_of_periods, settling_time=settling_time) w1_300 = logspace(log10(100),log10(300),20) G3 = fra(P, w1_300, amplitude=1.0, nbr_of_periods=nbr_of_periods, settling_time=settling_time) # Concatenate (overlapping) estimates to be used and sort based on freq G123 = sortfqs([G1; G2; G3]) #gr() #to plot in a window plotly() # to plot in a browser bopl(G123, m=:star) nypl(G123, m=:star) #@save "workspace.jld" #= ## Control ================================================================================== # Specify poles, zeros and gain for your feedback controller and feedforward compensator # The input `G` to `fbdesign` is the data you have measured, e.g, `G123` #@load "workspace.jld" polevect = [-10] zerovect = [] gain = 1 sysFBc,L,T,C = fbdesign(G123, polevect, zerovect, gain) polevect = [-10] zerovect = [] gain = 1 sysFFc,YR,FF = ffdesign(T, polevect, zerovect, gain)G1 = fra(P, w1_100, amplitude=1.0 nbr_of_periods=nbr_of_periods, settling_time=settling_time) #gr() #to plot in a window plotly() # to plot in a browser bopl(C, lab="Controller") bopl!(L, lab="Open-loop system r->y without FF") bopl!(FF, lab="Feedforward compensator") bopl!(YR, lab="Closed-loop system r->y") sysFB,sysFF = c2d(sysFBc,h)[1],c2d(sysFFc,h)[1] y,u,r = run_control_2DOF(P, sysFB, sysFF, duration=5, reference = t->2sign(sin(2π/3*t))) plot([y u r], lab = ["y" "u" "r"]) =# #= # If you have time, estimate ARX model ===================================================== prbs = PRBSGenerator() function prbs_experiment(;amplitude = 1, duration = 10) y = zeros(length(0:h:duration)) u = zeros(length(0:h:duration)) LabProcesses.initialize(P) for (i,t) = enumerate(0:h:duration) @periodically h begin y[i] = measure(P) u[i] = amplitude*(prbs()-0.5) #u[i] = amplitude*randn() #u[i] = amplitude*sum(sin(ω*(t-1)) for ω in logspace(-1,3,10)) control(P, u[i] + bias(P)) end end LabProcesses.finalize(P) y,u end y,u = prbs_experiment(duration = 10, amplitude = 1.0) plot([y u], lab=["y" "u"]) na = 6 # Order of A polynomial nb = 5 # Order of B polynomial arxtf, Σ = arx(h, y, u, na, nb) # Estimate transfer function with ARX method bopl(G123, lab="Measured transfer function") bodeconfidence(arxtf, Σ, ω = logspace(0,3,200), linecolor=:red, lab="ARX") =# # The exclamation mark (!) # adds to the previous plot. If there is no plot open, remove the !