CN105122004A - 预测向临界转变的系统轨迹 - Google Patents

Abstract

描述了一种用于预测向临界转变的系统轨迹的系统。该系统将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列。随后，确定成对时间序列的传递熵(TE)度量，其中，所述TE度量量化在复杂系统中从源传递至目的地的信息量。确定关联传递熵(ATE)度量，该关联传递熵(ATE)度量将TE的成对时间序列分解为不对称的定向信息流的相关状态，其中，ATE度量包括ATE+影响类和ATE-影响类。该系统随时间估计ATE+轨迹、TE轨迹和ATE-轨迹，并且使用ATE+轨迹、TE轨迹和ATE-轨迹中的至少一个以预测复杂系统中的临界转变。

Description

预测向临界转变的系统轨迹

相关申请的交叉引用

本申请为2013年5月29日在美国提交的标题为“利用动态频谱检测和识别方向影响”的美国非临时申请No.13/904,945的部分继续申请。本申请还是在2013年3月14日提交的标题为“向临界转变的预测系统轨迹”的美国临时申请No.61/784,167的非临时专利申请。

技术领域

本发明涉及一种用于自动检测临界转变的系统，更具体地，涉及一种利用预测的系统轨迹来自动检测临界转变的系统。

背景技术

复杂的社交-技术系统的一个核心特征是其自组织的新兴交互，自组织的新兴交互提供了有效地交换信息和资源的益处，但与此同时增加了传播攻击或故障的风险和速度。在高风险不稳定区域中操作的复杂系统上的一个小扰动可诱发导致不可阻挡的致命故障的临界转变。

最近有一些关于用于异构网动态系统的早期预警信号的工作。Kwon和Yang(见并入的参考文献列表，参考文献1)使用传递熵来分析金融市场数据。他们计算成对的传递熵以形成传递熵矩阵，并利用传递熵论证来自成熟市场的非对称影响。他们的方法估计全局传递熵，而不是随时间变化的局部传递熵；因此，他们的方法无法处理在定向影响下的动态和结构变化。

Staniek和Lehnertz(见参考文献4)介绍了一种利用符号技术来估计传递熵的鲁棒性和计算上有效的计算方法。他们对脑电活动数据的符号传递熵(STE)的论证表明STE能够检测不对称的依赖关系，并在不观察实际癫痫活动的情况下识别含有癫痫病灶的脑半球。然而，他们的重点也不是识别复杂系统结构的动态。

Scheffer等(见参考文献5)研究了起源于生态学领域的先进方法。这些方法只注重具有诸如种群动态的增加的时间相关性、偏态和空间相关性的信号的同构格，因此无法处理异构网动态系统。

Moon和Lu(见参考文献2和3)研发了一种频谱早期预警信号(EWS)理论，该理论检测临界转变的到来，并且利用协方差矩阵估计临界转变附近的系统结构和网络连接性。虽然频谱EWS量化了有多少系统实体一起移动，但是协方差频谱的对称性不允许在实体之间进行定向影响的分析。

因此，一直需要一种利用允许在实体间进行定向性影响的分析的信息动态频谱架构在同构网动态系统中动态检测临界转变的系统。

发明内容

本发明涉及一种用于预测向临界转变的系统轨迹的系统。该系统包括：一个或更多个处理器和存储器，所述存储器具有指令，使得当所述指令被执行时，所述一个或更多个处理器执行多个操作。将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列。确定成对序列的传递熵(TE)度量。确定关联传递熵(ATE)度量。ATE度量包括ATE+正影响类和ATE-负影响类。然后，随时间估计ATE+轨迹、TE轨迹和ATE-轨迹。最后，利用所述ATE+轨迹、所述TE轨迹和所述ATE-轨迹中的至少一个轨迹预测所述复杂系统中的临界转变。

在另一方面，所述系统根据下式g(t)＝aln(t)+c利用未知系数a和常数c通过自然对数ln估计所述TE轨迹，

其中，t表示时间。

在另一方面，所述系统通过对g(t)＝aln(t)+c求导以获得然后根据下式获得所述未知系数a的离散形式，来利用多个离散时间步长估计所述未知系数a。

在另一方面，所述系统根据以下确定所述未知系数a：

在窗口[T_start，T_end]内针对固定时间步长Δt，确定以下矩阵方程：

$\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ .\\ .\\ .\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ .\\ .\\ .\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right],$

其中，t₁，…，t_k+Δt∈[T_start，T_end]；以及

根据下式针对所述固定时间步长Δt确定a_Δt的近似值：

$a_{\mathrm{\Δ}t}={\mathrm{argmin}}_a{||\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ .\\ .\\ .\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ .\\ .\\ .\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right]||}_2,$

其中，argmin表示最小化，t_k表示第k个时间步长，其中，与该时间步长相对应的常数c为c_Δt＝g(T_end)-a_Δtln(T_end)。

在另一方面，所述系统根据下式针对未来时间T_end+t_d确定预测值：

F_Δt(t_d)＝a_Δtln(T_end+t_d)+c_Δt，其中，t_d是用于预测的期望时间长度。

附图说明

本发明的目的，特点和优点从以下参照以下附图对本发明各个方面的详细说明将是显而易见的，其中：

图1是描述根据本发明的原理的用于临界转变系统的自动检测的系统的组件的框图；

图2是根据本发明的原理的计算机程序产品的示意图；

图3A示出根据现有技术的电压随时间变化的非寄养电路；

图3B示出根据现有技术的随时间的传递熵(TE)/关联传递熵(ATE)+/ATE-曲线；

图4A示出根据本发明原理的与指数函数的图拟合的曲线；

图4B示出根据本发明原理的与对数函数的图拟合的曲线；

图5是将根据本发明原理的预测轨迹与实际轨迹进行比较的图；

图6A示出根据本发明原理产生的预测锥(predictioncone)和预测轨迹；

图6B示出根据本发明原理产生的预测锥和预测轨迹；以及

图7是根据本发明原理的用于临界转变的自动检测的系统的流程图。

具体实施方式

本发明涉及一种用于自动检测临界转变的系统，更具体地，涉及一种利用预测的系统轨迹自动检测临界转变的系统。

呈现下面描述以使本领域的普通技术人员能够构造和使用本发明并将它合并到特定应用的上下文中。不同应用中的各种修改以及多种用途对本领域的技术人员将是显而易见的，并且在此定义的一般原理可以应用到宽范围的多方面。因此，没有意图将本发明限于所呈现的方面，而是与本文公开的原理和新颖特征一致给予最宽范围。

在下面的详细描述中，阐述了许多具体细节以便提供本发明的更透彻理解。然而，对于本领域的技术人员而言明显的是，可以实施本发明而没有必要限于这些具体细节。在其它实例中，为了避免模糊本发明，公知的结构和设备以框图的形式被示出，而非详细示出。

读者的注意力致力于与本说明书同时提交且和本说明书开放给公众查阅的所有论文和文献，所有这样的论文和文献的内容合并于此以供参考。除非有明确说明，本说明书公开的所有特征(包括任何所附权利要求、摘要和附图)可以由用于相同、等同或类似目的另选特征来代替。因此，除非有明确说明，否则公开的每个特征仅是通用系列的等同或类似特征的一个示例。

此外，如在35U.S.C.第112条第6款中规定的，权利要求中的任何元素在没有明确说明用于执行特定功能的“装置”或用于执行特定功能的“”，不应理解为“装置”或“步骤”条款。具体地，本文权利要求中使用“的步骤”或“的动作”并不意味着调用(invoke)35U.S.C.第112条第6款的规定。

在详细描述本发明之前，首先给出引用的参考文献的列表。然后，提供本发明的多个主要方面的描述。随后，概述向读者提供本发明整体理解。最后，提供本发明的具体细节以便给出特定方面的理解。

(1)合并的引用参考文献的列表。

下面的参考文献在本申请中自始至终被引用。为了清楚和方便起见，参考文献在此列出作为读者的中心资源。以下参考文献通过引用并入本文，如同在此完全阐述。本申请中的参考文献通过相应的参考文献编号来引用。

1.O.Kwon和J.-S.Yang，InformationFlowbetweenStockIndices，2008EPL8268003.

2.H.Moon和T.-C.Lu.EarlyWarningSignalofComplexSystems：NetworkSpectrumandCriticalTransitions，WIN2010.

3.H.Moon和T.-C.Lu.NetworkCatastrophe：Self-OrganizedPatternsRevealboththeInstabilityandtheStructureofComplexNetworks，preprint，2012.

4.M.Staniek和K.Lehnertz.SymbolicTransferEntropy，PhysicalReviewLetters100，15801，2008.

5.M.Scheffer，J.Bascompte，W.A.Brock，V.Brovkin，S.R.Carpenter，VDakos，H.Held，E.H.vanNes，M.Rietkerk和G.Sugihara，Early-warningsignalsforcriticaltransitions.Nature，461，2009.

6.T.Schreiber，MeasuringinformationTransfer，Phys.Rev，Lett.85，461，2000.

(2)主要方面

本发明具有三个“主要方面”，第一个主要方面是用于自动检测临界转变的系统。该系统通常是计算机系统操作软件的形式或者是“硬编码”指令集的形式。该系统可并入到提供不同功能的多种设备中。第二个主要方面是通过是软件形式并且通过使用数据处理系统(计算机)来操作的一种方法。第三个主要方面是一种计算机程序产品。该计算机程序产品通常代表存储在非瞬时计算机可读介质(例如，诸如紧凑盘(CD)或数字多功能光盘(DVD)的光学存储设备，或诸如软盘或磁带的磁存储装置)上的计算机可读指令。计算机可读介质的其它的非限制性示例包括硬盘、只读存储器(ROM)以及闪存类型的存储器。这些方面将在下面更详细地描述。

图1中提供描述本发明的系统(即，计算机系统100)的示例的框图。计算机系统100被配置为执行与程序或算法相关联的计算、处理、操作和/或功能。在一方面，本文讨论的特定处理和步骤被实现为驻留在计算机可读存储单元内且由计算机系统100的一个或更多个处理器执行的一系列指令(例如，软件程序)。当所述指令被执行时，指令使计算机系统100执行特定动作，并表现出如本文所述的特定行为。

计算机系统100可包括配置为传输信息的地址/数据总线102。另外，一个或更多个数据处理单元(例如，一个处理器104(或多个处理器))与地址/数据总线102连接。处理器104被配置为处理信息和指令。在一方面，处理器104是微处理器。另选地，处理器104可以是不同类型的处理器，例如并行处理器或现场可编程门阵列。

计算机系统100被配置为利用一个或更多个数据存储单元。计算机系统100可以包括与地址/数据总线102连接的一个易失性存储单元106(例如，随机存取存储器(“RAM”)、静态RAM、动态RAM等)，其中，易失性存储单元106被配置成存储用于处理器104的信息和指令。计算机系统100还可以包括与地址/数据总线102连接的非易失性存储器单元108(例如，只读存储器(“ROM”)、可编程ROM(“PROM”)、可擦除可编程ROM(“EPROM”)、电可擦除可编程ROM(“EEPROM”)、闪存等)，其中，非易失性存储单元108被配置为用于处理器104的静态信息和指令。另选地，计算机系统100可以例如在“云”计算中执行从在线数据存储单元检索到的指令。在一方面，计算机系统100还可以包括与地址/数据总线102连接的一个或更多个接口(例如，接口110)。一个或更多个接口被配置为使计算机系统100与其它电子设备和计算机系统交互。由所述一个或更多个接口实现的通信接口可以包括有线(例如，串行电缆，调制解调器，网络适配器等)和/或无线(例如，无线调制解调器，无线网络适配器等)通信技术。

在一方面，计算机系统100可以包括与地址/数据总线102连接的输入设备112，其中，输入设备112被配置成向处理器100传送信息和命令选择。根据一方面，输入设备112是可包括字母数字和/或功能键的字母数字输入设备(例如键盘)。另选地，输入设备112可以是除了字母数字输入设备之外的输入设备。在一方面，计算机系统100可以包括与地址/数据总线102连接的光标控制设备114，其中，所述光标控制设备114被配置成向处理器100传送用户输入信息和/或命令选择。在一方面，使用例如鼠标、跟踪球、跟踪垫、光学跟踪设备或触摸屏来实现所述光标控制设备114。尽管前述如此，在一方面，例如响应于与输入设备112相关的特定键的使用和键序列命令，经由来自输入设备112的输入来指导和/或激活光标控制设备114。在另选方面，光标控制设备114被配置成通过语音命令而被指导或引导。

在一方面，计算机系统100还可以包括与地址/数据总线102连接的一个或更多个可选的计算机可用的数据存储设备(例如，存储设备116)。存储设备116被配置成存储信息和/或计算机可执行指令。在一方面，存储设备116是诸如磁盘或光盘驱动器(例如，硬盘驱动器(“HDD”)、软盘、紧凑盘只读存储器(“CD-ROM”)、数字多功能盘(“DVD”))的存储设备。按照一方面，显示设备118与地址/数据总线102连接，其中，显示设备118被配置为显示视频和/或图形。在一方面，显示设备118可以包括阴极射线管(“CRT”)、液晶显示器(“LCD”)、场发射显示器(“FED”)、等离子显示器、或任何其它适合于显示视频和/或图形图像和字母数字字符的用户可识别的显示装置。

这里所呈现的计算机系统100是根据一方面的示例计算环境。但是，计算机系统100的非限制性示例不是严格限定为计算机系统。例如，一方面提供了计算机系统100表示用于根据本文描述的多方面的一种数据处理分析。并且，其它计算系统也可以被实现。事实上，本技术的精神和范围不限于任何单一的数据处理环境。因此，在一方面，本技术的多个方面的一个或更多个操作通过使用由计算机执行的诸如程序模块的计算机可执行指令来控制或实现。在一个实现中，这样的程序模块包括被配置为执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件和/或数据结构。此外，一方面提供了本技术的一个或更多个方面通过利用一个或更多个分布式计算环境来实现，例如特定任务由通过通信网络链接的远程处理设备来执行，或者例如各种程序模块位于包括存储器存储设备的本地和远程计算机存储介质二者。

图2中描述了实施本发明的计算机程序产品(即，存储设备)的示图。所述计算机程序产品被描述为软盘200或诸如CD或DVD的光盘202。但是，如前所述，该计算机程序产品通常表示存储在任何兼容的非瞬时计算机可读介质中的计算机可读指令。针对本发明使用的术语“指令”通常表示要在计算机上执行的一组操作，并且可以表示整个程序的多个片段或单独的、可分离的、软件模块。“指令”的非限制示例包括计算机程序代码(源或目标代码)和“硬编码”电子(即编码到计算机芯片的计算机操作)。“指令”被存储在任何非瞬时计算机可读介质上，例如在计算机的存储器中、或软盘、CD-ROM以及闪存驱动器上。在上述任何一种情况下，指令被编码在非瞬时计算机可读介质上。

(3)介绍

复杂的社交技术系统的一个核心特征是自组织新兴交互，自组织新兴交互提供了有效地交换信息和资源的益处，但同时也增加了传播攻击或故障的风险和速度。在高风险不稳定区域中操作的复杂系统上的一个小扰动可诱发导致不可阻挡的致命故障的临界转变。在根据本发明原理的系统中，可通过检测临界转变的早期预警并且预测系统轨迹的似然(likelihood)和临界点的前置时间(leadtime)来避免致命故障。

本文所描述的发明是建立在标题为“利用动态频谱检测和识别定向影响(DetectionandIdentificationofDirectionalInfluencesUsingDynamicSpectrum)”美国申请No.13/904,945(以下简称为′945申请)中描述的发明的基础上，该美国申请No.13/904,945是一个信息动态频谱框架，其越过单向扩散动态以研究具有瞬时定向性影响的异构网动态系统。该′945申请通过引用并入本文，就如同在这里完全阐述。

该′945申请描述了一种新的关联传递熵(ATE)度量的发展，其将传递熵的成对定向性影响分解为非对称的、定向性的信息流相关联的状态。复杂系统的多元时间序列被变换成传递熵矩阵(TEM)和关联传递熵矩阵(ATEM)的频谱以捕获该系统的信息动态。TEM和ATEM的新型频谱半径测量被开发以检测源驱动的不稳定性的早期预警信号，并诱发定向影响结构识别源并揭示定向影响的动态。在′945申请中描述的框架使人们不仅可以准确地检测、预测和缓和复杂系统内突发行为改变，还使人们理解、设计和构建更加可持续的和适应性强的复杂系统。

关联传递熵的构思是通过限制两个处理的关联状态来分解传递熵(TE)。TE最简单的情况是两个影响类：一个正的ATE+和一个负的ATE-。在这种条件下，ATE能够区分信息传送的量相同但具有相反效果的两种情况：源以与ATE+相同的方向驱动目的地，和源以相反的ATE-定向驱动目的地。对于动态数据，这意味着信息流的量不是恒定的，计算本地时间窗内的时间序列的TE/ATE，使得TE/ATE变为时间的函数。这种分解进一步使人们能够以ATEM+和ATEM-的形式研究ATEM。

假定ATE+/TE的交叉签名为临界转变的早期预警信号，本发明预测(1)系统轨迹的似然和(2)临界点的前置时间。首先观察到ATE+和TE被最大化，同时在系统切换到另选稳定的机制之前ATE-变平坦。因此，基于模型的统计预测被应用以估计ATE+/TE/ATE-轨迹的似然。通过基于ATE-轨迹的预期平坦估计临界点来预测前置时间。根据本发明原理的预测方法将随着系统进展而连续地输出概率锥(probabilisticcone)的更新。确定所述似然和临界点的前置时间将使缓解策略进一步发展以避免临界转变。这些方面中的每一方面将在下面进一步详细地描述。

(4)发明的具体细节

描述了这样一种系统，该系统解决了自动检测出现的转变的需求，并从复杂系统的时间序列识别不稳定的来源。根据本发明原理的技术可用于检测除了固有动态之外的定向影响，从而得出封装的交互模块以预测行为轨迹。根据本发明原理的系统在分析复杂网络的动态中将得以改进。

本文所描述系统和方法的优点包括但不限于：(1)提供一种对关联的、非对称的定向性影响定量的有效测量方法，而其它方法仅可处理对称的或定向性的影响；(2)提供一种用于捕获系统级的定向影响的有效公式，而其它方法只能处理成对的定向影响或系统级的对称影响；(3)提供了对在具有定向影响动态的系统中检测不稳定性的有效测量，而其它方法假设非定向扩散动态；(4)提供了一种表征汇集的定向影响的不稳定性趋势的有效方法，而其它方法只能用简单阈值检测不稳定性；以及(5)提供了一种揭示改变定向影响的动态的源和结构的有效方法，而其它方法只能处理固定的动态。

(4，1)背景：信息动态频谱

非对称性测量，传递熵

$T_{x_j\→x_i}={\Σ}_{(x_{i,t+1},x_{t,t},x_{j,t})\∈D}p(x_{i,t+1},x_{i,t},x_{j,t})\log \frac{p\left(x_{i,t+1}\right|x_{i,t},x_{j,t})}{p\left(x_{i,t+1}\right|x_{i,t})},---\left(1\right)$

量化从源x_j至目的地x_i的信息传输的量。该公式被用于通过分析股市指数来论证从成熟市场到新兴市场的不对称影响，如参考文献1和6所述。然而，现有技术的当前状态只考虑平均成对信息传输和统计分析，而不包括用于解释定向影响的出现的信息动态频谱的更基本的对象(即，随时间变化的信息传递)。

在′945申请中描述的之前研发的信息动态频谱理论，复杂网络的动态表示为：dX＝F(X)dt+σdW，其中是解释m个元素X(t)＝[x₁(t)x₂(t)...x_m(t)]^T的系统的动态的可微函数，且通过F(X(t))＝[f₁(X(t))f₂(X(t))...f_m(X(t))]^T来定义。这捕获了单个实体的动态和实体间的网络相互作用(interaction)二者。例如，元素x_i(t)可以是第i个交易日的第i个股票指数。最后一项σdW与按照σ缩放的白噪声相对应。为了揭示未知的F(X)，对于邻近X₀的X，此问题通过线性近似F(X)≈F(X₀)+J(X₀)(X-X₀)被进一步简化，其中是雅可比矩阵。因此，局部时间动态可成对获得。对于i≠j，状态D_k的测量关联传递熵(ATE)通过以下来定义：

$T_{x_j\→x_i}^{D_k}={\Σ}_{(x_{i,t+1},x_{i,t},x_{j,t})\∈D_k}p(x_{i,t+1},x_{i,t},x_{j,t})\log \frac{p\left(x_{i,t+1}\right|x_{i,t},x_{j,t})}{p\left(x_{i,t+1}\right|x_{i,t})}---\left(2\right)$

其中D_k是表示x_i和x_j之间的特定关联状态的D的子集。ATE的目的是捕获针对特定状态关联的两个变量之间的信息传递。例如，TE可分解成两个影响类：一个正的和一个负的。

如果与节点之间的相互作用相比，每一节点内的内部动态可以忽略不计，则对于某线性函数g可推测使为关联传递熵矩阵(ATEM)：

$T^{D_k}=\left[\begin{array}{cccc}{T}_{x_1\→x_1}^{D_k}& T_{x_2\→x_1}^{D_k}& \·\·\·& T_{x_m\→x_1}^{D_k}\\ T_{x_1\→x_2}^{D_k}& T_{x_2\→x_2}^{D_k}& \·\·\·& T_{x_m\→x_2}^{D_x}\\ & \·& & \\ & \·& & \\ & \·& & \\ T_{x_1\→x_m}^{D_k}& T_{x_2\→x_m}^{D_k}& \·\·\·& T_{x_m\→x_m}^{D_k}\end{array}\right]---\left(3\right)$

为了适当地处理动态数据，这意味着从一个到另一个的信息流的量不是静态的，TE/ATE曲线作为时间的函数生成，并通过滑动窗口来计算。

由于TEM矩阵是非对称的，所以它的特征值是复数值。使用TEM矩阵的频谱半径以测量整个网络的信息流的总量。

为了有效地估计连续数据的传递熵，Staniek和Lehnertz(见参考文献4)提出了使用符号化技术的方法，该方法是鲁棒性的且计算速度快，引入排列熵的概念。在’945申请中，Staniek和Lehnertz的符号化技术适于计算ATE。具体地说，通过针对1，2，...，m}对{x(t)，x(t+1)，...x(t+m)}的值求一阶导并且表示＝在t处的符号x(t)的m阶的相关排列，来将连续值的时间序列符号化。然后，通过计算的ATE来估计的ATE。

图3A和图3B示出从′945申请中描述的发明的非寄养网络的ATE分析。电路在没有振荡的稳定区域中初始地操作。然后，添加小的扰动。不知道的是所述电路将变得振荡(不稳定)还是保持稳定。执行ATE分析以检测电路是否将变得同步(不稳定)。图3A中的绘图是随时间的电路中的电压，其中，Ant1表示天线1，Ant2表示天线2，i1表示端口1，i2表示端口2。图3B中的绘图示出了随时间的TE/ATE+/ATE-曲线。所述曲线分别是从TEM/ATEM+/ATEM的频谱的绝对值之和获得的。发现TEM/ATEM+/ATEM-的频谱半径也显示类似的结果，但是由于TEM是ATEM+和ATEM-的和，所以ATE+曲线决不会与TE曲线交叉。因此，频谱的绝对值之和提供了更多信息。如图3B所示，ATE+曲线在时间点800附近与TE曲线交叉，这表明同步性增加。此外，ATE+曲线恰在完全同步之前达到其波峰，而ATE-曲线变得平坦，这是由于不存在负相关性。

(4.2)用于轨迹预测的概率锥

在根据本发明原理的系统中，使用基于模型的统计预测以估计TE/ATE+/ATE-轨迹。分析TE/ATE曲线的优点是TE/ATE去除了原始数据的尖峰，如图3A中所示的电路数据的例子。观察到，可以根据下式由自然对数与未知系数a和常数c来很好地近似这种情况下的TE曲线：

g(t)＝aln(t)+c，(4)

图4A示出与指数函数的绘图拟合的曲线，图4B示出与对数函数的绘图拟合的曲线。粗体实曲线400是来自图3B的TE曲线。粗体虚曲线402是在间隔t∈[700，1100]内的拟合的曲线。非粗体虚曲线404是从拟合的曲线的预测。图4B示出对数适合于拟合TE曲线。

为了从TE/ATE曲线估计系数a，用多个离散时间步长来估计变化率，而不是用公式(4)中的对数函数确定性地对它们进行拟合。这就产生了预测锥，预测锥将在下面详细描述。从公式(4)的导数

$g^{,}\left(t\right)=\frac{a}{t},---\left(5\right)$

获得未知系数a的下面的离散形式：

$\frac{\mathrm{\Δ}g}{\mathrm{\Δ}t}=a\frac{1}{t},---\left(6\right)$

对于时间窗[T_start，T_end]内的固定时间步长Δt，使用最小二乘法来求解未知系数a。首先，写出如下矩阵方程：

$\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ .\\ .\\ .\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ .\\ .\\ .\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right],---\left(7\right)$

其中，t₁，…，t_k+Δt∈[T_start，T_end]，t_k表示第k时间点。

随后，对于固定时间步长Δt，a_Δt的近似值根据以式获得：

$a_{\mathrm{\Δ}t}={\mathrm{argmin}}_a\left|\right|\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ .\\ .\\ .\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ .\\ .\\ .\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right]|{|}_2---\left(8\right)$

然后，与该时间步长相对应的常数为c_Δt＝g(T_end)-a_Δtln(T_end)。因此，对于未来时间T_end+t_d的预测值为：

F_Δt(t_d)＝a_Δtln(T_end+t_d)+c_Δt(9)

t_d是针对预测期望的时间长度。

既然已经建立了用于估计常数c和a的方法，随着每个变化Δt，可以在每个时间点产生概率预测锥(或光锥(lightcone))，以获得c和a的多个估计值。因此，在给定时间的概率预测锥由从该点开始的自然对数曲线的集合组成。

为了针对紧接着的下一时间步长估计该预测的误差，进行如下计算：

error_Δt(i)＝g(t_i+Δt)-[a_Δtln(t_i+Δt)+c(i)]，(10)

其中，c(i)＝g(t_i-1)-a_Δtln(t_i-1)(11)

使E是所有时间步长Δt上的所有error_Δt的集合，并使σ＝E的标准偏差。使G为所有时间步长Δt上所有预测值G_Δt(t_d)的集合，且μ＝mean(G)。因此，对于未来时间T_end+t_d的95％置信区间是：

$CI=\[\mathrm{\μ}-1.96\frac{\mathrm{\σ}}{\sqrt{N}},\mathrm{\μ}+1.96\frac{\mathrm{\σ}}{\sqrt{N}}\]---\left(12\right)$

其中N是E的大小。

图5示出了非粗体实线500中的预测轨迹的95％置信区间。TE的预测轨迹由非粗体虚线502表示。粗体实线504描绘了TE的实际轨迹。如上所述，根据预测锥生成预测值和置信区间。实际轨迹的两个点在95％置信区间外，恰好在时间t＝1500附近非寄养电路完全同步之前。这些点由菱形图(diamonds)504表示，并表明了向临界转变的趋势的概率。

图6A和6B示出根据上述方法生成的预测锥600和预测轨迹值(由粗体虚曲线602所示)的两个时间快照。使用该预测锥以生成具有图5中的置信区间的预测值。粗体实曲线604代表实际的TE轨迹。

图7是描绘根据本发明的原理预测向临界转变的系统轨迹的系统的流程图。在第一步骤700，将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列。在第二步骤702，该系统确定成对时间序列的TE度量。在第三步骤704，确定ATE度量。在第四步骤706，该系统随时间估计ATE+、TE和ATE-轨迹。在第五步骤708，利用ATE+、TE和ATE-轨迹中的至少一个预测系统中的临界转变。

综上所述，本发明是一种用于分析复杂系统的时间序列以检测向临界转变的趋势的系统和方法，特别是用于具有定向影响动态的复杂系统。所述系统和方法可以被采用为嵌入式的决策支持模块，以提供复杂系统中的临界转变的早期指示，包括但不限于，车辆健康管理、对复杂电子的预报、供应链中的危机管理、社会动荡、金融市场不稳定性和大脑网络的癫痫发作。本文描述的技术将准确检测系统行为的即将到来的临界转变，准确定性向临界转变的不稳定趋势，准确识别不稳定的来源和传播结构的变化。本发明的应用包括但不限于，车辆健康管理、复杂电子、智能电网、社交网络分析、社会不稳定、供应链危机管理、网络安全和癫痫发作的检测。

Claims (15)

1.一种用于预测向临界转变的系统轨迹的系统，该系统包括：

一个或更多个处理器和非瞬时计算机可读介质，在所述非瞬时计算机可读介质上编码有可执行指令，使得当所述可执行指令被执行时，所述一个或更多个处理器执行以下操作：

将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列；

确定成对时间序列的传递熵TE度量；

确定关联传递熵ATE度量，其中，所述ATE度量包括ATE+正影响类和ATE-负影响类；

随时间估计ATE+轨迹、TE轨迹和ATE-轨迹；以及

利用所述ATE+轨迹、所述TE轨迹和所述ATE-轨迹中的至少一个轨迹预测所述复杂系统中的临界转变。

2.根据权利要求1所述的系统，其中，所述一个或更多个处理器还执行以下操作：根据下式g(t)＝aln(t)+c利用未知系数a和常数c通过自然对数ln估计所述TE轨迹，

其中，t表示时间。

3.根据权利要求2所述的系统，其中，所述一个或更多个处理器还执行以下操作：通过对g(t)＝aln(t)+c求导以获得然后根据下式获得所述未知系数a的离散形式，来利用多个离散时间步长估计所述未知系数a。

4.根据权利要求3所述的系统，其中，所述一个或更多个处理器还执行根据以下确定所述未知系数a的操作：

在窗口[T_start，T_end]内针对固定时间步长Δt，确定以下矩阵方程：

$\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right],$

其中，t₁，…，t_k+Δt∈[T_start，T_end]；以及

根据下式针对所述固定时间步长Δt确定a_Δt的近似值：

$a_{\mathrm{\Δ}t}={\mathrm{argmin}}_a\left|\right|\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right]|{|}_2,$

其中，argmin表示最小化，t_k表示第k个时间步长，其中，与该时间步长相对应的常数c为c_Δt＝g(T_end)-a_Δtln(T_end)。

5.根据权利要求4所述的系统，其中，所述一个或更多个处理器还执行根据下式针对未来时间T_end+t_d确定预测值的操作：

F_Δt(t_d)＝a_Δtln(T_end+t_d)+c_Δt，

其中，t_d是用于预测的期望时间长度。

6.一种用于预测向临界转变的系统轨迹的计算机实现的方法，所述方法包括以下动作：

使一个或更多个处理器执行存储在非瞬时存储器上的指令，使得当执行时，所述一个或更多个处理器执行以下操作：

将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列；

确定成对时间序列的传递熵TE度量；

确定关联传递熵ATE度量，其中，所述ATE度量包括ATE+正影响类和ATE-负影响类；

随时间估计ATE+轨迹、TE轨迹和ATE-轨迹；以及

利用所述ATE+轨迹、所述TE轨迹和所述ATE-轨迹中的至少一个轨迹预测所述复杂系统中的临界转变。

7.根据权利要求6所述的方法，其中，所述一个或更多个处理器还执行以下操作：根据下式g(t)＝aln(t)+c利用未知系数a和常数c通过自然对数ln估计所述TE轨迹，

其中，t表示时间。

8.根据权利要求7所述的方法，其中，所述数据处理器还执行以下操作：通过对g(t)＝aln(t)+c求导以获得然后根据下式获得所述未知系数a的离散形式，来利用多个离散时间步长估计所述未知系数a。

9.根据权利要求8所述的方法，其中，所述数据处理器还执行根据以下确定所述未知系数a的操作：

在窗口[T_start，T_end]内针对固定时间步长Δt，确定以下矩阵方程：

$\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right],$

其中，t₁，…，t_k+Δt∈[T_start，T_end]；以及

根据下式针对所述固定时间步长Δt确定a_Δt的近似值：

$a_{\mathrm{\Δ}t}={\mathrm{argmin}}_a\left|\right|\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right]|{|}_2,$

其中，argmin表示最小化，t_k表示第k个时间步长，其中，与该时间步长相对应的常数c为c_Δt＝g(T_end)-a_Δtln(T_end)。

10.根据权利要求9所述的方法，其中，所述数据处理器还执行根据下式针对未来时间T_end+t_d确定预测值的操作：

F_Δt(t_d)＝a_Δtln(T_end+t_d)+c_Δt，

其中，t_d是用于预测的期望时间长度。

11.一种用于预测向临界转变的系统轨迹的计算机程序产品，所述计算机程序产品包括存储在非瞬时计算机可读介质上的计算机可读指令，所述计算机可读指令能被具有处理器的计算机执行，以使得所述处理器执行如下操作：

将复杂系统的可观测的一组多元时间序列变换成一组符号多元时间序列；

确定成对时间序列的传递熵TE度量；

确定关联传递熵ATE度量，其中，所述ATE度量包括ATE+正影响类和ATE-负影响类；

随时间估计ATE+轨迹、TE轨迹和ATE-轨迹；以及

利用所述ATE+轨迹、所述TE轨迹和所述ATE-轨迹中的至少一个轨迹预测所述复杂系统中的临界转变。

12.根据权利要求11所述的计算机程序产品，所述计算机程序产品还包括用于使所述处理器执行以下操作的指令：根据下式g(t)＝aln(t)+c利用未知系数a和常数c通过自然对数ln估计所述TE轨迹，

其中，t表示时间。

13.根据权利要求12所述的计算机程序产品，所述计算机程序产品进一步包括使所述处理器执行以下操作的指令：通过对g(t)＝aln(t)+c求导以获得然后根据下式获得所述未知系数a的离散形式，来利用多个离散时间步长估计所述未知系数a。

14.根据权利要求13所述的计算机程序产品，所述计算机程序产品还包括使所述处理器执行根据以下确定所述未知系数a的操作的指令：

在窗口[T_start，T_end]内针对固定时间步长Δt，确定以下矩阵方程：

$\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right],$

其中，t₁，…，t_k+Δt∈[T_start，T_end]；以及

根据下式针对所述固定时间步长Δt确定a_Δt的近似值：

$a_{\mathrm{\Δ}t}={\mathrm{argmin}}_a\left|\right|\left[\begin{array}{c}\frac{1}{t_1+\mathrm{\Δ}t/2}\\ \frac{1}{t_2+\mathrm{\Δ}t/2}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{1}{t_k+\mathrm{\Δ}t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1+\mathrm{\Δ}t)-g\left(t_1\right)}{\mathrm{\Δ}t}\\ \frac{g(t_2+\mathrm{\Δ}t)-g\left(t_2\right)}{\mathrm{\Δ}t}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{g(t_k+\mathrm{\Δ}t)-g\left(t_k\right)}{\mathrm{\Δ}t}\end{array}\right]|{|}_2,$

其中，argmin表示最小化，t_k表示第k个时间点，其中，与该时间步长相对应的常数c为c_Δt＝g(T_end)-a_Δtln(T_end)。

15.根据权利要求14所述的计算机程序产品，所述计算机程序产品还包括使所述处理器执行根据下式针对未来时间T_end+t_d确定预测值的操作的指令：

F_Δt(t_d)＝a_Δtln(T_end+t_d)+c_Δt，

其中，t_d是用于预测的期望时间长度。