李聪，时宏伟

（四川大学计算机学院，成都610065）

0 引言

非线性动力系统辨识在现代控制领域内一直是致力于解决的难题之一。它基于这样的假设，系统能够由一组确定的常微分方程所确定，然而由于实际的建模过程中，这样的方程往往难以建模或者建立的模型非常复杂，于是就想通过机器学习的方法来解决该问题。文献[1]通过从数据中建模出动态化线性模型，对比自适应控制模型，该模型性能优异。文献[2]利用电压数据预测电路故障。文献[3]使用高斯过程回归建立预测模型，文献[4]从数据当中建立导弹退化故障模型，从而预测导弹的故障。这些模型都是从历史的数据中建模出系统动力模型，然后用该模型预测下一时刻的系统变化。即从系统采集到的数据发掘系统的动态变化规律。对于一般的动力系统，我们往往能获取系统中各个时刻的状态参数，想通过这些数据来揭示出系统模型。现在常用的方法有基于动态模式分解（Dynamic mode Decomposition，DMD）[4]，试图在局部时间关系建立线性关系，通过奇异值分解的方法求解伪逆的方法来求解动力系统中的参数。也有通过加入先验知识通过搜索的方法揭示出动力学模型，事先在一个备用的高次多项式的字典中搜索[5]，直到找到满足条件的动力学系统，但是文中所示实验仅仅以Lorenz 三维系统，方法新奇但不适合一般的较复杂的动力系统。也有基于神经网络结合Jacobian 正则的方法来建模动力系统[7]，能够很好的长时建模和预测动力系统，但是模型本身没有建模长时动力系统，随着迭代增长动力系统误差必然增大。文献[8]建模了多个时间步动力系统，但是该方法在Lorenz 系统上表现不佳，随着迭代时间步的增长，一点点误差导致系统状态参数误差增长不可控制。基于深度Koopman 算子建模非线性动力系统，然而系统仅建模了单步动力系统[11]，所以本文考虑在单步模型上考虑到了减少模型误差，所以引入了总均方误差（Total Least Square，TLS）[12]的方法求解单步动力模型。文献[15]首次使用TLS 在非线性动力系统辨识中，文中模型首次将输入数据和输出数据都考虑为含噪数据，且噪声性质相同，即输入数据和输出数据是对称的，此外在模型上扩展了均方误差，文中使用广义线性模型来建模数据。由于扩展了模型，使得广义线性模型具有更一般的表达能力。

非线性动力系统的数据也可以看作是时间序列数据，时间序列数据往往具有周期性，但往往不具有标准的周期，在一定的时间内具有某种非规范的周期性。为解决这个问题，文献[13]首次提出将时间序列数据分解为多个内蕴函数（Intrinsic Mode Function，IMF），非线性动力系统所表示的数据也可以表征为某种形式的内蕴函数的加权叠加，也有一些在这方面的尝试[20]。还有一些模型考虑从过去的多个时间步中提取信息，建模非线性动力学，借助长短时记忆网络（Long Short Memory Network，LSTM）[23]的记忆能力来表达系统的非线性特征。例如线性循环自编码网络[24]考虑了前面多个时刻的状态参数，将多个时刻的数据作为系统输入，再利用自编码机将数据映射到某一线性空间，相邻时间步的系统状态数据也建模为线性关系，最后在几个数据集上验证，获得不错的结果。

与之类似的方法在文献[25]提出，文章试图用自编码机从数据中学习出Koopman 不变子空间，文中方法在编码后的数据中建模时序依赖关系，这仍能建模出数据的局部线性关系。基于深度学习和Koopman 算子的结合的方法还有很多[26-27]，大致和文献[14]类似，都是试图学习出不变子空间，在线性空间上做状态参数演化，这种方法均是仅考虑状态参数在局部空间上的演化，没有考虑模型随时间演化误差较大的问题。

为了解决上述问题，本文在深度动态模式分解的算法框架之下，引入了Jacobian 正则化方法增加神经网络的鲁棒性；其次在线性模型的部分广义化了原模型，引入了Total Least Square 作为求解广义化模型的方法；最后在几个常用的测试数据集上进行实验验证，得到非常不错的效果。算法模型

本文使用模型包含三部分。第一，本文中的基础模型采用深度动态模式分解[14]，该模型在单步预测中取得较好效果，但在长期（时间步在1000～4000）预测上表现不佳，因此需要引入额外的约束。第二，在基础模型之上，通过在损失函数中加入Jacobian 正则，使得模型能够在保持单步预测不失真的条件下，又能让模型在循环预测中抑制误差的传播，相比较于基础模型，加入正则后的模型能够在理论上保证误差减小，得在后续的预测中更加准确。第三，总体最小二乘法（Total Least Square，TLS）方法减小基础模型中线性模型的误差，TLS 也是减小误差的方法，只不过在输入上也考虑到误差，且TLS 是广义的最小二乘法，如果输入数据中没有误差，TLS 会退化到普通最小二乘法（Ordinary Least Square，OLS）。

1 建模方法

1.1 非线性动力模型

本文中考虑的非线性动力系统的形式为：

其中系统状态向量xt∈Rn为系统t 时刻系统状态参数，xt+1∈Rn为下一时刻的系统状态参数，函数f ∈Rn表示当前时刻系统参数到下一时刻系统状态参数的映射，P 为其参数。在模型训练完成之后，可以使用单步预测。对于长时预测，使用模型需要循环预测，初始化系统参数x0已知，公式化表达如下：

1.2 深度动态模式分解

文献[14]基于扩展动态系统（Extend Dynamic Mode Decomposition）提出深度动态模式分解（Deep Dynamic Mode Decomposition，D-DMD），首次使用Autoencoder将系统状态向量映射到线性空间，模型假设经过编码后的数据，在时刻t 与t+1 时刻的zt，zt+1之间是线性关系。系统示意图如图1 所示。

图1 深度动态模式分解模型框架

需要注意的是在模型仅考虑单个时刻间的动力学性质时。在实际使用模型时，t 时刻的系统参数参数xt总是被假设为已知，且被用来预测下一时刻的系统参数xt+1，这样模型本身只具备有单步预测的能力，然而在实际应用中往往需要系统具备更长时间的预测能力。在文献[14]中，上述可以被公式化描述为如下：

其 中 xt,xt+1∈Rn，zt,zt+1∈Rl，其 中 系 统 参 数P={θenc,K,θdec} ，n 是系统原始参数，l 是隐空间参数的维度。 fenc，fdec自编码机的编解码器，参数为Θ={θenc,θdec}，可以是全连接神经网络，对图像数据也可以是卷积神经网络。记优化目标函数为J( K,θ )，上述算法的优化目标如下：

其中超参数λ1，λ2通常选取较小的值，例如0.01，0.001，||*||2表示l2范数，其中神经网络参数Θ 通过梯度下降法训练可以得到。矩阵参数K 通过求解伪逆得到[15]。为了叙述的连贯性，下面给出求解过程。对于经过神经网络编码的数据：Z=[z0,z1,z2,…,zT]，矩阵Zl×T表示隐空间中的数据集合。为了满足线性关系，不妨记: X=[z0,z1,z2,…,zT-1] ，矩 阵 形 式 如（5）所 示。Y=[ z1,z1,z2,…,zT]，如（6）所示，矩阵X，Y ∈Rl×T，参数矩阵K ∈Rl×l。

线性表达式：

则：

其中X†表示矩阵X 的Moore-Penrose 伪逆运算[12，15]。根据（7）式要求解矩阵K，则需要求解矩阵X的伪逆，一般使用奇异值分解[5，15，28，29]的方法来求解矩阵K，首先对矩阵Y 进行奇异值分解，得到X=UΣV*，其中U 为l×l 正交阵，V 为T×T 的正交阵，Σ 为l×T 的奇异矩阵。则有：

通常而言在DMD 算法中需要多次左乘矩阵K，所以文献[5，15，28，29]使特征值特征向量的方法计算。求解线性矩阵的过程中隐含了求解多变量最小均方误差模型[5，15，28，29]：

即优化目标：

整个实际优化过程中，神经网络参数的反向传播和矩阵求解（11）式交替进行。

虽然模型（3）能够很好地表达单步动力系统的关系，但是由于深度学习往往只能找到局部最优解，模型总存在一个较小的误差，导致模型在实际使用时每一步的预测都存在着较小的误差。假设上一步t=i 时刻预测的值为，t 时刻真实值xi，则t 时刻误差为。如果采用长时预测，则下一步t+1 时刻的模型计算输出为：

由于输出数据中误差εi的存在，导致模型输出的不可控制，随着迭代次数的增加，后续模型预测的值与真实值误差越来越大。是要想表达长期的关系，需要在迭代过程中控制单步的误差εi，以减小误差对后续迭代的影响，让单步误差在传递过程中被减小，因此引入Jacobian 正则化[7]的方法来减小在单步迭代中的误差。

1.3 Jacobian正则化方法

正则化方法在多个领域中有引入，在文献[16]中，作者首次使用Jacobian 正则化方法提高卷积神经网络识别物体的鲁棒性，使得模型在对抗攻击上表现良好，模型具备有抗噪的能力。在文献[7]中作者首次使用Jacobian 正则化在长时动力系统预测中，文中给出了误差在下一次预测中如何得到抑制，并给出了完整的分析过程。本文简述了Jacobian 正则如何在非线性动力预测系统中抑制噪声。

由（12）式可知，我们要想对含有噪声的输入进行预测进而获得输出：

于是对f( xi+ εi)在xi处进行泰勒展开，得到：

要想只保留f( xi)项，则需要在迭代的过程中使得除f( xi)以外的项尽可能的小，考虑到εi较小，高阶项忽略，所以有：

也即是Jacobian 正则项。通过（15）式可知，当Jacobian 矩阵的值为0 矩阵时，系动力系统中的误差则不被传递到下一次预测当中，所以在原模型的基础上加上Jacobian 正则。考虑到（16）是矩阵，通常使用Frobenius 范式，则优化目标函数J( K,θ )为：

1.4 Total Least-Square线性模型

为更进一步减小模型误差，需要对线性模型（10）进行扩展，文献[15]首次将线性模型进行广义化，文章指出公式（10）中的模型是不对称的，由于zt，zt+1,t=0,1,2…,T 在时间上具有连贯性，所以zt，zt+1都有相同性质的误差，且指出（10）式是（18）式的一种特殊情况，所以本文将原框架中线性模型（10）进行广义化：

注意到（18）式在输入数据X 上引入了误差∆X，这是因为在多步预测时，当前时刻的预测值（输出值）又作为下一时刻的输入，而当前时刻的数据又包含有误差，所以对于多步预测数据时，模型必须考虑到输入数据中包含误差。其次在单步模型中，输入也是包含有误差的，因为本文算法架构中，线性模型是对经过encoder 编码之后的数据来建模的，所以单步模型也能适用。对于求解模型（18）中的方程即求解总体均方误差（Total Least Square，TLS）的方法在文献[15，17]中给出详细介绍，利用奇异值分解的方法求解该问题。TLS模型作为最小均方误差的扩展模型在文献[9，10，11]中有广泛应用。在文献[15]中首次应用到动态模式求解算法中。对于求解（18）式，一般而言是先是对（18）式中右边展开移项后得到：

对增广矩阵Z ∈R2l×T奇异值分解，得到分解矩阵U2l×2l,Σ2l×T,VT×T，只保留矩阵V 的前l 个列向量Vl，计算得到投影后的时序数据：

再将（21）式的按照（8）式求解，转化成求解伪逆问题即可[15]。

1.5 算法框架

本文算法在实际计算步骤分为两部分，一是Autoencoder 的将数据的编解码，二是使用TLS 求压缩后求解压缩后数据的线性表达式。所以本文在文献[7]的基础上加上Jacobian 正则和TLS，提出精确深度模式分解算法（Exact Dynamic Mode Decomposition，ED-DMD）来建模非线性动力系统。

算法1：精确深度模式分解算法

输入：时序数据X=[x0,x1,x2,x3,…,xT]

输出：模型参数P={ }

θenc,K,θdec

初始化：模型参数θ={ }θenc,θdec，迭代次数Iter，学习率η，因子λ，编码维度l

重复直到迭代次数为Iter:

阶段1：固定参数θ，更新矩阵K

1） 使用参数θ，使用当前encoder 编码时序数据数据，Z=fenc( X )；则Z=[ z0,z1,z2,z3,…,zT]

2） 组织数据：

3） 对Zt奇异值分解，得到U,Σ,V

4） 更新矩阵K:K=Zt+1VΣ-1U*

阶段2：固定参数K，更新网络参数θ

1） 计算损失函数fl(x,θ)

2） 计算Jacobian 矩阵的F 范式JF(x,θ)

3） θ=θ-η∇θ(fl+λJF)

算法1 结束。

由整个算法的过程可知，更新神经网络参数θ 和K 是交替迭代的，更新其中一项参数时，另一项参数保持不变。参数θ 最先是给定初始值，算法交替更新参数时首先更新的是参数K，然后更新神经网络参数θ，交替迭代的过程也可以由图2 看出。为对比较时间复杂，假设输数据维度n，数据量T，m 层encoder 神经网络，每层神元数依次l1,l2…lm，且li≪T(i=1,2…m)，对应的decoder 的m 层网络神经元数为lm…l2,l1。在训练数据相同的情况下，对于原算法[14]提出模型，第i(i=1,2…m)层神经元计算复杂度为O(T*li-1*li)，对于求解线性部分是求解lm个最小二乘回归，其时间复杂度O(lm*lm2*T)，所以总的时间复杂度O(lm3*T)，其中lm个最小人可以并行求解。本文中模型计算量大的部分是求解总体最小二乘线性模型和计算Jacobian 矩阵。求解TLS 主要是对矩阵Z2lm×T求解SVD，选取l(l ≤lm)个特征向量，时间复杂度为O(T2*2l)，计算Jacobian 矩阵时间复杂度O(T*lm2)，总的时间复杂度为O（T2*2l）。

图2 算法流程图

2 实验结果

本文中采两个常用的常微分方程组，分别从这方程组中采取一定长度的轨迹数据，然后用数据来训练模型，直至模型收敛完成训练；在使用模型时只给出初始值x0，作为训练模型的输入数据，后续迭代按照（2）式所示进行。用表示模型循环预测的结果，xi表示采集第i 步的真实数据。最后对比训练好的模型预测数据与真实数据。特别地，在使用模型时，单个数据的误差为：

所有数据的总体误差为：

系统维度n，数据长度为T。在本文模型中，模型参数还需要满足总体误差Esum尽可能的小。由于仅从减小Jacobian 误差和总体均方误差最小化的角度来建模，以期望达到长时间精准建模动力系统。所以并没有明确的将长期误差考虑到优化目标函数（17）之内，且（23）式难以优化。于是采取在迭代过程中取出满足（23）式的模型参数作为最终模型，以达到最好的实验效果。

2.1 范德波（Van der Pol）方程

第一个实验采用的一个电路系统中常用的方程，可以使用如下方程组描述：

方程（24）描述的即是一个二维微分动力系统。参数u 取值u=4.0，在采集方程（24）的数据时，使用系统初始参数为[x0,y0]T=[-1,-4]T作为系统的初始状态。时间间隔Δt=0.01，从t=0 时刻开始采集数据，时间步为3000 步。神经网络采用自编码机，编码部分的神经网络采用三层神经网络，分别使用256，64，32 个神经元。解码器部分使用与之对应的神经元数目。实验结果如图3 所示。

图3

左侧系统状态参数随时间变化，右侧为系统的相图。可以看到模型预测得到的数据在整体上系统还是大致接近真实系统数据，作为第一个阐述性的例子，既描述出了系统参数图，系统参数随时间变化的值（真实值和预测值）；又描述出了系统相图。模型在整体上预测还是较好，但是在局部还是有误差，如图3。

2.2 Glycolysis动力系统

Glycolysis 动力系统方程组是在生物化学领域用于描述糖酵过程，该系统描述了7 个系统状态参数的演化过程。这一组方程常常被用来做基准测试动力模型算法。在文献[14，21，22]都用方程（25）来测试模型，所以：

本文也用它来测试模型的表达能力。数据采集时间间隔0.01，从t=0 时刻开始采集数据，采集4000步。神经网络编码器网络层数为：7，256，128，64，32，32，32。解码器有对应的网络层数。从实验结果可以看出，长期预测结果得到不错性能；整体上，本文模型更接近真实数据，由图4 的误差对比可以看出，本文模型误差在整体上还是比较小，模型性能有所提升。需要注意的是在算法在某些局部空间误差较大，如图4。

图4 模型预测结果误差对比

但是模型在整体上比原模型更接近真实数据如图5。以上两个实验结果初步证明的本文提出的改进的算法在两个微分方程组所确定的数据集上取得不错效果。尤其是对Glycolysis 动力系统的改进效果更加明显。原文算法[14]在仅仅能预测100～500 步，本文算法减小了前一步误差对后续的影响，所以如实验中的数据所示本文提出的算法达到了3000～4000 步的精确预测。

图5 两模型系统参数图和真实数据对比

在模型的训练过程中，模型的损失函数随着时间变化如图6 所示。然而本文算法对三维Lorenz 系统却表现不佳，本文模型能预测500 步左右的状态参数，有很多文献使用Lorenz 作为验证算法模型[6，8]，均不能取得较好的结果，这是由于该系统的Lyapunov 系数为正，微小的系统参数的误差将在后续模型预测中导致巨大误差，误差呈现指数级的增长。相比较于D-DMD算法，虽然本文提出的ED-DMD 算法收敛速度较慢，这是由于新增了Jacobian 正则和总体最小二乘求解运算，增加了计算时间，导致收敛速度较慢；但是随着迭代次数的增加，两个算法都趋于收敛。本文实验在NVIDIA GTX 1060 6GB 显卡和TensorFlow 框架完成。

图6 模型loss值随时间（单位：秒）的变化

本文算法在两个模型中取得较好的结果，误差曲线总体来说也比较低，分析了算法的时间复杂度和模型loss 值随时间的变化，由图6 可以看出本文算法时间消耗没有太大差距。实验结果也证明了算法的优越性。

3 结语

本文提出的算法旨在揭示一般系统的演化规律，即非线性动力学。这些一般的演化规律能够由一组确定的微分方程所表示，系统中的数据是由确定的但未知的规律所确定的。本文从系统数据入手，采用机器学习方法来揭示系统内在的规律，提出的方法能够在一定的时间步中学习到数据的内在演化规律，并在长时预测商保持较低的模型误差，实验结果也证明这点。然而模型还有不足的地方需要改进，例如模型没有明确的建模长期动力系统的关系[8]，模型没有学习到数据长时的内蕴模式[13]，模型没有使用循环神经网络来维持数据之间的时序关系[17-18]等；这些问题希望能在后续的工作中得到解决。