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

Открыть пример в Engee

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

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

using Plots
using MATLAB
using CSV
using Statistics
using ControlSystems
plotlyjs();

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

demoroot = @__DIR__

"/user/start/examples/base_simulation/juliacodeintegration"

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

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

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

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

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

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

$$I \cdot \frac{d^2\theta}{dt^2} + mgl \cdot \sin{\theta} = \tau \qquad (I = ml^2)$$

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

$$\frac{d^2\theta}{dt^2} + \frac{g}{l} \cdot \sin{\theta} = \frac{1}{ml^2} \cdot \tau $$

Для малых углов $\sin{\theta} \approx \theta$:

$$\frac{d^2\theta}{dt^2} + \frac{g}{l} \cdot \theta = \frac{1}{ml^2} \cdot \tau$$

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

$$\dot{x} = \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = Ax + Bu$$ $$y = Cx + Du$$

$$x_1 = \theta = y \qquad \dot{x}_1 = \dot{\theta} = x_2$$ $$x_2 = \dot{\theta} \qquad \dot{x}_2 = \ddot{\theta} = -\frac{g}{l} \cdot x_1 + \frac{1}{ml^2}\cdot u$$

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

$$A = \begin{bmatrix} 0 & 1 \\ -g/l & 0 \end{bmatrix}$$ $$B = \begin{bmatrix} 0 \\ 1/(ml^2) \end{bmatrix}$$ $$C = \begin{bmatrix} 1 & 0 \end{bmatrix}$$ $$D = 0$$

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

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

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

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

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

$$W = \frac{4}{s^2 + 19.62}$$

Блоками 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!("Угловое положение маятника")

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

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

Проанализируем модель 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))) * 100

3.9882799356848946e-14

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

plot(skf_time, broadcast(abs, tol_abs), legend = false)

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

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

Выводы

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

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

• Integrator
• Engee Function
• Transfer Fcn
• To CSV
• Step