Документация Engee
Notebook

Интеграция Julia кода в модели Engee с использованием блока Engee Function

В данном примере представлен способ интеграции Julia кода в модели Engee с использованием блока Engee Function.

Перед началом работы подключим необходимые библиотеки:

In [ ]: 
Pkg.add(["Statistics", "ControlSystems", "CSV"])
In [ ]: 
using Plots
using MATLAB
using CSV
using Statistics
using ControlSystems
plotlyjs();

Используем макрос @__DIR__ для того, чтобы узнать папку, в которой лежит интерактивный скрипт:

In [ ]: 
demoroot = @__DIR__
Out[0]:
"/user/start/examples/base_simulation/juliacodeintegration"

Постановка задачи

Предположим, что нас интересует фильтрация зашумлённых измерений углового положения маятника при помощи фильтра Калмана.

1.PNG

Определим коэффициент усиления фильтра Калмана с использованием библиотеки ControlSystems.jl и интегрируем Julia код, реализующий фильтр, в модель Engee с использованием блока Engee Function.

Упрощённая математическая модель маятника в пространстве состояний

Для реализации фильтра Калмана нам потребуется модель маятника в пространстве состояний.

Простая маятниковая система без трения может быть представлена уравнением:

Данное уравнение можно переписать следующим образом:

Для малых углов :

Заменим крутящий момент на входной вектор и представим линеаризованную систему моделью в пространстве состояний:

Тогда матрицы пространства состояний будут определяться следующим образом:

Фильтрация сигнала в Engee с использованием блока Engee Function

Анализ модели

Проанализируем модель juliaCodeIntegration.engee.

Верхний уровень модели:

ef_1.PNG

Маятник представлен передаточной функцией:

Блоками Band-Limited White Noise ProcessNoise и MeasurementNoise представлены внешнее воздействие и шум измерений.

Подсистема Kalman Filter:

ef_2.PNG

Настройки блока KF (Engee Function):

  • Вкладка Main:
ef_3_1.PNG

Количество входных и выходных портов определяется свойствами Number of input ports и Number of output ports.

  • Вкладка Ports:
ef_3_2.PNG

Размерность входных и выходных сигналов определяется свойством Size:

  • на вход блока подаётся два скалярных сигнала и вектор [2x1];
  • на выходе формируется один вектор [2x1].
  • Вкладка Parameters:
parameters.PNG

Параметрам KF_A, KF_B, KF_C, KF_D, KF_K сопоставляются переменные A, B, C, D и K из рабочей области Engee. Их значения будут определены в разделе "Запуск симуляции".

Содержимое вкладки ExeCode блока KF:

include("/user/start/examples/base_simulation/juliacodeintegration/filter.jl")

struct KF <: AbstractCausalComponent; end

function (kalman_filter::KF)(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

Структура KF наследуется от AbstractCausalComponent.

Строка 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:

In [ ]: 
juliaCodeIntegration = engee.load("$demoroot/juliaCodeIntegration.engee", force = true)
Out[0]:
Model(
	name: juliaCodeIntegration
	id: 053374c2-1fab-48ee-8dbe-e698974ce662
)

И инициализируем переменные, необходимые для её симуляции:

In [ ]: 
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:

In [ ]: 
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)
Out[0]:
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:

In [ ]: 
K = kalman(sys, Q, R)
Out[0]:
2×1 Matrix{Float64}:
 3.2413602377158526
 0.2532080953226874

Промоделируем систему:

In [ ]: 
simulationResults = engee.run(juliaCodeIntegration)
Out[0]:
Dict{String, DataFrame} with 2 entries:
  "filtered" => 40001×2 DataFrame…
  "noisy"    => 40001×2 DataFrame…

Закроем модель:

In [ ]: 
engee.close(juliaCodeIntegration,force = true);

Результаты моделирования

Импортируем результаты моделирования:

In [ ]: 
juliaCodeIntegration_noisy_t = simulationResults["noisy"].time;
juliaCodeIntegration_noisy_y = simulationResults["noisy"].value;
In [ ]: 
juliaCodeIntegration_filtered_t = simulationResults["filtered"].time;
juliaCodeIntegration_filtered_y = simulationResults["filtered"].value;

Построим график:

In [ ]: 
plot(juliaCodeIntegration_noisy_t, juliaCodeIntegration_noisy_y, label = "Исходный сигнал")
plot!(juliaCodeIntegration_filtered_t, juliaCodeIntegration_filtered_y, label = "Отфильтрованный сигнал")

title!("Фильтрация сигнала в Engee (Engee Function)")
xlabel!("Время, [с]")
ylabel!("Угловое положение маятника")
Out[0]:

Фильтрация сигнала в Simulink с использованием встроенного блока Kalman Filter

Анализ модели

Проанализируем модель simulinkKalmanFilter.slx:

skf_1.PNG

Настройки блока Kalman Filter:

skf_2.PNG

Во входной порт data модели simulinkKalmanFilter.slx передаются значения зашумлённого сигнала noisy из модели juliaCodeIntegration.engee, благодаря чему отфильтрованные сигналы, полученные в двух средах моделирования, можно сравнивать.

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

Извлечём переменные time и data, необходимые для запуска симуляции в Simulink, из файла juliaCodeIntegration_noisy.csv, полученного после завершения симуляции в Engee:

In [ ]: 
mat"""
skf_T = readtable('/user/CSVs/juliaCodeIntegration_noisy.csv', 'NumHeaderLines', 1);
time = skf_T.Var1;
data = skf_T.Var2;
"""

Промоделируем систему:

In [ ]: 
mdlSim = joinpath(demoroot,"simulinkKalmanFilter")
mat"""myMod = $mdlSim;
    out = sim(myMod);""";

Результаты моделирования

Импортируем результаты моделирования:

In [ ]: 
skf_data = mat"out.yout{1}.Values.Data";
skf_time = mat"out.tout";

Построим график:

In [ ]: 
plot(juliaCodeIntegration_noisy_t, juliaCodeIntegration_noisy_y, label = "Исходный сигнал")
plot!(skf_time, skf_data, label = "Отфильтрованный сигнал")

title!("Фильтрация сигнала в Simulink (Kalman Filter)")
xlabel!("Время, [с]")
ylabel!("Угловое положение маятника")
Out[0]:

Сравнение с результатом, полученным в Engee

Определим среднюю абсолютную погрешность вычислений:

In [ ]: 
tol_abs = skf_data - juliaCodeIntegration_filtered_y;
tol_abs_mean = abs(mean(tol_abs))
Out[0]:
8.708215107605803e-17

Определим среднюю относительную погрешность вычислений:

In [ ]: 
tol_rel = (skf_data - juliaCodeIntegration_filtered_y)./skf_data;
tol_rel_mean = abs(mean(filter(!isnan, tol_rel))) * 100
Out[0]:
2.312573232451167e-12

Построим график абсолютной погрешности вычислений:

In [ ]: 
plot(skf_time, broadcast(abs, tol_abs), legend = false)

title!("Сравнение результатов, полученных в Engee и Simulink")
xlabel!("Время, [с]")
ylabel!("Абсолютная погрешность вычислений")
Out[0]:

Исходя из полученных результатов можно сделать вывод, что реализация фильтра Калмана с использованием
блока Engee Function не уступает по точности встроенному в Simulink блоку Kalman Filter.

Выводы

Данный пример демонстрирует способ интеграции Julia кода в модели Engee и особенности использования блока Engee Function - работу с многомерными сигналами, передачу параметров, подключение файлов внешнего исходного кода.

Блоки, использованные в примере