diff --git a/src/BallAndBeam.jl b/src/BallAndBeam.jl
index 2d185fd95d9151d853c354674aa26dc53e641e91..2b399863a7c30642f63804e0c06ef2897c1d8522 100644
--- a/src/BallAndBeam.jl
+++ b/src/BallAndBeam.jl
@@ -109,30 +109,57 @@ function sortfqs(G)
G[perm,:]
end
-function bopl!(G; kwargs...)
- plot!(G[:,1], abs.(G[:,2]); ylabel="Magnitude", title="Bode plot", xscale=:log10, yscale=:log10, subplot=1, kwargs...)
+@userplot Bopl
+@userplot Nypl
+
+
+@recipe function Bopl(p::Bopl)
+ G = p.args[1]
+ link := :x
+ layout := (2,1)
+ @series begin
+ ylabel --> "Magnitude"
+ title --> "Bode plot"
+ xscale --> :log10
+ yscale --> :log10
+ subplot := 1
+ G[:,1], abs.(G[:,2])
+ end
phase = angle.(G[:,2]).*(180/π)
phase[phase .> 0] .-= 360
- plot!(G[:,1], phase; xlabel="Frequency [rad/s]", ylabel="Phase", xscale=:log10, subplot=2, kwargs...)
- plot!(link=:x)
+ @series begin
+ xlabel --> "Frequency [rad/s]"
+ ylabel --> "Phase"
+ xscale --> :log10
+ subplot := 2
+ G[:,1], phase
+ end
+ nothing
end
-nypl!(G; kwargs...) = plot!(real.(G[:,2]), imag.(G[:,2]);
- xlabel="\$Re(G)\$", ylabel="\$Im(G)\$", title="Nyquist plot", kwargs...)
+@recipe function Nypl(p::Nypl)
+ G = p.args[1]
+ xlabel --> "\$Re(G)\$"
+ ylabel --> "\$Im(G)\$"
+ title --> "Nyquist plot"
+ real.(G[:,2]), imag.(G[:,2])
+end
+
+
"""
- bopl(G)
+ bopl(G; kwargs...)
Plot a bodeplot when `G` is a matrix with frequencies in first column and complex magnitudes
in the second. Add to an existing plot with `bopl!(...)` See also [`nypl`](@ref)
"""
-bopl(G; kwargs...) = (plot(layout=(2,1),link=:x);bopl!(G;kwargs...))
+bopl
"""
- nypl(G)
+ nypl(G; kwargs...)
Plot a nyquistplot when `G` is a matrix with frequencies in first column and complex magnitudes
in the second. Add to an existing plot with `nypl!(...)` See also [`bopl`](@ref)
"""
-nypl(G; kwargs...) = (plot();nypl!(G;kwargs...))
+nypl
"""
diff --git a/src/FRTN35_lab1.jl b/src/FRTN35_lab1.jl
index 2d4860e3718eefe9278470ba9db62b3c652d6a08..75472ae1b56dc6cebdcbd09bcdd99fdab69188dd 100644
--- a/src/FRTN35_lab1.jl
+++ b/src/FRTN35_lab1.jl
@@ -1,8 +1,8 @@
-using BallAndBeam, ControlSystems, JLD, LabProcesses
+using BallAndBeam, ControlSystems, JLD, LabProcesses, Plots
# @load "workspace.jld" # Run this command to restore a saved workspace
Bias = 0.01 # Change this if your process drifts over time
-P = LabProcesses.Beamr(0.01, Bias) # Replace for BeamSimulator(0.01) to run simulations
+P = LabProcesses.Beam(0.01, Bias) # Replace for BeamSimulator(0.01) to run simulations
h = sampletime(P)
settling_time = 1
@@ -12,7 +12,7 @@ nbr_of_periods = 3
# 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)
+w1_100 = logspace(log10(1),log10(300),800)
G1 = fra(P, w1_100, amplitude=1, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
@save "workspace.jld"
@@ -62,21 +62,19 @@ function prbs_experiment(;amplitude = 1, duration = 10)
for (i,t) = enumerate(0:h:duration)
@periodically h begin
y[i] = measure(P)
- u[i] = amplitude*(prbs()-0.5)
+ u[i] = sum(sin(ω*(t-1)) for ω in logspace(-1,3,10))#amplitude*(prbs()-0.5)# randn()
control(P, u[i])
end
end
LabProcesses.finalize(P)
- y,u,times
+ y,u
end
y,u = prbs_experiment(duration = 10, amplitude = 1)
-plot([u y])
+plot([y u], lab=["y" "u"])
-na = 3 # Order of A polynomial
+na = 4 # Order of A polynomial
nb = 2 # Order of B polynomial
-arxtf = arx(h, y, u, na, nb) # Estimate trasfer function with ARX method
+arxtf, Σ = arx(h, y, u, na, nb) # Estimate trasfer function with ARX method
-mag, phase, ω = bode(arxtf, logspace(-1,3,200))
bopl(G123, lab="Measured transfer function")
-plot!(ω, mag[:], subplot=1, lab = "ARX estimate")
-plot!(ω, phase[:], subplot=2)
+bodeconfidence!(arxtf, Σ, ω = logspace(0,3,200))
diff --git a/src/arx.jl b/src/arx.jl
index 5eecff1bb446b67e8cb6300696eccd17358670f8..bb6dc0eca7418666104d669345bb96be3978b2f4 100644
--- a/src/arx.jl
+++ b/src/arx.jl
@@ -1,4 +1,4 @@
-export toeplitz, getARXregressor, find_na, arx
+export toeplitz, getARXregressor, find_na, arx, bopl_confidence, bopl_confidence!
## Helper functions
@@ -60,7 +60,7 @@ end
"""
find_na(y::AbstractVector,n::Int)
-Plots the RMSE and AIC for model orders up to `n`. Useful for model selection
+Plots the RMSE and AIC For model orders up to `n`. Useful for model selection
"""
function find_na(y::AbstractVector,n::Int)
error = zeros(n,2)
@@ -75,23 +75,108 @@ function find_na(y::AbstractVector,n::Int)
end
"""
- Gtf = arx(h,y, u, na, nb; λ = 0)
-Fit a transfer function to data using an ARX model.
+ Gtf, Σ = arx(h,y, u, na, nb; λ = 0)
+Fit a transfer Function to data using an ARX model.
`nb` and `na` are the orders of the numerator and denominator polynomials.
"""
function arx(h,y::AbstractVector{Float64}, u::AbstractVector{Float64}, na, nb; λ = 0)
na -= 1
y_train, A = getARXregressor(y,u, na, nb)
- n_points = size(A,1)
if λ == 0
w = A\y_train
else
w = (A'A + λ*eye(size(A,2)))\A'y_train
end
+ a,b = params2poly(w,na,nb)
+ model = tf(b,a,h)
+ Σ = parameter_covariance(y_train, A, w, λ)
+ return model, Σ
+end
+"""
+ a,b = params2poly(params,na,nb)
+"""
+function params2poly(w,na,nb)
a = [1; -w[1:na]]
b = w[na+1:end]
- model = tf(b,a,h)
- return model
+ a,b
+end
+
+"""
+ Σ = parameter_covariance(y_train, A, w, λ=0)
+"""
+function parameter_covariance(y_train, A, w, λ=0)
+ σ² = var(y_train .- A*w)
+ iATA = if λ == 0
+ inv(A'A)
+ else
+ ATA = A'A
+ ATAλ = factorize(ATA + λ*I)
+ ATAλ\ATA/ATAλ
+ end
+ iATA = (iATA+iATA')/2
+ Σ = σ²*iATA + sqrt(eps())*eye(iATA)
+end
+
+"""
+ bodeconfidence(arxtf::TransferFunction, Σ::Matrix; ω = logspace(0,3,200))
+Plot a bode diagram of a transfer function estimated with [`arx`](@ref) with confidence bounds on magnitude and phase.
+"""
+bodeconfidence
+
+@userplot BodeConfidence
+
+@recipe function BodeConfidence(p::BodeConfidence; ω = logspace(-2,3,200))
+ arxtfm = p.args[1]
+ Σ = p.args[2]
+ L = chol(Hermitian(Σ))
+ am, bm = -denpoly(arxtfm)[1].a[2:end], arxtfm.matrix[1].num.a
+ wm = [am; bm]
+ na,nb = length(am), length(bm)
+ mc = 100
+ res = map(1:mc) do _
+ w = L'randn(size(L,1)) .+ wm
+ a,b = params2poly(w,na,nb)
+ arxtf = tf(b,a,arxtfm.Ts)
+ mag, phase, _ = bode(arxtf, ω)
+ mag[:], phase[:]
+ end
+ magmc = hcat(getindex.(res,1)...)
+ phasemc = hcat(getindex.(res,2)...)
+ mag = mean(magmc,2)[:]
+ phase = mean(phasemc,2)[:]
+ uppermag = getpercentile(magmc,0.95)
+ lowermag = getpercentile(magmc,0.05)
+ upperphase = getpercentile(phasemc,0.95)
+ lowerphase = getpercentile(phasemc,0.05)
+
+ layout := (2,1)
+
+ @series begin
+ subplot := 1
+ title --> "ARX estimate"
+ ylabel --> "Magnitude"
+ ribbon := (lowermag, uppermag)
+ yscale --> :log10
+ xscale --> :log10
+ alpha --> 0.3
+ ω, mag
+ end
+ @series begin
+ subplot := 2
+ ribbon := (lowerphase, upperphase)
+ ylabel --> "Phase [deg]"
+ xlabel --> "Frequency [rad/s]"
+ xscale --> :log10
+ alpha --> 0.3
+ ω, phase
+ end
+ nothing
+end
+
+function getpercentile(mag,p)
+ uppermag = mapslices(mag, 2) do magω
+ sort(magω)[round(Int,endof(magω)*p)]
+ end
end