动态供应链网络中企业合作关系的链路预测

卢志刚，陈 倩

上海海事大学 经济管理学院，上海 201306

随着全球一体化的深入，供应链从单一的线性系统演变成了复杂的供应链网络。演变成的复杂供应链网络，由多个不同类型的企业组成，这些企业为了实现利益最大化，根据外部环境的变化随时不断地完善自身业务，调整与其他企业间的合作关系，使得供应链网络具有动态性[1]。预测动态供应链网络中企业合作关系可以有效地帮助企业储备潜在合作伙伴，解决动态供应链网络中企业节点或链路（企业与企业之间的联系，比如合作关系）出现的突发事故或中断，提高企业的互补性和集成能力。

企业间合作关系的预测可以通过数据包络分析法（DEA）[2]、聚类分析法[3]、层次分析法（AHP）[4]和人工神经网络法（ANN）[5]等来分析获取到的企业节点的属性信息来实现。然而，随着信息量的爆炸性增长，企业节点信息的获取成本越来越高，而且即使获得了企业节点的信息，也无法确保信息的真实性和准确性。相比获取供应链网络企业节点的属性信息，供应链网络结构信息更容易获得，也更加可靠。近年来，基于网络结构的链路预测方法受到了广泛的关注，给本文研究动态供应链网络中企业合作关系提供了新的思路。

链路预测作为数据挖掘领域的研究方向之一，在外部信息不确定的环境下能够有效地利用已知的自身网络结构信息，预测网络中当前没有产生连边的两个节点之间连接的可能性[6]。目前多数链路预测研究都集中在静态的简单网络，即不考虑时间的因素，基于当前节点的属性信息或网络结构对未知链路进行预测，然而，现实世界中的大部分网络都是动态复杂的。

在动态供应链网络链路预测中，网络结构的复杂性和随时间演变的动态性使得链路预测问题变得更加困难，但是它更加贴近真实的网络，所以充分考虑供应链网络的复杂结构信息和其动态性可以最大限度地保证信息的完整性，从而获得更加准确的预测效果。

本文致力于整合供应链网络结构的复杂性和随时间演变的动态性对动态供应链网络中企业合作关系（链路）进行预测，主要工作有：

（1）提出链路预测模型，该模型将动态供应链网络划分为由若干个时间片组成的网络快照，基于投影和时间事件计算潜在合作链路（potential cooperation link，PCL）集合里的每条链路的预测得分。

（2）基于链路预测模型，给出动态链路预测（dynamic link prediction，DLP）算法，该算法精度优于传统的静态链路预测算法，且发现其预测性能可通过改变供应链网络演化过程中网络快照的大小和时间事件变化率来提高。

1 相关工作

链路预测既包含了对网络中实际上存在但尚未被探测到的链路预测，即对静态网络中未知链路的预测;也包含了目前不存在和应该存在或者未来可能会存在的链路预测，即对动态网络中未来链路的预测。本文的主要相关研究包括静态网络中的链路预测和动态网络中的链路预测。

1.1 静态网络链路预测

早期的链路预测方法把重点放在了对复杂网络的静态预测。这些方法根据Lu等可以分为三类[7]：（1）基于局部信息相似性的链路预测[8-9]；（2）基于路径的链路预测[10-12]；（3）基于随机游走的链路预测[13-15]。除此之外，目前常见的一种方法是把静态网络划分为多个二部网络，然后根据二部网络的结构信息分别进行投影。Allali等[16]提出了用加权投影的方法来分析二部网络，通过该方法列出所有不存在的链路，然后将链路权重与预先设定的阈值进行比较，最后筛选出潜在链路。然而，列出所有不存在的链路会花费大量的计算时间，并且很难预先选择合适的阈值。Gao等[17]虽然在投影的基础上仅在候选节点对中寻找可能出现的节点对，弥补了文献[16]的不足，但未对其最终链路预测得分进行二次筛选以获得最优结果。本文旨在充分利用复杂网络结构信息以完善投影分析。

1.2 动态网络链路预测

在动态网络链路预测方面，Pujari等[18]提出了一个有监督的链路预测方法去预测动态网络中的潜在链路。Dunlavy等[19]通过整合基于矩阵和基于张量的方法进行动态链接预测。Vu等[20]在动态网络链路预测中加入了一个连续时间回归模型来考虑网络的动态演变性。Xu等[21]提出一种时序链路预测算法，该算法扩展了用于静态链路预测的潜在空间（或生成模型）算法，以考虑节点潜在位置的时间相关性。Lahiri等[22]首先提出用频繁子图集合序列来表示时间序列，然后估计频繁子图出现的概率，最后使用一种通用的自适应算法以预测未来任意节点处的链路。上述基于机器学习的动态网络链路预测方法尽管比静态网络链路预测方法具有更高的预测精度和准确性，但都没有结合动态网络演化过程中的时间事件。

时间事件是描述动态网络中链路从一个状态（连接或断开）过渡到另一状态的重要指标[23]。鉴于时间事件的重要性，Dong等[24]提出了一种基于时空特征的链接预测方法，用于社交媒体中的时间事件确定。Sett等[25]结合多域拓扑特征和时间事件，提出一种有效的时间感知多关系链路预测（TMLP）特征集，并在特征提取过程中开发时间序列模型，最终用于动态异构网络中的链路预测。Aslan等[26]基于链路上的时间事件提出了一种称为link-score基于路径的方法，该方法根据路径上链路的时间戳显示了路径接近度。Ahmed等[27]提出了一种基于相似性的方法来预测动态单峰网络中的潜在链路，该方法使用特殊的评分因子，以使其方法中新时间事件更具重要性。O’Madadhain等[28]基于目标时间段之前发生的时间事件构成的静态图，创建了不同的度量指标来预测在特定时间段内出现的链路。Ibrahim等[29]通过整合网络中的时间事件、社区结构和节点中心性，提出一种适用于动态网络的链路预测方法，为动态网络中的链路预测生成了一个更加贴近现实的模型。Tylenda等[30]分析了相关节点对的发展，并提出了一个包含链路权重的新网络模型，该模型中的链路权重代表节点之间发生的时间事件的年龄。然而上述方法都是基于网络的动态演化性提出的，并没有结合当前状态下的复杂网络结构。本文旨在整合网络结构的复杂性和随时间演变的动态性进行链路预测。

2 问题描述

动态供应链网络中企业合作关系的链路预测可以描述为：在动态供应链网络Gt=(S,M,Et)中，根据在t1时刻观察到的现有合作链路，预测tn时刻的未来合作链路。如图1所示，其中S={s1,s2,…,si}是彼此不相交的供应商节点集合，M={m1,m2,…,mj}是彼此不相交的制造商节点集合，Et是S和M之间随时间t变化的合作链路集合，实线代表已知的企业合作链路，虚线代表预测到的企业合作链路。

图1 动态供应链网络中企业合作关系的链路预测Fig.1 Link prediction of enterprise cooperation relationship in dynamic supply chain network

3 链路预测模型

为了在动态供应链网络中准确地预测企业合作关系，提出链路预测模型，该模型由三个阶段组成。首先，根据时间将动态供应链网络划分为一系列网络快照。其次，利用投影挖掘每一个网络快照中供应链网络的结构信息，找到供应链网络中的潜在合作链路（PCL）集合。最后，基于时间事件计算出PCL集合中每个元素的链路预测得分，比较得分，确定动态供应链网络中企业间未来最可能出现的合作关系。模型中所用主要参数如表1所示。

表1 模型中主要参数说明Table 1 Description of main parameters in this model

3.1 动态划分

为了充分考虑供应链网络的动态性，把t1~tn时间内可以观察到的动态供应链网络划分为由若干个时间片组成的网络快照，现有相关定义如下。

定义1（网络快照）在供应链网络G=(S,M,E)中，Gt1=(S,M,Et1)是在t1时刻记录的网络快照，Gtn=(S,M,Etn)则是在tn时刻记录下的网络快照，其单位是帧。

定义2（演化结构）记n为供应链网络中网络快照的数量，w是预测的时间范围，则以帧记录的演化结构定义为：

3.2 投影分析

为了挖掘动态供应链网络的结构信息，本模型采用投影的方法，具体定义如下。

定义3（投影）将t时刻的供应链网络快照Gt=(S,M,Et)投影为两个单模图，分别是供应商单模图Gts=(S,Ets)和制造商单模图Gtm=(M,Etm)，其中：

Γ(si)和Γ(sj)分别是si和sj在供应链网络快照中的合作伙伴集合。从上述定义可以看出，若节点si和sj在供应链网络快照中至少有一个共同制造商伙伴，投影后的供应商单模图将会增加一条链路，制造商单模图同理。如图2所示，供应链网络快照图2（a）经过投影后得到了供应商单模图2（b）和制造商单模图2（c）。

图2 供应链网络快照的投影Fig.2 Projection of supply chain network snapshot

基于投影后的单模图，给出潜在合作链路（PCL）的定义来过滤掉其他不存在的链路。

定义4（潜在合作链路）在t时刻的供应链网络快照Gt=(S,M,Et)添加一条链路(si,mu)∈S×M，构造出一个新的供应链网络快照，其中和分别是投影后的供应商单模图和制造商单模图，若，(si,mu)就被定义为t时刻供应链网络快照中的潜在合作链路（PCL），以供应链网络快照记录的演化结构为单位，确定潜在合作链路集合。

为了避免花费大量的时间考虑所有不存在的链路，除了潜在合作链路，t时刻供应链网络快照中其他不存在的链路不予考虑，这样一来，缩小了查找范围，大大减少了计算量，降低了链路预测的复杂度。同时，潜在合作链路的加入对原t时刻的供应链网络快照没有影响，确保了供应链网络快照中链路信息的完整性。

3.3 得分计算

为了对PCL集合里的每一条潜在合作链路进行评分，以表示动态供应链网络中该潜在合作链路出现的可能性，现给出如下定义。

定义5（潜在合作链路覆盖的模）假定(si,mu)是t时刻供应链网络快照中的潜在合作链路。对于任意的供应商节点sj，当满足时，(si,sj)∈Ets被定义为潜在合作链路(si,mu)覆盖的模，其本质是供应商单模图Gts中的一条链。其中Γs(si)是供应商节点si在供应商单模图Gts中的邻居集合，Γ(mu)是制造商节点mu在供应链网络快照中的合作伙伴集合。

在动态供应链网络中，模权重在从一帧过渡到其下一帧的过程中可能会增加或减少，模权重的描述如下。

定义6（模权重）假定(si,mu)是t时刻供应链网络快照中的潜在合作链路，(si,mu)覆盖的模(si,sj)∈Ets的权重即模权重（pattern weight）定义为：

其中，d(si)是供应商节点si的度，d(sj)是供应商节点sj的度。Γ(si)和Γ(sj)是节点si和sj在供应链网络快照中的合作伙伴集合，mv是他们的共同合作伙伴，d(mv)是节点si和sj共同合作的制造商节点mv的度。

根据模权重的变化，可以把供应链网络动态演化划分为连续时间事件、强化时间事件和弱化时间事件，具体描述如下。

（1）连续时间事件

当在第()ω-1帧网络快照中模权重增加或减少的幅度不超过第ω帧网络快照中的预定时间事件变化率时，连续时间事件（consistent temporal events）就会发生。对于模(si,sj)∈Ets的连续时间事件分值Cw(si,sj)定义为：

其中，PWw-1(si,sj)和PWw(si,sj)分别是网络快照Gw-1和Gw中节点si和sj的模权重。时间事件变化率p(0

（2）强化时间事件

当在第(ω-1)帧网络快照中的模权重增加的幅度超过第ω帧网络快照中的预定时间事件变化率时，就会发生强化时间事件（reinforced temporal events）。在处理从第(ω-1)帧到第ω帧的过渡期间的强化时间事件时，模(si,sj)∈Ets的强化时间事件分值Rw(si,sj)定义如下：

因为在强化时间事件中从一帧网络快照过渡到另一帧网络快照的过程中模权重会提高，所以常数r是介于(0

（3）弱化时间事件

当在第(ω-1)帧网络快照中的模权重减少的幅度超过第ω帧网络快照中的预定时间事件变化率时，则会发生弱化时间事件（descended temporal events）。在处理从第(ω-1)帧到第ω帧的过渡期间的弱化时间事件时，模(si,sj)∈Ets的弱化时间事件分值Dw(si,sj)定义如下：

因为在弱化时间事件中从一帧网络快照到另一帧网络快照的过渡中模权重会减小，所以常数d是介于(-1≤d<0)之间的负惩罚值。

时间事件可以用于分析动态供应链网络中每条潜在合作链路的演化发展，对于模(si,sj)∈Ets，可以将演化分为主演化、次演化和总演化，具体定义如下。

定义7（主演化）主演化（primary evolution）是指模(si,sj)∈Ets从第(ω-1)帧网络快照过渡到第ω帧网络快照中的演化，这一演化可以是连续时间事件，强化时间事件或弱化时间事件，P(si,sj,w)用于计算主演化得分。

定义8（次演化）次演化（secondary evolution）是指模(si,sj)∈Ets中节点si、sj与其合作伙伴连接的链路从第(ω-1)帧网络快照过渡到第ω帧网络快照中的演化。S(si,sj,w)用于计算次演化得分，其中Γ(si)和Γ(sj)分别是节点si和sj在供应链网络快照中的合作伙伴集合，mv是节点si和sj的共同合作伙伴。

定义9（总演化）总演化（total evolution）是指所有与模(si,sj)∈Ets相关演化的总和。T(si,sj)用于计算总演化得分，其中φ是一个预设的影响因子，它表示从第(ω-1)帧网络快照过渡到第ω帧网络快照中次演化对总演化的影响程度。

定义10（链路预测得分）假定(si,mu)是PCL集合里的潜在合作链路，(si,mu)覆盖的模(si,sj)∈Ets在每个过渡状态下的总演化之和是该潜在合作链路的最终链路预测得分（link prediction score）。

其中，λ(si,mu)是(si,mu)覆盖的模的集合。通过比较PCL集合里每个元素的链路预测得分，确定得分高于平均值的即为预测结果。

4 动态链路预测算法

基于所提模型，本章给出动态链路预测（DLP）算法用于预测动态供应链网络中的企业间合作关系。算法的框架描述如下：

输入：供应链网络G=(S,M,E)，用于训练的供应链网络：Gtrain=(S,M,Etrain)，预测时间间隔w，网络快照的个数n,PCL集合，以帧记录的Gtrain演化结构：

输出：PCL集合中每一个元素的最终链路预测得分

1./*构造供应链网络投影后的单模图*/

2./*计算潜在合作链路覆盖的模的权重*/

3./*计算潜在合作链路的最终链路预测得分*/

5 实验结果与分析

为了验证DLP算法的预测性能，本文从同花顺iFinD行业数据库采集数据，构建了2002年到2018年消费电子供应链网络Gt=(S,M,Et)，并用Matlab R2018a作为实验工具完成整个代码的编写，最后用Pajek实现数据可视化。

5.1 实验数据

实验以同花顺iFind智能投研给出的消费电子图谱为依据，从同花顺iFind金融数据终端采集整理了2002年到2018年消费电子供应商与制造商的合作数据。其中以[2002,2004)年的供应链网络快照为例，如图3所示，该网络快照由134个节点，941条链路构成，其中供应商节点个数为125，制造商节点个数为9。黄色节点代表供应商，绿色节点代表制造商，蓝色链路代表彼此之间的合作，其中9个制造商是由生产苹果、华为和小米品牌的智能穿戴、智能手机和平板电脑的企业组成，125个供应商是由提供指纹识别模组、摄像头模组、摄像头芯片、液晶面板模组、镜头、玻璃盖板、驱动IC、VCM马达、触控芯片和触摸屏的公司构成。

图3 2002—2004年的供应链网络快照Fig.3 Supply chain network snapshot of 2002—2004

5.2 实验说明

为了评估所给算法，本文将2002年到2018的消费电子供应链网络划分成两种网络快照大小不同的网络，并分别进行30次实验。在第一种供应链网络中，选取[2002,2016]的供应链网络作为训练网络Gtrain=(S,M,Etrain),[2016,2018]的供应链网络作为测试网络Gtest=(S,M,Etest)。在第二种供应链网络中，选取[2002,2014]的供应链网络作为训练网络Gtrain=(S,M,Etrain)，[2014,2018]的供应链网络作为测试网络Gtest=(S,M,Etest)。

在第一种供应链网络中，根据提出的DLP算法，将用于训练的供应链网络划分为n=7张网络快照，其中时间间隔w=2。A1={F1→F2},A2={F2→F3},A3={F3→F4},A4={F4→F5},A5={F5→F6}和A6={F7→F8}分别是供应链网络演化过程的6个过渡演化结构，F8是供应链网络的预测帧。第一种供应链网络划分后的8帧具体如下所示：

在第二种供应链网络中，将用于训练的供应链网络划分为n=3张网络快照，其中时间间隔w=4。A1={F1→F2}和A2={F2→F3}是供应链网络演化过程的2个过渡演化结构，F4是供应链网络的预测帧。第二种供应链网络划分后的4帧具体如下所示：

5.3 性能评估

在上述两种划分不同的供应链网络上，通过利用不同的时间事件变化率p(0

表2 传统的链路预测算法指标Table 2 Traditional link prediction algorithms measure

表2中kx和ky分别是供应链网络G=()S,M,E中节点x和y的度。和分别表示为矩阵l+行中第x行第x列、第y行第y列、第x行第y列的位置所对应的元素。Pxy和Pyx分别是一个粒子从节点sx出发走到节点my的概率和一个粒子从节点my出发走到节点sx的概率。Pxy(t)和Pyx(t)分别是一个粒子t时刻从节点sx出发，t+1时刻这个粒子正好走到节点my的概率和一个粒子t时刻从节点my出发，t+1时刻这个粒子正好走到节点sx的概率，qx和qy分别是节点sx和my的初始资源分布。在LRW的基础上，将t步及其以前的结果加总便得到SRW的值（PA：preferential attachment，ACT：average commute time，RWR：random walk with restart，LRW：local random walk，SRW：superposed random walk）。

5.4 实验结果

实验通过蛮力法评估和确定了连续时间事件分值，强化时间事件分值，弱化时间事件分值和影响因子的最佳值，分别为c=0.1，r=0.5，d=-0.1和φ=0.09。在划分后的两种供应链训练网络中，分别基于时间事件变化率p=40%和p=60%计算出DLP算法的精度。最后将该算法的精度与基于PA、ACT、RWR、LRW、SRW指标的链路预测算法分别进行对比，其中RWR的参数分别取0.85和0.95，LRW的步数分别取3和5，SRW的步数分别取3、4和5。

在第一种供应链网络中，分别基于时间事件变化率p=40%和p=60%计算出30次实验中DLP算法的精度，并将其平均值与基于传统链路预测指标的算法进行了比较。如图4所示，当供应链网络被划分成7张网络快照时，此时每一个网络快照的大小较小，DLP算法的精度优于其他静态链路预测算法。其中当时间事件变化率p=40%时，DLP算法的精度最高，当时间事件变化率p=60%时，DLP算法的精度排在第二，即当划分的供应链网络快照大小较小且时间事件变化率较低时，所提算法可以实现最佳的预测性能。这是因为当网络快照的大小较小时，动态供应链网络演化中的时间事件变化是由时间事件变化率来确定的。一方面，在事件变化率较低的情况下，观察到的事件变化更多。另一方面，在事件变化率较高的情况下，观察到的事件变化较少。上述情况直接影响到DLP算法的预测性能。

图4 第一种供应链网络中各链路预测算法精度Fig.4 Precision of link prediction algorithms in supply chain network 1

在第二种供应链网络中，如图5所示，当供应链网络被划分成3张网络快照时，此时每一个网络快照的大小较大，具有时间事件变化率p=60%的DLP算法的精度优于具有时间事件变化率p=40%的DLP算法和其他静态链路预测算法，即当划分的供应链网络快照大小较大且时间事件变化率较高时，所提算法可以实现最佳的预测性能。这是因为在网络快照大小较大的情况下，在事件变化率较高的网络中会观察到更多的事件变化。但是，在事件变化率较低的网络中，观察到的事件变化较少。

图5 第二种供应链网络中各链路预测算法精度Fig.5 Precision of link prediction algorithms in supply chain network 2

在上述两种供应链网络中，如图4和图5所示，DLP算法在精度方面比其他基于传统链路预测指标的算法表现出更好的性能，这是因为所提算法在计算潜在合作链路覆盖的模阶段考虑了供应链网络动态演化中三个关键的结构特征：节点的共同合作伙伴数、网络中节点的度以及共同合作伙伴节点的度。基于PA指标的链路预测算法仅仅考虑了当前网络中节点度的特性，因此具有较低的预测精度。基于ACT和RWR指标的链路预测算法只从概率的角度考虑了网络结构对链路形成的作用，完全忽略了网络的结构特征，其预测精度低于基于PA指标的链路预测算法。基于LRW和SRW指标的链路预测算法虽然和基于ACT和RWR指标的链路预测算法都隶属于基于随机游走相似性算法，但是由于前者考虑了节点的初始资源分布情况，所以它们的预测精度高于基于ACT和RWR指标的链路预测算法。

表3分别列举出了时间事件变化率p=40%和p=60%的DLP算法和基于PA、ACT、RWR（参数分别取0.85和0.95）、LRW（步数分别取3和5）、SRW（步数分别取3、4和5）指标的链路预测算法在选取[ ]2002,2016的供应链网络作为训练网络且被划分为7张网络快照，时间间隔为2的供应链网络1和选取[ ]2002,2014的供应链网络作为训练网络且被划分为3张网络快照，时间间隔为4的供应链网络2中的实验结果，通过对比发现当供应链训练网络被划分成7张网络快照时，时间事件变化率p=40%的DLP算法发现的潜在合作链路个数最多。如表3所示，基于PA、ACT、RWR、SRW指标的链路预测算法由于只考虑当前的供应链网络状态，忽略其动态演化性，所以它们在两种按时间划分具有不同网络快照大小的动态供应链网络中所发现的潜在合作链路个数是相同的。DLP算法之所以发现的潜在合作链路个数高于其他静态链路预测算法，是因为其不仅考虑了潜在合作链路覆盖的模的演化还考虑了其相关合作伙伴的演化，这使得动态供应链网络拓扑结构中的信息损失最小化。在供应链网络的动态演化中隐藏着很多有价值的信息，这些隐藏信息的背后包含着有关未来网络的重要提示，DLP算法克服了上述传统算法的限制，表现出了最佳的预测性能。为了使结果可视化，给出了时间事件变化率p=40%的DLP算法在第一种供应链网络中预测到的2016年—2018年的供应链网络快照图6。与2002年—2004年的供应链网络快照图3相比，图6新增了124条合作链路，节点与节点之间的合作链路更加复杂密集。

图6 2016—2018年供应链网络快照Fig.6 Supply chain network snapshot of 2016—2018

表3 实验结果算法对比Table 3 Algorithm comparison of experimental results

6 结束语

为了预测动态供应链网络中的企业合作关系，本文将动态供应链网络划分成一系列网络快照，提出了一种基于投影和时间事件的链路预测模型，并给出了动态链路预测算法。实验结果表明该算法与传统链路预测算法相比，在预测动态供应链网络中两个企业未来是否合作方面具有更好的预测效果，且发现其预测精度可以通过改变网络快照的大小和时间事件变化率来提高。真实世界的大多数网络不同于本文中的供应链网络，它们的节点种类繁多，链路有权有向，未来研究可以进一步优化本文算法，使其在更加复杂的大型动态网络中实现准确的预测。