随机梯度估计及其在系统管理中的应用

张公伯 ,李海东 ,彭一杰

(北京大学a.光华管理学院;b.工学院,北京 100871)

随着信息技术的快速发展,系统管理的技术方法正在发生深刻变革。系统仿真研究是一个经典的、传统的领域,其优点在于可以对一般的复杂系统做高置信度建模,受到学界与业界的广泛关注。随着我国产业升级不断推进,国际竞争程度日益加深,企业对系统管理精益化程度的要求不断提高,仿真技术在系统管理中应用的重要性日益显现。大数据、工业物联网、云计算、人工智能和数字孪生技术等作为实践系统管理的新兴技术在仿真分析中起着重要作用。复杂系统通常具有动态性和随机性,准确估计系统的平均表现需要做大量仿真,且仿真耗费的计算时间长。复杂随机系统仿真中通常包含影响系统结构的控制参数与不连续样本表现,导致系统表现的优化分析困难。因此,需要从理论与方法论层面研究如何将仿真与优化有效结合,应用于复杂随机系统的管理。有关仿真优化方法可参考文献[1-7]。

对于具有连续可控参数的优化问题,梯度蕴含了目标函数随参数变化的信息,在迭代搜索优化方法中有着重要的应用,使得结果逐步逼近局部或全局最优。复杂随机系统优化问题的目标函数的解析表达通常未知,只能通过蒙特卡洛仿真得到的系统样本表现进行估计。由于蒙特卡洛仿真存在噪声,因而观测到的系统样本表现具有随机性。通过蒙特卡洛仿真估计的随机系统中目标函数关于参数的梯度称为随机梯度估计。随机梯度估计在仿真分析和优化问题中有着广泛的应用,如敏感性分析、仿真模型的不确定性量化和随机优化及机器学习中的梯度下降算法。

随机梯度估计可应用于排队系统、库存管理系统、统计质量控制[8-14]、机器维护系统[15-18]和随机活动网络[19-21]等系统的优化和调控决策。现实中的系统通常结构复杂并具有随机性,例如,生产过程中产品订单的到达、机器加工时间及是否发生故障都是随机的。这些复杂系统通常不具有解析形式的输出分布,但具有部分结构信息可以分析,被称为灰箱系统。恰当地利用这部分结构信息可以在仿真估计系统表现的同时得到其梯度信息,从而实现求解优化问题时高效搜索的目标。例如,质量控制图中希望在真实生产环境下运行系统的同时实现控制边界的在线优化设计,但质量控制图的控制边界属于结构参数,它的表现关于控制边界不连续,使得质量控制图的结构分析非常困难;生产系统中通过排队系统对产品的加工流程进行建模,系统中通常存在重加工、批量加工和优先处理等结构,因而使得排队网络系统的结构分析非常困难。

采用随机梯度估计方法对分位数敏感性进行计算和分析,也是仿真和复杂系统管理领域关注的重要问题。相较于系统的期望表现,分位数能够提供更多随机系统表现尾部分布的特征,从而更好地刻画系统在极端环境下的表现,提高极端事件下系统的稳健性。例如,生产订单在排队系统中的等待时间过长就可能失效,因此,等待时间的分位数过大将有损订单的数量,降低企业的效益。

随机梯度估计方法主要分为间接和直接估计。间接梯度估计通常具有两个特征:①它只估计真实梯度的一个近似值,如在标量情况下通过割线近似;②它仅使用来自原始(未修改)系统的函数评估(绩效指标输出样本)[22-23]。直接梯度估计利用仿真模型中随机变量的分布信息或样本表现函数的结构信息(如输入的概率分布或系统表现与样本路径的函数关系)给出真实梯度的无偏梯度估计量。

间接梯度估计方法主要有有限差分法(Finite Differences,FD)和同时扰动法(Simultaneous Perturbation,SP)[22-23]两种。这两种方法均不需要了解仿真模型随机性的产生或样本表现函数,即仿真模型可以被视为黑箱。有限差分法通过仿真两个临近输入点的函数值并用微积分基本原理来近似导数,由此得到的估计量既有导数近似带来的偏差,也有函数评估噪声带来的方差。对于高维问题,有限差分法需要独立地估计每一个维度,因此,估计一次梯度所需要的仿真样本的数量与维度相关。同时扰动法的优点是一次仿真同时扰动自变量的多个维度,因此,估计一次梯度所需要的仿真样本的数量与自变量的维度无关。

尽管间接的梯度估计具有通用性,但直接的梯度估计具有如下优点:①它得到的梯度的估计量是无偏的,以无偏的梯度估计量作为迭代方向可以加快仿真优化算法(如随机逼近算法[24])的收敛速度;②它避免了间接的梯度估计方法中确定合适的扰动大小的需要,而扰动大小的选择直接影响梯度估计量的准确性;③所得的估计量在计算上更高效,在高维问题中直接的梯度估计方法达到给定的估计精度所需要的仿真样本的数量通常要比间接的梯度估计方法更少[22-23]。常用的无偏随机梯度估计方法有无穷小扰动分析法(Infinitesimal Perturbation Analysis,IPA)、似然比法(Likelihood Ratio,LR)(也称为分数函数法(Score Function,SF))、弱导数法(Weak Derivative,WD)(也称为测度求导法(Measured-Value Differentiation,MVD))和广义似然比法(Generalized Likelihood Ratio,GLR)。无穷小扰动分析方法可以用于解决带有结构参数(即出现在样本路径与系统表现度量的函数关系中的参数)和分布参数(即出现在输入随机变量的概率分布中的参数)的问题,但要求样本表现具有连续性;似然比方法虽然能够处理不连续的样本表现,但仅适用于系统表现度量关于分布参数的梯度估计而不适用于结构参数;弱导数方法能降低似然比方法的方差,但需要仿真额外的样本轨道。Gong 等[25]提出了平滑扰动分析法(Smoothed Perturbation Analysis,SPA),对无穷小扰动分析法进行了拓展,通过取条件期望来解决样本表现不连续的问题,但所选取的条件通常和问题本身相关,而且条件期望估计量的计算可能需要函数取逆和额外的仿真。Rubinstein等[19]提出了推出式似然比法(Pushout LR)来解决同时带有结构参数和不连续的样本表现的问题,但是一般需要解析的函数转换将结构参数从样本表现中移出。Peng等[11]提出的广义似然比方法既适用于带有结构参数和不连续的样本表现的问题,同时具有无偏性和单样本轨道性,也不需要取条件期望和函数转换。

系统管理领域信息化和数字化的转型,对基于梯度的仿真优化方法的设计提出了挑战。复杂系统的优化设计需充分考虑问题情境特点及其中的模型假设,以设计简洁有效的基于梯度的仿真优化方法。本文将首先介绍常用的随机梯度估计方法,进而考虑4个随机系统的梯度估计问题,给出其在系统管理中的应用场景。

1 随机梯度估计方法

其中,当Y是连续随机变量时,α=P（Y（θ）≤qα（θ））。随机梯度估计问题研究的是通过对随机模型进行仿真,得到系统表现度量关于参数的导数的估计。随机模型中的随机梯度估计方法可以用于仿真优化的输入和对参数进行敏感性分析。仿真优化中常用的随机逼近方法的参数迭代更新过程就采用梯度估计作为其迭代方向,第(n+1)步迭代的计算式为

1.1 间接的梯度估计方法

1.1.1有限差分法 有限差分法不需要获取复杂系统的信息,仅通过分析系统仿真得到的输出样本,就可得到有偏的梯度估计量。有限差分法将复杂系统视为“黑箱”,分别对参数向量θ中每一个元素进行摄动,可得到系统期望表现梯度估计量[22-23]。单边前向有限差分估计量是有限差分方法中常用的估计量,其第i个元素为:[Y（θ+ci ei）-Y（θ）]/ci。其中:ci为第i个方向的摄动值;ei为第i个方向的单位向量;有限差分法中另一种常用估计量是双边对称有限差分估计量,其第i个元素为

当参数向量的维度很高时,有限差分法需要的计算量很大。此外,有限差分法中需要选取摄动值的大小,如果摄动值的选取过小,得到的系统表现梯度估计量的偏差小、方差大。因此,有限差分法摄动值的选取需要平衡表现梯度估计量的偏差与方差。在实际应用中,双边对称有限差分估计量比单边前向有限差分估计量的精度高,但在每一次梯度估计时,双边对称有限差分估计量需要消耗更多的仿真资源(前者需要2d个仿真资源,而后者仅需要(d+1)个仿真资源)。对于二阶连续可微函数,双边对称有限差分估计量的均方误差的最优收敛阶为n-2/3,所对应的摄动值的阶为n-1/6,其中n为仿真资源的数量。得到最优摄动值(即包括n-1/6前的常数)对于小样本情况下双边对称有限差分估计量的表现十分重要,但往往非常困难。

由于不需要获取模型的信息,故有限差分法普遍适用于复杂系统表现函数的梯度估计。但采用有限差分法估计梯度需要多次重复仿真,运算成本大,因此不适用于高维参数问题。在基于梯度的优化方法中,有限差分法比直接的梯度估计法的表现差。

1.1.2同时扰动法 同时扰动法由Spall[26]提出。同时扰动估计量消耗的仿真资源与参数向量的维度无关,适用于求解高维参数问题,它的第i个元素为:［Y（θ+CΔ）-Y（θ-CΔ）］/（2ciΔi）。其中:Δ=(Δ1,Δ2,…,Δd)是d-维扰动向量,通常假设各元素Δi独立;C是对角线为摄动值ci的对角矩阵。在每一次梯度估计时,同时扰动法需要生成d个随机数来构造扰动向量Δ。与双边对称有限差分估计式(1)相比,同时扰动估计量的每个元素的分子都相同,只需要两个样本表现函数扰动估计,与参数向量的维度d无关。由于生成样本表现函数扰动估计的运算成本大,因而同时扰动法比有限差分法的计算效率高。为保证同时扰动法几乎处处收敛,要求扰动向量中每一项独立,均值为0且逆二阶矩有限。均匀分布和正态分布不满足上述条件,常用的Δi的分布是成功概率为1/2的伯努利分布。类似于拟蒙特卡洛方法的核心思想,文献[27-28]中提出了采用确定的扰动向量代替随机扰动向量Δ的方法。

文献[29]中介绍了采用同时扰动法作为迭代方向的随机逼近方法。同时,扰动法求解高维参数问题的优越性使其在排队系统[30]、交通系统[31-35]和统计质量控制[36]等系统的优化和调控决策中有广泛的应用。Fu等[30]采用同时扰动法估计排队系统的系统停留时间关于服务台服务速率的梯度;朱承元等[34]优化设计珠三角地区多机场系统的航班时刻,采用同时扰动法估计系统中航班总延误时间的梯度,降低了航班总延误成本;齐心[35]优化设计城市交叉口排队系统的信号灯时长,采用同时扰动法估计车辆在交叉口的通行时间的梯度,提升了城市交叉口的通行能力。

1.2 直接的梯度估计方法

1.2.1扰动分析法 扰动分析法(Perturbation Analysis,PA)是一种通过分析单个样本路径对随机模型的期望表现梯度进行估计的方法。文献[37-38]中最早对该方法进行了研究。这种方法通过样本路径对表现度量的影响以及系统参数对样本路径产生的扰动来衡量系统表现梯度。扰动分析法中最常用的是无穷小扰动分析法,该方法通过样本表现关于参数的导数来估计系统期望表现梯度。系统期望表现为

如果期望与导数是可以交换的,则系统期望表现的梯度为

无穷小扰动分析法可以解决样本表现带有结构参数的随机梯度估计问题,但不能解决样本表现不连续的问题。为克服这一理论缺陷,文献[25,42-43]中提出利用取条件期望处理样本表现不连续的平滑扰动分析法,但该方法需要解析地计算一个条件期望。在实际应用中,该方法需要视实际问题情况而选取适当的条件。Hong[44]首次提出利用无穷小扰动分析法解决分位数敏感性估计问题,但得到的梯度估计量是渐近无偏,且需要对数据进行批处理才可以保证估计量的弱相合性。Jiang等[45]在文献[44]的基础上提出一种不需要进行批处理的无穷小扰动分析估计量,但该估计量仅在特定的条件下具有强相合性。有关无穷小扰动分析法和平滑扰动分析法的详细介绍,可以参考文献[46]。在复杂系统的优化设计中,基于扰动分析法的随机逼近方法可以实现生产参数的估计。排队系统是对生产工艺流程建模的一种常用方法,可以刻画生产订单到达和装配生产线加工等生产流程。优化工艺流程的资源配置(如订单到达速率、生产线机器数量和缓冲区数量等)可以提高生产效率、降低生产成本。Cheng[47]考虑多阶段存货式生产系统,采用无穷小扰动分析法估计样本路径的梯度,实现减少订单等待时间和提高生产加工效率的权衡。Fu等[48]采用平滑扰动分析法估计多服务台单队列的排队系统的稳态系统时间的梯度。有关扰动分析法在排队系统优化设计中的更多应用,可以参考文献[49-53]。企业需要管理生产原料和成品的库存,根据客户的需求确定补充库存的时间和订货数量。优化设计企业的库存管理策略,可以避免企业超储或缺货,实现降低库存成本的目标。Fu[54]考虑带有随机客户需求和延迟交货的(s,S)策略库存系统,应用扰动分析法估计表现函数的样本路径梯度。Bashyam 等[55]采用扰动分析法估计每周期的期望费用关于参数s和S的梯度。扰动分析法在其他复杂系统管理中也有广泛的应用。Heidergott[16]采用平滑扰动分析法估计维护系统修复时间关于系统阈值的梯度,优化设计维护系统的修复时间。Donohue等[56]应用扰动分析法得到收益最大化的生产线产能配置方案。Yan等[57]基于扰动分析法提出了一种求解近似最优阈值的算法,用于控制带有两台串联机器的不可靠柔性制造系统。Liberopoulos等[58]采用无穷小扰动分析法估计生产流的一阶和二阶梯度,用于生产流的设计和控制。Brooks等[59]采用无穷小扰动分析法估计异步传输网络最小平均延迟的渐近无偏梯度。何宁等[60]采用扰动分析法估计卫星网络信息流的延时关于带宽的梯度,实现最优的信息流带宽分配。

1.2.2似然比法 似然比法,也称为分数函数法。该方法将系统期望表现对参数的依赖全部转移至输入分布的密度函数中,即参数θ是输入分布的密度函数中的参数(分布参数)。似然比法最早在文献[61]中进行讨论,随后,文献[62-65]也独立地发现并对该方法进行了研究。系统期望表现为

如果期望与导数可以交换,则系统期望表现的梯度为

称h(x)为似然比估计量。似然比法可以解决样本表现不连续的问题,但不能解决样本表现带有结构参数的问题。文献[46,66]中声称在同一模型表示下,无穷小扰动分析法比似然比法方差小,该结论在仿真分析中被视为经验准则。Cui等[67]给出了该经验准则成立的充分条件,也给出了该经验准则不成立的反例。

为克服似然比方法无法解决样本表现带有结构参数问题的方法局限,文献[19]中提出了推出似然比法(Push-out LR),该方法利用变量替换将不连续的样本表现中的结构参数推出至输入分布,再采用似然比法;但采用变量替换推出的似然比法依赖于解析形式的变量替换,在实际应用中这种方法需要视实际问题情况而定,这给该方法的理论分析和实际应用带来了困难。与无穷小扰动分析法相比,似然比法的无偏性条件通常更容易满足。L’Ecuyer[68]在样本表现带有结构参数且连续的情况下,采用导数中的乘积法则,将无穷小扰动分析法与似然比法相结合,提出了IPA-LR 估计量。该估计量具有单样本轨道的性质,但在处理样本表现不连续的问题时是有偏的。

在复杂系统的优化设计中,似然比法不适用于估计系统表现关于系统结构参数的梯度。例如,在生产订单随机到达的(s,S)策略库存系统中,似然比法不适用于估计库存量关于参数s和S的梯度,但可以估计库存量关于订单随机到达时间分布中的参数的梯度。Fu[1]采用似然比法估计了GI/G/1排队系统停留时间关于随机服务时间分布中的参数的梯度。Fu等[8]采用似然比法估计质量控制图的平均运行长度关于控制边界的梯度,用于控制图控制边界的敏感性分析,更加准确地判定系统是否处于受控状态。

1.2.3弱导数法 弱导数法,也称为测度求导法。该方法最早由Pflug[69-70]进行讨论。与似然比法相似,系统的期望表现梯度可以表示为

称c（θ）（h（X(+)（θ））-h（X(-)（θ）））为弱导数估计量,其中,c（θ）是常数,随机变量X(+)（θ）的概率密度函数为,随机变量X(-)（θ）的概率密度函数为由于概率密度(或概率质量)函数的分解不是唯一的,因而弱导数估计量不是唯一的。

与无穷小扰动分析法和似然比法相比,弱导数法的无偏性条件通常更容易满足,但要仿真额外的样本轨道。当参数向量的维度很高时,弱导数法需要的计算成本大。Heidergott等[71]采用弱导数法研究了分位数敏感性估计问题,解决了文献[44]中估计量无法处理样本表现不连续的问题,并将其应用于排队系统中,优化了设计考虑系统停留时间分位数的系统服务速率。

在复杂系统的优化设计中,弱导数法不适用于估计系统表现关于系统结构参数的梯度,且不适用于高维生产参数的估计问题,因而有关弱导数法在复杂系统中的应用的文献较少。Heidergott[17]考虑带有老化更换策略的多组件维护系统的优化设计问题,应用弱导数法估计维护成本关于阈值参数的梯度。

1.2.4广义似然比估计法 复杂随机系统中经常存在不连续样本表现与结构参数,例如,质量管理中控制图的控制边界属于结构参数,它的表现函数关于控制边界不连续。存在不连续样本表现与结构参数时,经典梯度估计方法的无偏性无法得到保障,如何系统地解决该问题是仿真领域的重要难题之一,很多学者对这类问题进行了研究。Wang等[72]提出一种混合无穷小扰动分析与似然比的方法,但该方法处理不连续样本表现的手段非常依赖于问题的特定结构,适用范围有限(主要在金融领域)。Liu等[73]提出了一种核函数方法来处理不连续样本表现,但该方法得到的随机梯度估计是有偏的。文献[71,74-75]也为解决该问题做了大量研究。这些方法在应用时需要视问题的特定情况分析,并且通常不具有单样本轨道。

Peng等[11]首次提出了广义似然比估计法,该方法是目前仿真优化领域解决复杂随机系统包含结构参数和不连续样本表现的单样本轨道无偏随机梯度估计问题中适应性最强的方法之一。扰动分析法与似然比法均可以放到广义似然比法的理论框架中分析。“广义”的含义是当样本表现没有结构参数时,广义似然比估计量可以退化为似然比估计量。系统的期望表现梯度可以表示为

称φ（h（X;θ））w（X;θ）为广义似然比估计量,其中,φ是几乎处处连续的可测函数,

d（x;θ）是与函数h、函数h的雅可比矩阵和样本表现函数有关的量。广义似然比估计法不但适用于绝大多数带有结构参数的问题,而且适用于现有方法通常很难处理的不连续且非线性的样本表现问题。该方法同时解决了无穷小扰动分析法不能处理不连续样本表现,与似然比法不能处理结构参数的问题。该方法为无偏随机梯度估计提出了新的估计方法或途径,提供了扰动分析法与似然比法由于样本表现对于结构参数不连续导致的偏差的解析表示,并提出了可积性条件用来补偿光滑性条件从而得到无偏梯度估计的方法创新。

广义似然比法在复杂系统的优化设计中具有广泛的应用前景。广义似然比法可以解决排队系统分析、统计过程控制、供应链管理、制造工厂和服务系统等一系列梯度估计问题。Peng等[11-12,76]讨论了广义似然比法在质量控制图、随机活动网络、维护系统和具有积压的单一产品、离散时间及多期库存模型中的应用。Zhang等[13]和张公伯等[14]提出了一种质量控制图的经济优化设计方法,采用广义似然比法估计平均运行成本关于控制图控制边界的梯度。Zhou等[77]考虑在给定当前时刻信息的条件下,随机系统在未来的预期性能,采用广义似然比法估计了样本表现条件期望的梯度。Peng等[78]采用广义似然比法研究了复杂随机模型的分布敏感性,导出了仿真模型输出分布函数关于参数与变量各阶导数的估计,可以解决生产服务与管理等复杂系统数据驱动的随机建模,减小模型结构假设不准确导致的误导性输出结果。Glynn 等[79]采用广义似然比法研究了基于分布敏感性的分位数敏感性估计。但在采用广义似然比法进行分位数敏感性估计的问题中,仿真样本有限时估计量存在偏差,因而在梯度下降算法中,计算估计量需要的仿真样本数量随着迭代次数增加,从而消耗大量的仿真资源。如何在分位数梯度优化中逐步消除分位数敏感性估计量的偏差值得进一步研究。

表1比较了几种梯度估计方法及其在复杂系统管理中的应用场景。

表1 随机梯度估计方法的比较

2 系统管理中的应用

随机梯度估计在复杂系统管理中有着广泛的应用。复杂系统的输入和输出的关系未知,输入通常受到测量误差、信息缺失等不确定性因素的影响。敏感性分析可以通过识别模型输入来减少不确定性,提高系统的稳健性。随机梯度估计方法可以估计复杂系统表现函数关于生产参数的敏感性,实现优化复杂系统结构、控制风险的目标。由于间接的梯度估计方法在实际问题总是可用,本节主要介绍直接的随机梯度估计方法在随机活动网络、质量控制图、排队系统和(s,S)策略库存系统4个系统中的应用。

2.1 随机活动网络

网络计划技术是制定项目进度计划的一种常用技术。在项目实施过程中,受天气、资金和材料供应等随机因素的影响,生产活动持续时间具有随机性。项目管理人员采用随机活动网络表示项目中计划完成的各项生产活动的先后顺序和逻辑关系,进而可以预测项目的工期和费用,优化资源配置,实现总工期最短的目标。应用随机梯度估计方法,可以估计项目的工期和费用关于某项生产活动的参数的梯度。如下考虑一种简单的随机活动网络模型[23],并介绍随机梯度估计方法在该模型中的应用。

随机活动网络由具有M个节点、N条边的有向非循环图表示,其中,每条有向边表示生产活动,生产活动时间为随机变量Xi,其分布为Fi,i=1,2,…,N,X1,X2,…,XN相互独立。假定节点1为任务起始点,节点M为任务终止点,P为构成起始点至终止点路径的有向边的集合,为所有路径的集合,P*∈为最优路径(活动时间最短或最长),最优路径的时间可以表示为:Y(X)=问题的目标是估计d E[Y]/dθ,其中,θ为Xi的分布参数。应用随机梯度估计方法,可得无穷小扰动分析估计量为

其中,1{·}是示性函数。似然比估计量和弱导数估计量分别为:

数值算例。如图1所示,考虑包含5个节点、6条边的随机活动网络。假定节点5为任务终止点,所有路径的集合为假设生产活动时间Xi服从均值为i的指数分布,即令θ=θ1=1为X1中的分布参数,P*为使生产活动时间最短的最优路径。应用随机梯度估计方法,可得无穷小扰动分析估计量、似然比估计量和弱导数估计量分别为:

图1 随机活动网络

实验分别在指数分布和正态分布的设定下,比较3种无偏随机梯度估计量和单边前向有限差分估计量FD(c)。基于106次独立的仿真实验,得到各估计量的数值(平均值±标准误差),其中,标准误差由估计的标准偏差除以样本数的平方根计算得到。

由表2可以观察到,当摄动值为0.1时,有限差分估计量的偏差较大;当摄动值为0.01时,有限差分估计量的方差较大。无穷小扰动分析估计量和弱导数估计量的方差比似然比估计量的方差小。基于108次独立的仿真实验,在指数分布的设定下,FD(0.001)的梯度估计量的数值为0.585±0.103 8;在正态分布的设定下,FD(0.001)的梯度估计量的数值为0.843±0.126 2(平均值±标准误差)。这表明,无穷小扰动分析估计量,似然比估计量和弱导数估计量具有无偏性。

表2 随机活动网络最优路径时间的梯度估计量(平均值±标准误差)

2.2 质量控制图

在生产过程中,受生产人员、机械设备、材料、工艺方法和环境等因素的影响,生产过程可能处于失控状态,导致产品的质量特性值偏离规定的质量特性值。质量控制图是统计过程控制的一种重要方法,它利用控制图提供的信息判定生产过程是否处于统计控制状态,有助于将生产过程维持在受控状态,减少生产阶段因产品异常带来的损失。生产过程输出的质量特性值有随机性,应用随机梯度估计方法,可以估计生产过程平均运行长度关于质量控制图的控制边界的梯度。如下考虑指数加权移动平均控制图 (Exponentially Weighted Moving Average,EWMA)[13-14],并介绍随机梯度估计方法在该模型中的应用。此外,通过一个考虑最小化平均运行成本的控制边界优化设计的数值算例,以期为随机梯度估计方法在复杂系统管理中的具体应用、改善的指标和得到的结论提供范例。

在统计过程控制中,生产过程的平均运行长度被定义为

其中:N为系统失控的时间,即N=min{i:Yi∉[θ1,i,θ2,i]};Yi为在得到第i个样本后观测到的检验统计量,对于指数加权移动平均控制图:Yi=αXi+(1-α)Yi-1,i＞1,Y1=X1。Xi为第i个样本输出,θ是一类参数,可以是θ1,i或θ2,i,分别为第i个检验统计量的控制下界和控制上界。当系统处于受控状态和非受控状态时,系统将输出具有不同分布的样本。取条件于Z=z,其中,Z=R/Δ,R为系统受控的随机持续时间(不可观测),服从密度为q0(·)的分布,Δ为采样间隔,则Xi的条件密度为

其中,f1(·)和f2(·)分别为系统受控和不受控时输出Xi的密度函数。假定系统失控时间服从参数为λ的指数分布,受控样本输出服从均值为μ0、方差为σ2的正态分布,失控样本输出服从均值为μ1、方差为σ2的正态分布。令θ1,i=1,θ2,i=2,问题的目标是估计d ARL（θ）/d2。由于结构参数和不连续样本表现的存在,无穷小扰动分析估计量,似然比估计量和弱导数估计量均不存在,由此可得广义似然比估计量为

2.2.1数值算例 考虑如下参数设定的指数加权移动平均控制图:α=0.5;系统失控时间服从均值为20的指数分布,即λ=0.05;采样间隔Δ=1;控制下界1=-1,控制上界2=1;受控或失控样本方差σ=1;受控样本均值μ0=0,失控样本均值μ1=1,μ1=3。通过实验比较广义似然比估计量和单边前向有限差分估计量FD(c)。分别基于104和106次独立的仿真实验,得到各估计量的数值(平均值±标准误差),其中,标准误差由估计的标准偏差除以样本数的平方根计算得到。

由表3可观察到,当摄动值为0.1时,有限差分估计量的偏差较大;当摄动值为0.01时,有限差分估计量的方差较大。基于108次独立的仿真实验,在失控样本均值μ1=1的设定下,FD(0.001)的梯度估计量的数值为10.2±1.12;在失控样本均值μ1=3的设定下,FD(0.001)的梯度估计量的数值为7.9±1.02(平均值±标准误差)。这表明,广义似然比估计量具有无偏性。

表3 指数加权移动平均控制图的梯度估计量(平均值±标准误差)

2.2.2控制上界的优化设计 生产过程的单位时间内的平均运行成本被定义为

考虑一个指数加权移动平均控制图,假设失控样本均值μ1=3,方差σ=2,其余参数设定与例2相同。令损耗成本c=2,系统维修成本C=10,随机逼近方法中控制图的初始控制上界为θ(0)=1,迭代步长an=1/n,第n步迭代使用(1 000+n)次独立的仿真实验估计梯度ΔE[N]和ΔE[(N -R/Δ)+]。采用随机逼近方法迭代10 000次得到的控制上界的样本轨道如图2(a)所示。基于108次独立的仿真实验估计系统的平均运行成本,得到的平均运行成本的样本轨道如图2(b)所示。

求解得到的控制图最优控制上界的值为θ*=1.139。由图2(b)观察得到结论:通过控制上界的最优设计,生产过程的单位时间内的平均运行成本从1.52降低为1.32。采用基于广义似然比梯度估计的随机逼近方法,解决了控制图控制上界的优化设计问题,降低了生产过程的单位时间内的平均运行成本,有助于提升企业经济效益。

图2 样本轨道:(a)控制上界(b)平均运行成本

2.3 单服务台排队系统

排队论又称随机服务系统理论,是通过分析随机服务系统中排队等待现象的稳态概率特征进行系统优化设计和控制的理论。排队论广泛应用于制造业产品生产、运输、库存等资源分配管理的随机服务系统。例如,在产品加工生产线流程设计问题中,可以构建加工设备服务待加工制品的排队系统;在企业设备维修问题中,可以构建维修员工服务故障设备的排队系统;在产品运输的集装箱港口资源配置问题中,可以构建堆场服务集装箱的排队系统。排队系统中存在着大量随机性,如产品到达时间间隔、工作站服务时间、系统等待时间和系统停留时间等,应用随机梯度估计方法,可以估计排队模型表现函数关于模型参数的梯度。如下考虑一种简单的单服务台排队系统[23],并介绍随机梯度估计方法在该模型中的应用。

考虑一个先到先服务的G/G/1 的排队网络。记Ai为生产订单的间隔到达时间,Xi为生产加工时间,Ti为生产订单在排队系统的停留时间。前N个订单在排队系统中的平均停留时间可表示为:令θ为Xi的分布Fi中的参数,问题的目标是估计假定系统初始状态是空的,生产订单到达过程与服务过程独立,且到达过程互相是独立的。根据Lindley方程,一个生产订单在排队系统的停留时间可表示为:Ti+1=Xi+1+(Ti-Ai+1)+。应用随机梯度估计方法,可得无穷小扰动分析估计量为

其中:M为观测到的忙期数;nm为在第m个忙期中的最后一个订单;似然比估计量和弱导数估计量分别为:

2.4 (s,S)策略库存系统

库存管理是企业物流和供应链管理中重要的研究问题。企业的库存是为了保持持续生产、满足客户需求而存储的原材料、在制品和制成品。库存管理是通过对客户需求的预测,在保证企业生产活动的正常进行和产品供应的前提下,加速物资流动速度,减少库存量,实现降低库存成本、提升企业经济效益和提高生产柔性的目标。供应链的库存管理通过企业间的协调配合和信息共享,减少牛鞭效应带来的需求误差,实现按需生产和降低库存总成本的目标。库存问题的随机性来自供应者、生产者和客户的不确定性。应用随机梯度估计方法,可以估计库存量和库存成本关于库存策略中的参数的梯度。如下考虑一种简单的(s,S)策略库存系统[22],并介绍随机梯度估计方法在该模型中的应用。

考虑一种生产原材料的定期检查(s,S)策略库存问题,在每个周期检查库存量,确定是否需要补货。(s,S)订购策略为如果当前原材料库存量低于s,则补货到库存量S。假设所有的超额需求被积压,会在下一周期被满足,原材料的订单没有延迟,在每一周期结束时决定是否补货。令Dn为第n周期的需求量,,Vn为第n周期结束时原材料的库存量,可以表示为

其中:E[D]为平均需求量;ξn=Yn-s;Yn为满足需求前的原材料库存量。

3 结论

随着新一代信息技术与先进仿真技术的集群式创新、融合发展与突破,大规模复杂系统的仿真数据日益可测可获,不仅有效提升企业生产质量、效率和效益,也孕育着企业管理决策理论与方法的重大变革,推动管理决策研究向仿真优化领域转变。仿真技术为企业提供了能够覆盖产品论证、研发、设计、采购生产、销售和服务等供应链各阶段的全生命周期解决方案,是企业数字化转型升级的重要技术手段。复杂随机系统的梯度估计包含比系统表现估计更多有用的信息,梯度估计的准确性常常影响着系统表现的优化分析及后续管理决策的进行。为提升仿真及优化效率,需结合问题情境特点选取高效准确的随机梯度估计方法。

在实际应用中,根据复杂系统的特定情况选择合适的随机梯度估计方法对后续优化分析及管理决策的制定至关重要。在进行方法选择时,可从以下几方面考虑。首先,如果没有获取到复杂系统的信息,优先考虑间接的梯度估计方法。进一步,如果获取到复杂系统的一些信息,如系统输入的分布或样本表现的特征等信息,则可以考虑应用直接的梯度估计方法。对于不连续型样本表现,如示性函数或非平滑函数,无穷小扰动分析法不适用。如果样本表现带有结构参数,则似然比法和弱导数法不适用。平滑扰动分析法需要视实际问题情况选取适当的条件计算条件期望,会消耗额外的仿真量。推出似然比法依赖于解析形式的变量替换来推出结构参数。对于样本表现不连续且带有结构参数的问题,优先考虑广义似然比法。其次,选择不同随机梯度估计方法的另一影响因素是参数向量的维度。有限差分法和弱导数法在处理高维参数向量的梯度时会消耗很多额外的仿真量。在基于梯度的参数估计方法中,选择合适的随机梯度估计方法需要同时考虑估计量的偏差和方差。有限差分法得到的梯度估计量是有偏差的。一般情况下,无穷小扰动分析法比似然比法方差小。

基于梯度的参数估计方法的效果提升及其在复杂系统管理实践中的应用是仿真优化领域关注的方向。

(1) 随机梯度估计方法与强化学习和人工智能算法的结合可以改善大型装配线在复杂环境中的质量控制任务,提高系统智能决策的抗干扰能力,实现人机共融的协同决策。深度学习与强化学习中训练人工神经网络所采用的反向传递算法属于仿真优化中研究的随机梯度估计算法。开发AlphaGo的谷歌公司Deep Mind 研究团队在文献[80]中对随机梯度估计方法在机器学习中的应用做了长篇综述。Peters等[81]讨论了有限差分法和似然比法在强化学习策略梯度估计中的应用。Peng等[82]采用推出似然比方法训练人工神经网络,该方法通过在神经元传递信号中加入噪声,不但可以处理离散型激活函数和损失函数的问题,也提高了人工神经网络在对抗攻击时的鲁棒性。

(2) 随机梯度估计通常可以用于参数的敏感性分析和基于梯度的仿真优化算法的输入,实现优化复杂系统结构的目标,进而度量和控制复杂系统的风险。随机梯度估计方法在项目计划制定、产品生产、运输、库存、设备维修和生产过程控制等系统管理领域有广泛的应用。经典的随机梯度估计的目标是随机模型样本表现的期望,而期望刻画的是系统在随机环境下的平均表现,未反映系统在极端环境下的表现。考虑随机模型样本表现的分位数的随机梯度估计可以得到在极端环境下更稳健的系统管理策略。在供应链系统管理中,考虑报酬收益分位数敏感性的报童模型可以得到在极端市场需求出现时表现更稳健的订购策略。Hu等[83]研究了以分位数为目标的仿真优化问题,基于广义似然比法提出了单样本轨道多尺度的随机逼近方法,解决了分位数敏感性估计中仿真样本轨道的不连续性以及由比值引起偏差的困难,建立了算法的收敛性并分析了其渐进收敛到平衡点的稳定性。Jiang等[84]研究了累计回报的分位数敏感性的强化学习问题,基于似然比法提出了两尺度的随机逼近方法,得到了在极端事件下表现更稳健的策略。近两年疫情导致的极端事件对全球供应链造成了巨大影响,极端事件下的系统管理将会是后疫情时代的重点研究方向之一。

综上所述,本文在系列相关研究的基础上对系统管理领域中随机梯度估计的应用进行了梳理,同时还介绍了常用的随机梯度估计方法和系统管理模型。针对常用的随机梯度估计方法,分别总结了其研究状态和应用场景,并通过模型介绍了这些方法在企业项目进度计划制定、生产过程控制、生产调度流程优化和供应链库存管理领域的应用。后续研究进一步探索随机梯度估计方法的求解路径和优化设计,如有限差分法最优摄动值的选取,并进一步思考相关方法在复杂系统管理实践中的具体应用。