Интеграция Julia кода в модели Engee с использованием блока Engee Function
В данном примере представлен способ интеграции Julia кода в модели Engee с использованием блока Engee Function.
Перед началом работы подключим необходимые библиотеки:
using Plots
using MATLAB
using CSV
using Statistics
using ControlSystems
plotlyjs();Используем макрос @DIR для того, чтобы узнать папку, в которой лежит интерактивный скрипт:
demoroot = @__DIR__"/user/start/examples/base_simulation/juliacodeintegration"Постановка задачи
Предположим, что нас интересует фильтрация зашумлённых измерений углового положения маятника при помощи фильтра Калмана.
Определим коэффициент усиления фильтра Калмана с использованием библиотеки ControlSystems.jl и интегрируем Julia код, реализующий фильтр, в модель Engee с использованием блока Engee Function.
Упрощённая математическая модель маятника в пространстве состояний
Для реализации фильтра Калмана нам потребуется модель маятника в пространстве состояний.
Простая маятниковая система без трения может быть представлена уравнением:
Данное уравнение можно переписать следующим образом:
Для малых углов :
Заменим крутящий момент на входной вектор и представим линеаризованную систему моделью в пространстве состояний:
Тогда матрицы пространства состояний будут определяться следующим образом:
Анализ модели
Проанализируем модель juliaCodeIntegration.engee.
Верхний уровень модели:
Маятник представлен передаточной функцией:
Блоками Band-Limited White Noise ProcessNoise и MeasurementNoise представлены внешнее воздействие и шум измерений.
Подсистема Kalman Filter:
Настройки блока KF (Engee Function):
-
Вкладка
Main:
Количество входных и выходных портов определяется свойствами Number of input ports и Number of output ports.
-
Вкладка
Ports:
Размерность входных и выходных сигналов определяется свойством Size:
-
на вход блока подаётся два скалярных сигнала и вектор
[2x1]; -
на выходе формируется один вектор
[2x1]. -
Вкладка
Parameters:
Параметрам parameter1, parameter2, parameter3, parameter4, parameter5 сопоставляются переменные A, B, C, D и K из рабочей области Engee. Их значения будут определены в разделе "Запуск симуляции".
Содержимое вкладки ExeCode блока KF:
include("/user/start/examples/base_simulation/juliacodeintegration/filter.jl")
struct KalmanFilter <: AbstractCausalComponent
A::Matrix{Real}
B::Vector{Real}
C::Matrix{Real}
D::Real
K::Matrix{Real}
function KalmanFilter()
new(parameter1, parameter2, parameter3, parameter4, parameter5)
end
end
function (kf::KalmanFilter)(t::Real, u, y, x)
xhat = xhat_predict(u, x, kf.A, kf.B) + xhat_update(u, y, x, kf.C, kf.D, kf.K)
return xhat
end
update!(kf::KalmanFilter, t::Real, u, y, x) = kfСтруктура KalmanFilter наследуется от AbstractCausalComponent, определяет свойства блока KF - A, B, C, D, K и их типы данных.
Функция function KalmanFilter() является внутренним конструктором, инициализирующим свойства структуры A, B, C, D, K параметрами parameter1, parameter2, parameter3, parameter4 и parameter5.
Строка include("/user/start/examples/base_simulation/juliacodeintegration/filter.jl") необходима для подключения внешнего файла исходного кода filter.jl.
В файле filter.jl определяются функции xhat_predict и xhat_update, предназначенные для формирования выходного сигнала xhat блока KF:
function xhat_predict(u, x, A, B)
return (B * u) + (A * x)
end
function xhat_update(u, y, x, C, D, K)
return (K * (y .- ((C * x) .+ (D * u))))
endМетод update! необходим для обновления внутренних состояний блока. В данном примере внутренние состояния отсутствуют, поэтому он не оказывает никакого влияния на результат.
Запуск симуляции
Инициализируем переменные, необходимые для симуляции модели juliaCodeIntegration.engee.
g = 9.81;
m = 1;
L = 0.5;
Ts = 0.001;
Q = 1e-3;
R = 1e-4;
processNoisePower = Q * Ts;
measurementNoisePower = R * Ts;Сформируем модель маятника в пространстве состояний с использованием функции ss из библиотеки ControlSystems.jl:
A = [0 1; -g/L 0];
B = [0; 1/(m*L^2)];
C = [1.0 0.0];
D = 0.0;
sys = ss(A, B, C, D)StateSpace{Continuous, Float64}
A =
0.0 1.0
-19.62 0.0
B =
0.0
4.0
C =
1.0 0.0
D =
0.0
Continuous-time state-space modelОпределим коэффициент усиления фильтра Калмана с использованием функции kalman из библиотеки ControlSystems.jl:
K = kalman(sys, Q, R)2×1 Matrix{Float64}:
3.2413602377158526
0.2532080953226874Промоделируем систему:
juliaCodeIntegration = engee.load("$demoroot/juliaCodeIntegration.engee", force = true)
simulationResults = engee.run(juliaCodeIntegration)Dict{String, DataFrames.DataFrame} with 2 entries:
"filtered" => 40001×2 DataFrame…
"noisy" => 40001×2 DataFrame…Закроем модель:
engee.close(juliaCodeIntegration,force = true);Результаты моделирования
Импортируем результаты моделирования:
juliaCodeIntegration_noisy_t = simulationResults["noisy"].time;
juliaCodeIntegration_noisy_y = simulationResults["noisy"].value;juliaCodeIntegration_filtered_t = simulationResults["filtered"].time;
juliaCodeIntegration_filtered_y = simulationResults["filtered"].value;Построим график:
plot(juliaCodeIntegration_noisy_t, juliaCodeIntegration_noisy_y, label = "Исходный сигнал")
plot!(juliaCodeIntegration_filtered_t, juliaCodeIntegration_filtered_y, label = "Отфильтрованный сигнал")
title!("Фильтрация сигнала в Engee (Engee Function)")
xlabel!("Время, [с]")
ylabel!("Угловое положение маятника")
Анализ модели
Проанализируем модель simulinkKalmanFilter.slx:
Настройки блока Kalman Filter:
Во входной порт data модели simulinkKalmanFilter.slx передаются значения зашумлённого сигнала noisy из модели juliaCodeIntegration.engee, благодаря чему отфильтрованные сигналы, полученные в двух средах моделирования, можно сравнивать.
Запуск симуляции
Извлечём переменные time и data, необходимые для запуска симуляции в Simulink, из файла juliaCodeIntegration_noisy.csv, полученного после завершения симуляции в Engee:
mat"""
skf_T = readtable('/user/CSVs/juliaCodeIntegration_noisy.csv', 'NumHeaderLines', 1);
time = skf_T.Var1;
data = skf_T.Var2;
"""Промоделируем систему:
mdlSim = joinpath(demoroot,"simulinkKalmanFilter")
mat"""myMod = $mdlSim;
out = sim(myMod);""";Результаты моделирования
Импортируем результаты моделирования:
skf_data = mat"out.yout{1}.Values.Data";
skf_time = mat"out.tout";Построим график:
plot(juliaCodeIntegration_noisy_t, juliaCodeIntegration_noisy_y, label = "Исходный сигнал")
plot!(skf_time, skf_data, label = "Отфильтрованный сигнал")
title!("Фильтрация сигнала в Simulink (Kalman Filter)")
xlabel!("Время, [с]")
ylabel!("Угловое положение маятника")
Сравнение с результатом, полученным в Engee
Определим среднюю абсолютную погрешность вычислений:
tol_abs = skf_data - juliaCodeIntegration_filtered_y;
tol_abs_mean = abs(mean(tol_abs))1.0422585954599654e-17Определим среднюю относительную погрешность вычислений:
tol_rel = (skf_data - juliaCodeIntegration_filtered_y)./skf_data;
tol_rel_mean = abs(mean(filter(!isnan, tol_rel))) * 1003.9882799356848946e-14Построим график абсолютной погрешности вычислений:
plot(skf_time, broadcast(abs, tol_abs), legend = false)
title!("Сравнение результатов, полученных в Engee и Simulink")
xlabel!("Время, [с]")
ylabel!("Абсолютная погрешность вычислений")
Исходя из полученных результатов можно сделать вывод, что реализация фильтра Калмана с использованием блока Engee Function не уступает по точности встроенному в Simulink блоку Kalman Filter.