From 9f416c0ece74c8ef655d32100a06345224e4eb6a Mon Sep 17 00:00:00 2001
From: baggepinnen <cont-frb@ulund.org>
Date: Fri, 7 Dec 2018 16:55:34 +0100
Subject: [PATCH] updates to julia1
---
REQUIRE | 1 +
src/BallAndBeam.jl | 31 +++++++++++++++++++------------
src/FRTN35_lab1.jl | 18 +++++++++---------
src/arx.jl | 30 +++++++++++++++---------------
4 files changed, 44 insertions(+), 36 deletions(-)
diff --git a/REQUIRE b/REQUIRE
index ac86564..bc087dc 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -5,3 +5,4 @@ Polynomials
ProgressMeter
JLD
Documenter
+DSP
diff --git a/src/BallAndBeam.jl b/src/BallAndBeam.jl
index 0cfa11a..b668b9a 100644
--- a/src/BallAndBeam.jl
+++ b/src/BallAndBeam.jl
@@ -1,3 +1,4 @@
+using Pkg
installed_packages = Pkg.installed()
if "LabProcesses" ∉ keys(installed_packages)
Pkg.clone("https://gitlab.control.lth.se/processes/LabProcesses.jl.git")
@@ -19,7 +20,7 @@ bopl, bopl!, nypl, nypl!, plot, fbdesign, ffdesign, opendoc
export SysFilter, run_control_2DOF # For documentation
-using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter
+using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter, DSP, LinearAlgebra, Statistics
include("arx.jl")
@@ -30,18 +31,18 @@ Perform fra-experiemnt For a single frequency `ω`. Called from inside `fra`
"""
function run_experiment(P::Beam, ω, duration, settling_time, amplitude)
h = sampletime(P)
- data = zeros(0:h:duration)
+ data = zeros(length(0:h:duration))
LabProcesses.initialize(P)
for t = 0:h:settling_time
@periodically h begin
- u = amplitude*sin(ω*t)
+ u = amplitude*sin(ω*t)
control(P, u + bias(P))
end
end
for (i,t) = enumerate(0:h:duration)
@periodically h begin
- data[i] = measure(P)
- u = amplitude*sin(ω*t)
+ data[i] = measure(P)
+ u = amplitude*sin(ω*t)
control(P, u + bias(P))
end
end
@@ -51,7 +52,7 @@ end
function run_experiment(P::BeamSimulator, ω, duration, settling_time, amplitude)
h = sampletime(P)
- data = zeros(0:h:duration)
+ data = zeros(length(0:h:duration))
LabProcesses.initialize(P)
for t = 0:h:settling_time
u = amplitude*sin(ω*t)
@@ -86,10 +87,10 @@ function fra(P::AbstractBeam, Ω::AbstractVector;
nbr_of_periods = 10,
amplitude = 1)
h = sampletime(P)
- G = zeros(Complex, size(Ω)) # Transfer function
+ G = zeros(Complex, size(Ω)) # Transfer function
@showprogress 1 "Running experiments" for (i,ω) = enumerate(Ω)
- T = 2π/ω # Period time
+ T = 2π/ω # Period time
settling_time_i = ceil(settling_time/T)*T # Settling time even number of periods
duration = nbr_of_periods*T
data = run_experiment(P,ω,duration,settling_time_i, amplitude)
@@ -139,10 +140,16 @@ end
@recipe function Nypl(p::Nypl)
G = p.args[1]
+ @series begin
+ linestyle := :dash
+ color := :black
+ α = LinRange(0,2π,100)
+ cos.(α),sin.(α)
+ end
xlabel --> "\$Re(G)\$"
ylabel --> "\$Im(G)\$"
title --> "Nyquist plot"
- real.(G[:,2]), imag.(G[:,2])
+ @series real.(G[:,2]), imag.(G[:,2])
end
@@ -175,7 +182,7 @@ function fbdesign(G::AbstractMatrix, polevect, zerovect, gain)
sysFB = ss(zpk(zerovect,polevect,gain*prod(abs.(polevect))/pzv))
C = Number[ω freqresp(sysFB, ω)]
L = Number[ω G[:,2].*C[:,2]]
- T = Number[ω L[:,2]./(1+L[:,2])]
+ T = Number[ω L[:,2]./(1 .+L[:,2])]
sysFB,L,T,C
end
@@ -201,8 +208,8 @@ poles2poly(poles) = reduce((r,l) -> conv(r,[1, l]), [1], poles) #|> reverse |> P
Build and launch the documentation website.
"""
function opendoc()
- include(Pkg.dir("BallAndBeam", "docs","make.jl"))
- docpath = Pkg.dir("BallAndBeam", "docs","build","index.html")
+ include(joinpath(dirname(pathof("BallAndBeam")), "docs","make.jl"))
+ docpath = joinpath(dirname(pathof("BallAndBeam")), "docs","build","index.html")
run(`xdg-open $docpath`)
end
diff --git a/src/FRTN35_lab1.jl b/src/FRTN35_lab1.jl
index 7d0d793..2eb3cb0 100644
--- a/src/FRTN35_lab1.jl
+++ b/src/FRTN35_lab1.jl
@@ -1,4 +1,4 @@
-using BallAndBeam, ControlSystems, JLD, LabProcesses, Plots
+using BallAndBeam, ControlSystems, LabProcesses, Plots, JLD
# @load "workspace.jld" # Run this command to restore a saved workspace
P = LabProcesses.Beam(bias = 0.) # Replace for BeamSimulator() to run simulations
@@ -11,15 +11,15 @@ 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 = exp10.(LinRange(log10(1),log10(300),8)) # exp10.(LinRange = logspace
G1 = fra(P, w1_100, amplitude=1, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
@save "workspace.jld"
-w1_200 = logspace(log10(50),log10(100),10)
+w1_200 = exp10.(LinRange(log10(50),log10(100),10))
G2 = fra(P, w1_200, amplitude=2, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
@save "workspace.jld"
-w1_300 = logspace(log10(100),log10(300),20)
+w1_300 = exp10.(LinRange(log10(100),log10(300),20))
G3 = fra(P, w1_300, amplitude=2, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
@save "workspace.jld"
@@ -35,7 +35,7 @@ nypl(G123, m=:star)
polevect = [-10]
zerovect = []
gain = 1
-sysFBc,L,T,C = fbdesign(G, polevect, zerovect, gain)
+sysFBc,L,T,C = fbdesign(G123, polevect, zerovect, gain)
polevect = [-10]
zerovect = []
@@ -55,8 +55,8 @@ 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(0:h:duration)
- u = zeros(0:h:duration)
+ 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
@@ -71,12 +71,12 @@ function prbs_experiment(;amplitude = 1, duration = 10)
y,u
end
y,u = prbs_experiment(duration = 10, amplitude = 1)
-plot([y u], lab=["y" "u"])
+plot([u y], lab=["u" "y"])
na = 4 # Order of A polynomial
nb = 2 # 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) # The exclamation mark (!)
+bodeconfidence!(arxtf, Σ, ω = exp10.(LinRange(0,3,200)), linecolor=:red) # The exclamation mark (!)
# adds to the previous plot. If there is no plot open, remove the !
diff --git a/src/arx.jl b/src/arx.jl
index 13017c2..f430a73 100644
--- a/src/arx.jl
+++ b/src/arx.jl
@@ -8,10 +8,10 @@ fit(y,yh) = 100 * (1-rms(y-yh)./rms(y-mean(y)));
aic(x,d) = log(sse(x)) + 2d/length(x)
"""
- toeplitz{T}(c::AbstractArray{T},r::AbstractArray{T})
+ toeplitz(c::AbstractArray,r::AbstractArray)
Returns a Toeplitz matrix where `c` is the first column and `r` is the first row.
"""
-function toeplitz{T}(c::AbstractVector{T},r::AbstractVector{T})
+function toeplitz(c::AbstractVector{T},r::AbstractVector{T}) where T
nc = length(c)
nr = length(r)
A = zeros(T, nc, nr)
@@ -38,13 +38,13 @@ B = [10,5] # B(z) coeffs
u = randn(100) # Simulate 100 time steps with Gaussian input
y = filt(B,A,u)
yr,A = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
-x = A\yr # Estimate model polynomials
+x = A\\yr # Estimate model polynomials
plot([yr A*x], lab=["Signal" "Prediction"])
```
For nonlinear ARX-models, see [BasisFunctionExpansions.jl](https://github.com/baggepinnen/BasisFunctionExpansions.jl/)
"""
-function getARXregressor(y::AbstractVector{Float64},u::AbstractVecOrMat{Float64}, na, nb)
- assert(length(nb) == size(u,2))
+function getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
+ @assert(length(nb) == size(u,2))
m = max(na+1,maximum(nb))
n = length(y) - m+1
offs = m-na-1
@@ -86,7 +86,7 @@ function arx(h,y::AbstractVector{Float64}, u::AbstractVector{Float64}, na, nb;
if λ == 0
w = A\y_train
else
- w = (A'A + λ*eye(size(A,2)))\A'y_train
+ w = (A'A + λ*I)\A'y_train
end
a,b = params2poly(w,na,nb)
model = tf(b,a,h)
@@ -116,27 +116,27 @@ function parameter_covariance(y_train, A, w, λ=0)
ATAλ\ATA/ATAλ
end
iATA = (iATA+iATA')/2
- Σ = σ²*iATA + sqrt(eps())*eye(iATA)
+ Σ = σ²*iATA + sqrt(eps())*I
end
"""
- bodeconfidence(arxtf::TransferFunction, Σ::Matrix; ω = logspace(0,3,200))
+ bodeconfidence(arxtf::TransferFunction, Σ::Matrix; ω = exp10.(LinRange(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))
+@recipe function BodeConfidence(p::BodeConfidence; ω = exp10.(LinRange(-2,3,200)))
arxtfm = p.args[1]
Σ = p.args[2]
- L = chol(Hermitian(Σ))
+ L = cholesky(Hermitian(Σ))
am, bm = -reverse(denpoly(arxtfm)[1].a[1:end-1]), reverse(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
+ w = L.L*randn(size(L,1)) .+ wm
a,b = params2poly(w,na,nb)
arxtf = tf(b,a,arxtfm.Ts)
mag, phase, _ = bode(arxtf, ω)
@@ -144,8 +144,8 @@ bodeconfidence
end
magmc = hcat(getindex.(res,1)...)
phasemc = hcat(getindex.(res,2)...)
- mag = mean(magmc,2)[:]
- phase = mean(phasemc,2)[:]
+ mag = mean(magmc,dims=2)[:]
+ phase = mean(phasemc,dims=2)[:]
uppermag = getpercentile(magmc,0.95)
lowermag = getpercentile(magmc,0.05)
upperphase = getpercentile(phasemc,0.95)
@@ -176,7 +176,7 @@ bodeconfidence
end
function getpercentile(mag,p)
- uppermag = mapslices(mag, 2) do magω
- sort(magω)[round(Int,endof(magω)*p)]
+ uppermag = mapslices(mag, dims=2) do magω
+ sort(magω)[round(Int,length(magω)*p)]
end
end
--
GitLab