数据驱动识别偏微分方程的序列奇异值过滤方法

1.本发明属于计算力学领域，具体涉及一种数据驱动识别偏微分方程的序列奇异值过滤方法。


背景技术：

2.近年来，微型机械和纳米结构被应用于医学诊断、药物研发、电子器件、能源催化等各种研究领域。微型机械可以完成常规方法难以完成的微加工、微操作，纳米材料具有高强度、高韧性等优异的力学性能，因此微纳米材料广泛渗透到物理学、化学、材料学、生物医学等领域，微纳米力学目前已成为交叉性的前沿学科，对微纳米材料力学性质的研究具有重要意义。
3.实验表明，微纳米尺度下材料的力学性能与响应与宏观尺度下有着较大差异，宏观的经典弹性模型不具备描述微纳米材料力学行为的能力。针对微纳米尺度效应，目前提出的各种理论模型是通过引入不同的微观变形机理，从而在经典的二阶微分理论中引入不同的高阶项来预测微纳米材料的力学响应。然而，大量的实验和分子模拟结果表明，这些理论大多只能描述某一种变形机制，并不能准确地预测材料的所有力学行为。建立完整统一的微纳米力学模型，是固体力学领域的重要挑战。
4.研究表明，数据驱动方法是一种从已知数据中建立模型的有效方法。数据驱动直接从观测的数据出发，为系统建立模型，跳过了人为建模再验证的过程。一些具有代表性的系统识别方法近年来被相继提出，如符号回归、深度学习等，然而将上述方法应用于微纳米尺度力学方程的识别时存在许多问题：它们需要人为地预先提供结构化的先验知识，才能进行有监督地学习，而对于微纳米尺度的控制方程，其个数和具体形式(如待定的左端项)实际上都是未知的；此外，符号回归在处理较多的变量数目和函数库数目时，会陷入局部最优以及面临过拟合的问题；以神经网络为代表的深度学习技术在拟合变量之间非线性关系的能力上占据优势，然而深度学习的可解释性依旧是目前学术界没有攻克的难题，数据驱动得到的结果难以被显式地表达，从而限制了人们的理解和进一步利用。上述问题的存在，使得发展适用于微纳米力学系统的数据驱动识别方法目前仍然存在很大的挑战，具有重要的研究价值。
5.最近的研究表明，稀疏学习是一种非线性系统建模的有效方法，其通过将稀疏向量的范数加入到优化目标函数中，从而能解决构建非线性函数库面临的维数暴增所带来的矩阵病态问题。brunton等人于2016年首次将稀疏学习引入到动力系统的方程识别当中(sindy)：对于各个观测量的时间序列数据，使用给定的非线性基函数组成数据矩阵，并使用最小绝对收缩选择算子、阈值最小二乘法等稀疏学习算法寻找稀疏解，以平衡拟合精度与方程复杂程度之间的矛盾。然而该方法也必须预先人为设定控制方程的个数以及为每个方程指定左端项，从而无法用于识别微纳米力学方程。本发明将借鉴其构建非线性函数库并识别稀疏模型的思想，提出奇异值过滤方法，把依赖于先验知识的有监督学习范式转变为不依靠先验知识、自动化的无监督学习，将方法的应用场景扩展至控制方程个数和方程
各项完全未知的偏微分方程组的识别问题，从而使得方法能够识别微纳米材料的力学控制方程。
6.为了能够从数据集中准确地识别微分方程的个数和方程组中存在的所有非线性项，本发明创新性地设计了一套奇异值过滤的算法流程，通过迭代实施奇异值分解方法来识别隐藏在数据集中的本质属性。奇异值分解(svd)是一种矩阵分解技术，主要用于矩阵的特征分解和降维。其作为很多机器学习算法的基石，能够处理模式识别、数据压缩、信号降噪等问题，并且应用于推荐系统、自然语言处理等机器学习领域。奇异值分解是谱分析理论在任意矩阵上的推广，其主要功能是在平方意义下寻找最为线性无关的一组正交向量，从而张成初始矩阵的原象。因此，奇异值分解具有寻找大数据中隐藏的本征线性结构和计算矩阵的低秩近似的能力。依据各个奇异值的大小关系，算法可以自动判断矩阵中各个模态所占的比重，从而识别控制方程的个数，而不依赖人为的预先设定。
7.为了确定每个方程的最简形式，使其在拟合精度与方程复杂度之间做出平衡，本发明拟采用处理列子集选择问题的代表性算法：强秩揭示正交三角分解(srrqr)方法。由于假设存在的非线性函数库会存在大量冗余，随意选取方程的左端项进行回归将无法得到正确的结果。要设计方法使其能够实现自动识别，则这一组左端项应该满足两个条件：一是它们之间不应存在线性相关关系；二是它们应是某种范数意义下最能代表剩余列的一组基向量。这对应线性代数中的列子集选择问题(cssp)，或称为内插分解(id)。其中，srrqr在满足上述两个条件的同时，具有较低的计算复杂度，能够保证方法的效率。
8.目前，对于非线性微分方程组的无监督数据驱动识别的研究还处于发展阶段，应用于微纳米力学方程识别所面临的依赖先验知识、可解释性差等问题尚未被解决。因此，本发明将提出一种数据驱动识别偏微分方程的序列奇异值过滤方法，结合快速高效的先进数值方法，设计鲁棒性强、适用广泛的算法流程，为微纳米力学系统的数据驱动识别提供一种行之有效的途径。


技术实现要素：

9.本发明要解决的技术问题：本发明基于矩阵的奇异值分解技术，并结合强秩揭示正交三角分解方法，针对微纳米力学系统的数据驱动建模设计算法流程，创新性地提出了一种序列奇异值过滤方法(sequential singular value filtering method，简称seq-svf)，其目的在于解决现有技术存在的以下问题：采用奇异值分解克服现有方法需要预先人为设定方程个数的缺点；采用强秩揭示正交三角分解克服现有方法对控制方程结构化先验知识的依赖，使得该方法能够突破其只能应用于动力系统的局限，从而将应用场景扩展到形式完全未知的微纳米力学微分方程组的数据驱动识别问题中；设计自动化的奇异值筛选流程来滤去非线性函数库中的冗余项，减少了矩阵分解的计算复杂度，提高了识别结果的泛化能力。
10.本发明的技术方案：
11.数据驱动识别偏微分方程的序列奇异值过滤方法(seq-svf)，
12.为了从微纳米尺度材料的位移等数据中提取高阶控制微分方程组，提高计算结果的效率与精度，本发明将设计一套基于奇异值分解方法和强秩揭示正交三角分解方法的算法流程，从测量数据中无监督地提取非线性控制方程组，提供了一种数据驱动识别偏微分
方程的序列奇异值过滤方法，具体步骤如下：
13.步骤1、预处理编码器：设置候选函数库，根据质点位移数据和候选函数库计算数据矩阵，对数据矩阵的每列实施正则化，使用奇异值分解方法识别偏微分方程的个数；
14.首先，给出本方法识别微纳米力学系统的基本格式；设该微纳米力学系统由一组偏微分方程控制：
[0015][0016]
其中，x1～xn代表所有假设存在的项，包括质点位置对于空间坐标的二阶以及更高阶导数项和它们通过相乘得到的非线性项；根据在各个采样点处测得的质点位移，通过数值计算得到数据矩阵a所有元素的大小；
[0017]
接着，为了防止各项的数量级相差过大对计算结果造成影响，对数据矩阵a的每一列实施l2范数的归一化；最后，对数据矩阵a实施奇异值分解，通过判断矩阵的小奇异值的数量来确定偏微分方程的个数；
[0018]
步骤2、函数过滤器：删除候选函数库中冗余的项；
[0019]
通过删除数据矩阵a中的各列，根据奇异值分解的结果判断该列是否存在于偏微分方程中；具体实现过程如下：
[0020]
(1)初始化变量，令i＝1，k为空向量；p为数据矩阵a中小奇异值的个数；
[0021]
(2)若i≤数据矩阵a的列数n，开始循环；
[0022]
(2.1)删掉数据矩阵a的第i列，做奇异值分解，计小奇异值的个数为pi；
[0023]
(2.2)如果pi＝p-1，令k＝(k i)；否则不实施操作；令i＝i+1，返回步骤(2)继续循环；
[0024]
(3)a
′
＝a(：，k)，即为过滤掉无关项后的数据矩阵，以供后续步骤使用；
[0025]
步骤3、线性解耦器：使用强秩揭示正交三角分解方法识别所有左端项在通过步骤2得到的过滤后数据矩阵(a
′
)中的位置，并更新数据矩阵；
[0026]
步骤4、稀疏解码器：确定各个偏微分方程中的项并计算各项系数；
[0027]
对每个偏微分方程各实施一次步骤2，即可在候选项中确定每个偏微分方程的稀疏模型；之后使用最小二乘法计算所有系数；其具体实施过程如下：
[0028]
(1)初始化变量，令i＝1，j＝1，k
p
×
(n
′‑
p)
＝0；其中n
′
为过滤掉无关项后的数据矩阵a
′
的列数；k中的各个元素，用于记录该列对应的项是否存在于该行对应的偏微分方程中；
[0029]
(2)若i≤p，开始循环；
[0030]
(2.1)将第i个左端项和所有右端项组成矩阵，即a
′i＝a
′
(：，(i p+1：n
′
))；
[0031]
(2.2)若j≤n
′‑
p，开始循环；
[0032]
(2.2.1)删掉矩阵a
′i的第j+1列，做奇异值分解，计小奇异值的个数为p
i，j
；
[0033]
(2.2.2)如果p
i，j
＝0，令k(i，j)＝1；否则不实施操作；令j＝j+1，返回步骤(2.2)继续循环；
[0034]
(2.3)令i＝i+1，j＝1；返回步骤(2)继续循环；
[0035]
(3)使用最小二乘法计算偏微分方程组所有项的系数。
[0036]
进一步，在步骤1中，为了防止各项的数量级相差过大对计算结果造成影响，对数据矩阵a的每一列实施l2范数的归一化；具体操作方法为将各列的所有分量平方求和后开根号，即为各列的l2范数，之后令矩阵每个元素除以该列的l2范数即可；归一化后的数据矩阵a每列的l2范数都等于1，保证各列的数量级一致。
[0037]
进一步，数据驱动识别偏微分方程的序列奇异值过滤方法seq-svf的具体实现过程如下：
[0038]
(1)对各质点的原始质点位移数据进行高斯滤波处理，并计算所有导出项的值，将其组成数据矩阵a，对各列实施l2范数的归一化；对数据矩阵a实施奇异值分解，小奇异值的个数计为p；
[0039]
(2)初始化变量，令i＝1，k为空向量；
[0040]
(2.1)若i≤矩阵a的列数n，开始循环；
[0041]
(2.1.1)删掉矩阵a的第i列，做奇异值分解，计小奇异值的个数为pi；
[0042]
(2.1.2)如果pi＝p-1，令k＝(k i)；否则不实施操作；令i＝i+1，返回步骤(2)继续循环；
[0043]
(2.2)a
′
＝a(：，k)，即为过滤掉无关项后的数据矩阵；
[0044]
(3)对a
′
实施强秩揭示正交三角分解方法，得到置换矩阵n，其后n
′‑
q列中单位元素所处的行数即为左端项在a
′
中的列数，将a
′
中的左端项移至最左侧；
[0045]
(4)初始化变量，令i＝1，j＝1，k
p
×
(n
′‑
p)
＝0；其中n
′
为a
′
的列数；k中的各个元素用于记录该列对应的项是否存在于该行对应的偏微分方程中；
[0046]
(4.1)若i≤p，开始循环；
[0047]
(4.1.1)将第i个左端项和所有右端项组成矩阵，即a
′i＝a
′
(：，(i p+1：n
′
))；
[0048]
(4.1.2)若j≤n
′‑
p，开始循环；
[0049]
(4.1.2.1)删掉矩阵a
′i的第j+1列，做奇异值分解，计小奇异值的个数为p
i，j
；
[0050]
(4.1.2.2)如果p
i，j
＝0，令k(i，j)＝1；否则不实施操作；令j＝j+1，返回步骤(3.1.2)继续循环；
[0051]
(4.1.3)根据k中记录的各偏微分方程对应项的信息，令(＝i+1，j＝1，返回步骤(3.1)继续循环；
[0052]
(4.2)使用最小二乘法计算偏微分方程组所有项的系数。
[0053]
本发明的有益效果：
[0054]
采用本发明提供的技术方案与现有技术相比，具有如下显著效果：
[0055]
(1)本发明提供的数据驱动识别偏微分方程的序列奇异值过滤方法，为微纳米力学系统的数据驱动建模提供了一套简便的算法流程。将所有可能存在的质点位移对空间坐标的高阶导数项和非线性项补充至数据集中，有效克服了传统建模方法只能识别线性系统的不足。通过采用奇异值分解、强秩揭示正交三角分解等先进的数值计算方法，保证了计算效率和结果的稳定性。与现有的基于优化范式的系统识别方法相比，能够避免结果陷入局部最优的问题。
[0056]
(2)本发明提供的数据驱动识别偏微分方程的序列奇异值过滤方法，克服了现有方法过于依赖预先人为地对系统的控制方程个数和方程形式加以限制的不足，将有监督的学习范式转变为无监督的算法流程，将数据驱动识别方法的应用场景从动力系统拓展至由
一般非线性偏微分方程组所控制的系统。使用奇异值分解技术自动识别隐藏在数据矩阵中的本征低秩结构，可以避免算法对关于控制方程个数的先验知识的依赖。使用强秩揭示正交三角分解方法可以从函数库中自动提取最具代表性的一组基底作为控制方程的右端项，避免了目前的有监督学习方法对于方程结构的先验知识的依赖。
[0057]
(3)本发明提供的数据驱动识别偏微分方程的序列奇异值过滤方法，为显式提取微纳米力学控制方程提供了一套自动化的算法流程，克服了现有的深度学习方法可解释性差的不足。创新性地设计了一套奇异值过滤流程，通过矩阵秩的变化情况判断各项与其余项的线性相关关系，能够自动排除与系统无关的干扰变量的影响，有效地确保了模型的稀疏结构，能够确保本发明的鲁棒性和泛化能力，并且为以该结果为基础的后续工作，如分析微纳米力学方程的物理意义和进一步的优化与数值模拟奠定了坚实基础。
附图说明
[0058]
图1为本发明的一种数据驱动识别偏微分方程的序列奇异值过滤方法(seq-svf)的操作流程图；
[0059]
图2为本发明的实施例1使用seq-svf识别三维梁弹性变形控制方程的流程示意图和计算结果；
[0060]
图3为本发明的实施例2使用seq-svf识别纳米铜位移控制方程的流程示意图和计算结果；
[0061]
图4为本发明的实施例3热力强耦合系统的(a)位移u和(b)温度t关于时间t和坐标x的变化图；
[0062]
图5为本发明的实施例4使用seq-svf识别lorenz系统控制方程的流程示意图和去噪前后计算结果的对比；
[0063]
图6为本发明的实施例5二维不可压缩流场在(a)t＝0，(b)t＝10和(c)t＝20时位移u的原始数据云图。
具体实施方式
[0064]
下面结合附图和实施例对本发明的性能做出进一步详细说明。以下实施例用于说明本发明，但不能用来限制本发明的适用范围。
[0065]
一种数据驱动识别偏微分方程的序列奇异值过滤方法，步骤如下：
[0066]
首先，给出序列奇异值过滤方法识别偏微分方程组的基本格式；对于一个在n维的时空坐标(t x
1 x2ꢀ…ꢀ
x
n-1
)下，由m个变量(u
1 u2ꢀ…ꢀ
um)组成的微纳米力学系统，设该系统由一组偏微分方程控制：
[0067][0068]
将偏微分方程中所有可能存在的项罗列出来，如式(3)所示：
[0069][0070]
其中，为待定系数；设采样点的数量为m，式(3)应在所有采样点上处处成立，即：
[0071][0072]
其中，x1～xn代表所有假设存在的项，包括高阶导数和非线性项，因此根据在各个采样点处测得的变量值，通过数值计算得到数据矩阵a所有元素的大小；序列奇异值过滤方法的目标是求解式(4)中的系数矩阵λ：识别偏微分方程的个数p，并计算待定系数的值，同时使每一列系数向量尽可能稀疏；
[0073]
接着，使用奇异值分解方法来确定偏微分方程的个数，并过滤掉函数库中冗余的项；式(4)表明，由于偏微分方程的存在，使得数据矩阵a分布在低维的子空间上，而不是整个n维空间：对于p个偏微分方程控制的系统，偏微分方程会对数据集产生p维的约束，使得数据矩阵a分布在n-p维的子空间上；使用奇异值分解方法提取出p的值；奇异值分解是将m
×
n的数据矩阵a分解为正交矩阵u，v和对角矩阵∑的乘积：
[0074]
a＝u∑v
t
ꢀꢀ
(5)
[0075]
由于采样点个数m远大于函数库的规模n，将上式展开有：
[0076][0077]
其中，σ1，σ2，...，σn称为a的奇异值，满足σ1≥σ2≥
…
≥σn≥0；当a不是满秩矩阵时，即rank(a)《n，此时a的奇异值只有前rank(a)个非零，即：σ
rank(a)+1
＝σ
rank(a)+2
＝
…
＝σn＝0；此时得到约化的奇异值分解：
[0078][0079]
测量数据会受到噪声的影响，并且数值微分也会带来难以避免的数值误差，导致实际的数据矩阵a并不是完全亏秩的，是一个近似奇异的病态矩阵；此时后n-rank(a)个奇异值的大小不严格等于0，是与前面奇异值相比的小量，它们的大小关系取决于噪声的量级大小；当噪声控制在一定范围时，a的后n-rank(a)个奇异值会远小于前面的奇异值，即σ
rank(a)+1
，σ
rank(a)+2
，...，σn＜＜σ
rank(a)
；基于上述原理，通过判断矩阵的小奇异值的数量来确定偏微分方程的个数，n-rank(a)即为p的大小；此外，为了防止各项的数量级相差过大对计算结果造成影响，在预处理阶段对数据矩阵a的每一列实施l2范数的归一化；
[0080]
为了提高后续流程的计算效率和精度，通过设计迭代流程，反复对数据矩阵进行缩减并实施奇异值分解，来筛选所有与偏微分方程相关的项；具体过程如下：依次删掉数据矩阵a中的每一列进行奇异值分解，删掉某列后如果仍然识别出了p个偏微分方程，说明这一列对应的项不包含在偏微分方程组里面；如果识别出了p-1个偏微分方程，说明这项存在
于偏微分方程组里面；计算n次后可筛选出n
′
个与偏微分方程有关的项，将这些项对应的列保留下来得到重构的数据矩阵a
′
，以供后续步骤使用；
[0081]
之后，使用强秩揭示正交三角分解方法srrqr提取每个偏微分方程的左端项；给出srrqr的计算格式，对于任一数据矩阵a，srrqr将其分解成式(8)形式：
[0082][0083]
其中，q为正交矩阵，r为上三角矩阵，q1，q2，b，c，d分别为q，r的子块；并满足：
[0084][0085]
σ(d)≤σ
q+1
(a)
·
q1(q，n)
ꢀꢀꢀꢀꢀꢀ
(10)
[0086]
|(b-1
c)
1～q，1～n-q
|≤q2(q，n)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0087]
其中σ(b)是矩阵b的奇异值，q是需要预先设定的值，q1(q，n)，q2(q，n)具有多项式级的上界；式(11)指矩阵b-1
c每个元素的大小都被限制不大于q2(q，n)；式(8)中的π是一个置换矩阵；将前面通过奇异值分解确定的矩阵的秩rank(a)代入为q的值，此时q＝rank(a)＝n
′‑
p，所以有n
′‑
q＝p，对缩减后的数据矩阵a
′
使用srrqr方法；根据计算得到的n来提取a
′
π的后n
′‑
q列在原矩阵a
′
中的位置信息，即找出了n
′‑
q个即p个左端项在a
′
中的位置；
[0088]
最后，使用奇异值分解方法和最小二乘法确定稀疏形式的偏微分方程组；确定偏微分方程稀疏形式的方法与用奇异值分解方法过滤函数库中冗余项的方法相同：将p个左端项分别逐个与n
′‑
p个候选右端项组成矩阵，此时的矩阵受一个偏微分方程控制，即这个左端项对应的偏微分方程；依次删掉每个候选右端项进行奇异值分解，删掉某项后如果仍然识别出了1个偏微分方程，说明这一列对应的项不包含在此偏微分方程里面；如果识别出了0个偏微分方程，说明这项存在于此偏微分方程中；进行p(n
′‑
p)次奇异值分解后，所有偏微分方程的左端项和右端项都已确定；使用最小二乘法计算对应系数，结果如式(12)所示：
[0089][0090][0091]
根据上述理论推导得出的奇异值分解和强秩揭示正交三角分解方法的基本格式，嵌入到所述的迭代流程，即实现序列奇异值过滤方法，具体步骤如下，
[0092]
步骤1：根据具体实例的需要设置候选函数库，包括设定导函数的最高阶次，乘积的最高幂次，以及其它非线性函数的形式，根据初始数据和候选函数库计算数据矩阵，对数据矩阵的每列实施正则化，使用奇异值分解方法识别偏微分方程的个数；
[0093]
步骤2：根据删掉某一项后识别出的偏微分方程个数是否减少判断该项是否存在于偏微分方程组中，并更新数据矩阵；
[0094]
步骤3：使用强秩揭示正交三角分解方法识别所有左端项在更新后数据矩阵中的位置，并进一步更新数据矩阵；
[0095]
步骤4：依次将每个左端项与所有右端项组成矩阵，并依次删除每个右端项，根据
偏微分方程的个数是0还是1确定该项是否存在于当前左端项对应的偏微分方程中；确定各个偏微分方程中的项之后使用最小二乘法计算各项系数。
[0096]
为了使本发明的目的、技术方案和具体实施效果展示更加清晰明了，下面通过五个具体实施例结合附图2～6对本发明提出的seq-svf方法的准确性和有效性作进一步的详细说明。首先，通过一个三维的经典弹性力学算例验证本发明所提出方法对识别一般微分方程组的准确性和有效性。其次，通过一个纳米铜变形算例验证本发明所提出方法对识别微纳米力学高阶微分方程组的适用性。然后，通过采用热力耦合算例说明本发明提出的seq-svf方法识别复杂耦合问题的能力。最后，通过采用非线性动力系统和非线性微分方程控制的流场算例说明本发明提出方法识别非线性系统的能力。
[0097]
实施例1：三维弹性梁复合加载下的受力变形(附图2)
[0098]
本实施例为三维弹性体受力变形算例，通过在仿真软件中对弹性梁施加弯扭剪拉多种载荷，使其产生弯曲、扭转、剪切和轴向拉伸多种变形模式，如图2所示，并收集若干离散点的初始位置(x，y，z)和变形后的位置(u，v，w)作为原始数据集。弹性范围内固体的位移受一组微分方程控制：
[0099][0100]
其中为拉普拉斯算子，λ，μ称为拉梅常数，在本实例中将其代入式(13)并展开，为一组包含三个方程的二阶微分方程组：
[0101][0102][0103]
seq-svf假设该系统包含所有0、1、2阶的导数作为函数库，即：
[0104][0105]
在添加噪声的情况下，seq-svf成功从30个候选函数中提取出了15个正确的项，同时将其余无关的项过滤掉，计算结果如(16)所示，各项系数的相对误差在1％以内：
[0106][0107]
此实施例有效地说明了本发明提出的序列奇异值过滤方法对识别一般偏微分方程组的有效性和准确性。
[0108]
实施例2：纳米铜梁复合加载下的受力变形(附图3)
[0109]
本实施例为纳米铜梁受力变形算例，通过分子动力学模拟软件对纳米铜梁施加弯扭剪拉多种载荷，使其产生弯曲、扭转、剪切和轴向拉伸多种变形模式，如图3所示，并收集
若干铜原子的初始位置(x，y，z)和变形后的位置(u，v，w)作为原始数据集。
[0110]
seq-svf假设该系统包含位移对空间坐标的所有2阶和4阶导数作为函数库，即：
[0111][0112]
seq-svf成功从60个候选函数中提取出了30个项，同时将其余无关的项过滤掉，计算结果如(18)所示：
[0113]uxx
+0.3770u
yy
+0.3770u
zz
+1.0554v
xy
+1.0554w
xz
+2.5390u
xxxx
+0.8263u
yyyy
+0.8263u
zzzz
+2.0818v
xxxy
+2.0818w
xxxz
＝0
[0114]vyy
+0.3770v
xx
+0.3770v
zz
+1.0554u
xy
+1.0554w
yz
+2.5390v
yyyy
+0.8263v
xxxx
+0.8263v
zzzz
+2.0818u
xyyy
+2.0818w
yyyz
＝0
[0115]wzz
+0.3770w
xx
+0.3770w
yy
+1.0554u
xz
+1.0554v
yz
+2.5390w
zzzz
+0.8263w
xxxx
+0.8263w
yyyy
+2.0818u
xzzz
+2.0818v
yzzz
＝0
[0116]
(18)
[0117]
此实施例有效地说明了本发明提出的序列奇异值过滤方法对识别微纳米力学高阶偏微分方程组的适用性。
[0118]
实施例3：热力强耦合系统的数据驱动控制方程识别(附图4)
[0119]
本实施例为物体的变形与温度互相影响的耦合算例。seq-svf作为一种无监督的偏微分方程识别方法，具有从多个变量中提取方程的能力，在识别复杂耦合微分方程问题方面具有显著优势。下面通过一个热力强耦合算例来说明序列奇异值过滤方法在处理多场问题方面的优势和必要性。如图4所示，考虑一维情况下物体位移(u)和温度的传播(t)。物体在变形过程中，弹性势能不断积累和释放，其中一部分转化为热能，会导致温度的变化；温度的变化会使得应力重新分布，从而反过来影响位移。热力强耦合的控制方程如式(19)所示，其中各物理量经过了无量纲化处理：
[0120][0121][0122]
seq-svf将所有0、1、2阶的导数作为函数库，在添加噪声的情况下，识别出了正确形式的控制方程，并且各项系数的相对误差在0.5％以内，识别结果如(20)所示：
[0123][0124][0125]
此实例强调了本发明提出的序列奇异值过滤方法对识别耦合多变量系统的有效性。
[0126]
实施例4：lorenz系统的数据驱动控制方程识别(附图5)
[0127]
本实施例为非线性动力系统的验证算例。lorenz系统是一种典型的非线性系统，在相空间内的轨迹如图5所示，其混沌吸引子具有对数值微扰的极端敏感性，因此成为非线性系统识别的一个经典算例。系统的控制方程如式(21)所示：
[0128][0129]
本发明提出的序列奇异值过滤方法将非线性函数库设为乘积形式，并将最高阶数假设至2阶，即数据矩阵为seq-svf识别控制方程的流程如图5所示，可以看出导数数据的计算对噪声十分敏感，未经去噪处理的数据会导致方法得不到正确的结果。因此，此实例证明了对原始数据进行高斯滤波处理是必要并且有效的。计算得到的直接结果如(22)所示：
[0130][0131]
将各变量的一阶导数项移至左侧，以便对比计算结果与真实系数：
[0132][0133]
其中多出来的一项是由于对矩阵做线性变换所致，其系数可以小到忽略不计。该结果进一步证明了seq-svf方法的准确性与鲁棒性。
[0134]
实施例5：由navier-stokes方程控制的二维不可压缩流场(附图6)
[0135]
最后一个实施例将考察序列奇异值过滤方法对识别一般非线性偏微分方程的有效性。不可压缩流体的行为可以由navier-stokes方程来描述。二维流场的n-s方程如式(24)所示：(雷诺数为100)
[0136]ut
+uu
x
+vuy+p
x-0.01u
xx-0.01u
yy
＝0
[0137]vt
+uv
x
+vvy+p
y-0.01v
xx-0.01v
yy
＝0
ꢀꢀ
(24)
[0138]
u，v为位移，p为压强。当seq-svf没有假设非线性项在函数库中时，数据矩阵如下所示：
[0139]
a＝[u，v，p，ut，ux，uy，v
t
，v
x
，vy，p
t
，p
x
，py]
ꢀꢀꢀꢀ
(25)
[0140]
由于数据矩阵中没有非线性项，所以没有识别出n-s方程，而是识别出了不可压缩条件：
[0141]ux
+1.0064vy＝0
ꢀꢀꢀꢀ
(26)
[0142]
其理论方程中两项的系数为1。接着，在函数库中添加非线性项：
[0143]
a＝[u，v，p，u
t
，uut
x
，vuy，p
x
，u
xx
，u
yy
，v
t
，uv
x
，vvy，py，v
xx
，v
yy
]
ꢀꢀꢀꢀ
(27)
[0144]
seq-svf识别出了正确形式的控制方程，结果如(28)所示：
[0145]
0.9928u
t
+0.9916uu
x
+1.0042vuy+p
x-0.0104u
xx-0.0102u
yy
＝0
[0146]
0.9977v
t
+0.9962uv
x
+vvy+1.0009p
y-0.0100v
xx-0.0099v
yy
＝0
ꢀꢀ
(28)
[0147]
本实施例进一步说明了本发明提出的序列奇异值过滤方法的通用性和准确性，可
以通过实际应用场景的需要合理设计非线性函数库，以求通过序列奇异值过滤方法得到预期的结果。
[0148]
综上所述，我们首先通过实施三维弹性体位移驱动的偏微分方程的系统识别验证了本发明所提出的序列奇异值过滤方法seq-svf在一般微分方程组的系统识别中的重要性和准确性，说明了seq-svf从构建数据集到过滤有关项并识别左端项，最后计算方程系数的整体流程。随后，将seq-svf应用于一个纳米铜受载变形算例，从铜原子的位移数据中识别出了一组四阶偏微分控制方程组，验证了seq-svf对识别微纳米力学高阶控制微分方程组的适用性。之后，考察固体的热力耦合算例，数据驱动地提取出了物体位移与温度变化的显式关系，有效说明了seq-svf面对多场耦合系统识别问题的正确性和必要性。接着，通过对非线性lorenz系统的识别，有效说明了seq-svf在识别非线性动力系统问题中的准确性，并通过对比高斯滤波前后的结果说明了本方法对原始数据实施去噪处理的必要性。最后，一个由非线性微分方程组控制的二维流场系统的识别问题被考察，说明了seq-svf在识别一般形式的多变量非线性偏微分方程组的通用性和有效性。因此，本发明所提出的序列奇异值过滤方法seq-svf是一种极具发展前景的偏微分方程数据驱动识别方法。
[0149]
本发明的实施例是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显而易见的。选择和描述实施例是为了更好的说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

技术特征：
1.一种数据驱动识别偏微分方程的序列奇异值过滤方法，其特征在于，步骤如下：首先，给出序列奇异值过滤方法识别偏微分方程组的基本格式；对于一个在n维的时空坐标(t x
1 x2ꢀ…ꢀ
x
n-1
)下，由m个变量(u
1 u2ꢀ…ꢀ
u
m
)组成的微纳米力学系统，设该系统由一组偏微分方程控制：将偏微分方程中所有可能存在的项罗列出来，如式(2)所示：其中，为待定系数；设采样点的数量为m，式(2)应在所有采样点上处处成立，即：其中，x1～x
n
代表所有假设存在的项，包括高阶导数和非线性项，因此根据在各个采样点处测得的变量值，通过数值计算得到数据矩阵a所有元素的大小；序列奇异值过滤方法的目标是求解式(3)中的系数矩阵λ：识别偏微分方程的个数p，并计算待定系数的值，同时使每一列系数向量尽可能稀疏；接着，使用奇异值分解方法来确定偏微分方程的个数，并过滤掉函数库中冗余的项；式(3)表明，由于偏微分方程的存在，使得数据矩阵a分布在低维的子空间上，而不是整个n维空间：对于p个偏微分方程控制的系统，偏微分方程会对数据集产生p维的约束，使得数据矩阵a分布在n-p维的子空间上；使用奇异值分解方法提取出p的值；奇异值分解是将m
×
n的数据矩阵a分解为正交矩阵u,v和对角矩阵∑的乘积：a＝u∑v
t
ꢀꢀꢀꢀꢀꢀꢀ
(4)由于采样点个数m远大于函数库的规模n，将上式展开有：其中，σ1，σ2，...，σ
n
称为a的奇异值，满足σ1≥σ2≥
…
≥σ
n
≥0；当a不是满秩矩阵时，即rank(a)＜n，此时a的奇异值只有前rank(a)个非零，即：σ
rank(a)+1
＝σ
rank(a)+2
＝
…
＝σ
n
＝0；此时得到约化的奇异值分解：测量数据会受到噪声的影响，并且数值微分也会带来难以避免的数值误差，导致实际
的数据矩阵a并不是完全亏秩的，是一个近似奇异的病态矩阵；此时后n-rank(a)个奇异值的大小不严格等于0，是与前面奇异值相比的小量，它们的大小关系取决于噪声的量级大小；当噪声控制在一定范围时，a的后n-rank(a)个奇异值会远小于前面的奇异值，即σ
rank(a)+1
，σ
rank(a)+2
，...，σ
n
＜＜σ
rank(a)
；基于上述原理，通过判断矩阵的小奇异值的数量来确定偏微分方程的个数，n-rank(a)即为p的大小；此外，为了防止各项的数量级相差过大对计算结果造成影响，在预处理阶段对数据矩阵a的每一列实施l2范数的归一化；为了提高后续流程的计算效率和精度，通过设计迭代流程，反复对数据矩阵进行缩减并实施奇异值分解，来筛选所有与偏微分方程相关的项；具体过程如下：依次删掉数据矩阵a中的每一列进行奇异值分解，删掉某列后如果仍然识别出了p个偏微分方程，说明这一列对应的项不包含在偏微分方程组里面；如果识别出了p-1个偏微分方程，说明这项存在于偏微分方程组里面；计算n次后可筛选出n
′
个与偏微分方程有关的项，将这些项对应的列保留下来得到重构的数据矩阵a
′
，以供后续步骤使用；之后，使用强秩揭示正交三角分解方法srrqr提取每个偏微分方程的左端项；给出srrqr的计算格式，对于任一数据矩阵a，srrqr将其分解成式(7)形式：其中，q为正交矩阵，r为上三角矩阵，q1，q2，b，c，d分别为q，r的子块；并满足：σ(d)≤σ
q+1
(a)
·
q1(q，n)
ꢀꢀꢀꢀ
(9)|(b-1
c)
1～q，1～n-q
|≤q2(q，n)
ꢀꢀꢀꢀ
(10)其中σ(b)是矩阵b的奇异值，q是需要预先设定的值，q1(q，n)，q2(q，n)具有多项式级的上界；式(10)指矩阵b-1
c每个元素的大小都被限制不大于q2(q，n)；式(7)中的π是一个置换矩阵；将前面通过奇异值分解确定的矩阵的秩rank(a)代入为q的值，此时q＝rank(a)＝n
′‑
p，所以有n
′‑
q＝p，对缩减后的数据矩阵a
′
使用srrqr方法；根据计算得到的π来提取a
′
π的后n
′‑
q列在原矩阵a
′
中的位置信息，即找出了n
′‑
q个即p个左端项在a
′
中的位置；最后，使用奇异值分解方法和最小二乘法确定稀疏形式的偏微分方程组；确定偏微分方程稀疏形式的方法与用奇异值分解方法过滤函数库中冗余项的方法相同：将p个左端项分别逐个与n
′‑
p个候选右端项组成矩阵，此时的矩阵受一个偏微分方程控制，即这个左端项对应的偏微分方程；依次删掉每个候选右端项进行奇异值分解，删掉某项后如果仍然识别出了1个偏微分方程，说明这一列对应的项不包含在此偏微分方程里面；如果识别出了0个偏微分方程，说明这项存在于此偏微分方程中；进行p(n
′‑
p)次奇异值分解后，所有偏微分方程的左端项和右端项都已确定；使用最小二乘法计算对应系数，结果如式(11)所示：
根据上述理论推导得出的奇异值分解和强秩揭示正交三角分解方法的基本格式，嵌入到所述的迭代流程，即实现序列奇异值过滤方法，具体步骤如下，步骤1：根据具体实例的需要设置候选函数库，包括设定导函数的最高阶次，乘积的最高幂次，以及其它非线性函数的形式，根据初始数据和候选函数库计算数据矩阵，对数据矩阵的每列实施正则化，使用奇异值分解方法识别偏微分方程的个数；步骤2：根据删掉某一项后识别出的偏微分方程个数是否减少判断该项是否存在于偏微分方程组中，并更新数据矩阵；步骤3：使用强秩揭示正交三角分解方法识别所有左端项在更新后数据矩阵中的位置，并进一步更新数据矩阵；步骤4：依次将每个左端项与所有右端项组成矩阵，并依次删除每个右端项，根据偏微分方程的个数是0还是1确定该项是否存在于当前左端项对应的偏微分方程中；确定各个偏微分方程中的项之后使用最小二乘法计算各项系数。2.根据权利要求1所述的数据驱动识别偏微分方程的序列奇异值过滤方法，其特征在于，步骤如下：步骤1、预处理编码器：设置候选函数库，根据质点位移数据和候选函数库计算数据矩阵，对数据矩阵的每列实施正则化，使用奇异值分解方法识别偏微分方程的个数；首先，给出本方法识别微纳米力学系统的基本格式；设该微纳米力学系统由一组偏微分方程控制：其中，x1～x
n
代表所有假设存在的项，包括质点位置对于空间坐标的二阶以及更高阶导数项和它们通过相乘得到的非线性项；根据在各个采样点处测得的质点位移，通过数值计算得到数据矩阵a所有元素的大小；接着，为了防止各项的数量级相差过大对计算结果造成影响，对数据矩阵a的每一列实施l2范数的归一化；最后，对数据矩阵a实施奇异值分解，通过判断矩阵的小奇异值的数量来确定偏微分方程的个数；步骤2、函数过滤器：删除候选函数库中冗余的项；通过删除数据矩阵a中的各列，根据奇异值分解的结果判断该列是否存在于偏微分方程中；具体实现过程如下：(1)初始化变量，令i＝1，k为空向量；p为数据矩阵a中小奇异值的个数；(2)若i≤数据矩阵a的列数n，开始循环；(2.1)删掉数据矩阵a的第i列，做奇异值分解，计小奇异值的个数为p
i
；(2.2)如果p
i
＝p-1，令k＝(k i)；否则不实施操作；令i＝i+1，返回步骤(2)继续循环；(3)a
′
＝a(：，k)，即为过滤掉无关项后的数据矩阵，以供后续步骤使用；步骤3、线性解耦器：使用强秩揭示正交三角分解方法识别所有左端项在通过步骤2得到的过滤后数据矩阵(a
′
)中的位置，并更新数据矩阵；步骤4、稀疏解码器：确定各个偏微分方程中的项并计算各项系数；
对每个偏微分方程各实施一次步骤2，即可在候选项中确定每个偏微分方程的稀疏模型；之后使用最小二乘法计算所有系数；其具体实施过程如下：(1)初始化变量，令i＝1，j＝1，k
p
×
(n
′‑
p)
＝0；其中n
′
为过滤掉无关项后的数据矩阵a
′
的列数；k中的各个元素，用于记录该列对应的项是否存在于该行对应的偏微分方程中；(2)若i≤p，开始循环；(2.1)将第i个左端项和所有右端项组成矩阵，即a
′
i
＝a
′
(：，(i p+1：n
′
))；(2.2)若j≤n
′‑
p，开始循环；(2.2.1)删掉矩阵a
′
i
的第j+1列，做奇异值分解，计小奇异值的个数为p
i，j
；(2.2.2)如果p
i，j
＝0，令k(i，j)＝1；否则不实施操作；令j＝j+1，返回步骤(2.2)继续循环；(2.3)令i＝i+1，j＝1；返回步骤(2)继续循环；(3)使用最小二乘法计算偏微分方程组所有项的系数。3.根据权利要求1所述的数据驱动识别偏微分方程的序列奇异值过滤方法，其特征在于，在步骤1中，为了防止各项的数量级相差过大对计算结果造成影响，对数据矩阵a的每一列实施l2范数的归一化；具体操作方法为将各列的所有分量平方求和后开根号，即为各列的l2范数，之后令矩阵每个元素除以该列的l2范数即可；归一化后的数据矩阵a每列的l2范数都等于1，保证各列的数量级一致。4.根据权利要求1所述的数据驱动识别偏微分方程的序列奇异值过滤方法，其特征在于，数据驱动识别偏微分方程的序列奇异值过滤方法seq-svf的具体实现过程如下：(1)对各质点的原始质点位移数据进行高斯滤波处理，并计算所有导出项的值，将其组成数据矩阵a，对各列实施l2范数的归一化；对数据矩阵a实施奇异值分解，小奇异值的个数计为p；(2)初始化变量，令i＝1，k为空向量；(2.1)若i≤矩阵a的列数n，开始循环；(2.1.1)删掉矩阵a的第i列，做奇异值分解，计小奇异值的个数为p
i
；(2.1.2)如果p
i
＝p-1，令k＝(k i)；否则不实施操作；令i＝i+1，返回步骤(2)继续循环；(2.2)a
′
＝a(：，k)，即为过滤掉无关项后的数据矩阵；(3)对a
′
实施强秩揭示正交三角分解方法，得到置换矩阵n，其后n
′‑
q列中单位元素所处的行数即为左端项在a
′
中的列数，将a
′
中的左端项移至最左侧；(4)初始化变量，令i＝1，j＝1，k
p
×
(n
′‑
p)
＝0；其中n
′
为a
′
的列数；k中的各个元素用于记录该列对应的项是否存在于该行对应的偏微分方程中；(4.1)若i≤p，开始循环；(4.1.1)将第i个左端项和所有右端项组成矩阵，即a
′
i
＝a
′
(：，(i p+1：n
′
))；(4.1.2)若j≤n
′‑
p，开始循环；(4.1.2.1)删掉矩阵a
′
i
的第j+1列，做奇异值分解，计小奇异值的个数为p
i，j
；(4.1.2.2)如果p
i，j
＝0，令k(i，j)＝1；否则不实施操作；令j＝j+1，返回步骤(3.1.2)继续循环；
(4.1.3)根据k中记录的各偏微分方程对应项的信息，令i＝i+1，j＝1，返回步骤(3.1)继续循环；(4.2)使用最小二乘法计算偏微分方程组所有项的系数。

技术总结
本发明属于计算力学领域，提供一种数据驱动识别偏微分方程的序列奇异值过滤方法，该方法使用奇异值分解方法和强秩揭示正交三角分解方法提取隐藏在数据集中的本征线性结构，提高了从观测数据提取非线性偏微分方程算法的计算效率、准确度和鲁棒性。本发明通过设计奇异值过滤的算法，对数据矩阵进行反复缩减并计算矩阵的奇异值，可以有效识别控制方程的个数并提取有关项。采用强秩揭示正交三角分解方法为每个控制方程设置左端项，并通过奇异值分解和最小二乘法有效识别控制方程组的最简形式。本发明所提出的方法作为一种新的控制方程数据驱动识别方法，借助高效的矩阵分解技术，快速准确，并且可以通过编写简单循环程序实现，降低了数值实施复杂度。降低了数值实施复杂度。降低了数值实施复杂度。


技术研发人员：郑勇刚 吴哲同 刘振海 叶宏飞
受保护的技术使用者：大连理工大学
技术研发日：2023.04.21
技术公布日：2023/8/1