Engee 文档
Notebook

使用Engee功能块将Julia代码集成到Engee模型中

此示例显示了使用Engee函数块将Julia代码集成到Engee模型中的方法。

在开始工作之前，我们将连接必要的库。:

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]（https://juliacontrol.github.io/ControlSystems.jl/stable /）并使用Engee函数块将实现过滤器的Julia代码集成到Engee模型中。

状态空间中钟摆的简化数学模型

为了实现卡尔曼滤波器，我们需要一个状态空间中的钟摆模型。

一个简单的无摩擦摆系统可以用方程表示:

这个方程可以改写如下:

对于小角度 :

$ $

更换扭矩 到输入向量 让我们将线性化系统表示为状态空间中的模型。:

然后状态空间矩阵将定义如下:

使用Engee功能块在Engee中进行信号滤波

模型分析

我们来分析一下模型 juliaCodeIntegration.engee.

模型的上层:

ef_1.PNG

钟摆由传递函数表示:

块带限白噪声 ProcessNoise 和 MeasurementNoise 给出了外部影响和测量噪声.

子系统 Kalman Filter:

ef_2.PNG

块设置 KF (工程功能):

*标签 Main:

ef_3_1.PNG

输入和输出端口数由属性输入端口数和输出端口数确定。

*标签 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;结束

函数(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）

返回xhat

结束

结构 KF 继承自 AbstractCausalComponent.

行 include("/user/start/examples/base_simulation/juliacodeintegration/filter.jl") 需要连接外部源代码文件 filter.jl.

在文件中 filter.jl 函数被定义 xhat_predict 和 xhat_update 设计用于产生输出信号 xhat 街区 KF:

``'茱莉亚
函数xhat_predict(u,x,A,B)
返回(Bu)+(Ax)
结束

函数xhat_update(u,y,x,C,D,K)
返回(K*(y。-((Cx)。+(Du)))))
结束

在这个例子中，没有内部状态，所以方法 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]（https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.ss ）从图书馆[ControlSystems.jl]（https://juliacontrol.github.io/ControlSystems.jl/stable/):

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
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中使用内置卡尔曼滤波单元进行信号滤波

模型分析

我们来分析一下模型 simulinkKalmanFilter.slx:

skf_1.PNG

块设置卡尔曼滤波器:

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函数块在精度上并不逊色于Simulink内置的卡尔曼滤波器块。

结论

此示例演示如何将Julia代码集成到Engee模型中，以及使用Engee函数块处理多维信号、传输参数和连接外部源代码文件的功能。

示例中使用的块