Engee documentation
Notebook

Integrating Julia code into Engee models using the Engee Function block

This example shows a way to integrate Julia code into the Engee model using the Engee Function block.

Before starting work, we will connect the necessary libraries.:

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

Using a macro @__DIR__ in order to find out the folder where the interactive script is located:

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

Setting the task

Let's assume that we are interested in filtering noisy angular position measurements. a pendulum using a Kalman filter.

1.PNG

Let's determine the gain of the Kalman filter using the library [ControlSystems.jl](https://juliacontrol.github.io/ControlSystems.jl/stable /) and integrate the Julia code implementing the filter into the Engee model using the Engee Function block.

Simplified mathematical model of a pendulum in the state space

To implement the Kalman filter, we need a model of a pendulum in the state space.

A simple frictionless pendulum system can be represented by the equation:

This equation can be rewritten as follows:

For small angles :

Replace the torque to the input vector and let's represent the linearized system as a model in the state space.:

Then the state space matrices will be defined as follows:

Signal filtering in Engee using the Engee Function block

Model analysis

Let's analyze the model juliaCodeIntegration.engee.

The upper level of the model:

ef_1.PNG

The pendulum is represented by a transfer function:

In Blocks Band-Limited White Noise ProcessNoise and MeasurementNoise The external impact and measurement noise are presented.

Subsystem Kalman Filter:

ef_2.PNG

Block Settings KF (Engee Function):

  • Tab Main:
ef_3_1.PNG

The number of input and output ports is determined by the properties Number of input ports and Number of output ports.

  • Tab Ports:
ef_3_2.PNG

The dimension of the input and output signals is determined by the Size property:

  • Two scalar signals and a vector are applied to the input of the unit. [2x1];
  • one vector is formed at the output [2x1].
  • Tab Parameters:
parameters.PNG

Parameters KF_A, KF_B, KF_C, KF_D, KF_K variables are compared A, B, C, D and K from the Engee workspace. Their values will be defined in the section "Starting the simulation".

Contents of the ExeCode block tab 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

Structure KF inherited from AbstractCausalComponent.

Line include("/user/start/examples/base_simulation/juliacodeintegration/filter.jl") required to connect an external source code file filter.jl.

In the file filter.jl The functions are defined xhat_predict and xhat_update designed to generate an output signal xhat the block 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

In this example, there are no internal states, so the method update!, designed to update them, is missing.

Running the simulation

Loading the model juliaCodeIntegration.engee:

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

And initialize the variables necessary to simulate it.:

In [ ]: 
g = 9.81;
m = 1;
L = 0.5;

Ts = 0.001;
Q = 1e-3;
R = 1e-4;

processNoisePower = Q * Ts;
measurementNoisePower = R * Ts;

Let's form a model of a pendulum in the state space using the function ss from the library 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

Let's determine the gain of the Kalman filter using the function kalman from the library ControlSystems.jl:

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

Let's model the system:

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

Closing the model:

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

Simulation results

Importing the simulation results:

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;

Let's build a graph:

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

Signal filtering in Simulink using the built-in Kalman Filter unit

Model analysis

Let's analyze the model simulinkKalmanFilter.slx:

skf_1.PNG

Block Settings Kalman Filter:

skf_2.PNG

To the input port data models simulinkKalmanFilter.slx the values of the noisy signal are transmitted noisy from the model juliaCodeIntegration.engee due to this, the filtered signals obtained in the two simulation environments can be compared.

Running the simulation

Extracting the variables time and data required to run the simulation in Simulink, from the file juliaCodeIntegration_noisy.csv received after completion of the simulation in Engee:

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

Let's model the system:

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

Simulation results

Importing the simulation results:

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

Let's build a graph:

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

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

Comparison with the result obtained in Engee

Let's determine the average absolute error of calculations:

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

Let's determine the average relative error of calculations:

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

Let's plot the absolute error of calculations.:

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

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

Based on the results obtained, it can be concluded that the implementation of the Kalman filter using
The Engee Function block is not inferior in accuracy to the Kalman Filter block built into Simulink.

Conclusions

This example demonstrates how to integrate Julia code into Engee models and the features of using the Engee Function block - working with multidimensional signals, transferring parameters, and connecting external source code files.

Blocks used in example