From c7d40de88815c75cf731350735c4a996ecf00e3d Mon Sep 17 00:00:00 2001 From: Mladen Gibanica Date: Sat, 13 Feb 2021 16:36:17 +0100 Subject: [PATCH 1/3] doc space fix --- julia/src/estimate_d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/julia/src/estimate_d.jl b/julia/src/estimate_d.jl index aae5a62..42b70a1 100644 --- a/julia/src/estimate_d.jl +++ b/julia/src/estimate_d.jl @@ -33,7 +33,7 @@ Parameters: * `type`: data type of model either `Real` or `Complex` * `estimd`: if set to false no d matrix is esimated and a zero d matrix is returned -Returns : +Returns: * `c`: LS-optimal c matrix * `d`: LS-optimal d matrix """ From d4fd9a6be0d796001038e66b4596b13e55e40010 Mon Sep 17 00:00:00 2001 From: Mladen Gibanica Date: Sat, 13 Feb 2021 18:04:45 +0100 Subject: [PATCH 2/3] add special matrices --- julia/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/julia/Project.toml b/julia/Project.toml index 76fa440..5a8e429 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -8,6 +8,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +SpecialMatrices = "928aab9d-ef52-54ac-8ca1-acd7ca42c160" [compat] julia = "1" From 6b738f1cb4cd3f990db65abfdbb7b08555343c9e Mon Sep 17 00:00:00 2001 From: Mladen Gibanica Date: Mon, 12 Jul 2021 00:04:57 +0200 Subject: [PATCH 3/3] Add lrm, moebius and markov_kung algorithms There are some numerical probelms with lrm. Many samples are needed to pass tests. --- julia/examples/examples.jl | 198 +++++++++++---------- julia/src/fsid.jl | 7 +- julia/src/math.jl | 19 ++ julia/src/misc.jl | 343 ++++++++++++++++++++++++++++++++++--- julia/test/runtests.jl | 210 ++++++++++++++++------- 5 files changed, 586 insertions(+), 191 deletions(-) diff --git a/julia/examples/examples.jl b/julia/examples/examples.jl index 15e91a6..38ca347 100644 --- a/julia/examples/examples.jl +++ b/julia/examples/examples.jl @@ -1,7 +1,7 @@ using LinearAlgebra, FFTW, Plots -using fsid: +using fsid: estimate_cd, estimate_bd, transpose_ffdata, @@ -13,6 +13,7 @@ using fsid: ffsid, bilinear_c2d, bilinear_d2c, + lrm, lsim, fdsim, cf2df, @@ -48,14 +49,14 @@ function example_ffsid_rank() fde = fsid.fresp(z, Ae, Be, Ce, De) err = computerr(ffdata, fde) println("||H - He||/||H|| = ", err) - + z2 = vcat(z, conj.(z)) ffdata2 = vcat(ffdata, conj.(ffdata)) Ae1, Be1, Ce1, De1, s = fsid.gffsid(z2, ffdata2, n, 2*n, type = Complex, estimd = true) fde = fsid.fresp(z2, Ae1, Be1, Ce1, De1) err = computerr(ffdata2, fde) println("||H - He||/||H|| = ", err) - + #return nothing A, B, C, D = buildss(n, m, p) @@ -65,7 +66,7 @@ function example_ffsid_rank() w = 2pi * fset wexp = ztrans(fset) fd = fsid.fresp(wexp, A, B, C, D) - + # estimate ss model from ffdata Ae, Be, Ce, De, s = fsid.ffsid(w, fd, n, n*2) @@ -95,7 +96,7 @@ function example_ffsid() w = 2pi * fset wexp = ztrans(fset) fd = fsid.fresp(wexp, A, B, C, D) - + # estimate ss model from ffdata Ae, Be, Ce, De, s = fsid.ffsid(w, fd, n, n*2) @@ -116,8 +117,8 @@ function example_fdsid() p = 4 N = 100 - A, B, C, D = buildss(n, m, p) - + A, B, C, D = buildss(n, m, p) + # create a random stable DT system lam = eigen(A).values rho = maximum(abs.(lam)) @@ -125,13 +126,13 @@ function example_fdsid() # random excitation signal u = randn(N, m) - - # time domain simulation + + # time domain simulation y = fsid.lsim((A, B, C, D), u) - + #plot(y) - ## Crfeate the N point DFT of the signals + ## Crfeate the N point DFT of the signals yf = fft(y[1], 1) uf = fft(u, 1) fset = buildfset(N) @@ -141,7 +142,7 @@ function example_fdsid() fddata = (w, yf, uf) Ae, Be, Ce, De, xt, s = fsid.fdsid(fddata, n, 2*n, estTrans = true) - println("Singular vales = $s") + # println("Singular vales = $s") fde = fsid.fresp(wexp, Ae, Be, Ce, De) fd = fsid.fresp(wexp, A, B, C, D) @@ -149,7 +150,7 @@ function example_fdsid() println("||G - Ge||/||G|| = ", computerr(fd, fde)) Ae1, Be1, Ce1, De1, xt1, s = fsid.fdsid(fddata, n, 2*n, estTrans = false) - println("Singular vales = $s, model order n = $n") + # println("Singular vales = $s, model order n = $n") fde = fsid.fresp(wexp, Ae, Be, Ce, De) fd = fsid.fresp(wexp, A, B, C, D) @@ -168,22 +169,22 @@ function example_1() p = 4 N = 100 - A, B, C, D = buildss(n, m, p) - + A, B, C, D = buildss(n, m, p) + # create a random stable DT system lam = eigen(A).values rho = maximum(abs.(lam)) A = A/rho/1.01 - + # random excitation signal u = randn(N, m) - # time domain simulation + # time domain simulation y = fsid.lsim((A, B, C, D), u) - + #plt.plot(y) - ## Crfeate the N point DFT of the signals + ## Crfeate the N point DFT of the signals yf = fft(y[1], 1) uf = fft(u, 1) fset = buildfset(N) @@ -194,7 +195,7 @@ function example_1() # Test estimation of B D and Xt Be, De, xt, resid2 = fsid.fdestim_bd(wexp, yf, uf, A, C, estTrans = true, type = Real) - println("fdestim_bd and xt residual = ", resid2) + # println("fdestim_bd and xt residual = ", resid2) ## Check that the frequency functions coincide fde = fsid.fresp(wexp, A, Be, C, De) @@ -206,19 +207,19 @@ function example_1() yy = fsid.fdsim((A,B,C,D), uf, wexp, xt) Be, De, xte, resid1 = fsid.fdestim_bd(wexp, yy[1], uf, A, C; estTrans = true) - println("fdestim_bd residual = ", resid1) + # println("fdestim_bd residual = ", resid1) println("fdestim_bd residual xt = ", computerr(xt, xte)) #Ce, De, resid1 = fsid.fdestim_cd(wexp, yy[1], uf, A, B, xt; estTrans = true) Ce, De, resid1 = fsid.fdestim_cd(wexp, yy[1], uf, A, B, xt) # Be, De, resid1 = fdestim_bd(wexp, yy, uf, A, C, estTrans = false) - println("fdestim_cd residual = ", resid1) + # println("fdestim_cd residual = ", resid1) fddata = (w, yf, uf) Ae, Be, Ce, De, xt, s = fsid.fdsid(fddata, n, 2n, estTrans = true) # Ae, Be, Ce, De, xt, s = fsid.fdsid(fddata, n, 2*n, estTrans = false) - println("Singular vales = ", s) + # println("Singular vales = ", s) fde = fsid.fresp(wexp, Ae, Be, Ce, De) fd = fsid.fresp(wexp, A, B, C, D) @@ -237,8 +238,8 @@ function example_2() p = 4 N = 100 - A, B, C, D = buildss(n, m, p) - + A, B, C, D = buildss(n, m, p) + # create a random stable DT system lam = eigen(A).values rho = maximum(abs.(lam)) @@ -247,7 +248,7 @@ function example_2() u = randn(N, m) y = fsid.lsim((A, B, C, D), u) - + plot(y) yf = fft(y[1], 1) @@ -257,10 +258,10 @@ function example_2() wexp = ztrans(fset) fddata = (w, yf, uf) - + Be, De, xt, resid2 = fsid.fdestim_bd(wexp, yf, uf, A, C; estTrans = true, type = Real) - println("fdestim_bd residual = ", resid2) - + # println("fdestim_bd residual = ", resid2) + fde = fsid.fresp(wexp, A, Be, C, De) fd = fsid.fresp(wexp, A, B, C, D) println("||H - He||/||H|| = ", computerr(fd, fde)) @@ -270,18 +271,18 @@ function example_2() yy = fsid.fdsim((A,B,C,D), uf, wexp, xt) Be, De, xt, resid1 = fsid.fdestim_bd(wexp, yy[1], uf, A, C; estTrans = true) # Be, De, resid1 = fdestim_bd(wexp, yy, uf, A, C, estTrans = false) - println("fdestim_bd residual = ", resid1) + # println("fdestim_bd residual = ", resid1) #Ce, De, resid1 = fsid.fdestim_cd(wexp, yy, uf, A, B, xt; estTrans = true) Ce, De, resid1 = fsid.fdestim_cd(wexp, yy[1], uf, A, B, xt) # Be, De, resid1 = fdestim_bd(wexp, yy, uf, A, C, estTrans = false) - println("fdestim_cd residual = ", resid1) + # println("fdestim_cd residual = ", resid1) fddata = (w, yf, uf) Ae, Be, Ce, De, xt, s = fsid.fdsid(fddata, n, 2*n, estTrans = true) # Ae, Be, Ce, De, s = fdsid(fddata, n, 2*n, estTrans = false) - println("Singular vales = ", s) + # println("Singular vales = ", s) fde = fsid.fresp(wexp, Ae, Be, Ce, De) fd = fsid.fresp(wexp, A, B, C, D) @@ -309,23 +310,23 @@ function example_ct_ffsid() w = 2pi * fset # wexp = ztrans(fset) fd = fsid.fresp(im*w, A, B, C, D) - + wd = fsid.cf2df(w,T) - + # estimate ss model from ffdata Ae, Be, Ce, De, s = fsid.ffsid(wd, fd, n, n*2) # convert to CT Aec, Bec, Cec, Dec = fsid.bilinear_d2c((Ae, Be, Ce, De), T) - + # frequency response of estimated model fde = fsid.fresp(im*w, Aec, Bec, Cec, Dec) println("CT FF: ||G - Ge||/||G|| = ", computerr(fd, fde)) - + Ae, Be, Ce, De, s = fsid.ffsid(w, fd, n, n*2, CT = true, T = T, estimd = false) fde2 = fsid.fresp(im*w, Ae, Be, Ce, De) println("CT2 FF: ||G - Ge||/||G|| = ", computerr(fd, fde2)) - println(De) + # println(De) end @@ -343,33 +344,89 @@ function example_ct_fdsid() A, B, C, D = buildss(n, m, p) D = D*0 - println(D) + # println(D) # Create frequency function fset = buildfset(N) w = 2pi * fset # wexp = ztrans(fset) fd = fsid.fresp(im*w, A, B, C, D) uf = randn(N, m) + im*randn(N, m) - + yf = fsid.fdsim((A ,B ,C ,D), uf, im*w) - + wd = fsid.cf2df(w,T) - + # estimate ss model from ffdata Ae, Be, Ce, De, xt, s = fsid.fdsid((wd, yf[1], uf), n, 2n; estTrans = false, estimd = true) # convert to CT Aec, Bec, Cec, Dec = fsid.bilinear_d2c((Ae, Be, Ce, De), T) - + # frequency response of estimated model fde = fsid.fresp(im*w, Aec, Bec, Cec, Dec) println("CT FD: ||G - Ge||/||G|| = ", computerr(fd, fde)) - println(Dec) - + # println(Dec) + Ae, Be, Ce, De, xt, s = fsid.fdsid((w, yf[1], uf), n, 2n; CT = true, T = T, estimd = false) fde2 = fsid.fresp(im*w, Ae, Be, Ce, De) println("CT2 FD: ||G - Ge||/||G|| = ", computerr(fd, fde2)) - println(De) + # println(De) +end + + +""" + example_ffsid_lrm() + + +""" +function example_ffsid_lrm() + N = 2000 + nmpset = [(4, 1, 1), (1, 1, 1), (2, 4, 12)] + nmpset = [(10, 5, 5)] + for (n, m, p) in nmpset + A, B, C, D = buildss(n, m, p) + + # create a random stable DT system + lam = eigen(A).values + rho = maximum(abs.(lam)) + A = A/rho/1.01 + + fset = buildfset(N) + z = ztrans(fset) + fd = fresp(z, A, B, C, D) + u = randn(N, m) + y, x = lsim((A, B, C, D), u, type = Real) + ff = lrm(u, y, n = 2, Nw = 30) + # yf = np.fft.fft(y, axis=0) + # uf = np.fft.fft(u, axis=0) + + println("||fd - ff||/||fd|", computerr(fd, ff)) + #plot(20log10.(abs.(fd[:, 1, 1]))) + #plot(20log10.(abs.(ff[:, 1, 1]))) + #plot(20log10.(abs.(ff[:, 1, 1] - fd[:, 1, 1]))) + end + + N = 2000 + A = [0.9999 1; 0 .9] + B = [0; 1][:, :] + C = [1 0] + D = [0][:, :] + u = randn(N, 1) + y, x = lsim((A, B, C, D), u; type = Real) + + fset = buildfset(N) + z = ztrans(fset) + fd = fresp(z, A, B, C, D) + u = randn(N, 1) + y, x = lsim((A, B, C, D), u; type = Real) + ff = lrm(u, y; n=1, Nw=3) + # yf = np.fft.fft(y, axis=0) + # uf = np.fft.fft(u, axis=0) + + println("||fd - ff||/||fd|", computerr(fd, ff)) + #plot(20log10.(abs.(fd[1:20, 1, 1]))) + #plot(20log10.(abs.(ff[1:20, 1, 1]))) + #plot(20log10.(abs.(ff[1:20, 1, 1] - fd[1:20, 1, 1]))) end @@ -386,54 +443,5 @@ function examples() example_ct_ffsid() example_ct_fdsid() example_ffsid_rank() - - - # need to add lrm for this - # N = 1000 - # nmpset = [(4, 1, 1), (1, 1, 1), (2, 4, 12)] - # nmpset = [(10, 5, 5)] - # for (n, m, p) in nmpset - # A, B, C, D = buildss(n, m, p) - - # # create a random stable DT system - # lam = eigen(A).values - # rho = maximum(abs.(lam)) - # A = A/rho/1.01 - - # fset = buildfset(N) - # z = ztrans(fset) - # fd = fsid.fresp(z, A, B, C, D) - # u = randn(N, m) - # y = fsid.lsim((A, B, C, D), u, type = Real) - # ff = fsid.lrm(u, y, n = 2, Nw = 30) - # # yf = np.fft.fft(y, axis=0) - # # uf = np.fft.fft(u, axis=0) - - # println("||fd - ff||/||fd|", computerr(fd, ff)) - # plot(20log10(abs.(fd[:, 1, 1]))) - # plot(20log10(abs.(ff[:, 1, 1]))) - # plot(20log10(abs.(ff[:, 1, 1] - fd[:, 1, 1]))) - # end - - # N = 100 - # A = [0.9999, 1; 0, .9] - # B = [0; 1] - # C = [1, 0] - # D = [0] - # u = randn(N, 1) - # y = fsid.lsim((A, B, C, D), u, type = Real) - - # fset = buildfset(N) - # z = ztrans(fset) - # fd = fsid.fresp(z, A, B, C, D) - # u = randn(N, 1) - # y = fsid.lsim((A, B, C, D), u, type = Real) - # ff = fsid.lrm(u, y, n=1,Nw=3) - # # yf = np.fft.fft(y, axis=0) - # # uf = np.fft.fft(u, axis=0) - - # println("||fd - ff||/||fd|", computerr(fd, ff)) - # plot(20log10(abs.(fd[1:20, 1, 1]))) - # plot(20log10(abs.(ff[1:20, 1, 1]))) - # plot(20log10(abs.(ff[1:20, 1, 1] - fd[1:20, 1, 1]))) + example_ffsid_lrm() end \ No newline at end of file diff --git a/julia/src/fsid.jl b/julia/src/fsid.jl index 7127673..983d3c0 100644 --- a/julia/src/fsid.jl +++ b/julia/src/fsid.jl @@ -1,7 +1,6 @@ module fsid -using Documenter -using LinearAlgebra +using LinearAlgebra, FFTW, SpecialMatrices, Documenter include("types.jl") include("math.jl") @@ -10,6 +9,4 @@ include("estimate_d.jl") include("misc.jl") include("n4sid.jl") -end # module - - +end diff --git a/julia/src/math.jl b/julia/src/math.jl index c2cdb1a..6dcf1a0 100644 --- a/julia/src/math.jl +++ b/julia/src/math.jl @@ -62,4 +62,23 @@ function optimal_print(s, i, j) optimal_print(s, s[i, j] + 1, j) print(")") end +end + + +@doc raw""" + vander(x::AbstractVector{T}, n = length(x), rev = false) + +Create Vandermonde matrix. +""" +function vander(x::AbstractVector{T}, n = length(x), rev = false) where T<:Number + m = length(x) + V = Matrix{T}(undef, m, n) + + for i in 1:m + for j in 1:n + V[i, j] = x[i]^(j-1) + end + end + + return rev ? reverse(V, dims=2) : V end \ No newline at end of file diff --git a/julia/src/misc.jl b/julia/src/misc.jl index dc6f098..25f702e 100644 --- a/julia/src/misc.jl +++ b/julia/src/misc.jl @@ -1,3 +1,294 @@ +@doc raw""" + lrm(u, y, n = 2, Nw = nothing) + + +""" +function lrm(u::AbstractArray, y::AbstractArray; n = 2, Nw = nothing) + us = size(u) + if length(us) == 1 + m = 1 + u = reshape(u, (us[1], 1)) + else + m = us[2] + end + Nu = us[1] + + ys = size(y) + if length(ys) == 1 + p = 1 + y = reshape(y, (ys[1], 1)) + else + p = ys[2] + end + Ny = ys[1] + + if Nu != Ny + println("lrm: Incorrect number of rows in u and y.") + return false + end + + if Nw === nothing + Nw = 2Int(ceil(((m + 2)*n + 1)/2)) + end + + if Nw < Int(ceil(((m + 2)*n + 1)/2)) + println("Error: Nw too small.") + return false + end + + yf = fft(y, 1) + uf = fft(u, 1) + ff = zeros(Complex{typeof(first(u))}, Ny, p, m) + iset = -Nw:Nw + R = vander(iset, n+1, true) + yfe = vcat(yf[end-Nw+1:end, :], yf, yf[1:Nw, :]) + ufe = vcat(uf[end-Nw+1:end, :], uf, uf[1:Nw, :]) + + for i = 1:Ny + for pidx = 1:p + yy = yfe[Nw .+ i .+ iset, pidx] + RR = hcat(R, diagm(-yy) * R[:, 1:n]) + + for midx = m:-1:1 + RR = hcat(RR, diagm(ufe[Nw .+ i .+ iset, midx]) * R) + end + + ht = RR\yy + ht = reverse(ht; dims = 1) + ff[i, pidx, :] = ht[1:(n + 1):(n + 1)*m] + end + end + + return ff +end + +@doc raw""" + kung_realization(mp::AbstractArray{T}, n::Int, q::Int = 0) where T <: Number + +Calulates the state-space realization `(a, b, c)` from Markov parameters using Kung's relization algorithm. + +Parameters: +* `m`: array of Marov parameters `m[i, j - 1, k - 1]` hold row `j` and column `k` of the Markov parameter of sample `i` +* `n`: model order state-space system +Optional: +* `q`: number of rows in Hankel matrix. The default is `n + 1` + +Returns: +* `a`: system matrix +* `b`: input matrix +* `c`: output matrix +""" +function kung_realization(mp::AbstractArray{T}, n::Int, q::Int = n + 1) where T <: Number + N = size(mp, 1) + p = size(mp, 2) + m = size(mp, 3) + ncol = N - q + 1 + + if ncol < n + println("n, q and N are not compatible. N >= q + n - 1 must be satisfied") + return false + end + + H = Matrix{T}(undef, p*q, ncol*m) + for j = 1:ncol + for i = 1:q + H[p*(i - 1) + 1:i*p, m*(j - 1) + 1:j*m] = mp[i + j - 1, :, :] + end + end + + u, s, vh = svd(H) + c = u[1:p, 1:n] + lh = u[1:p*(q-1), 1:n] + rh = u[p+1:end, 1:n] + a = lh\rh + b = diagm(s[1:n]) * vh[1:m, 1:n]' + + return a, b, c, H +end + +@doc raw""" + markov(sys::Tuple, N::Int) + +Calculate markov parameters from `sys = (a,b,c)`. + +Parameters: +* `sys`: `(a, b, c)` state-space matrices +* `N`: numer of Markov parameters to generate + +Returns: +* `mp`: `mp[i, :, :]` is Markov parameter `C(A^i)B` +""" +function markov(sys::Tuple, N::Int) + a, b, c = sys + n = size(a, 1) + m = size(b, 2) + p = size(c, 1) + mp = Array{typeof(first(a))}(undef, N, p, m) # change this to function acting on type + + aa = Matrix{typeof(first(a))}(I, n, n) # and this too + for i in 1:N + mp[i, :, :] = c * (aa * b) + aa = aa * a + end + + return mp +end + + +@doc raw""" + make_sys_real(sys::Tuple) + +Convert realization sys into a real-valued realization. + +Parameters: +* `sys`: `(a, b, c)` or `(a, b, c, d)` + +Returns: +* `sysr`: `(ar, br, cr, dr)` the realization with real valued matrices +""" +function make_sys_real(sys::Tuple) + a, b, c = sys[1:3] + n = size(a, 1) + mp = markov(sys, 2n) + mpr = real.(mp) + a, b, c = kung_realization(mpr, n) + + if length(sys) == 3 + return (a, b, c) + else + return (a, b, c, real.(sys[4])) + end +end + + +@doc raw""" + make_obs_real(a::Matrix{<:Number}, c::Matrix{<:Number}) + +Convert `(a, c)`` into real-valued matrices by approximating a real valued observability range space to the original one. + +Parameters: +* `a`: complex matrix +* `b`: complex matrix + +Returns: +* `ar`: real valued matrix +* `br`: real valued matrix +""" +function make_obs_real(a::Matrix{<:Number}, c::Matrix{<:Number}) + + n = size(a, 1) + p = size(c, 1) + obs = Matrix{Complex}(undef, p*(n + 1), n) + obs[1:p, :] = c + + for i in 1:n + obs[p*(i + 1):p*(i + 2), :] = obs[p*i:p*(i+1), :] * a + end + + obsr = hcat(real.(obs), imag.(obs)) + u, s, vh = svd(obsr) + c = u[1:p, 1:n] + lh = u[1:p*n, 1:n] + rh = u[p:end, 1:n] + lsres = lh\rh + a = lsres[1] + + return a, c +end + + +""" + moebius(sys; alpha = 1, beta = 0, gamma = 0, delta = 1) + +Calculates the bilinear transformation D->C for ss-system sys. +""" +function moebius(sys; alpha = 1, beta = 0, gamma = 0, delta = 1) + a, b, c, d = sys + n = size(a, 1) + ainv = inv(alpha * I(n) - gamma*a) + ac = (delta * a - beta * I(n)) * ainv + bc = (alpha * delta - gamma * beta) * ainv * b + cc = c * ainv + dc = d + gamma * multidot(c, ainv, b) + + return ac, bc, cc, dc +end + + +""" + moebius_arg(z; alpha = 1, beta = 0, gamma = 0, delta = 1) + + +""" +function moebius_arg(z; alpha = 1, beta = 0, gamma = 0, delta = 1) + nz = size(z, 1) + s = Vector{contype(Complex)}(undef, nz) + + for idx = 1:nz + s[idx] = (alpha*z[idx] + beta)/(gamma*z[idx] + delta) + end + + return s +end + + +""" + moebius_inv(sys; alpha = 1, beta = 0, gamma = 0, delta = 1) + +Calculates the bilinear transformation D->C for ss-system sys. +""" +function moebius_inv(sys; alpha = 1, beta = 0, gamma = 0, delta = 1) + a, b, c, d = sys + n = size(a, 1) + ainv = inv(delta * I(n) + gamma*a) + ac = (alpha * a +beta * I(n)) * ainv + bc = ainv * b + cc = -(gamma * beta - alpha * delta) * (c * ainv) + dc = d - gamma * multidot(c, ainv, b) + + return ac, bc, cc, dc +end + + +""" + moebius_arg_inv(s; alpha = 1, beta = 0, gamma = 0, delta = 1) + + +""" +function moebius_arg_inv(s; alpha = 1, beta = 0, gamma = 0, delta = 1) + nz = size(s, 1) + z = Vector{contype(Complex)}(undef, nz) + + for idx = 1:nz + z[idx]= (beta - delta*s[idx])/(gamma*s[idx] - alpha) + end + + return z +end + + +""" + uq_cond(z, q) + +""" +function uq_cond(z, q) + m = 1 + nw = size(z, 1) + u = Matrix{contype(Complex)}(undef, m*q, nw*m) + + for widx = 1:nw + u[1:m, (widx - 1)*m + 1:widx*m] = Matrix(I, m, m) + zx = z[widx] + for qidx = 1:q + u[(qidx-1)*m + 1:qidx*m, (widx-1)*m + 1:widx*m] = zx*Matrix(I, m ,m) + zx *= z[widx] + end + end + + return cond(u) +end + + @doc raw""" function fresp_slow!( frsp::AbstractArray{<:Number, 3}, @@ -84,7 +375,7 @@ Frequency response of ss-model `(a,b,c,d)` a rational matrix function ``fresp[i,:,:] = d+c*inv(z[i]*I-a)*b`` `fresp[i,:,:]` is the function value of the rational matrix function evaluated at `z[i]` - + Parameters: * `z`: vector with samples of the function argument * `a`: matrix @@ -143,7 +434,7 @@ function ltifr_slow!( ) n, m = size(b) nw = length(z) - + for widx in 1:nw fkern[widx, :, :] = (I * z[widx] - a)\b end @@ -178,18 +469,18 @@ function ltifr_fast!( bb = eig.vectors\b for widx in 1:nw da = o * z[widx] - eig.values - fkern[widx, :, :] = multidot(eig.vectors, diagm(0 => (1 ./ da)), bb) + fkern[widx, :, :] = multidot(eig.vectors, diagm(0 => (1 ./ da)), bb) end return nothing end -""" +""" function ltifr( a::AbstractMatrix{<:Number}, b::AbstractMatrix{<:Number}, z::AbstractVector{<:Number}, - noWarning::Bool = true + noWarning::Bool = true ) Calculates the frequency kernel @@ -203,20 +494,20 @@ Parameters: * `noWarning`: if true suppress information message Returns: -* `fkern`: frequency response +* `fkern`: frequency response """ function ltifr( a::AbstractMatrix{<:Number}, b::AbstractMatrix{<:Number}, z::AbstractVector{<:Number}, - noWarning::Bool = true + noWarning::Bool = true ) n, m = size(b) nw = length(z) fkern = zeros(ComplexF64, nw, n, m) eig = eigen(copy(a)) - + if rank(eig.vectors) < n if !noWarning println("Matrix a defective. Eigenvectors do not form a basis. Using slow mode.") @@ -262,7 +553,7 @@ function ltifd_slow!( eyen = Matrix{Float64}(I, n, n) for widx in 1:nw - fkern[:, widx] = (eyen * z[widx] - a)\(b * u[widx, :]) + fkern[:, widx] = (eyen * z[widx] - a)\(b * u[widx, :]) end return nothing @@ -307,7 +598,7 @@ function ltifd_fast!( for widx in 1:nw da = o * z[widx] - eig.values fkern[:, widx] = multidot(eig.vectors, diagm(0 => (1 ./ da)), bb * u[widx, :]) - end + end return nothing end @@ -351,7 +642,7 @@ function ltifd( if rank(eig.vectors) < n if !noWarning - print("Matrix a defective. Eigenvectors do not form a basis. Using slow mode.") + println("Matrix a defective. Eigenvectors do not form a basis. Using slow mode.") end ltifd_slow!(fkern, a, b, u, z) else @@ -365,13 +656,13 @@ end ffdata2fddata(ffdata::AbstractArray{<:Number, 3}, z::AbstractVector{<:Number}) Converts `ffdata` to `fddata`. - + Converts frequency function data `ffdata` to input/output data format `fddata`. Parameters: * `ffdata`: frequency function data in the format such that `ffdata[i,:,:]` corresponds to the frequency function matrix of size `(p,m)` at frequency index `i` corresponding to function argument `z[i]` at a total number of samples `nz` * `z`: array with the corresponding frequency function argument of size `nz` - + Returns: * `u`: Fourier transform of input of size `(m*nz, m)` * `y`: Fourier transform of output of size `(m*nz, p)` @@ -406,7 +697,7 @@ end ) Calculates the time domain input to state respone. - + Calculates the time domain state response ``x[i+1,:] = a*x[i,:]) + b*u[i,:]`` @@ -416,7 +707,7 @@ Parmeters: * `b`: a matrix of size `(n, m)` * `u`: an array of input vectors such that `u[i,:]` is the input vector of size `m` at time index `i`. The array has size `(N, m)` * `x0`: intial vector of size `n`, i.e. `x[0,:] = x0`. Default value is the zero vector - + Returns: * `x`: the resulting state-sequence of size `(N, n)` * `x[k,:]`: is the state at sample `k` @@ -462,7 +753,7 @@ Parameters: * `u`: an array of input vectors such that `u[i,:]` is the input vector of size m at time index `i`. The array has size `(N, m)` *Optional:* * `x0`: intial vector of size n, i.e. `x[0,:] = x0`. Default value is the zero vector - + Returns: * `y`: the resulting output sequence of size (N, p) * `x`: the resulting state sequence of size (N, n) @@ -497,7 +788,7 @@ function lsim( for idx in 1:nu y[idx, :] = c * x[idx, :] + d * u[idx, :] end - + return y, x end @@ -519,8 +810,8 @@ Parameters: * `sys`: typle `sys = (a, b, c, d)` or `sys = (a, b, c)` where `a` is a square matrix of size (n,n), `b` is a matrix of size (n, m), `c` is a matrix of size (p, n) and (optionally) `d` is a matrix of size (p, m). * `u`: an array of input vectors such that `u[i,:]` is the input vector of size m at sample index `i`. The array has size (N, m) * `z`: vector with the samples of the frequency function argument\\ -* `xt`: transient vector of size n, Default value is the zero vector. - +* `xt`: transient vector of size n, Default value is the zero vector. + Returns: * `y`: the resulting output sequence of size (N,p) * `x`: the resulting state sequence of size (N,p) @@ -538,17 +829,17 @@ function fdsim( p, nc = size(c) nr, m = size(b) d = zeros(p, m) - elseif nn == 4 + elseif nn == 4 a, b, c, d = sys p, nc = size(c) nr, m = size(b) else - print("fdsim: Incorrect number of matrices in sys.") + println("fdsim: Incorrect number of matrices in sys.") return false end y = Matrix{ComplexF64}(undef, nwu, p) - + if length(xt) > 0 ue = hcat(u, reshape(z, nwu, 1)) be = hcat(b, reshape(xt, nr, 1)) @@ -586,7 +877,7 @@ function bilinear_d2c(sys::Tuple, T::Real = 1) ac = (ainv * (a - I)) * 2/T bc = (ainv * b) * 2/sqrt(T) cc = (c * ainv) * 2/sqrt(T) - dc = d .- multidot(c, ainv, b) + dc = d .- multidot(c, ainv, b) return ac, bc, cc, dc end @@ -610,7 +901,7 @@ function bilinear_c2d(sys::Tuple, T::Real = 1) a, b, c, d = sys n = size(a, 1) ainv = inv(I * 2/T - a) - ad = (I * 2/T + a) * ainv + ad = (I * 2/T + a) * ainv bd = (ainv * b) * 2/sqrt(T) cd = (c * ainv) * 2/sqrt(T) dd = d .+ multidot(c, ainv, b) @@ -628,7 +919,7 @@ Parameters: * `T`: frequency scaling factor Returns: -* `wd`: vector of transformed frequencies +* `wd`: vector of transformed frequencies """ cf2df(wc::AbstractVector{<:Number}, T::Real) = 2*atan.(wc*T/2) @@ -642,7 +933,7 @@ Parameters: * `T`: frequency scaling factor Returns: -* `wc`: vector of transformed frequencies +* `wc`: vector of transformed frequencies """ df2cf(wd::AbstractVector{<:Number}, T::Real) = 2*tan.(wd/2)/T diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 4400edf..47b449d 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -18,7 +18,16 @@ using fsid: fdsim, cf2df, df2cf, - ltifr + ltifr, + lrm, + kung_realization, + markov, + make_sys_real, + make_obs_real, + moebius, + moebius_inv, + moebius_arg_inv, + moebius_arg # Helper functions @@ -28,6 +37,9 @@ ztrans(fset) = exp.(2pi * im * fset) computerr(fd, fde) = norm(fd-fde)/norm(fd) +const TOL = 1e-8 + + # Unit test functions function unit_test_estimate_cd_1(fset, m, n, p) z = ztrans(fset) @@ -36,7 +48,7 @@ function unit_test_estimate_cd_1(fset, m, n, p) Ce, De = estimate_cd(fd, z, A, B) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -48,7 +60,7 @@ function unit_test_estimate_cd_2(fset, m, n, p) Ce, De = estimate_cd(fd, z, A, B; estimd = false) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -59,7 +71,7 @@ function unit_test_estimate_cd_3(fset, m, n, p) Ce, De = estimate_cd(fd, z, A, B; type = Complex) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -70,7 +82,7 @@ function unit_test_estimate_bd(fset, m, n, p) Be, De = estimate_bd(fd, z, A, C) fde = fresp(z, A, Be, C, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -82,7 +94,7 @@ function unit_test_estimate_bd_2(fset, m, n, p) Be, De = estimate_bd(fd, z, A, C, estimd=false) fde = fresp(z, A, Be, C, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -93,7 +105,7 @@ function unit_test_estimate_bd_3(fset, m, n, p) Be, De = estimate_bd(fd, z, A, C, type=Complex) fde = fresp(z, A, Be, C, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -103,8 +115,7 @@ function unit_test_transpose_ffdata(fset, m, n, p) fd = fresp(z, A, B, C, D) fde = transpose_ffdata(transpose_ffdata(fd)) - - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -124,7 +135,7 @@ function unit_test_fresp(fset, m, n, p) for fidx in 1:nw emad = z[fidx] * ones(n) for didx in 1:n - dd[didx] = 1 / (emad[didx] - adiag[didx]) + dd[didx] = 1 / (emad[didx] - adiag[didx]) end fd[fidx, :, :] = ((C * diagm(0 => dd)) * B) + D end @@ -135,19 +146,8 @@ function unit_test_fresp(fset, m, n, p) De = D fde = fresp(z, Ae, Be, Ce, De) #SLOW? fdef = fresp(z, Ae, Be, Ce, De) - err = computerr(fd, fde) - errf = computerr(fd, fdef) - - return max(err, errf) # compare largest error with TOL - - #if err > 1e-8: - # print('Unit test "fresp_slow" failed') - # return false - #if errf > 1e-8: - # print('Unit test "fresp" failed', errf) - # return false - # print('Unit test "fresp" "fresp_slow" passed') - # return true + @test computerr(fd, fde) < TOL + @test computerr(fd, fdef) < TOL end function unit_test_fresp_def(fset, m, n, p) @@ -157,7 +157,7 @@ function unit_test_fresp_def(fset, m, n, p) frsp = fresp(z, A, B, C, D) #SLOW? frspf = fresp(z, A, B, C, D) - return computerr(frspf, frsp) + @test computerr(frspf, frsp) < TOL end @@ -171,7 +171,7 @@ function unit_test_fdestim_bd(fset, m, n, p) Be, De, resid = fdestim_bd(zn, Y, U, A, C) fde = fresp(z, A, Be, C, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -186,8 +186,7 @@ function unit_test_fdestim_bd_no_d(fset, m, n, p) Be, De, resid = fdestim_bd(zn, Y, U, A, C, estimd=false) fde = fresp(z, A, Be, C, De) - - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -203,7 +202,7 @@ function unit_test_fdestim_bd_cmplx(fset, m, n, p) Be, De = fdestim_bd(zn, Y, U, A, C; type = Complex) fde = fresp(z, A, Be, C, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -217,7 +216,7 @@ function unit_test_fdestim_cd(fset, m, n, p) Ce, De, resid = fdestim_cd(zn, Y, U, A, B) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -232,12 +231,12 @@ function unit_test_fdestim_cd_no_d(fset, m, n, p) Ce, De, resid = fdestim_cd(zn, Y, U, A, B, estimd=false) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end function unit_test_fdestim_cd_cmplx(fset, m, n, p) - A, B, C, D = buildss(n, m, p) .+ im .* buildss(n, m, p) + A, B, C, D = buildss(n, m, p) .+ im .* buildss(n, m, p) w = 2pi * fset z = ztrans(fset) fd = fresp(z, A, B, C, D) @@ -246,9 +245,9 @@ function unit_test_fdestim_cd_cmplx(fset, m, n, p) Ce, De, resid = fdestim_cd(zn, Y, U, A, B, type=Complex) fde = fresp(z, A, B, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end - + function unit_test_ltifr_slow(fset, m, n, p) z = ztrans(fset) @@ -276,21 +275,19 @@ function unit_test_ltifr_slow(fset, m, n, p) Be = T * B fkern = ltifr(Ae, Be, z) #SLOW? fkernf = ltifr(Ae, Be, z) - err = computerr(fkern1, fkern) - errf = computerr(fkern1, fkernf) - - return max(err, errf) # compare largest error with TOL + @test computerr(fkern1, fkern) < TOL + @test computerr(fkern1, fkernf) < TOL end - + function unit_test_ltifr_def(fset, m, n, p) A = [0 0; 1 0] B = randn(n, m) z = ztrans(fset .+ 1.0e-5) - fkern = ltifr(A, B, z) #SLOW? + fkern = ltifr(A, B, z) #SLOW? fkernf = ltifr(A, B, z) - return computerr(fkernf, fkern) + @test computerr(fkernf, fkern) < TOL end @@ -302,7 +299,7 @@ function unit_test_fconv(fset, m, n, p) wd = cf2df(wc, T) w = df2cf(wd, T) - return computerr(wc, w) + @test computerr(wc, w) < TOL end @@ -310,9 +307,9 @@ function unit_test_fdsid(fset, m, n ,p) N = length(fset) A, B, C, D = buildss(n, m, p) lam, = eigen(A) - rho = maximum(abs.(lam)) - - ## Here we create a random stable DT system + rho = maximum(abs.(lam)) + + # create a random stable DT system A = A/rho/1.01 w = 2pi * fset z = ztrans(fset) @@ -327,7 +324,7 @@ function unit_test_fdsid(fset, m, n ,p) Ae, Be, Ce, De, xt, s = fdsid(fddata, n, 2n; estTrans = true, type = Real) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -336,9 +333,9 @@ function unit_test_fdsid_no_d(fset, m, n ,p) A, B, C, D = buildss(n, m, p) D = zeros(p,m) lam, = eigen(A) - rho = maximum(abs.(lam)) - - ## Here we create a random stable DT system + rho = maximum(abs.(lam)) + + # create a random stable DT system A = A/rho/1.01 w = 2pi * fset z = ztrans(fset) @@ -353,7 +350,7 @@ function unit_test_fdsid_no_d(fset, m, n ,p) Ae, Be, Ce, De, xt, s = fdsid(fddata, n, 2n; estTrans = true, type = Real, estimd = false) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -361,9 +358,9 @@ function unit_test_fdsid_cmplx(fset, m, n ,p) N = length(fset) A, B, C, D = buildss(n, m, p) lam, = eigen(A) - rho = maximum(abs.(lam)) + rho = maximum(abs.(lam)) - # Here we create a random stable DT system + # create a random stable DT system A = A/rho/1.01 B = B + im*randn(n, m) D = D + im*randn(p, m) @@ -378,7 +375,7 @@ function unit_test_fdsid_cmplx(fset, m, n ,p) Ae, Be, Ce, De, xt, s = fdsid(fddata, n, 2*n; estTrans = true, type = Complex) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -390,7 +387,7 @@ function unit_test_ffsid(fset, m, n, p) Ae, Be, Ce, De, s = ffsid(w, fd, n, n+1; type = Real, estimd = true) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -403,7 +400,7 @@ function unit_test_ffsid_no_d(fset, m, n, p) Ae, Be, Ce, De, s = ffsid(w, fd, n, n+1; type = Real, estimd = false) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -419,7 +416,7 @@ function unit_test_ffsid_complex(fset, m, n, p) Ae, Be, Ce, De, s = ffsid(w, fd, n, n+1; type = Complex, estimd = true) fde = fresp(z, Ae, Be, Ce, De) - return computerr(fd, fde) + @test computerr(fd, fde) < TOL end @@ -428,37 +425,117 @@ function unit_test_bilinear(fset, m, n, p) A, B, C, D = buildss(n, m, p) a, b, c, d = bilinear_c2d((A, B, C, D), 2) Ae, Be, Ce, De = bilinear_d2c((a, b, c, d), 2) - - return computerr(A, Ae) + computerr(B, Be) + computerr(C, Ce) + computerr(D, De) -end + + @test (computerr(A, Ae) + computerr(B, Be) + computerr(C, Ce) + computerr(D, De)) < TOL +end + + +function unit_test_lrm(fset, m, n, p) + N = length(fset) + + A, B, C, D = buildss(n, m, p) + + # create a random stable DT system + lam = eigen(A).values + rho = maximum(abs.(lam)) + A = A/rho/1.01 + + z = ztrans(fset) + + fd = fresp(z, A, B, C, D) + u = randn(N, m) + y, x = lsim((A, B, C, D), u; type = Real) + + ff = lrm(u, y) + + @test computerr(fd, ff) < TOL +end + + +function unit_test_moebius(fset, m, n, p) + N = length(fset) + + sset = randn(N) + im*randn(N) + alpha, beta, gamma, delta = randn(4) + im*randn(4) + zset = moebius_arg_inv(sset; alpha, beta, gamma, delta) + sset1 = moebius_arg(zset; alpha, beta, gamma, delta) + + @test computerr(sset, sset1) < TOL + + A, B, C, D = buildss(n, m, p) + zset = randn(N) + im*randn(N) + alpha, beta, gamma, delta = randn(4) + im*randn(4) + sset = moebius_arg(zset; alpha, beta, gamma, delta) + fd = fresp(sset, A, B, C, D) + Ae, Be, Ce, De = moebius((A, B, C, D); alpha, beta, gamma, delta) + fde = fresp(zset, Ae, Be, Ce, De) + + @test computerr(fd, fde) < TOL + + sset = randn(N) + im*randn(N) + zset = moebius_arg_inv(sset; alpha, beta, gamma, delta) + fd = fresp(zset, A, B, C, D) + Ae, Be, Ce, De = moebius_inv((A, B, C, D); alpha, beta, gamma, delta) + fde = fresp(sset, Ae, Be, Ce, De) + + @test computerr(fd, fde) < TOL + + Ae, Be, Ce, De = moebius((A, B, C, D); alpha, beta, gamma, delta) + a, b, c, d = moebius_inv((Ae, Be, Ce, De); alpha, beta, gamma, delta) + sys = (a, b, c, d) + sys0 = (A, B, C, D) + + err = 0 + for idx = 1:4 + err += computerr(sys[idx], sys0[idx]) + end + @test err < TOL + +end + + +function unit_test_markov_kung(fset, m, n, p) + a, b, c, d = buildss(n, m, p) + mp = markov((a, b, c), 2n) + ae, be, ce = kung_realization(mp, n) + me = markov((ae, be, ce), 2n) + @test computerr(me, mp) < TOL +end -const TOL = 1e-8 # Main test function function unit_test(testfunc) N = 100 - if string(testfunc) == "unit_test_fresp_def" || - string(testfunc) == "unit_test_ltifr_def" + if string(testfunc) == "unit_test_fresp_def" || string(testfunc) == "unit_test_ltifr_def" nmpset = [(2, 4, 12), (2, 3, 6)] + elseif string(testfunc) == "unit_test_ltifr_slow" nmpset = [(4, 1, 1), (1,1,1), (2, 1, 1), (2, 4, 12), (2, 3, 12)] # change from (1,1,1) to (2,1,1) + elseif string(testfunc) == "unit_test_fconv" nmpset = [(1, 1, 1)] + + elseif string(testfunc) == "unit_test_lrm" + N = 3000 + nmpset = [(4, 1, 1), (1, 1, 1), (2, 4, 12)] + + elseif string(testfunc) == "unit_test_markov_kung" + nmpset = [(4, 1, 1), (1, 1, 1), (4, 2, 1), (2, 4, 12)] + else nmpset = [(4, 1, 1), (1, 1, 1), (2, 4, 12)] + end fset = buildfset(N) for (n, m, p) in nmpset - error = testfunc(fset, m, n, p) - - @test error .< TOL + testfunc(fset, m, n, p) end return nothing end - + # Run the unit tests @testset "fsid.jl" begin unit_test(unit_test_estimate_cd_1) @@ -486,4 +563,7 @@ end unit_test(unit_test_ffsid_no_d) unit_test(unit_test_ffsid_complex) unit_test(unit_test_bilinear) + unit_test(unit_test_lrm) + unit_test(unit_test_moebius) + unit_test(unit_test_markov_kung) end