Engee 文档
Notebook

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

本例展示了如何使用 ** Engee Function** 块将 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"

任务设置

假设我们希望使用卡尔曼滤波器过滤$\theta$ 摆角位置的噪声测量值。

1.PNG

我们使用 ControlSystems.jl库确定卡尔曼滤波器的增益，并使用 Engee Function 块将实现滤波器的 Julia 代码集成到 Engee 模型中。

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

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

一个简单的无摩擦钟摆系统可以用公式表示：

$$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$$.

我们用输入向量$u$ 代替扭矩$\tau$ ，并用状态空间模型表示线性化系统：

$$\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 中进行信号滤波

模型分析

让我们对模型juliaCodeIntegration.engee 进行分析。

模型的上层：

ef_1.PNG

摆锤由传递函数表示：

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

块 带限白噪声ProcessNoise 和MeasurementNoise 代表外部影响和测量噪声。

子系统Kalman Filter ：

ef_2.PNG

单位设置KF （Engee 功能）：

  • TabMain ：

ef_3_1.PNG

输入和输出端口的数量由 ** 输入端口数** 和 ** 输出端口数** 属性定义。

  • TabPorts ：

ef_3_2.PNG

输入和输出信号的维度由 Size 属性定义：

  • 两个标量信号和一个矢量[2x1] 输入到程序块；
  • 输出端生成一个矢量[2x1] 。
  • TabParameters ：

parameters.PNG

参数KF_A,KF_B,KF_C,KF_D,KF_K 被映射到Engee工作区的变量A,B,C,D 和K 。它们的值将在模拟运行部分定义。

KF 块的ExeCode选项卡内容：

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)

    返回 xhat
结束

结构KF 继承自AbstractCausalComponent 。

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

文件filter.jl 定义了函数xhat_predict 和xhat_update ，用于生成模块KF 的输出信号xhat ：

函数 xhat_predict(u, x, A, B)
    返回 (B * u) + (A * x)
结束

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

在这个例子中，没有内部状态，因此没有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;

使用库 ControlSystems.jl 中的函数 ss 在状态空间中形成摆的模型：

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

使用库 ControlSystems.jl 中的函数 kalman 确定卡尔曼滤波增益：

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

模型分析

让我们对模型simulinkKalmanFilter.slx 进行分析：

skf_1.PNG

卡尔曼滤波单元设置：

skf_2.PNG

模型simulinkKalmanFilter.slx 的输入端口data 输入了来自模型juliaCodeIntegration.engee 的噪声信号noisy 的值，这样就可以比较在两个仿真环境中获得的滤波信号。

运行仿真

从在Engee中模拟后得到的文件juliaCodeIntegration_noisy.csv 中提取在Simulink中运行模拟所需的变量time 和data ：

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块卡尔曼滤波器。 块 **Engee 函数的精度并不比内置Simulink块 Kalman 滤波器低。

结论

本例演示了如何将Julia代码集成到Engee模型中，以及使用Engee Function块的特殊性--处理多维信号、传递参数、连接外部源代码文件。

示例中使用的块