Describing function analysis

describing_function.jl

# ref: https://www.control.lth.se/fileadmin/control/Education/EngineeringProgram/FRTN05/2019/lec06eight.pdf
# ref: https://www.kth.se/social/upload/4fd0a7fdf276547736000003/lec08.pdf
using ControlSystemsBase
using ControlSystems
using ControlSystemIdentification

s = tf("s")
G = 4/(s*(s+1)^2) # A test system

nl = ControlSystemsBase.saturation(1)
# nl = ControlSystemsBase.deadzone(1)
# nl = ControlSystemsBase.nonlinearity(x->abs(x)<1 ? 0 : sign(x))

e(t, A, w) = A*sin(w*t)
integrate(fun,data,ω,h) = h*sum(fun(ω*(i-1)*h)*data[i] for i = eachindex(data))


W = exp10.(LinRange(-1, 2, 30)) # A frequency vector
As = 0.1:0.1:10 # A vector of amplitudes
h = 0.001 # Discretization time
# res = lsim(nl, (x,t)->e(t, 1, 1), 0:h:10)
# plot(res)
Gs = map(eachindex(As)) do ai # Compute descibing function
    A = As[ai]
    Nw = ComplexF64[]
    for i in eachindex(W)
        w = W[i]
        t = 0:h:(2pi/w)
        u = nl.f[1].(e.(t, A, w)) # NOTE: this only works for simple static nonlinearities. Use lsim for more complex stuff like backlash
        a = w/(A*pi) * integrate(cos, u, w, h)
        b = w/(A*pi) * integrate(sin, u, w, h)
        push!(Nw, complex(b, a))
    end
    FRD(W, Nw)
end

NA = map(Gs) do G # For odd nonlinearities, compute N(A,w) = N(A)
    median(real.(G.r)) # NOTE: this only holds for odd nonlinearities
end
plot(As, NA, title="N(A)", xlabel="A")
##

@userplot Describingfunctionplot2
@recipe function nyquist_df(p::Describingfunctionplot2)
    A, NA = p.args[1:2]
    iNA = -1 ./ NA
    anninds = 1:length(A)÷8:length(A)
    anns = [(real(iNA[i]), imag(iNA[i])+0.3, text("A = $(A[i])", 10)) for i in anninds]
    @series begin
        hover := A
        linewidth --> 2
        lab --> "\$N(A)\$"
        real(iNA), imag(iNA)
    end
    @series begin
        seriestype := :scatter
        primary := false
        annotations := anns
        real(iNA[anninds]), imag(iNA[anninds])
    end
end

# Plot describing function on top of Nyquist curve. The point where the Nyquist curve crosses the describing function predicts the amplitude and frequency of the limit cycle (use hover information with the plotly backend to see amplitude and frequency).
# plotly()
nyquistplot(G)
describingfunctionplot2!(As, NA)

##
# plot()
# for G in Gs
#     plot!(G)
# end
# display(current())

## Simulate the nonlinear system to see if the analysis was accurate
plot(step(feedback(G*nl), 40))