ЛАФЧХ из модели с помощью PRBS
Автор
svetlanaСоавторы
andrewraysОпределение частотной характеристики с помощью блока генератора сигнала псевдослучайной двоичной последовательности
Для построения частотных характеристик системы, необходимо получить частотный отклик модели: подать тестовый сигнал на вход системы и измерить ее выход. Чтобы сократить время эксперимента, без потери качества результатов можно подавать не гармонический сигнал, а сигнал псевдослучайной двоичной последовательности (PRBS). PRBS — это периодический сигнал, который может принимать только два определенных значения.
В этом примере в модели model_analysis_with_prbs.engee с помощью специальных блоков подается PRBS сигнал и записывается отклик системы. В этом скрипте данные обрабатываются и строится ЛАФЧХ.
Pkg.add("ControlSystems")
Pkg.add("DSP")
In [ ]:
using ControlSystems
using DSP
using EngeeDSP
Период дискретизации алгоритма
In [ ]:
Ts = 0.001 # @param {type:"number"}
;
Массив из частот, на которых будет собираться частотный отклик в логарифмическом масштабе.
wmin - минимальная частота, на которой строить ЛАФЧХ
wmax - максимальная частота, на которой строить ЛАФЧХ
n_ponts - количество точек на графике
In [ ]:
wmin = 0.01 # @param {type:"number"}
wmax = 100 # @param {type:"number"}
n_points = 20 # @param {type:"number"}
w = collect(logrange(wmin, wmax, length=n_points));
Запуск симуляции
In [ ]:
modelName = "model_analysis_with_prbs" # @param {type:"string"}
;
In [ ]:
if modelName ∉ getfield.(engee.get_all_models(), :name)
engee.load("$(@__DIR__)/$(modelName).engee");
end
engee.run(modelName);
Запись данных симуляции в переменные
In [ ]:
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
Обрабатываем сигналы перед вычислением частотного отклика системы
In [ ]:
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
In [ ]:
# Подключаем файл со служебной функцией
include("$(@__DIR__)/interpfreqresp.jl")
plant_response = interpfreqresp(freq, u2y[1:length(freq)], w);
Вычисляем ЛАФЧХ системы по частотному отклику
In [ ]:
magnitude_estimated = 20 * log10.(abs.(plant_response));
phase_estimated = rad2deg.(DSP.unwrap((atan.(imag.(plant_response), real.(plant_response)))));
Строим ЛАФЧХ системы по частотному отклику
In [ ]:
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)
Out[0]:
Вычисляем ЛАФЧХ передаточной функции системы.
Сравним экспериментальную ЛАФЧХ с характеристикой, построенной по передаточной функции.
In [ ]:
# Объект управления
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];
sys = tf(tfnum, tfden)
W = sys
Out[0]:
In [ ]:
W = c2d(sys, Ts) # Дискретизируем передаточную функцию с периодом квантования, равным периоду дискретизации алгоритма.
magnitude_ideal, phase_ideal, w_ideal = bode(W);
magnitude_ideal = reshape(magnitude_ideal, size(magnitude_ideal, 3),)
phase_ideal = reshape(phase_ideal, size(phase_ideal, 3),);
Сравнение ЛАФЧХ
In [ ]:
p1 = plot(w, magnitude_estimated, linealpha = 0.0, markershape = :circle, markersize = 4, ylabel = "L, дБ", legend = :none)
plot!(w_ideal, 20 * log10.(magnitude_ideal), xlims = (wmin, wmax))
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)
Out[0]:
