Школа системного моделирования

Определение частотной характеристики с помощью блока генератора сигнала псевдослучайной двоичной последовательности

Pkg.add("ControlSystems")
Pkg.add("DSP")

using ControlSystems
using DSP
using EngeeDSP

Период дискретизации алгоритма

Ts = 0.001 # @param {type:"number"}
;

Массив из частот, на которых будет собираться частотный отклик в логарифмическом масштабе.

wmin - минимальная частота, на которой строить ЛАФЧХ
wmax - максимальная частота, на которой строить ЛАФЧХ
n_ponts - количество точек на графике

wmin = 0.1 # @param {type:"number"}
wmax = 100 # @param {type:"number"}
n_points = 20 # @param {type:"number"}
w = collect(logrange(wmin, wmax, length=n_points));

Запуск симуляции

modelName = "prbs" # @param {type:"string"}
;

if modelName ∉ getfield.(engee.get_all_models(), :name)
    engee.load("$(@__DIR__)/$(modelName).engee");
end
engee.run(modelName);

Запись данных симуляции в переменные

model_data = collect(PRBSSigGenData);
t = model_data[:,1];
data = model_data[:,2]

nd = size(data, 1)
nw = length(w)
    
ready = zeros(nd)
perturbation = zeros(nd)
input = zeros(nd)
output = zeros(nd)

for i = 1:nd
    ready[i] = data[i][1]
    perturbation[i] = data[i][2]
    input[i] = data[i][3]     
    output[i] = data[i][4]
end

Обрабатываем сигналы перед вычислением частотного отклика системы

ts = t[2] - t[1]
dupidx = findall(diff(t) .< ts/2)
deleteat!(ready, dupidx)
deleteat!(perturbation, dupidx)
deleteat!(input, dupidx)

idx = findall(ready .== 1.)
d = perturbation[idx]
u = input[idx]
y = output[idx]

nfft = sum(idx .> 0.)

if !iszero(rem(nfft, 2))
    freq = (2*pi/ts) .* collect(range(0., 1., nfft + 1))
    lastel = Int64((nfft + 1) / 2)
    freq = freq[1:lastel]
else
    freq = (pi/ts) .* collect(range(0., 1., Int64(nfft/2) + 1))
end

pd = 2*pi / (nfft - 1)

filterwindow = [0.5 .* ones(nfft-1) .- 0.5 .* cos.(pd .* collect(0:nfft-2)); 0.]

fd = EngeeDSP.Functions.fft(d .* filterwindow);
fu = EngeeDSP.Functions.fft(u .* filterwindow);
fy = EngeeDSP.Functions.fft(y .* filterwindow);

d2u = fu ./ fd;
d2y = fy ./ fd;
u2y = d2y ./ d2u;

Построение частотной характеристики

Получаем частотный отклик системы

Обратите внимание, что в папке должен находиться файл с функцией interpfreqresp.jl

# Подключаем файл со служебной функцией
include("$(@__DIR__)/interpfreqresp.jl")
plant_response = interpfreqresp(freq, u2y[1:length(freq)], w);

Вычисляем ЛАФЧХ системы по частотному отклику

magnitude_estimated = 20 * log10.(abs.(plant_response));
phase_estimated = rad2deg.(DSP.unwrap((atan.(imag.(plant_response), real.(plant_response)))));

Строим ЛАФЧХ системы по частотному отклику

gr()
p1 = plot(w, magnitude_estimated, linealpha = 0.0, markershape = :circle, markersize = 4, ylabel = "L, дБ", legend = :none)
p2 = plot(w, phase_estimated, linealpha = 0.0, markershape = :circle, markersize = 4, label = "Модель", xlabel = "ω, рад/с", ylabel = "Ψ, град", legend = :bottomleft)
plot(p1, p2, layout = (2,1), grid = false, framestyle = :box, xscale=:log10)

Вычисляем ЛАФЧХ передаточной функции системы.

Сравним экспериментальную ЛАФЧХ с характеристикой, построенной по передаточной функции.

# Объект управления
L = 0.006   # Индуктивность Гн
k = 0.10    # Коэффициент противо-ЭДС (В⋅с)/рад
R = 0.25    # Cопротивление Ом
J = 0.015   # Момент инерции ротора кг⋅м^2
K = 0.10    # Постоянная крутящего момента Н⋅м/А
tfnum = [K];
tfden = [L*J, R*J, K*k];

# ПИД регулятор
kp = 0.1322157366217838
ki = 0.8546367917451264
C = pid(kp, ki, 0; form=:parallel)

# Разомкнутая система
sys = C*tf(tfnum, tfden)

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
         0.01322157366217838s + 0.08546367917451264
------------------------------------------------------------
8.999999999999999e-5s^3 + 0.00375s^2 + 0.010000000000000002s

Continuous-time transfer function model

W = c2d(sys, Ts) # Дискретизируем передаточную функцию с периодом квантования, равным периоду дискретизации алгоритма.
magnitude_ideal, phase_ideal, w_ideal = bode(W);
magnitude_ideal = 20 * log10.(reshape(magnitude_ideal, size(magnitude_ideal, 3),))
phase_ideal = reshape(phase_ideal, size(phase_ideal, 3),);

Сравнение ЛАФЧХ

gr()
p1 = plot(w, magnitude_estimated, linealpha = 0.0, markershape = :circle, markersize = 4, ylabel = "L, дБ", legend = :none)
plot!(w_ideal, magnitude_ideal, xlims = (wmin, wmax))#, ylabel = "L, дБ", legend = :none)
p2 = plot(w, phase_estimated, linealpha = 0.0, markershape = :circle, markersize = 4, label = "Модель", xlabel = "ω, рад/с", ylabel = "Ψ, град", legend = :bottomleft)
plot!(w_ideal, phase_ideal, xlims = (wmin, wmax), label = "Передаточная функция")
plot(p1, p2, layout = (2,1), grid = false, framestyle = :box, xscale=:log10)

Модуль 4. Лекция 12 - Школа системного моделирования