Skip to content

Commit

Permalink
Some corrections.
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasvarga
andreasvarga committed Mar 11, 2022
1 parent dd18805 commit ac9a8e6
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 43 deletions.
12 changes: 6 additions & 6 deletions src/FDIsynthesis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -386,7 +386,7 @@ function efdsyn(sysf::FDIModel{T}; rdim::Union{Int,Missing} = missing, poles::Un
h = [ hbase zeros(T, rdim, nvec-nout) ]
h = h[:,ip] # permute columns to match unpermuted QR
end
else
else
if rdim == nout
if emptyHD
h = Htemp[ip[1:rdim],:]
Expand All @@ -401,7 +401,7 @@ function efdsyn(sysf::FDIModel{T}; rdim::Union{Int,Missing} = missing, poles::Un
h = [ hbase zeros(T, rdim, nvec-nout) ]
h = h[:,ip] # permute columns to match unpermuted QR
end
end
end
else
if minimal
if rdim == nout
Expand Down Expand Up @@ -618,7 +618,7 @@ function efdbasesel(S::BitArray, degs::Vector{Int}, rdim::Int, nout::Int, simple
(rdim >=1 && rdim <= nvec) || error("ndim must have a positive value not exceeding $nvec")
(nout >=rdim && nout <= nvec) || error("nout must have a value at least $rdim and at most $nvec")

nvec == 1 && (return [1], nodegs ? Int[] : degs )
nvec == 1 && (return [[1]], nodegs ? Int[] : degs )

# find rdim combinations of nout vectors which solve the EFDP
seli = collect(combinations(Vector(1:nvec),nout))
Expand Down Expand Up @@ -1080,7 +1080,7 @@ function afdbasesel(S::BitArray, rwgain::Matrix, degs::Vector{Int}, rdim::Int, n
(rdim >=1 && rdim <= nvec) || error("ndim must have a positive value not exceeding $nvec")
(nout >=rdim && nout <= nvec) || error("nout must have a value at least $rdim and at most $nvec")

nvec == 1 && (return [1], nodegs ? Int[] : degs )
nvec == 1 && (return [[1]], nodegs ? Int[] : degs )

# determine rank of rwgain
if size(rwgain,2) > 0
Expand Down Expand Up @@ -2942,7 +2942,7 @@ function emmbasesel(rgain::Matrix, degs::Vector{Int}, nout::Int, simple::Bool, a
(rdim >=1 && rdim <= nvec) || error("mf must have a positive value not exceeding $nvec")
(nout >=rdim && nout <= nvec) || error("nout must have a value at least $rdim and at most $nvec")

nvec == 1 && (return [1], nodegs ? Int[] : degs )
nvec == 1 && (return [[1]], nodegs ? Int[] : degs )


# find rdim combinations of nout vectors which solve the AFDP
Expand Down Expand Up @@ -3834,7 +3834,7 @@ function ammbasesel(rgain::Matrix, degs::Vector{Int}, nout::Int, simple::Bool, a
nodegs || length(degs) == nvec || error("the dimension of degs must be equal to the number of rows of rgain")
nout <= min(mr,nvec) || error("nout must have a value at most $(min(mr,nvec))")

nvec == 1 && (return [1], nodegs ? Int[] : degs )
nvec == 1 && (return [[1]], nodegs ? Int[] : degs )


# find nout combinations of nvec vectors which solve the AMMP
Expand Down
33 changes: 3 additions & 30 deletions src/MDsynthesis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,20 +391,16 @@ function amdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
# check j-th model detectability
if isel != j
temp = gir(Htemp*qtemp*[sysm.sys[j]; eye(T1,mu,m[j])]; atol)
#i == 1 && (println(size(temp)))
if strongMD
if emdtest
St = fdisspec_(temp[:,1:mu+md[j]], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
else
St = fdisspec_(temp[:,1:mu], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
end
#i == 1 && (println("i= $i j = $j St = "); display(St))
Sj = trues(size(St,1),lfreq)
for k = 1:lfreq
Sj[:,k] = any(view(St,:,:,k),dims=2)
end
#i == 1 && (println("j = $j Sj = "); display(Sj))
#println("i = $i j = $j Sj = "); display(Sj)
if any(maximum(Sj,dims=1))
S[:,:,j] = Sj
else
Expand Down Expand Up @@ -449,7 +445,6 @@ function amdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
else
rdim = min(rdimtarget[i],nvec)
end
#i == 1 && println("rdim = $rdim")

# adjust rmin to be at most the rank of Rwii
rwi > 0 && (rdim = min(rdim,rwi))
Expand All @@ -472,8 +467,6 @@ function amdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
# filter with rdim outputs:
# basesel(i,:) contains the indices of candidate basis vectors;
# ordsel(i) contains the presumably achievable least orders
#println("i = $i degs = $degs rdim = $rdim nout = $nout simple = $simple S = ")
#i == 1 && display(S)
basesel, ordsel = emdbasesel(S, degs, rdim, nout, simple)
#
# update the synthesis using the selections of candidate vector(s),
Expand Down Expand Up @@ -547,7 +540,6 @@ function amdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
@warn "amdsyn: expected reduced order not achieved"
end
h = Htemp[ip[1:rdim],:]
#i == 1 && println("h = $h")
else
qtest, _, info2 = glmcover1([Htemp; eye(nvec)]*qtemp[ip,:], rdim; atol1, atol2, rtol)
end
Expand Down Expand Up @@ -593,15 +585,12 @@ function amdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
Sj[:,k] = any(view(St,:,:,k),dims=2)
end
notOK = !any(maximum(Sj,dims=1))
#i == 1 && println("Sj = $Sj")
#notOK = lfreq == 1 ? !any(St) : !all([any(view(St,:,:,k)) for k = 1:lfreq])
else
if emdtest
Sj = fditspec_(temp[:,1:mu+md[j]]; atol1, atol2, rtol, FDtol = MDtol)
else
Sj = fditspec_(temp[:,1:mu]; atol1, atol2, rtol, FDtol = MDtol)
end
#println("j = $j Sj = $Sj")
notOK = !any(Sj)
end
end
Expand Down Expand Up @@ -1145,7 +1134,6 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
end



# decode input information
mu = sysm.mu
inpu = 1:mu
Expand Down Expand Up @@ -1223,24 +1211,20 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
end
#strongMD && (qtemp = glcf(qtemp; atol1, atol2, rtol, fast, sdeg = sdegNS)[1])
defer = !emptyHDi && minimal
for j = 1:N
for j = 1:N
# check j-th model detectability
if isel != j
temp = gir(Htemp*qtemp*[sysm.sys[j]; eye(T1,mu,m[j])]; atol)
#i == 1 && (println(size(temp)))
if strongMD
if emdtest
St = fdisspec_(temp[:,1:mu+md[j]], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
else
St = fdisspec_(temp[:,1:mu], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
end
#i == 1 && (println("i= $i j = $j St = "); display(St))
Sj = trues(size(St,1),lfreq)
for k = 1:lfreq
Sj[:,k] = any(view(St,:,:,k),dims=2)
end
#i == 1 && (println("j = $j Sj = "); display(Sj))
#println("i = $i j = $j Sj = "); display(Sj)
if any(maximum(Sj,dims=1))
S[:,:,j] = Sj
else
Expand Down Expand Up @@ -1270,10 +1254,6 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
end
end
end
#i == 1 && println("i = $i S = ")
#display(S)
#i == 1 && println("rdimtarget = $rdimtarget")
# setup the number of filter outputs
if ismissing(rdimtarget)
if emptyHDi
rdim = minimal ? 1 : nvec
Expand All @@ -1283,7 +1263,6 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
else
rdim = min(rdimtarget[i],nvec)
end
#i == 1 && println("rdim = $rdim")


# Step 2): compute admissible Q2 to reduce the order of Q2*Q;
Expand All @@ -1307,8 +1286,6 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
# filter with rdim outputs:
# basesel(i,:) contains the indices of candidate basis vectors;
# ordsel(i) contains the presumably achievable least orders
#println("i = $i degs = $degs rdim = $rdim nout = $nout simple = $simple S = ")
#i == 1 && display(S)
basesel, ordsel = emdbasesel(S, degs, rdim, nout, simple)
#
# update the synthesis using the selections of candidate vector(s),
Expand Down Expand Up @@ -1382,7 +1359,6 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
@warn "emdsyn: expected reduced order not achieved"
end
h = Htemp[ip[1:rdim],:]
#i == 1 && println("h = $h")
else
qtest, _, info2 = glmcover1([Htemp; eye(nvec)]*qtemp[ip,:], rdim; atol1, atol2, rtol)
end
Expand Down Expand Up @@ -1428,15 +1404,12 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
Sj[:,k] = any(view(St,:,:,k),dims=2)
end
notOK = !any(maximum(Sj,dims=1))
#i == 1 && println("Sj = $Sj")
#notOK = lfreq == 1 ? !any(St) : !all([any(view(St,:,:,k)) for k = 1:lfreq])
else
if emdtest
Sj = fditspec_(temp[:,1:mu+md[j]]; atol1, atol2, rtol, FDtol = MDtol)
else
Sj = fditspec_(temp[:,1:mu]; atol1, atol2, rtol, FDtol = MDtol)
end
#println("j = $j Sj = $Sj")
notOK = !any(Sj)
end
end
Expand Down Expand Up @@ -1494,7 +1467,7 @@ function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing,
Qt[i] = Htemp*qtemp
end
end

# compute Q3i such that Q3i*Q[i] has a desired stability degree;
# update Q[i] <- Q3i*Qi[i]
k = 1;
Expand Down Expand Up @@ -1619,7 +1592,7 @@ function emdbasesel(S::BitArray, degs::Vector{Int}, rdim::Int, nout::Int, simple
(rdim >=1 && rdim <= nvec) || error("ndim must have a positive value not exceeding $nvec")
(nout >=rdim && nout <= nvec) || error("nout must have a value at least $rdim and at most $nvec")

nvec == 1 && (return [1], nodegs ? Int[] : degs )
nvec == 1 && (return [[1]], nodegs ? Int[] : degs )

# find rdim combinations of nout vectors which solve the EFDP
seli = collect(combinations(Vector(1:nvec),nout))
Expand Down
26 changes: 19 additions & 7 deletions test/test_emdsyn.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,9 @@ using Random
using Test

# Testing examples for EMDSYN
println("Test_emdsyn")
@testset "emdsyn" begin
#rand()
rand()

## Actuator faults
s = rtf('s');
Expand Down Expand Up @@ -60,15 +61,12 @@ gpole(Q)
@test mind == 2 && argmin(mdgain) == mind

@time Q, R, info = emdsyn(sysm; sdeg = -15, poles = [-20], MDSelect = Int[]); info.MDperf
display(Q)
display(R)
@time Q, R, info = emdsyn(sysm; sdeg = -15, poles = [-20], MDSelect = [1]); info.MDperf
@time Q, R, info = emdsyn(sysm; sdeg = -15, poles = [-20], MDSelect = [2]); info.MDperf
@time Q, R, info = emdsyn(sysm; sdeg = -15, poles = [-20], normalize = true); info.MDperf
@test info.MDperf[1,:] == info.MDperf[:,1]

sysc1=MDModel(rss(3,2,6,stable = true); mu = 2, md = 1,mw = 2, ma = 1)
display(sysc1)
sysc2=MDModel(rss(3,2,6,stable = true); mu = 2, md = 2, mw = 1)
sysm = mdmodset([sysc1,sysc2])
@test mddist2c([sysc1,sysc2],[sysc1,sysc2])[1] == mddist2c([sysc1,sysc2],sysm)[1] == mddist2c(sysm,[sysc1,sysc2])[1]
Expand All @@ -88,6 +86,23 @@ sysm1 = mdmodset([sys1, sys2, sys3, sys4], controls = 1:mu,
disturbances = [mu .+ (1:md[1]), mu .+ (1:md[2]),mu .+ (1:md[3]),mu .+ (1:md[4])]);
sysm2 = MDMModel([sys1, sys2, sys3, sys4]; mu, md);

@time distgap, fpeak = mddist(sysm, distance = "Inf")
@time distgap1, fpeak1 = mddist(sysm, distance = "Inf", MDfreq = fpeak[:])
@test distgap ≈ distgap1 && fpeak ≈ fpeak1

@time distgapd, fpeakd = mddist(sysm, distance = "Inf", cdinp = true)
@time distgapd1, fpeakd1 = mddist(sysm, distance = "Inf", cdinp = true, MDfreq = fpeakd[:])
@test distgapd ≈ distgapd1 && fpeakd ≈ fpeakd1 && all(abs.(distgapd-distgap+1.e-7*I) .>= 0)

@time distgapc, fpeakc = mddist2c(sysm,sysm, distance = "Inf")
@time distgapc1, fpeakc1 = mddist2c(sysm, sysm, distance = "Inf", MDfreq = fpeakc[:])
@test norm(distgapc-distgapc1,Inf) < 1.e-7

@time distgapcd, fpeakcd = mddist2c(sysm,sysm, cdinp = true, distance = "Inf")
@time distgapcd1, fpeakcd1 = mddist2c(sysm, sysm, cdinp = true, distance = "Inf", MDfreq = fpeakcd[:])
@test norm(distgapcd-distgapcd1,Inf) < 1.e-7 && all(abs.(distgapcd-distgapc+1.e-7*I) .>= 0)


@test all(iszero.(sysm.sys .- sysm1.sys,atol=1.e-7)) && sysm.mu == sysm1.mu &&
sysm.md == sysm1.md && sysm.mw == sysm1.mw && sysm.ma == sysm1.ma

Expand Down Expand Up @@ -161,13 +176,11 @@ end
# setup synthesis model
sysm = mdmodset(sysu, controls = 1:mu);


# call of EMDSYN with the options for stability degree -1 and pole -1 for
# the filters, tolerance and a design matrix H to form a linear combination
# of the left nullspace basis vectorsH = [ 0.7645 0.8848 0.5778 0.9026 ];
H = [ 0.7645 0.8848 0.5778 0.9026 ];
@time Q, R, info = emdsyn(sysm, sdeg = -1, poles = [-1], HDesign = H);
display(Q)
R1 = mdIFeval(Q, sysm, atol=1.e-7, minimal = true)
@test norm(diag(info.MDperf)) < 1.e-7 && all(iszero.(R.sys .- R1.sys,atol=1.e-7))

Expand All @@ -176,7 +189,6 @@ R11 = mdIFeval(Q1, sysm, atol=1.e-7, minimal = true)
@test norm(info.MDperf-info1.MDperf) < 1.e-7 && all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) &&
all(iszero.(R.sys .- R11.sys,atol=1.e-7))


@time Q, R, info = emdsyn(sysm, nullspace = true, sdeg = -1, poles = [-1], HDesign = H);
R1 = mdIFeval(Q, sysm, atol=1.e-7, minimal = true)
@test norm(diag(info.MDperf)) < 1.e-7 && all(iszero.(R.sys .- R1.sys,atol=1.e-7))
Expand Down

2 comments on commit ac9a8e6

@andreasvarga
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/56451

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.0.1 -m "<description of version>" ac9a8e642c75697eaf473017ecf2544b2e23cce4
git push origin v1.0.1

Please sign in to comment.