几种模型降阶方法的仿真对比研究

① 系统仿真技术主要有哪几种

1.基本概念 所谓系统仿真，就是根据系统分析的目的，在分析系统各要素性质及其相互关系的基础上，建立能描述系统结构或行为过程的、且具有一定逻辑关系或数量关系的仿真模型，据此进行试验或定量分析，以获得正确决策所需的各种信息。 2、系统仿真的实质 (1)它是一种对系统问题求数值解的计算技术。尤其当系统无法通过建立数学模型求解时，仿真技术能有效地来处理。 (2)仿真是一种人为的试验手段。它和现实系统实验的差别在于，仿真实验不是依据实际环境，而是作为实际系统映象的系统模型以及相应的“人造”环境下进行的。这是仿真的主要功能。 (3)仿真可以比较真实地描述系统的运行、演变及其发展过程。 3、系统仿真的作用 (1)仿真的过程也是实验的过程，而且还是系统地收集和积累信息的过程。尤其是对一些复杂的随机问题，应用仿真技术是提供所需信息的唯一令人满意的方法。 (2)对一些难以建立物理模型和数学模型的对象系统，可通过仿真模型来顺利地解决预测、分析和评价等系统问题。 (3)通过系统仿真，可以把一个复杂系统降阶成若干子系统以便于分析。 (4)通过系统仿真，能启发新的思想或产生新的策略，还能暴露出原系统中隐藏着的一些问题，以便及时解决。 编辑本段4、系统仿真方法 系统仿真的基本方法是建立系统的结构模型和量化分析模型，并将其转换为适合在计算机上编程的仿真模型，然后对模型进行仿真实验。 由于连续系统和离散(事件)系统的数学模型有很大差别，所以系统仿真方法基本上分为两大类，即连续系统仿真方法和离散系统仿真方法。 在以上两类基本方法的基础上，还有一些用于系统(特别是社会经济和管理系统)仿真的特殊而有效的方法，如系统动力学方法、蒙特卡洛法等。 系统动力学方法通过建立系统动力学模型(流图等)、利用DYNAMO仿真语言在计算机上实现对真实系统的仿真实验，从而研究系统结构、功能和行为之间的动态关系。

② 汽车发动机开环和闭环的区别

动态系统建模被各领域广泛应用，例如电动汽车，能源系统，航空航天。我们本文提到动态系统主要是被控对象，对被控对象进行建模是因为我们希望了解这个系统（被控对象）的物理特性以及接受一些外部输入（力，扭矩，电流等等）时会有什么样的动态响应，基于此从而可以更好的给出控制输入得到我们期望的系统的输出，以及理解系统的退化或最大化提升系统效率。

这些动态系统的行为是由多物理场复杂的交互作用决定的，因此系统行为和系统响应建模通常需要复杂的第一原理支撑，仿真时也需要大量的计算(例如有限元模型)。

这也是本文的出发点，提供数据驱动（主要介绍深度学习和系统辨识）的模型降阶（Reced Order Modeling）提速的方法，通过数据得到具有一定保真度的数据模型，在捕捉到系统动态特性的同时也提升仿真速度。

本文中将涉及多个 demo，数据以及脚本文件，若您感兴趣进一步获取这些链接，可以在文末填写反馈问卷，获取这些链接。

动态系统

动态系统包含状态空间 S, 时间集 T 和一个映射（规则）来描述状态随时间的演变规则 R:S×T→S。例如给定一个时刻 t 的状态 st ,通过这个规则可以计算后面一个或几个时刻状态 st+1=R(st)，st+2=R(R(st)) 等等。动态系统通常可以用随时间变化的方程或方程组来描述。尤其对于连续时间系统，可以通过微分方程来表示。

我们先看一个简单的常微分系统（ODE）, [链接1]

其中 y(t) 是系统状态。例如：一个简单二自由度线性系统，

其中 A 是一个 2x2 的矩阵。初始条件 x0 = [2; 0], 可以通过求解 ODE 方程得到的相应的二维输出 x(t)，包含两个状态，时序和动态图如下：

图表 1 系统输出x(t):

x0 = [2; 0];

A = [-0.1 -1; 1 -0.1];

trueModel = @(t,y) A*y; % 定义系统函数，此处就是一个状态空间方程

numTimeSteps = 2000;

T = 15;

odeOptions = odeset(RelTol=1.e-7);

t = linspace(0, T, numTimeSteps);

[~, xTrain] = ode45(trueModel, t, x0, odeOptions);

后面（在介绍 Neural ODE 部分）我们会尝试利用这个系统的数据 xTrain 进行深度学习模型的训练来得到这个系统的数据代理模型（Surrogate Model），这种思路也可以同样用于复杂系统。

既然数据驱动，有很多机器学习和深度学习算法可以用，那动态系统建模有什么特殊性呢，不是都适用吗？

目前工程中已经用到很多稳态（静态）模型。例如在发动机排放标定，通过 DoE 试验时我们会将发动机维持在不同的稳态工况（转速恒定，扭矩恒定等等），通过试验数据建模得到用于标定的稳态数据模型[链接2]。

稳态工况下，对于方程（1）这样一个简单系统，其中 y(t)' 可以看作 0，到达平衡点，于是 y(t) 和 u(t) 关系恒定，不再在时间维度上与历史状态 y(t-1)，y(t-2) 等等有关，因此稳态模型针对稳态工况是非常准的。

而在瞬态工况下通常 y(t)' 非零，因此方程（1）在求解系统输出 y(t) 时不仅由当前时刻的输入 u(t) 决定，还取决于 y(t-1)，y(t-2)，u(t-1), u(t-2) 等等，这就是动态系统的特殊性，当前输出不仅依赖于当前的输入，还依赖于系统过去的行为（历史输入和历史输出）。我们在下一篇关于系统辨识的文章中会基于示例详细说明这一点。

不考虑动态系统，单纯从系统建模来说，通常有如下的两个方向：基于第一原理的和基于数据驱动的。

第一原理建模是领域工程师都比较熟悉的，例如可以使用 M 脚本语言，Simulink 或 Simscape 建模语言从物理原理进行系统模型的创建。

那什么情况下搭建系统会考虑使用或部分使用数据代理模型（Surrogate Model）？

• 物理系统原理比较复杂或者不够清晰，无法构建第一原理模型

• 数据获取相对简单

• 第一原理物理模型过于复杂，求解耗时，在控制开发时有时效要求，需要加速仿真计算

可以看到右半部分主要是基于数据驱动的建模手段，其中针对动态系统的建模主要是系统辨识和神经网络/深度学习。也是本系列两个方向。

本文介绍神经网络的几种用于动态系统建模的模型，下一篇文章会介绍系统辨识的几种模型。

前馈神经网络（Feedforward Neural Network）

前神经网络大家都不陌生，各神经元分层排列。如图2

图2 简单前馈神经网络示意图


每层神经元只接受上层输出，结果只传给下一层，没有反馈。稍微复杂点的如图3。

图3 squeeze net MATLAB 示例

前馈神经网络是相对于循环神经网络（Recurrent Neural Network）而言，后者具有反馈，后面我们也会介绍。


那么问题来了，前馈神经网络因为没有反馈，如何表达动态系统（Dynamic System）的时间状态依赖呢？

通常主要手段使用不同尺度的滑动窗口来构建衍生特征从而表征系统在时间上的动态。

电池 SoC 预测示例

我们以电池管理系统为例，通过使用深度学习来估计电池的荷电状态 SoC（State of Charge）。

本示例主要介绍 SoC 数据驱动的建模方法，在即使不清楚电池电化学模型以及物理非线性特性的情况下，依然可以进行 SoC 估计。我们可以通过实验室中收集到的实测数据进行一个前馈神经网络代理模型训练实现 SoC 估计[1]。

脚本和数据[链接3]

图表4 数据集预览和模型的5个输入1个输出

数据准备

我们利用实测数据通过预处理得到训练数据集 cdsTrain，我们看其中一条数据记录：

preview(cdsTrain)

ans = 1×2 cell array

{5×669956 double} {1×669956 double}

数据包含五个特征，分别是电压 V、电流 I、温度 T、滑动平均电压 V_avg、滑动平均电流 I _avg,其中后面两个衍生特征用于表征动态信息。

搭建模型

接下来我们构建神经网络模型，结构比较简单，三层全连接网络，一个输出的 Feedforward Neural Network模型。

layers =[sequenceInputLayer(numFeatures,"Normalization","zerocenter")

fullyConnectedLayer(numHiddenUnits)

tanhLayer % HyperbolicTangent

fullyConnectedLayer(numHiddenUnits)

leakyReluLayer(0.3) % 激活函数

fullyConnectedLayer(numResponses)

clippedReluLayer(1) % 激活函数

regressionLayer];

设置训练选项并训练

options =trainingOptions('adam', ... % Adam optimizer

'MaxEpochs',Epochs,...

'ExecutionEnvironment','cpu',...%可以选择GPU

'InitialLearnRate',InitialLearnRate, ...

'LearnRateSchele','piecewise', ...

'LearnRateDropPeriod',LearnRateDropPeriod,...

'LearnRateDropFactor',LearnRateDropFactor,...

'ValidationData', {X,Y}, ...

'ValidationFrequency',validationFrequency,...

'MiniBatchSize',miniBatchSize, ...

'CheckpointPath', NET_Path);

图表 5 训练过程 Loss 变化

导入测试数据验证模型在测试集上的准确度

图表 6 深度学习模型的预测值与实测值比较

仿真测试以及代码生成

Simulink 中的深度学习推断模块[链接8]支持将我们训练好的模型作为 block 参数, 一起作为被控对象集成到整个电池管理系统中。

图表 7 电池管理系统和电池的系统模型

图中所示的 BMS 的 Simulink 模型可以监控电池状态，确保运行安全，还有一个电池模型用于仿真电池的动态和负荷。上面训练的深度学习 SoC 预测器和其他电池平衡逻辑 Block 一样嵌入在 BMS 中可以闭环仿真，以及后面做代码生成与硬件在环。

图表 8 Simulink 中原生的用于深度学习推断的 Block 和 BMS 中的闭环测试 SoC 预测效果

图表 9 deep learning 模块的 C 代码生成

上面这个示例我们利用衍生特征（通过时间滑窗构建时域依赖的特征，也可以通过不同长度的滑窗构建多尺度的时域特征），将前馈神经网络用于动态系统的建模。

Temporal convolutional network (TCN)

TCN 的主要构成是一个扩展因果卷积层。任何一个时刻的计算输出都是基于多个历史时刻的输入。

它用于构建动态系统的逻辑和前面介绍的衍生特征是类似的，都是考虑了历史的多个时间步的输入，因此也可以用于建模动态系统。

TCN 从之前的时间步构建依赖，通常需要将多个卷积层叠加在一起。为了获得更长期的依赖关系，卷积层的膨胀因子呈指数级增加，如下图所示。

假设第 k 个卷积层的膨胀因子为2⁽ᵏ⁻¹⁾ ,步长为 1，则该网络的考虑到的依赖的时间窗的大小可计算为 R=(f-1)(2ᵏ-1)+1,其中 f 为过滤器大小，K 为卷积层数。图中对应的 f=2，K=4, 于是 R=16, 也就是当前时刻输出可以考虑到前面 16 个时刻步长输入。

与循环网络（RNN）相比，TCN 的缺点之一是它在推理期间占用更大的内存。计算下一个时间步需要整个原始序列。下图是一个经典的 TCN 模型结构（结合残差网络）：[链接4]

for i = 1:numBlocks

dilationFactor = 2^(i-1);

layers = [

convolution1dLayer(filterSize,numFilters,DilationFactor=dilationFactor,Padding="causal",Name="conv1_"+i)

layerNormalizationLayer

spatialDropoutLayer(dropoutFactor)

convolution1dLayer(filterSize,numFilters,DilationFactor=dilationFactor,Padding="causal")

layerNormalizationLayer

reluLayer

spatialDropoutLayer(dropoutFactor)

additionLayer(2,Name="add_"+i)];

% Add and connect layers.

lgraph =addLayers(lgraph,layers);

lgraph =connectLayers(lgraph,outputName,"conv1_"+i);

一般的 TCN 架构(如[1]中所述)由多个残差块组成，每个残差块包含两组具有相同扩张因子的扩张因果卷积层，然后是归一化层、ReLU 激活层和空间 dropout 层。

网络将每个块的输入加到该块的输出上(当输入和输出之间的通道数量不匹配时，对输入进行 1 × 1 的卷积)，并应用最终的激活函数。

循环神经网络: LSTM/Gru

循环网络的结构与前馈神经网络不同，它带有一个指向自身的环，用来表示它可以传递当前时刻处理的信息给下一时刻使用，我们选取 LSTM 来介绍，结构如下，

LSTM 之所以可以用于动态系统建模，是因为 LSTM 和动态系统有类似的特性：对于时刻 t, LSTM 使用当前网络的状态 (ht-1,ct-1) 和当前的输入 xt来计算网络输出 ht，同时更新当前网络的状态 (ht,ct)，ht 也叫输出状态，就是当前时刻 t 的 LSTM 网络的输出，ct 称为 cell state,包含学习来的历史时刻的状态信息。每个时间步 LSTM 都会对 ct进行更新：添加信息或移除信息，主要通过四个门函数(f,g,i,o)来实现，将前面的计算过程示意性的描述一下：

其中，上面方程中的 Wf,g,i,o, Rf,g,i,o, bf,g,i,o是这些门函数各自的可学习参数，主要包括针对输入 xt 的权重，针对上一时刻输出状态（同时也是当前时刻的输入）ht-1 的权重，以及偏置，这些可学习参数本身是无状态的，被所有时刻共享。训练的过程就是优化这些学习参数的过程。


总结一下：当前时刻的输出 ht不仅依赖于当前的输入xt，还依赖于当前的状态 (ht-1,ct-1)。它可以根据训练选择性的记住每一时刻的“重要”信息，并且利用这个信息作为状态，结合当前输入 xt 进行输出状态 ht预测。

电机温度预测示例

接下来我们使用数据结合 LSTM 模型来搭建永磁同步电机（PMSM）的一个代理模型（Surrogate Modeling），用于电机不同位置的温度预测。

类似 Demo [链接5]

理解数据集

数据集来自多片时长不同的数据。从数据中可以看到 PMSM 不同位置温度与电气系统和热系统对应的工况参数的相互影响。

同时环境温度的变化也会对电机不同位置温度以及相应需求扭矩所需的电流电压有影响。

我们使用环境温度、冷却液温度、电压、电流作为输入，输出为不同位置 PMSM 的温度。

数据预处理与特征工程

同样作为一个动态系统，我们通过对初始数据再进行不同尺度滑窗实现衍生特征生成，所有特征结合在一起作为 LSTM 的输入（尽管 LSTM 本身也具有考虑时间依赖关系的特性）。

% create derived features using raw voltages and currents derivedinputs =computedrivedfeatures(tt_data); % check the noise in the data tt_data=[tt_data derivedinputs]; Vnames=tt_data.Properties.VariableNames; s1=620;s2=2915;s3=4487;s4=8825; % preprocess exponentially weighted moving average [t_s1,t_s2,t_s3,t_s4]=preprocmovavg(tt_data,s1,s2,s3,s4,Vnames); % preprocess exponentially weighted moving variance [t_v1,t_v2,t_v3,t_v4]=preprocmovvar(tt_data,s1,s2,s3,s4,Vnames); % attach features to the original table predictors=[tt_data,t_s1,t_s2,t_s3,t_s4,t_v1,t_v2,t_v3,t_v4,tt_profileid]; responses=[tt(:,9:12) tt_profileid]; VResponsenames=responses.Properties.VariableNames;

准备训练数据和验证数据

holdOuts =[657258]; % define profiles that are withheld from training. [xtrain,ytrain]= prepareDataTrain(predictors,responses,holdOuts);

我们将使用第 58 条记录作为验证集，其中包括 4.64 小时的数据。

validationdata =58; [xvalidation, yvalidation]= prepareDataValidation(predictors,responses,validationdata); numResponses = size(ytrain{1},1); featureDimension = size(xtrain{1},1); numHiddenUnits=125;

上面的 DAG 网络可以对许多复杂的系统进行建模。通过上面的结构（左右两支）帮助我们对依赖于时间历史的物理行为动态行为以及含状态的行为进行建模。

查看预测结果

如图所示，红色和蓝色分别代表了实测数据和模型预测结果，右侧图像显示的是他们的残差，可以看到大部分误差在 1% 以内，预测效果比较理想。模型在瞬态变化较快和较慢的工况下都能和实测数据保持一致，说明模型也保持了一定的保真度。

将模型导出 Simulink

我们将训练好的模型保存为 .mat 文件，并将其导入 Simulink Deep Neural Network Predict 模块[链接9]，这样我们就有了一个只有 50Kbyte 大小可以预测温度的代理模型（Surrogate Model）用于仿真。

Neural Ordinary Differential Equations：神经网络 ODE

这是 NIPS 2018 年最佳论文[2]提出的一种新的网络层。当然这个要理论上介绍还是比较复杂，我们可以通过应用场景直观的解释一些 Neural ODE 如何实现动态系统建模。

试想我们有一个动态系统，因为系统动力学过于复杂，我们没有真正的物理模型，但我们可以不断地通过测量得到系统的初始状态 y(t0)，动态输 入 u(t0), u(t1),…, 与动态输出 y(t1),y(t2),…,y(tn)。接下来我们想是否可以实现这样一个微分方程：

使得这个微分方程正好代表了我们的系统，也就是说在系统输入 y(t0) 的初始条件下通过求解(例如使用 ode45)这个微分方程得出的解 yp(t) 和我们实测结果是吻合的。但如何基于已有的数据 y(t0),序列 u(t) 和序列 y(t) 得到 f 呢？如果我们将 f 用一个神经网络 F(θ) 替代，即

我们现在有数据 u(t), y(t)。我们不断地利用数据训练参数θ，使得上述方程的解 yp(t1), yp(t2),...,yp(tn) 与实测 y(t)是吻合的，那就可以得到这个动态系统的微分方程模型了，从而可以用于后续系统仿真与预测。如何理解 yp(t) 与实测 y(t) 的接近度,也就是损失函数? 我们简单介绍一下训练时 Loss 函数计算。对于方程，我们在知道系统初始状态 y(t0)，可以通过很多数值积分求解器（例如常用的 ode45）得到任何时刻的推断输出 yp(t) (当然前提是系统的输入 u(t)也是已知).

θ 是神经网络 F 的静态参数。对于所有时刻都是不变的。我们就可以方便得到损失函数的值

其中 L 可以是任何自定义的损失函数。于是我们可以进行参数 θ 的训练。当然关于梯度计算与反向传播会有相对复杂的数学推导，论文提出了伴随方法（Adjoint Method）来实现这一过程，此处不做详细论述。

对应于上述过程，MATLAB中提供了dlode45 [链接10]，用于建模方程右侧非线性函数F的同时，计算 ODE 的时序解。

即 dlode45 接收一个含参神经网络 F(θ)、需要计算输出结果的时刻序列 [t0, t1,…,tN]、系统的 t0 时刻的初始状态 yt0、神经网络参数的一组值，就可以计算出时刻 [t1,…,tN] 所对应的输出状态。

使用 Neural ODE 为系统建模示例

我们通过一个示例介绍如何使用Neural ODE为动态系统建模。[链接6]

我们就借用文章刚开始的简单二自由度线性系统，x' (t)=Ax(t), 其中A是一个 2x2 的矩阵。我们用这个已知的系统产生一些数据，利用这些数据来训练一个 Neural ODE 的方程，使得这个基于数据训练好的系统（Neural ODE 方程）能够接近已知的这个动态系统。

生成物理系统的数据 xTrain 作为真值

x0 = [2; 0];

A = [-0.1 -1; 1-0.1];

trueModel = @(t,y)A*y;

[~, xTrain] =ode45(trueModel, t, x0, odeOptions);

xTrain 两个自由度的可视化


定义和初始化神经网络 F(t,x(t),θ) 的参数 θ

neuralOdeParameters.fc1= struct;

sz = [hiddenSizestateSize];

neuralOdeParameters.fc1.Weights= initializeGlorot(sz, hiddenSize, stateSize);

neuralOdeParameters.fc1.Bias = initializeZeros([hiddenSize 1]);

neuralOdeParameters.fc1

ans = 包含以下字段的 struct:


Weights: [20×2 dlarray]

Bias: [20×1 dlarray]

同样

neuralOdeParameters.fc2

ans = 包含以下字段的 struct:

Weights: [2×20 dlarray]

Bias: [2×1 dlarray]

定义神经网络模型 F(t,x(t),θ) 函数

function y = odeModel(~,y,theta)

y =tanh(theta.fc1.Weights*y + theta.fc1.Bias);

y =theta.fc2.Weights*y + theta.fc2.Bias;

end

结合定义好的 F(t,x(t),θ) 作为 dlode45 的输入来构建代理模型函数

function X =model(tspan,X0,neuralOdeParameters)

X = dlode45(@odeModel,tspan,X0,neuralOdeParameters,DataFormat="CB");

end

定义模型梯度函数

主要用于训练过程计算损失以及对应待训练参数的梯度

function [gradients,loss] =modelGradients(tspan,dlX0,neuralOdeParameters,targets)

% Compute predictions.

dlX = model(tspan,dlX0,neuralOdeParameters);

% Compute L1 loss.

loss =l1loss(dlX,targets,NormalizationFactor="all-elements",DataFormat="CBT");

% Compute gradients.

gradients =dlgradient(loss,neuralOdeParameters);

end

训练模型

不断地迭代训练，创建 miniBatch,并进行损失函数计算和自动微分，通过调用 adam 求解器进行参数学习

for iter=1:numIter

% Create batch

[dlx0, targets] = createMiniBatch(numTrainingTimesteps,neuralOdeTimesteps, miniBatchSize, xTrain);

% Evaluatenetwork and compute gradients

[grads,loss] = dlfeval(@modelGradients,timesteps,dlx0,neuralOdeParameters,targets);

% Update network

[neuralOdeParameters,averageGrad,averageSqGrad] =adamupdate(neuralOdeParameters,grads,averageGrad,averageSqGrad,iter,...

learnRate,gradDecay,sqGradDecay);

% Plot loss

currentLoss =double(extractdata(loss));

测试模型

选取新的初始条件作为训练好的模型的输入，来进行和物理系统输出的对比

x0Pred1 =sqrt([2;2]);

x0Pred2 =[-1;-1.5];

x0Pred3 = [0;2];

x0Pred4 = [-2;0];

可以看到模型对于新的初始条件依然表现优异。因此神经网络 ODE 在构建动态系统上很有潜力，目前在发动机建模上也有一些示例应用。

NARX（nonlinear autoregressive network with exogenous inputs）反馈神经网络

在深度学习网络爆发之前，在浅层神经网络的场景中，NARX 反馈神经网络是经常用于动态系统建模的，它具有反馈连接，即输出 y(t) 依赖于系统之前时刻的输出 y(t-1)，y(t-2) 等等和输入。

其中网络 F 主要是通过一个前馈神经网络实现。

如下图，其中F即为第一层隐含层前馈网络，其中输入和输出可以是多维的。

跟其他反馈神经网络类似（前面介绍的 LSTM）,它的训练过程与推断过程有些区别。网络本身的输出需要被反馈到前馈神经网络的输入。在训练时，因为我们能够拿到整段输出真实的数据，因此，我们会用当前时刻真实的输出值作为训练时模型输入而不是反馈预测的输出，换句话说，在训练时我们会把网络作为开环去训练。

这有两个好处。一是前馈网络的输入更准确。第二，生成的网络具有纯前馈架构，静态反向传播变得可用。当我们进行未来多步推断时，因为这种情况我们只能用推断的数据进行下一时刻预测，所以这次我们才把网络闭环，用于推断。

磁悬浮系统的 NARX 建模示例

接下来我们使用 narx 神经网络来对动态系统进行建模的示例[链接7]。示例系统是一个磁悬浮系统。目标是控制悬浮在电磁铁上方的磁铁的位置，在电磁铁的位置上，磁铁受到限制，只能在垂直方向上移动，如下图所示，

系统的运动方程为：

其中 y(t) 是磁铁在电磁铁上方的距离，i(t) 是经过电磁铁的电流，M 是磁铁的质量，g 是重力常数。其中 β 为粘性摩擦系数，由磁体运动材料决定；α 为场强常数，由电磁铁上导线匝数和磁体强度决定。我们搜集了系统输入 u(t)-施加在电磁铁上的电压和系统输出 y(t)-永磁体的位置，对应两个时间序列。

搭建网络和准备数据

d1 = [1:2];

d2 = [1:2];

narx_net =narxnet(d1,d2,10); % 使用narxnet功能创建NARX开环串联网络，10个隐藏层神经元

[p,Pi,Ai,t] =preparets(narx_net,u,{},y); % 用preparets准备数据

其中 y(t) 是一个反馈信号，既是输入，也是输出，训练时我们既可以拿到当前时刻数据，也可以拿到后面时刻的数据，所以可以用于开环训练，当推断时我们会将输出接到输入作为闭环进行推断。


训练网络

narx_net =train(narx_net,p,t,Pi); % 训练网络得到训练好的网络

验证开环网络推断效果

yp =sim(narx_net,p,Pi);

e =cell2mat(yp)-cell2mat(t);

可以看到误差很小。因为我们用的开环训练，所以推断结果是用前面时刻的真实输出数据（而非推断输出反馈），所以这里的误差是 one-step-ahead 推断误差。


测试闭环推断效果

如果要看网络真实准确度的表现，需要将开环的输出作为反馈接到输入，然后进行多步预测。

narx_net_closed =closeloop(narx_net);

现在可以使用闭环执行 900 个时间步的迭代预测。在这个闭环网络中只需要两个初始输入和两个初始输出作为初始条件。

y1 = y(1700:2600);

u1 = u(1700:2600);

[p1,Pi1,Ai1,t1] = preparets(narx_net_closed,u1,{},y1);

yp1 =narx_net_closed(p1,Pi1,Ai1);

从闭环预测的结果看，蓝线是磁铁的实际位置，红线是 NARX 神经网络预测的位置。即使网络预测的时间步预测了 900 步，预测依然是非常准确的。闭环多步预测准确的前提就是开环下单步推断误差要小。


总结

本文主要介绍了动态系统的特性和用于动态系统建模的神经网络模型从而可以实现模型降阶（Reced Order Modeling），包括前馈神经网络，TCN，循环神经网络，神经网络 ODE, NARX 网络和相应的一些手段，并结合示例与场景进行了说明这些手段的有效性。

在后续文章我会来介绍系统辨识（System Identification）的一些示例,欢迎继续关注。

本文中涉及到的多个 demo，数据以及脚本文件，若您感兴趣，可以通过扫描填写下面这个反馈问卷，或点击”阅读原文“进一步获取这些链接。

获取文中示例链接

参考文献


[1] Vidal, C., Kollmeyer, P., Naguib, M., Malysz, P. et al., “Robust xEV Battery State-of-Charge Estimator Design Using a Feedforward Deep Neural Network,” SAE Technical Paper 2020-01-1181, 2020, doi:10.4271/2020-01-1181.

[2] Ricky T. Q. Chen*, Yulia Rubanova*, Jesse Bettencourt*, David Duvenaud University of Toronto, Vector Institute “Neural Ordinary Differential Equations”

编辑：谢雅洁 校对 ：向映姣

③ 降阶模型的介绍

降阶模型是指对于用状态空间方法表达的模型，采用模型集结的方法降低状态空间模型的阶数，所获得的低阶模型，或对于用微分方程、差分方程或时间序列分析等方法建立的模型，忽略其高阶项而获得的低阶模型。

④ Krylov子空间迭代法二十世界比较经典的算法，但是还是不了解

摘要:算法比较研究，比较几种主要模型降阶方法的优缺点，为给工程应用提供方法参考。利用奇异值分解的模型降阶方法具有较好的理论性质，能够保持降阶系统结构特性，但计算成本较高故不适合大规模动态系统的降阶;采用矩匹配的模型降阶方法计算简便，适合大规模系统降阶，但无法保证降阶系统稳定性，也很难求得降阶误差界。最小二乘降阶法同时利用了系统的Gramian矩阵和Krylov子空间理论，结合了二者的优点，使得降阶过程计算简化，保持了降阶系统的结构特性，而且降阶误差进一步减小。仿真算例证明了最小二乘法较前两者具有优越性。

⑤ 查晓明的科研课题

主要项目如下：
1.国家自然科学基金面上项目(批准号：51177113)：基于微分几何同调的微电网等值建模理论与模型降阶方法
2. 国家自然科学基金面上项目(批准号：51277137)：微电网逆变器交互作用分析及建模方法研究
3. 国家自然科学基金青年项目(批准号：51107091)：一种基于部分单元能量回馈的级联多电平逆变器研究
4. 国家自然科学基金青年项目(批准号：51207115)：基于脉冲微分动力系统理论的风电功率预测的极大误差评估方法研究
5. 国家自然科学基金青年项目(批准号：51307126)：一种新型级联多电平变换器的高压大功率电机准矢量控制方法研究
6. 高等学校博士学科点专项科研基金(课题编号：20120141110074)：基于多自由度超越摄动的统一接口微电网复杂系统建模理论与等效方法
7. 国防973项目课题1项
8. 湖北省自然科学基金一般项目(批准号：2011CDB263)：基于部分单元能量回馈的级联多电平变换器研究
9. 2012年湖北省科技攻关计划项目：220V直流微电网多样性配电平台产业化研究
主要参与的项目如下：
1. 国家自然科学基金重大项目课题(批准号：51190102)：大规模风电场多时空尺度聚合建模理论和方法
2. 科技部973项目课题1项(批准号：2012CB215101)：大规模高集中度风电场出力多时空尺度爬坡特征分析、预测与控制
其他还主持承担了863重大项目子题，湖北省科技攻关项目，武汉市科技攻关项目，国家电力公司青年科技促进费项目，武汉市青年晨光计划项目，以及各级电网公司、电力科学研究院和电力电子产品生产制造公司多项委托项目。在国内外重要会议和期刊上发表教学和科研论文70多篇，其中SCI收录6篇，EI收录30余篇;获发明专利14项、实用新型专利11项。在编《电能质量控制技术》专着一部。 1、期刊文章
1) 查晓明、孙建军、陈允平，并联型有源电力滤波器的重复学习Boost变换控制策略，电工技术学报，2005年
2）查晓明、陈允平，无源与有源结合电力滤波器抑制谐波振荡的重复学习Boost变换控制，电工技术学报，2004.6
3) 孙建军、查晓明、丁凯、余利生，一种低耐压中性线谐波治理电路原理及实现，电力系统自动化，18期，2003年，pp34-37
4) 邓恒、殷波、查晓明、陈允平，一种改进的PWM方法在APF中的应用，电力系统自动化，14期，2002年，pp51-54
5) 黄敏、查晓明、陈允平，并联型电能质量调节器的模糊变结构控制，电网技术，26卷7期，2002年，pp11-14；
6) 查晓明、陈卫勇、陈允平，三相交流系统综合补偿的补偿分量检测研究，电网技术，23卷7期，1999年，pp24-26
7) 查晓明、陈一尧、周巍，LMS自适应信号处理在快速多频率阻抗测量中的应用，仪器仪表学报，19卷2期，1998年，pp211-214
8) 殷波、邓恒、陈允平、查晓明，a-b坐标系下瞬时无功功率理论与传统功率理论的统一数学描述及物理意义，电工技术学报，2003.5
9) 查晓明、王瑾，基于PWM控制的并联型有源电力滤波器的MATLAB仿真研究, 电力系统及其自动化学报, 2001年
10）查晓明、王瑾，串联有源滤波器拓扑结构的MATLAB仿真，电力系统及其自动化学报，2001年
11）殷波、查晓明，综合电能质量控制器，电网技术，2002
12）殷波、查晓明，一种新的降低三相有源滤波器开关损耗的方法，电力电子技术，2002
13）孙建军、查晓明，基于双DSP结构的有源滤波器检测及控制系统，武汉大学学报，2001
14）湛秀平、查晓明，基于AT89C2051的步进电动机驱动电路，电气传动自动化，2001
2、会议论文
1）查晓明、陈允平, An integrated learning control method for PWM-VSI based hybrid filter to dampen harmonic resonance in instrial power system, Proceedings of IEEE Annual Conference on Instrial Electronics (IECON’04), 2004
2) 查晓明、孙建军、陈允平，2003年，Application of iterative learning control to active power filter and its robustness design, Proceedings of IEEE power electronics specialist conference (PESC’03), pp785-790；
3）查晓明、陶骞、孙建军、陈允平，2002年，Development of iterative learning control strategy for active power filter, Proceedings of IEEE Canadian conference on electrical & computer engineering (CCECE’02), pp245-250;
4) 查晓明、刘开培、陶桂洪，2002年，An active AC bridge with adaptive current orthogonal decomposition for precision measurement of dielectric loss angle, Proceedings of IEEE conference on precision electromagnetic measurement (CPEM’02), pp182-183；
5) 陈允平、查晓明、王瑾、孙建军、汤洪海、刘会金，2000年，Unified power quality conditioner (UPQC): the theory, modeling and application, Proceedings of power system technology (PowerCon 2000), pp1329-1333;
6) 查晓明、陈允平, The iterative learning control strategy for hybrid filter to dampen harmonics resonance in instrial power system, Proceedings of IEEE International Symposium on Instrial Electronics (ISIE’03);
7) 陈允平、查晓明、赵磊The Measurement of Impulse Grounding Resistance and Inctance of Grounding Network, PP352-355(VOL.1), Proceeding of Energy Management and Power Delivery (EMPD’98), Singapore, Mar. 1998

⑥ 怎样验证仿真模型的可靠性

三套方法来验证。
第一种，用不同的软件（例如分别用ansys和abaqus）计算同一个问题，模拟得到的结果非常近似（不太可能保证完全一样）。这样就是让人信服的。
第二种，先证明你的方法是正确的，比如用你的方法去做一个别人已经做过的结果（例如已发布的文献中提到的模型和结果，你的模型和他近似，过程自己来操作，得到的结果和他一致，这样就可以说明你的方法是可信的）。然后就用这个证明过的操作方法，去完成需要完成的项目，得出的结果也较为可信。
第三种，试验模型验证。做一个实际比例模型，约束和加载与数值仿真的完全一致，然后比较试验结果和数值结果，曲线走向和趋势基本一致，数量级一致等等等等，就是可信的。

三种方法比较。
第三种最好，但难以实现，一般在研究所里有条件才采用，耗时耗力；第二种其次，最为简单，并且与权威杂志的结果有个比较，比较有说服力；第三种也可以，但一般是分别分给两个人用不同软件计算，或者两个人在没有交流的情况下用同一软件计算，这样的结果才具有一定可信性，在正式项目中一般不会采用，或在团队中具有较高资质的仿真工程师的情况下采用，毕业生或实习生的结果是不可信的。

⑦ POD流场降阶的意义在于什么

加速流场计算。CFD这块，计算量庞大是大家公认的，无论是用RANS计算整个飞机，还是用DNS计算湍流，计算耗时都不可忽略。发展降阶模型，一个重要的原因是，我们可以通过计算CFD的少数状态，就能够用降阶这种数据处理手段，把其他状态的流场特征外推或者内插出，而且这个过程的耗时相比CFD重新计算而言，绝对是可以忽略的。尤其是气动弹性问题，相信做过的都有体会，在没有先验知识的前提下，要准确捕捉颤振边界，基本都是需要不断试错（试动压），观察响应情况。要注意气动弹性仿真是需要非定常CFD求解器和CSD（计算结构力学）相互迭代的，所以这个计算量必然大于单独的CFD求解。有了降阶模型之后，我们只需要CFD仿真一段样本，之后的CFD和CSD迭代过程，可以转化为ROM和CSD求解器迭代（对于线性结构，可以直接根据耦合的线性系统特征值得到稳定性），这个效率提高确实很可观。另外，随着流场仿真的需求增加，这个效率的提高会越来越明显。在别的领域有叫PCA的，有叫KL变换的，其实是一个东西。个人感觉这种pod类的模态分解方法，其本质是提供一组低维的坐标系，在这组新的坐标系下，我们可以更加简洁的表达流场。至于说用cfd几个参数的，我想要实现上面的两个意义都是不大可能的吧。总结一点，降阶模型的核心其实是抓主要矛盾，主要矛盾抓好了，既能代替原来的复杂，高阶系统，实现仿真，预测，控制等功能，还能帮助我们了解物理问题的主要特征，因此还是很有价值的。最近机器学习的火也烧进了这一领域，凭借大数据处理能力，来优化流场控制过程。大导磨刀霍霍，我等只能慷慨赴之。这一领域，大牛太多，山头林立，我这样的杂鱼出头之日遥遥无期呐。

⑧ 高阶系统降阶的两个主要方法

高阶系统降阶的两个主要方法：用主导极点的方法。

主导极点就是离虚轴最近的一个或者一组极点，而且其余极点距离虚轴的距离远大于这组极点(大一个数量级)。另外主导极点不可以离零点太近。选取出主导极点后，就可以用这组极点近似分析系统，其余极点的影响可以忽略，相当于把系统降阶。

在静态条件

变量各阶导数为零——下描述变量之间关系的代数方程叫静态数学模型；描述变量各阶导数之间的关系的微分方程叫动态数学模型，这些方程就是运动方程。运动方程的阶次是指运动方程中出现的最高阶导数的阶次，如N阶运动方程就是指方程中出现的最高阶导数是N阶的。

⑨ 有人用机器学习的方法解决流体力学问题吗

可以做，但有几个困难：
1，如果训练样本是RANS的结果，最终预测精度不会超过RANS。
2，如果自变量是边界条件（包含几何形状），维度太高，样本不够

比较可行的方案：
1，不代替而是扩展CFD的能力
2，不做可以泛化的模型，缩减状态空间维度，只做具体问题的简单变化

可惜这个方案早有人做了，不过没打着机器学习的招牌，而是叫surrogate-based optimization。一般是在设计方案基础上解放几个参数的自由度，采一些样本（CFD或者实验），做一下回归模型，然后在参数空间里找最优，还是挺有趣的方向。

那机器学习研究湍流这样的美好想法，首先是要建立在对问题理解深刻的基础上。

⑩ 四阶魔方的还原方法之一：降阶法

关于四阶魔方的玩法，用的较多的方法是降阶法。关于降阶法，简单解释就是把高阶魔方设法降解成低阶魔方的模型，再利用低阶魔方的解法来完成。目前魔方界对高阶魔方的解法使用降阶法的较多。 复原第二层：这一步中先复原翼，再复原点。这里是变通地使用三阶魔方第二层棱块的复原公式。
2.1复原第二层的翼
公式2-1：URU’R’U’F’UF
公式2-1’：U’F’UFURU’R’
这个公式与三阶的复原第二层棱块公式完全一样。如果目标块不在顶层，那它一定在第三层，只要将目标块转到FR棱的第三层处，做公式2-1即可将其转到顶层。
2.2 复原第二层的点
公式2-2：Ru’R’U’f’Uf
公式2-2’：Ru’R’U’f’Uf
这个公式与三阶的复原第二层棱块公式类似，不同之处是有一部分转动了中间层。
如果目标块不在顶层，那它一定在第三层，转动第三层，只要将目标块转到F面近R面处，做公式2-2即可将目标块转到顶层（或转动第三层，将目标块放到R面近F面处，做公式2-2'）。
至此，下面两层全部完成。 复原第三、四层：这里是先用换角公式复原顶层角，之后用换棱公式复原顶层翼，再用换棱公式复原第三层翼，最后变化使用换棱公式复原剩下的点块。
3.1复原顶层角块位置
公式3-1：R’ L U L’ U’ R U L U’ L’
公式3-1’：L R’ U’ R UL’ U’ R’ U R
若有两个相邻角块位置正确。将正确的两个角块放在右手，做公式3-1。
有两个不相邻角块位置正确。直接做公式3-1或公式3-1’，便可以得到第2种情况
这与三阶魔方第六步一致。考虑到后续公式的一致性，这里介绍一个新公式，其实这个公式完全可以替代三阶换棱公式（R2 D2 R' U' R D2 R'U R'）。它的变换形式可以一直用到本方法最后一步。
3.2复原顶层角块方向
公式3-2：R’ U’ R U’ R’ U2 R
公式3-2’：R U R’ U R U2 R’
这与三阶魔方第五步顶面全部翻正的小鱼公式完全一样。从三阶魔方考虑，是很容易理解的。
3.3复原三、四层的翼
公式3-3：R’ LuL’u’ R uLu’L’
公式3-3’：L R’u’Ru L’ u’R’uR
这个公式与3-1极为相似只有U大写改为小写。
由于有十二个翼需要替换，建议先复原顶层的八个。找第三层翼含有顶面颜色(黄色)的，将黄色面放在右侧或左侧，做公式3-3(在左侧)或3-3'(在右侧)。
顶层八个翼复原后，看第三层翼，这时第三层翼已经不含顶层颜色，我们找在同一面上两个翼是否有相同的颜色。这有三种情况：
1．四个面的翼颜色都相同，说明第三层已经翼已经完成复原；
2．找到有一对翼在同一面颜色相同，将这对翼放在R面。做（公式3-3→ yU’→公式3-3）。就得到了第1种情况。
注意：这里是把同面同色的翼放在右手，拆一个顶面翼；拆完之后，用同一公式去复原顶面翼，这样就实现了第三层三个翼的旋转互换。
没有找到在同
3.4复原三、四层的点
这时只剩下几个点要对调了，用下面两点和三点对换公式可以很轻松完成。
公式3-6：r U l' U' r' U l U'
公式3-6'：l' U' r U l U' r’ U
做该公式时注意目标块必须放在R面和F面下半区域的对应位置。如果目标块不在下方，转动该面就可以了。做完公式再恢复回去即可。
三点轮换：
公式3-4：r’ L u L u' r u L u' L’
公式3-4’：l R' u’ R u l' u’ R' u R
公式3-5：r’ L d’ L d r d’ L d L’
公式3-5’：l R' d R d’ l' d R' d’ R
我的建议，如果你是初学者，能用两点对换就不要用三点轮换；如果你已经是高手了，随心所欲吧。