CN113485273B - 一种动态系统时延计算方法及系统 - Google Patents

Abstract

本发明涉及一种动态系统时延计算方法及系统，先根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列。然后根据输入序列计算目标加权邻接矩阵，最后根据目标加权邻接矩阵和输出序列计算图的全局光滑性度量，并对图的全局光滑性度量进行最小化，得到动态系统的时延。本发明所提供的计算方法及系统，能够对未知系统模型的动态系统进行时延计算，应用更为广泛，且能够降低时延计算的复杂度，提高时延计算结果对数据噪声的鲁棒性，进而提高时延计算结果的精度。

Description

一种动态系统时延计算方法及系统

技术领域

本发明涉及系统辨识技术领域，特别是涉及一种基于图Laplacian的动态系统时延计算方法及系统。

背景技术

在大多数过程系统中，时间延迟是很常见的，特别是由于运输时间滞后或由于大容量而滞后或由于传感器分析而滞后(例如，在使用气相色谱仪进行成分分析时)。时间延迟显著影响相关控制系统的性能，因为它们减缓了过程中的响应。因此，时延估计(TDE)及其后续的补偿或处理是控制工程师们面临的一个重要问题。一般情况下，工程师使用时滞模型来模拟实际过程，时延作为时滞模型的关键参数，其取值与时滞模型的性能密切相关。因此，时延估计在系统识别中受到了广泛的关注。近年来有很多专家学者根据时间迟延对象自身的特点提出了许多时间迟延估计方法，如基于常规相关性的方法，基于最小二乘法的方法，基于内插的方法，基于神经网络的方法，基于线性矩阵不等式的方法，基于非交换环理论的方法等。这些时间迟延估计方法可以被分类为两个更宏观的类别，即基于系统模型的方法和基于输入输出数据的方法。

基于系统模型的方法将动态系统看作灰盒子或白盒子，并在时延估计过程中考虑系统建模。如对线性系统进行时延估计时，系统时延可以作为传递函数中的一个参数，根据系统的输入输出信号，便可以通过如最小二乘法等回归方法估计出系统的时延参数。对于非线性系统，利用增量谐波平衡方法估计离散非线性系统的时延，该方法能有效应对数据噪声。针对多入多出时滞系统，可以首先把多入多出的线性时不变系统解耦成多个单入单出的子系统，然后基于Hilbert变换估计子系统的时延。或可以通过连续小波变换与交叉相关相结合的方法实现线性系统时延矩阵的估计。但基于系统模型的方法过度依赖系统结构信息，不适用于未知结构时滞系统的时延估计任务。而现实中的工业系统往往比较复杂，基于系统模型的方法不适用于估计这类系统的时间迟延。

基于输入输出数据的方法把系统看作黑箱，独立于系统模型的方法不需要系统结构先验知识，摆脱了系统建模的束缚。但由于这类方法不需要系统结构作为先验知识，虽然其应用范围更广，应用的对象系统也更复杂，但估计系统时延的难度较大，利用现有的基于输入输出数据的方法对时延进行估计时，估计精度低。

发明内容

本发明的目的是提供一种动态系统时延计算方法及系统，提高动态系统时延计算结果的精度。

为实现上述目的，本发明提供了如下方案：

一种动态系统时延计算方法，所述计算方法包括：

根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

根据所述输入序列计算目标加权邻接矩阵；

根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延。

一种动态系统时延计算系统，所述计算系统包括：

采样模块，用于根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

矩阵计算模块，用于根据所述输入序列计算目标加权邻接矩阵；

时延计算模块，用于根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明所提供的一种动态系统时延计算方法及系统，先根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列。然后根据输入序列计算目标加权邻接矩阵，最后根据目标加权邻接矩阵和输出序列计算图的全局光滑性度量，并对图的全局光滑性度量进行最小化，得到动态系统的时延。本发明所提供的计算方法及系统，能够对未知系统模型的动态系统进行时延计算，应用更为广泛，且能够降低时延计算的复杂度，提高时延计算结果对数据噪声的鲁棒性，进而提高时延计算结果的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例1所提供的计算方法的方法流程图；

图2为本发明实施例1所提供的计算原始加权邻接矩阵所用方法的方法流程图；

图3为本发明实施例1所提供的时延计算的方法流程图；

图4为本发明实施例1所提供的单入单出的线性惯性迟延系统的示意图；

图5为本发明实施例1所提供的固定迟延下线性惯性迟延系统输入输出之间的图的全局光滑性度量；

图6为本发明实施例1所提供的单入单出的线性迟延系统的示意图；

图7为本发明实施例1所提供的固定迟延下线性迟延系统输入输出之间的图的全局光滑性度量；

图8为本发明实施例2所提供的计算系统的系统框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种动态系统时延计算方法及系统，能够在降低计算复杂度的同时提高动态系统时延计算结果的精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例1：

现有的基于输入输出数据的方法包括相关法，当输入为白噪声时，输入输出之间的相关性为脉冲响应，其上升位置表示时间延迟。此外，还有使用小波预滤波器和经典互相关、应用交叉小波变换通过多项式插值、用Copula理论中的Schweizer和Wolff'sσ(SWσ)来测量输入-输出数据的依赖关系等方法得到时延估计结果。但这些方法的计算复杂度高，进而造成计算效率低，且计算精度低。

为了解决上述问题，本实施例用于提供一种动态系统时延计算方法，如图1所示，所述计算方法包括：

S1：根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

具体的，分别对输入信号u(t)和输出信号y(t)以预设周期Δt进行等间隔采样，即每隔固定时间(预设周期Δt)对动态系统的输入信号和输出信号同时进行采样，得到输入序列u＝{u_i}和输出序列y＝{y_i}。需要说明的是，输入序列和输出序列中相邻两个采样点之间的时间间隔相等，均为预设周期Δt。

在对动态系统的输入信号和输出信号进行均匀采样后，本实施例的计算方法还包括：分别对输入序列和输出序列进行归一化，得到归一化后输入序列和归一化后输出序列，并将归一化后输入序列作为新的输入序列，将归一化后输出序列作为新的输出序列，执行步骤S2。对输入序列和输出序列进行归一化即为对输入序列和输出序列的采样值进行归一化。以输入序列为例，确定输入序列所包括的所有采样点的采样值的最大值，对于每一采样点，将其采样值除以最大值，以将每一采样点的采样值归一化到0-1之间。进而通过对输入序列和输出序列进行归一化，能够消除输入输出之间的数值分布范围对计算结果的影响，提高时延计算结果的精度。

S2：根据所述输入序列计算目标加权邻接矩阵；

如图2所示，S2可以包括：

S21：以所述输入序列中的每一个采样点作为一个结点；即采样点与结点一一对应，所有结点组成的集合为结点集合V。结点的个数与采样点的个数相等，均为M。

S22：根据结点位置函数确定每一所述结点的位置；

要计算任意两个结点之间的距离，需要给每个结点赋予一个位置。为此，需要预先定义结点位置函数，结点位置函数可以是一个函数，也可简单表示为一个值，在本实施例中结点位置函数即为输入序列采样点的采样值。具体的，根据结点位置函数，将采样点对应的采样值作为与采样点相对应的结点的位置。

S23：根据所有所述结点的位置，计算任意两个所述结点之间的距离；

S24：根据所述距离，计算所述距离对应的两个所述结点的连接边的权重；

具体的，利用高斯函数计算连接边的权重，计算公式如下：

式1中，w_i，j为第i个结点和第j个结点的连接边的权重，连接边指第i个结点和第j个结点连接时形成的无向边；i＝1，2，...，M；j＝1，2，...，M；M为输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数，σ对加权无向图的几何形式有重要影响，经过计算测试，σ的值定为所有的任意两个结点间距离的均值的3倍。

作为一种可选的实施方式，考虑到距离比较远的两个结点之间的边接边的权重较小，因此对图的全局光滑性度量的贡献也较小，为了减小计算量，在计算连接边的权重时加入距离阈值约束，此时，边接边的权重计算公式如下：

式2中，w_i，j为第i个结点和第j个结点的连接边的权重，连接边指第i个结点和第j个结点连接时形成的无向边；i＝1，2，...，M；j＝1，2，...，M；M为输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数，σ对加权无向图的几何形式有重要影响，经过计算测试，σ的值定为所有的任意两个结点间距离的均值的3倍；δ为距离阈值。

S25：根据所有所述连接边的权重，构建原始加权邻接矩阵。

所构建的原始加权邻接矩阵为M×M矩阵，第i个结点和第j个结点的连接边的权重w_i，j即为加权邻接矩阵中第i行第j列的元素的值，进而根据所有连接边的权重，构建加权邻接矩阵。

此时，可直接将S25所获得的原始加权邻接矩阵记为目标加权邻接矩阵，执行步骤S3。

为了进一步提高计算效率，在S25构建原始加权邻接矩阵之后，S2还包括：根据原始加权邻接矩阵构造N-连接图，并利用N-连接图对原始加权邻接矩阵进行更新，得到更新后加权邻接矩阵，此时，将原始加权邻接矩阵和更新后加权邻接矩阵记为目标加权邻接矩阵，即目标加权邻接矩阵包括原始加权邻接矩阵和更新后加权邻接矩阵，然后再执行步骤S3。

具体的，N-连接图定义为对任意的结点v_i∈V，v_i与距离v_i最近的N个结点均有边相连，即每个结点只保留与距离其最近的N个结点之间的连接边的权重，进而由原始加权邻接矩阵进一步构造得到N-连接图。由于对任意结点v_i来说，离v_i距离越远的结点与v_i之间的连接边的权重也就相对越小，对计算结果的影响也就越小，故本实施例利用N-连接图对原始加权邻接矩阵进行更新，得到更新后加权邻接矩阵，仅保留距离结点v_i最近的N个结点与结点v_i的连接边的权重，将结点v_i与其他结点的连接边的权重置为0，以对原始加权邻接矩阵进行更新，得到更新后加权邻接矩阵，再次减小计算负荷，减小计算复杂度，提高计算效率。

S3：根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延。

如图3所示，S3可以包括：

S31：在周期系数的预设搜索范围内随机确定多个周期系数；

本实施例中的时延是整数倍的预设周期，即：

τ＝kΔt； (3)

式3中，τ为时延；k为周期系数；Δt为预设周期。

故为了计算时延，需要先计算得到周期系数。本实施例可以预先确定周期系数的预设搜索范围为(0，40)，在预设搜索范围内随机确定任意个周期系数。

S32：根据所述输出序列，计算每一所述周期系数对应的计算序列；

计算序列与输出序列之间的时延为k′Δt，k′为S31所得到的周期系数。计算序列可表示为：

即相当于让输出序列y＝{y_i}的第i个采样点的值y_i更改为第i+k′个采样点的值y_i+k′，得到计算序列。

S33：根据所述目标加权邻接矩阵和所述计算序列计算每一所述周期系数对应的图的全局光滑性度量；

图的全局光滑性度量定义为：动态系统的鲁棒性反映在对应的流形上即为流形的光滑性，图的全局光滑性度量则用来刻画流形的光滑性，当信号光滑时，则图的全局光滑性度量很小，反之，则很大。图的全局光滑性度量的计算公式如下：

式4中，S(y)为输入序列y对应的图的全局光滑性度量；

为输入序列y在第i个结点v_i处的梯度；y_j为输入序列y中第j个采样点的采样值；y_i为输入序列y中第i个采样点的采样值；L为拉谱拉斯矩阵；L的计算公式如下：

L＝D-W； (5)

式5中，D为对角矩阵；W为对称矩阵；D和W均根据目标加权邻接矩阵所确定；

当目标加权邻接矩阵仅包括原始加权邻接矩阵时，即未对原始加权邻接矩阵进行更新时，D由原始加权邻接矩阵所确定；W与原始加权邻接矩阵相同；

当目标加权邻接矩阵同时包括原始加权邻接矩阵和更新后加权邻接矩阵时，此时，D仍然由原始加权邻接矩阵所确定，其表达式与式6相同，但W不再与原始加权邻接矩阵相同，而是与更新后加权邻接矩阵相同。构建N-连接图的目的即为只保留每个结点与距离其最近的N个结点之间的连接边的权重，此时的更新后加权邻接矩阵是只保留了每个结点与距离其最近的N个结点之间的连接边的权重，所以此时W并不是每个位置都有元素值，计算量也是在此处减少的。

将S32得到的每一周期系数对应的计算序列替换式4中的输入序列，即可计算得到每一周期系数对应的图的全局光滑性度量。

S34：选取所有所述图的全局光滑性度量中的最小值所对应的周期系数作为实际系数；

当动态系统的输入序列和输出序列不存在时延时，图的全局光滑性度量最小，因此，本实施例选取S33中计算得到的所有图的全局光滑性度量的最小值所对应的周期系数作为实际系数。

实际系数

S35：根据所述实际系数和所述预设周期计算所述动态系统的时延。

实际系数和预设周期的乘积即为动态系统的时延，即以实际系数作为周期系数k，利用式3得到动态系统的时延。

以下，通过仿真实验证明本实施例所用计算方法的优势，图4所示为一个单入单出的线性惯性迟延系统，从系统的传递函数可以看到，该系统有20秒的时延。图5为应用本实施例的方法对该系统进行时延计算的时延估计结果，上图为输入的1000个采样点，中图为输出的1020个采样点，预设周期为1秒，下图为图的全局光滑性度量。可以看到，在时延为20秒时，图的全局光滑性度量非常小，进而通过求解图的全局光滑性度量的最小值即可得到该系统的时延为20秒，与真实时延相同，计算精度高。图6所示为一个单入单出的线性迟延系统，该系统有30秒的迟延。图7为应用本实施例方法对该系统进行时延计算的时延估计结果，上图为输入的2000个采样点，中图为输出的2030个采样点，采样周期为1秒，下图为图的全局光滑性度量。可以看到，在时延为30秒时，图的全局光滑性度量非常小，进而通过求解图的全局光滑性度量的最小值即可得到该系统的时延为30秒，与真实时延相同，计算精度高。

本实施例涉及一种基于图Laplacian的动态系统时延计算方法，属于系统辨识技术领域。计算方法包括：先均匀采样动态系统的输入输出信号并归一化采样值，根据输入序列计算系统输入的原始加权邻接矩阵，利用原始加权邻接矩阵构造N-连接图，计算输入与输出之间的图的全局光滑性度量，通过最小化图的全局光滑性度量得到动态系统的时延估计值。本实施例的时延计算方法可以用于工业系统建模、工业系统过程控制以及过程仿真等领域。如可将本方法应用于风速预测中，以提高风能利用率，降低碳排放。

本实施例所用计算方法针对基于模型的迟延估计方法须在已知系统模型的前提下对系统迟延计算的缺点，提出了对未知模型的系统迟延进行计算的方法，使得算法实用性更加广泛，适用于大部分复杂系统。针对传统的基于输入输出的迟延估计方法计算复杂度高，估计结果对数据噪声敏感的不足，利用基于系统输入输出构建出的图的全局光滑性度量，并通过最小化图的全局光滑性度量来确定系统的时延，降低了时延计算的复杂度，提高了时延计算结果对数据噪声的鲁棒性，进而提高计算精度。

实施例2：

本实施例用于提供一种动态系统时延计算系统，如图8所示，所述计算系统包括：

采样模块M1，用于根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

矩阵计算模块M2，用于根据所述输入序列计算目标加权邻接矩阵；

时延计算模块M3，用于根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延。

本实施例认为任意动态系统的输入与输出之间存在函数关系，并将动态系统的输入输出看作随机变量，计算输入序列各结点之间的距离，得到输入序列的原始加权邻接矩阵，进而得到N-连接图，再由N-连接图计算图的全局光滑性度量，最后对图的全局光滑性度量作最小化处理，得到系统时延，能够解决已有迟延估计系统过分依赖系统结构，不能应用于复杂系统的问题，同时解决迟延估计结果精度不高的问题。

本说明书中每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (5)

1.一种动态系统时延计算方法，其特征在于，所述计算方法包括：

根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

根据所述输入序列计算目标加权邻接矩阵；

根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延；

所述根据所述输入序列计算目标加权邻接矩阵具体包括：

以所述输入序列中的每一个采样点作为一个结点；

根据结点位置函数确定每一所述结点的位置；

根据所有所述结点的位置，计算任意两个所述结点之间的距离；

根据所述距离，计算所述距离对应的两个所述结点的连接边的权重；

根据所有所述连接边的权重，构建原始加权邻接矩阵；

将所述原始加权邻接矩阵记为目标加权邻接矩阵；

所述根据所述距离，计算所述距离对应的两个所述结点的连接边的权重具体包括：

其中，w_i，j为第i个结点和第j个结点的连接边的权重；i＝1，2，...，M；j＝1，2，...，M；M为所述输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数；

或者，

其中，w_i，j为第i个结点和第j个结点的连接边的权重；i＝1，2，...，M；j＝1，2，...，M；M为所述输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数；δ为距离阈值；

所述根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延具体包括：

在周期系数的预设搜索范围内随机确定多个周期系数；

根据所述输出序列，计算每一所述周期系数对应的计算序列；

根据所述目标加权邻接矩阵和所述计算序列计算每一所述周期系数对应的图的全局光滑性度量；

选取所有所述图的全局光滑性度量中的最小值所对应的周期系数作为实际系数；所述实际系数和所述预设周期的乘积即为所述动态系统的时延。

2.根据权利要求1所述的计算方法，其特征在于，在根据所述输入序列计算目标加权邻接矩阵之前，所述计算方法还包括：分别对所述输入序列和所述输出序列进行归一化，得到归一化后输入序列和归一化后输出序列；并将所述归一化后输入序列作为新的输入序列，将所述归一化后输出序列作为新的输出序列。

3.根据权利要求1所述的计算方法，其特征在于，在构建原始加权邻接矩阵之后，所述根据所述输入序列计算目标加权邻接矩阵还包括：

根据所述原始加权邻接矩阵构造N-连接图；

利用所述N-连接图对所述原始加权邻接矩阵进行更新，得到更新后加权邻接矩阵；

将所述原始加权邻接矩阵和所述更新后加权邻接矩阵记为目标加权邻接矩阵。

4.根据权利要求1或3所述的计算方法，其特征在于，所述根据结点位置函数确定每一所述结点的位置具体包括：根据结点位置函数，将所述采样点对应的采样值作为与所述采样点相对应的所述结点的位置。

5.一种动态系统时延计算系统，其特征在于，所述计算系统包括：

采样模块，用于根据预设周期，分别对动态系统的输入信号和输出信号进行采样，得到输入序列和输出序列；

矩阵计算模块，用于根据所述输入序列计算目标加权邻接矩阵；

时延计算模块，用于根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延；

所述根据所述输入序列计算目标加权邻接矩阵具体包括：

以所述输入序列中的每一个采样点作为一个结点；

根据结点位置函数确定每一所述结点的位置；

根据所有所述结点的位置，计算任意两个所述结点之间的距离；

根据所述距离，计算所述距离对应的两个所述结点的连接边的权重；

根据所有所述连接边的权重，构建原始加权邻接矩阵；

将所述原始加权邻接矩阵记为目标加权邻接矩阵；

所述根据所述距离，计算所述距离对应的两个所述结点的连接边的权重具体包括：

其中，w_i，j为第i个结点和第j个结点的连接边的权重；i＝1，2，...，M；j＝1，2，...，M；M为所述输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数；

或者，

其中，w_i，j为第i个结点和第j个结点的连接边的权重；i＝1，2，...，M；j＝1，2，...，M；M为所述输入序列中的采样点的个数；l(v_i)为第i个结点的位置；l(v_j)为第j个结点的位置；dist(l(v_i)，l(v_j))为第i个结点和第j个结点之间的距离；σ为影响参数；δ为距离阈值；

所述根据所述目标加权邻接矩阵和所述输出序列计算图的全局光滑性度量，并对所述图的全局光滑性度量进行最小化，得到所述动态系统的时延具体包括：

在周期系数的预设搜索范围内随机确定多个周期系数；

根据所述输出序列，计算每一所述周期系数对应的计算序列；

根据所述目标加权邻接矩阵和所述计算序列计算每一所述周期系数对应的图的全局光滑性度量；

选取所有所述图的全局光滑性度量中的最小值所对应的周期系数作为实际系数；所述实际系数和所述预设周期的乘积即为所述动态系统的时延。

Images