随机多目标

随机多目标（通用7篇）

随机多目标 篇1

0 引言

多目标性能评价指标在跟踪系统设计与评价方面具有非常重要的意义,同时,多目标跟踪性能评价也是一个非常困难的问题。多目标跟踪性能指标作为评价多目标跟踪算法性能好坏的一个量值,对于算法选择、分析和评估起着非常重要的作用。一般来说,算法的性能指标主要包括:精度、复杂度、运行时间、收敛性、鲁棒性、实时性、一致性等等。在本文中,我们重点考虑是集值估计算法精度指标问题。

对于多目标跟踪性能的评价指标,最为常见的就是均方误差(MSE)。MSE指的是真值和估计值之间的期望差值。实际上,由于期望通常很难获得。因此,直接计算MSE指标就非常困难。因此,经常使用的一个指标是均方根误差(RMSE),它利用Monte Carlo仿真的采样值来统计逼近该期望值。RMSE在多目标跟踪领域是一个最为常用的指标,但是,RMSE指标有几个不足:首先,它不是欧氏空间上的距离概念;其次,当目标个数很大的情况下,例如上百个批次,用RMSE作为多目标跟踪评价指标就显得过于冗余。实际上在大批次目标情况下,单个目标的跟踪性能指标越来越被弱化。换句话说,在这种情况下,人们更多地关注对整体目标群的跟踪性能的评价,而不再关注单个目标的跟踪性能。然而,RMSE的定义是建立在各个目标真值和其估计值之间存在明确对应关系的基础之上的。或者说,需要首先考虑关联的问题,否则,RMSE是没法直接使用的。因此,RMSE在大批次目标情况下应用是很繁琐的,也是不现实的。在这种情况下,自然会考虑目标集的跟踪评价指标,这种目标集的评价指标实际上可以抽象为集合之间的距离。Wasserstein距离就是建立在评价集合间差异的一种度量指标[1,2],并被进一步引入到多目标跟踪领域[3]。然而,Wasserstein距离对于集合元素个数的差别惩罚过重,因此,一致性比较差。为了克服这个问题,Schuhmacher等引入Optimal Subpattern Assignment(OSPA)距离[4],它是建立在Wasserstein距离基础之上,并且修正了Wasserstein距离中的不足。Circular position error probability(CPEP)把目标跟踪位置正则化到某个阈值半径的圆上[5],判断目标是否丢失,统计出目标的跟踪丢失率。本文将对这些指标进行比较分析,同时给出算例进行分析说明,指出各个指标的特点、算法实现和应用方法等。

1 随机集多目标跟踪性能指标

1.1 均方误差(MSE)和均方根误差(RMSE)

均方误差的定义如下:

其中:xk是k时刻目标真实状态值,是目标状态估计值,如前所述,上式中期望值计算比较困难。因此,通常利用RMSE来代替:

其中:M是总的Monte Carlo仿真次数,xk,i是第i次仿真数据,是相应的状态估计值,它利用样本值统计值代替真实的期望值。RMSE误差不是欧式空间上的距离,因此,其度量指标和常规距离概念有所不同。此外,在多目标情况下,RMSE需要考虑目标真实状态值和估计值之间的对应关系,也就是说首先需要先解决关联问题。否则,两个不同目标状态之间相减没有任何意义,在这种情况下计算RMSE误差就失去了参考意义,这一点也是RMSE不能用于集值估计算法的主要障碍。

1.2 圆丢失概率(CPEP)

该指标用于评价目标跟踪丢失率的情况,在最初的文献[5]中定义如下

这里是2范数,r表示圆阈值。这里“所有的”指的是所有的估计值,一般包括两类:第一类是指Monte Carlo仿真值;第二类是多个估计器的估计值。CPEP取值在[0,1]区间之间,值越大,表示跟踪丢失率越大。从式(3)可以看出,它同样需要考虑目标真实状态和估计状态之间的对应关系。文献[5]中通过Monte Carlo仿真来计算CPEP丢失率。Vo等在文献[6-7]中引入该指标,并给出了如下的公式:

其中:Hk=[I2,02],r是位置误差半径。不过,这个式子的表达容易给人造成误解:Xk是所有目标的估计状态集,并不是(3)中的所有某个目标估计值,二者具有不同含义。因此,本文建议该公式修改为

上式表示对于所有估计值并且是对应状态为xk的目标估计值。实际上,如果利用Mont Carlo仿真计算该丢失率,我们可以给出如下的计算式:

式中:M是Monte Carlo次数或者估计器个数,是k时刻所有Monte Carlo估计集值,即,表示没有丢失的目标标签,是集合中的元素个数。上式可以直接计算,不需要考虑估计和真实值之间的对应关系。

1.3 Wasserstein距离

Wasserstein距离是统计理论上的一个指标,用于度量两个密度函数之间的距离[8],p阶Wasserstein距离定义如下:

其中:µ,ν是边缘概率分布,其对应的随机变量分别是xµ,yν;G(x,y)是所有xµ,yν的联合概率分布;g是其中任意一个联合分布。Wasserstein开始主要应用在度量数字图像某些特征的相近程度,例如纹理和颜色的相近程度[9]。Oliver Drummond等把该概念用于多目标误差判断[1,2],Hoffman等给出了Lp意义上的Wasserstein距离[3],即:

其中:C={Ci,j}是一个n×m矩阵,其中的元素Ci,j和满足如下的条件

因此∑ni=1∑mj=1Ci,j=1,其中,Wasserstein距离用于评价两个集合之间的距离,两个集中任意一个为空集时,Wasserstein距离不存在。该距离不仅评价集合元素之间的差异,并且也评价集合间元素个数(集合势)的差异。

1.4 OSPA距离

Wasserstein距离对于集合之间势的差别惩罚过重,并且当集合元素和势之间同时出现误差的条件下。即使对于集合元素数值误差很小,但是由于集合势存在误差,使得评价数值出现很大的误差,这和实际的直观理解不一致。因此,Dominic Schuhmacher等提出了修正的OSPA距离[4],定义如下。

分别对应p阶和无穷大评价指标。其中的参数c是一个水平参数,其物理含义是目标状态估计误差阈值,用于调节集合势的估计误差比重。

2 误差解释及性能分析

我们主要从以下几个方面分析这几个指标数学和物理上的含义及各自的特点。

2.1 距离含义

MSE本身就是距离的概念,虽然RMSE本身和欧氏空间上的距离概念有所不同,也是向量空间中的一种范数概念。CPEP是一种概率测度上的度量指标,指的是目标在某个球面邻域内的丢失率。Wasserstein是一种集合之间距离的度量指标。它本质上可以解释为一个如下的优化问题:

即满足上述约束条件的集合X和X之间的最短距离。当两个集合势不同时,两个集合之间的Wasserstein距离会增加。这时候,由于包含集合势的误差,在这种条件下,比较Wasserstein距离是没有意义的。当两个集合势相同时,Wasserstein可以解释为两个集合之间的最小平均距离。

和Wasserstein距离有所不同,OSPA距离也可以解释成为如下的一个优化问题:

OSPA距离基本含义是单个目标的平均误差,也可以被解释为位置和集合势两部分距离:

OSPA有两个参数p和c,p是距离敏感性参数,对于任意的参数1≤p1≤p2≤∞,OSPA距离满足下面的关系:

c是一个水平调节数,用与调节集合势误差的影响。也就是说,如果我们认为距离误差比势误差更重要,那么水平数c可以选择小点。相反,c取大一些的值。对于参数1≤c1≤c2≤∞,OSPA距离满足如下的关系:

从上式(18),(19)可以看出,OSPA距离对于参数p,c是一致的。

2.2 计算实现

RMSE计算简单,这里我们重点考虑其他三个指标的计算问题。CPEP是概率测度的计算,分别根据目标的真实状态和估计状态之间的差值和阈值r之间的关系来计算,主要通过Monte Carlo仿真来实现。Wasserstein距离是一个线性优化问题,因此,我们可以通过优化程序处理包来实现。OSPA距离可以利用匈牙利算法来解决。

从上述分析也可以看出,对于基于随机集的误差评估,或者集值估计的误差评判,传统的RMSE,欧式距离,范数等等的概念很难应用。主要原因是传统的算法是一种点和点之间误差评判准则,因此,可以直接应用的范数概念来解释评价。而基于随机集的估计理论则需要对集合估计性能进行评判,它属于点集-点集之间的度量,评判起来更为复杂,点集-点集之间的评价指标至少应该满足如下的条件:

1)集合中元素的顺序不应该影响最终的评判结果,或者说,集合元素不变,顺序变化,那么评价指标结果应该一致,不发生变化。显然,传统的范数是不能直接应用的。这和范数的交换率类似,即对于任意的两个点x,y,范数d(x,y)=d(y,x)。而对于任意的两个集合X,Y,集值评价指标满足如下的两类交换性条件:

上式第一个条件对应集合元素的次序变化,其中,Li,Lj,Jm,Jn分别表示不同的排序;第二个条件对应两个评价集合位置交换,这个也是集值评价指标的基本要求。

2)从1)的要求也可以看出,集值估计的评价指标必须满足元素位置的不变性,也必须满足集合位置的不变性。这个问题似乎只能转化成为一个优化算法的最优解问题。这一点也可以从Wasserstein距离和OSPA距离定义看出来,这种最优解常常是一个最优分配问题,目前基于随机集的多目标跟踪性能评价指标都是基于这种思路的。

3)RMSE很难作为随机集多目标跟踪算法的评价指标,主要原因是RMSE要求确定目标估计状态和实际状态之间的对应关系。如前所说,在大批次目标情况下或者集值估计算法中无法直接使用。

3 算例

在这部分我们分别给出线性和非线性两个算例,来比较分析这些指标的评价效果。由于RMSE需要考虑真值和估计值之间的对应关系,我们重点分析集值估计算法的评价,因此,这里暂不考虑RMSE,只分析CPEP,Wasserstein距离和OSPA距离三个指标。

3.1 算例1:线性跟踪系统

考虑三个目标的运动,目标系统运动方程和量测方程如下

目标初始位置分别为x01,=(500,20,-500,10),x0,2=(500,25,500,-20),Qk=cov(wk,wk)=diag(25,25),Rk=cov(vk,vk)=diag(100,100)。我们利用GM-PHD滤波器跟踪目标[6]。假设采样时间为T=1 s,目标的出生强度为υγ(x)=.01N(x,x01,,P0)+.01N(x,x0,2,P0)。其中矩阵Fk-1,Gk-1定义如下:

图1是目标状态估计过程,图2分别是x,y方向目标状态的估计,图3是基于1 000次蒙特卡洛仿真(MC)的目标个数估计。可以看出,GM-PHD滤波器能很好地跟踪三个产生于不同时刻的目标,估计目标个数也能很好刻画了实际目标个数的变化过程。图4分别用三个不同指标对系统的跟踪过程进行评价,通过1 000MC仿真,可以看出,在目标个数变化的时刻,Wasserstein距离会出现很大的误差,这个误差大小不受控制。即如果目标位置都很接近,但是因为估计个数变化,会导致评价指标出现很大的误差,这个和人们直观上的感觉是不一致的。

另外一种情况是目标位置估计误差相差很大,而个数估计准确,这时候Wasserstein距离可能并不大,这个也和实际中的直觉不一致。而OSPA指标通过引入水平参数c来解决这个问题。可以看出,如果目标位置估计比较重要,而目标个数出现误差相对不太重要,那么取比较小的水平参数c,在图4中对应c取100对应的值。可以看出目标个数出现变化的地方波峰不明显,OSPA曲线整体比较平缓,主要反映的是跟踪位置误差的变化。如果目标个数估计误差比较重要,那么我们可以进一步增加c的取值,那么相应的OSPA距离也会增加。随着水平参数c值变大,评价指标也增加。在目标个数估计不正确的地方,OSPA距离会进一步增加,甚至出现一个波峰。并且c值越大,波峰越明显。因此,OSPA参数能更全面地描述集合误差的估计。这种基于集合的误差评判指标是随机集理论中必须要考虑的问题。

图4中,在第6∼9时刻之间,OSPA指标出现一个波峰,主要原因是由于目标漏检的存在,使得目标个数为零,这意味着出现空集,在这种情况下。按照OSPA定义,取最大的水平值,而这个水平值一般远大于正常的跟踪误差。这样使得第10时刻反而出现峰谷。目标出生和死亡时间如下表1所示。

图5给出的是目标丢失率曲线,可以看出目标丢失的地方都集中在目标个数发生变化的地方。分别对应目标出生,死亡的时刻。例如第10时刻目标2产生,第20个时刻目标3产生,40个时刻目标1死亡。

3.2 算例2:非线性跟踪系统

为了比较方便,假设目标的运动方程不变,量测方程为如下的角度、距离量测方程:

其中:量测协方差阵Rk=cov(vk,vk)=diag(σ2k,θ,σ2k,r),σk,θ=2π/180rad,σk,r=10005m,杂波个数还是40个,其它参数和上面的线性方程一致。

我们通过Particle-PHD滤波器获得PHD粒子分布,通过MCMC方法估计多目标状态[10]。比较线性的情况(虽然线性和非线性很难直接比较,但是从相同的运动方程,相同的杂波个数,相同的观测区域上来说,可以提供一个比较一致的参考),总体说来,非线性情况下跟踪性能有所降低。这可以从图6,图7中看出。目标个数估计也偏大,主要原因是非线性作用后,目标的状态分布变得复杂,在杂波作用下,使得虚警个数增加。同时,Wasserstein距离,OSPA距离相对比线性情况有所增加。在第10个时刻目标2产生,三个指标都出现波峰,主要原因是在非线性条件下,虽然从图8看,第10个时刻目标个数估计的平均值比较精确,但由于目标个数估计值的方差变大,使得三个指标都出现波峰。在其他目标个数变化的地方,OSPA波峰不明显,主要原因是虚警的增加,使得整个估计误差变大。在丢失率方面,在目标个数发生变化的地方,丢失率明显增加。比较图5和图10,可以看出非线性量测的丢失率要高于线性量测系统。

4 结论

本文着重分析了三种常用于随机集多目标算法性能评价指标:CPEP,Wasserstein距离和OSPA距离。可以看出CPEP主要用于判断目标丢失率,针对CPEP计算中的一些不足,本文进行了分析和修正;Wasserstein距离可以用于评价随机集估计性能,但是它在目标个数估计不正确的情况下往往会出现很大的误差,这一点和直观上的理解不一致。而OSPA距离弥补了Wasserstein距离在这方面的不足,通过一个水平参数来弥补这个不足,值越小,说明目标位置的估计误差越重要;相反,值越大,说明目标个数估计出现的误差越重要。

参考文献

[1]Drummond O E,Fridling B E.Ambiguities in evaluating performance of multiple target tracking algorithms[J].Proc.Signal Data Processing Small Targets(S0277-786X),1992,1698:326-337.

[2]Drummond O E.Performance evaluation of single target tracking in clutter[J].Proc.Signal Data Processing Small Targets(S0277-786X),1995,2468:92-104.

[3]Hoffman J,Mahler R.Multitarget miss distance via optimal assignment[J].IEEE Trans.Sys.,Man,and Cybernetics Part A(S1083-4427),2004,34(3):327-336.

[4]Schuhmacher Dominic,Vo Ba-Tuong,Vo Ba-Ngu.A Consistent Metric for Performance Evaluation of Multi-Object Filters[J].IEEE Transactions on Signal Processing(S1053-587X),2004,56(8):3447-3457.

[5]Ruan Y,Willett P.The turbo PMHT[J].IEEE Transactions on Aerospace and Electronic Systems(S0018-9251),2004,40(4):1388-1398.

[6]Vo Ba-Ngu,Ma Wing-Kin.The Gaussian Mixture Probability Hypothesis Density Filter[J].IEEE Transactions on Signal Processing(S1053-587X),2006,54(11):4091-4104.

[7]Ba-Tuong Vo,Ba-Ngu Vo,Antonio Cantoni.Analytic Implementations of the Cardinalized Probability Hypothesis Density Filter[J].IEEE Transactions on Signal Processing(S1053-587X),2007,55(7):3553-3567.

[8]Jordan Richard,Kinderlehrer David,Otto Felix.The variational formulation of the Fokker-Planck equation[J].SIAM J.Math.Anal(S0030-1410),1998,29(1):1-17.

[9]Bissacco A,Chiuso A,Soatto S.Classification and Recognition of Dynamical Models:The Role of Phase,Independent Components,Kernels and Optimal Transport[J].IEEE Transactions on Pattern Analysis and Machine Intelligence(S0162-8828),2007,29(11):1958-1972.

[10]LIU Wei-feng,HAN Chong-zhao,LIAN Feng,et al.Multitarget State Extraction for the Probability Hypotheses Density Using Markov Chain Monte Carlo[J].IEEE Transactions on Aerospace and Electronic Systems(S0018-9251),2010,46(2):864-883.

随机多目标 篇2

关键词：风电场,随机动态经济调度,多目标,场景法,非线性原对偶内点法,对角加边结构,场景解耦,异步块迭代法

0 引言

风电场接入电力系统后的动态经济调度问题实际上是一个随机优化问题,目前国内外对此进行了相关研究。文献[1]考虑到风电功率场景选取与系统有功调节能力之间的内在联系,提出一种风电功率场景自适应选取方法。在此基础上,建立了基于场景分析法的风电并网系统有功调度模型。 文献[2]在考虑风电功率预测误差发生概率及其所导致的经济补偿费用基础上,得到统计意义上更为合理的经济调度模型。文献[3]在考虑风电并网功率的典型波动方式下,研究了计及风电接纳能力的电网调度模型,该模型使系统总煤耗和电网风电接纳能力得到了有效协调。但文中没有给出火电机组运行点发生改变时,电网可接纳风电功率的波动范围。文献[4]采用两阶段随机凸规划求解含风电接入的电力系统经济调度问题。求解过程中采用L型截面法将大规模随机规划问题分解成若干子问题,以提高求解效率。文献[5]运用时间序列法进行风速预测,通过随机机会约束将含风电场随机动态经济调度模型中的约束条件表现为概率的形式,并采用综合随机模拟、神经元网络和遗传算法的混合智能算法求解。文献[6]提出不确定运行条件下电力系统鲁棒调度的概念,采用场景树的方法生成不同的风电功率预测场景。文中指出鲁棒调度的特点是能够根据当前已有信息做出对系统运行中不确定因素具有一定免疫力的调度决策;通过鲁棒调度一般模型的构建,指出场景树约束使鲁棒调度区别于确定性调度方式的组合,使系统各种可能存在的运行场景之间存在了联系。文献[7-10]以风电功率预测为基础,采用拉丁超立方抽样和场景消除技术构建误差场景,将随机性机组组合问题转化为确定性模型。上述方法在应用于求解大系统、多场景、多目标动态经济调度问题时还存在困难。

目前国内电网有功调度执行的是节能发电调度模式,目标函数只有一个或是几个目标函数的加权和。本文尝试同时优化两个目标:发电总燃料耗量最小和购电费用最小,最终找到帕累托最优解集。决策者可以根据实际情况,从这个解集中选择一个折中解安排发电计划。其基本思路是:利用场景法将多目标随机动态经济调度问题转化为大规模多目标确定性动态经济调度问题,再借助法线边界交叉(normal boundary intersection,NBI)法[11,12,13,14,15]将其转化为一系列大规模单目标非线性规划问题,并用非线性原对偶内点法求解[16]。在应用非线性原对偶内点法求解这些大规模单目标非线性规划问题过程中,按照场景顺序排列的简化修正方程的系数矩阵具有对角加边结构[16,17,18,19,20]。因此可对其实施解耦,并采用异步块迭代法对解耦后的低维修正方程组进行求解[16],并应用于求解某省级电力系统的多目标随机动态调度问题。

1 多目标随机动态经济调度模型

由于风电出力具有随机性,而现有的预测技术并不能保证预测值完全准确,故导致风电场出力实际上是在预测值附近随机波动,因此,需要在风电预测的基础上考虑其随机性。在随机模型中,风电出力可以通过一些可能出现的场景进行模拟。本文将风电预测出力定义为预测场景,将考虑风电出力预测误差所生成的场景定义为误差场景。误差场景生成是借助拉丁超立方抽样、Cholesky分解和场景消除[21,22]3个过程来实现的。假设风电出力的预测误差服从正态分布,数学期望即为风电场在各时段出力的预测值,标准方差为与风电波动有关的参数。

1)目标函数。其中f1为发电总燃料耗量,f2为购电费用,即

式中:Pi(t)为机组i在时段t的发电功率;T为调度周期总的时段数,在本文中取96;N为常规发电机组数目;对于燃煤机组,Ai,2,Ai,1,Ai,0分别为常规机组i的耗量特性系数,对于燃气机组常采用线性拟合,有Ai,2=Ai,0=0,对于水电机组,有Ai,2=Ai,1=Ai,0=0;Ci为电网向常规机组i的购电电价,不同机组的上网电价不同,并假定同一机组在不同时刻的上网电价不变;CWind为电网向风电场购电电价,本文假设风电场属于电网公司,即CWind=0;PWi(t)为风电场i在时段t的有功出力;NW为风电场数目。

2)预测场景约束。其中第1个式子为功率平衡方程;第2个式子为常规机组上下限约束;第3和第4个式子分别为常规机组的滑坡和爬坡约束,即

式中:t=1,2,…,T;PLoad(t)为系统在时段t的负荷预测值为机组i的有功出力上下限;rdi和-rui分别为机组i的滑坡和爬坡率;T15为一个运行时段,即15min。

3)误差场景约束。其中第1个式子为场景s的功率平衡方程;第2个式子为常规机组在场景s的上下限约束;第3和第4个式子分别为常规机组在场景s的滑坡和爬坡约束;第5个式子为误差场景和预测场景之间的功率增减速度约束,即

式中:t=1,2,…,T;s=1,2,…,S,其中S为误差场景数目;Pis(t)为常规机组i在时段t场景s的有功出力;PsWi(t)为风电场i在时段t场景s的有功出力;Δi为发电机组i在15min内有功出力可以迅速调节的水平。

4)动态经济调度模型的紧凑形式

将式(1)至式(3)描述的多目标优化模型写成式(4)所示的紧凑形式。 其中,第2 个式子代表式(2)中的第1个式子,即预测场景下的功率平衡方程;第3个式子代表式(3)中的第1个式子,即误差场景下的功率平衡方程;第4个式子代表式(2)中的不等式约束和式(3)中的第5个式子,即预测场景下的不等式约束及误差场景和预测场景之间的功率增减速度约束;第5 个式子代表式(3)中的第2 至第4个式子,即仅与误差场景相关的不等式约束。

式中:x0为预测场景下所有常规机组在各个时段的有功出力向量;xs(s=1,2,…,S)为第s个误差场景下,所有常规机组在各时段的有功出力向量;f1(x0)为总燃料耗量;f2(x0)为购电费用。

2 多目标非线性规划问题的求解

仅考虑总燃料耗量f1(x0)最小进行单目标优化,得到最优解x0(1)*,对应于平面几何坐标系下的点f1*;同理可得考虑f2(x0)最小时的最优解x0(2)*,对应于点f2*。在由两个目标函数构成的坐标平面中,点f1*和f2*构成帕累托前沿的端点,连接它们之间的直线称为乌托邦线,如图1所示。

生成帕累托前沿曲线的步骤如下[14]。

步骤1:规格化目标函数。由于两个目标函数具有不同的量纲,为避免产生数值问题,采用式(5)所示的映射对帕累托曲线上的点进行规格化。

步骤2:生成乌托邦线上均匀分布的点。

步骤3:对应于乌托邦线上的每个点,将多目标优化转化为单目标优化问题。

步骤4:求每个单目标优化问题。

步骤5:上述获得的每个点构成了帕累托最优解,将其光滑连接就构成帕累托前沿曲线。

如图1所示,帕累托前沿越靠近乌托邦点(f2(x1*),f2(x2*))越好,因此A和B两点之间的距离越远越好,即λ越大越好。因此,经过步骤1至步骤3,多目标优化问题(式(4))可转换为与乌托邦线上的点相对应的帕累托前沿上的点之间的距离最大为目标的单目标优化问题,即

式中:λ为乌托邦线上的点和相对应的帕累托前沿上的点之间的距离;ω=[ω1ω2]T,ω 的取值决定了乌托邦线上点的分布;n为乌托邦线的法线向量;F-(x0)为通过对F(x0)进行规格化后得到的目标函数向量[14];-Φ为支付矩阵[12],第2个式子用来表示乌托邦线与法线向量之间的垂直关系。

为了简化问题的讨论,将式(6)中的第2 个和第3个式子合并写成:

下面简要说明如何用非线性原对偶内点法求解由式(6)描述的单目标非线性规划问题,并得到简化修正方程组。引入松弛向量su1,ssu2(s=1,2,…,S)将不等式约束变为等式约束,再引入拉格朗日乘子向量y3,y2s(s=1,2,…,S),yu1,ysu2(s=1,2,…,S),并引入对数壁垒函数消去松弛向量的非负性约束,从而构成增广拉格朗日函数如下:

式中:Nu1和Nu2分别为h1(x0,x1,…,xS)和h2s(xs)的维数。

根据KKT最优性条件,对增广拉格朗日函数求偏导,得到一组非线性方程组,再用牛顿法求解可得到简化修正方程组。按场景对简化后的修正方程和变量进行排序,可得到如下简化修正方程组:

式中:ΔZ0为与预测场景相关的变量增量组成的向量,ΔZ0=[Δx0Δy3Δλ]T;ΔZs(s=1,2,…,S)为与误差场景相关的变量增量组成的向量,即ΔZs=[ΔxsΔy21Δy22… Δy2S]T;L0,Ls(s=1,2,…,S),Ms(s=1,2,…,S)均为对称稀疏矩阵,其维数分别为96(N+1)+3,96(N+1),96(N+1)。

3 简化修正方程的场景解耦

对大系统而言,方程(9)为一个稀疏高维线性方程组,对其实施有效求解是多目标随机优化动态经济调度问题的核心。方程(9)的系数矩阵维数为96(N+1)(S+1)+3。在第4节将会看到,对于大系统、多场景的情况,可能无法对其求解。方程(9)中的系数矩阵具有对角加边结构,类似的结构在潮流和最优潮流计算[15,16,17,18,19]中均有应用,因此,本文利用这种结构特点对其进行解耦后再求解。

将式(9)展开,得

采用异步块迭代法[16,17]求解式(10)和式(11)。假设已完成k次迭代,在进行第k+1次迭代时,根据式(10)可得到:

将式(13)代入式(12)中可得:

通过式(12)和式(13)的交替求解,最终可求得ΔZs和ΔZ0。

式(12)和式(13)的系数矩阵维数分别为96(N+1)+3和96(N+1)。在下节将会看到,对于大系统、多场景的情况,可对其实施有效求解。

4 算例分析

4.1 试验系统简介

以某省级电力系统的数据为例验证本文算法的有效性,动态调度周期取1d,等分为96个时段。该电网的装机容量为23 595 MW,包括:燃煤机组40台,容量为18 740 MW;水电机组10台,容量为4 855 MW;风电场1座。最大负荷为21 112 MW,负荷预测数据如图2所示。本文对风电接入最大功率分别取777 MW和1 716 MW两种情况进行研究,即对应风电占最大负荷百分比分别为3.68%和8.13%,其风电出力预测曲线分别对应图3中的“风电1”和 “风电2”。发电机组的燃料特性系数及上网电价见附录A表A1。

假设风电出力预测误差服从正态分布,数学期望为预测场景下的风电出力,标准方差设定为0.2。采用拉丁超立方抽样方法进行抽样,分别对图3中的风电预测曲线“风电1”和“风电2”产生原始样本200个,形成两个计算案例,以下简称案例1和案例2。并对每个案例经过场景削减分别生成10,50,100,200个误差场景的情况。

4.2 案例1

4.2.1 计算规模和计算性能对比分析

针对含上述确定的误差场景情况,采用通用代数建模软件GAMS框架下的CONOPT求解器[23](算法1)、非线性原对偶内点法(算法2)和本文提出的解耦算法(算法3),分别计算多目标动态经济问题的帕累托前沿。在将多目标问题转化为单目标问题时,将乌托邦线10等分,则在帕累托前沿上与这些等分点对应的有11个。非线性原对偶内点法和场景解耦法均用MATLAB工具实现,所用计算机为Intel Core i5CPU M480 2.67GHz/4GB内存。

表1和表2分别列出了多目标随机经济调度问题的规模和3种算法计算时间的比较。其中,算法1的收敛精度设定为:目标函数梯度不大于10-8;算法2和算法3的收敛精度设定为:补偿间隙不大于10-6,等式方程的最大残差不大于10-4。由表2可以看出,算法2和算法3的计算速度比算法1快得多。误差场景数为10时,算法3的计算速度略快于算法2。随着场景数的增加,算法2计算时提示内存不足,无法计算。

4.2.2 帕累托前沿对比分析

由于算法1和算法2只能求解10个误差场景的情况,故表3仅列出了3种算法求解含10个误差场景的动态经济调度问题的帕累托前沿的对比。

在帕累托前沿的两个端点上3种算法的结果存在一定差异:①购电费用最小点,3种方法计算的购电费用相同;②煤耗最小点,算法2和算法3计算的煤耗相同,与算法1 相差0.006%。在帕累托前沿的其他点上两种算法的结果也存在一定差异,造成这种差异的原因有:①两个端点的误差;②计算帕累托前沿时解耦本身造成的误差。3种算法购电费用和煤耗的误差范围为0~0.02%。

由燃料耗量和燃料价格可确定发电生产成本。取煤价800元/t,以表3中的算法3为例,列出帕累托前沿的两个端点和中间各折中解处以特殊的总成本形式表示的各目标函数值。对比发现,折中最优解(第7个点)对应的调度方案的总成本最低,与两个端点相比分别节约710万元和239万元以上,综合经济效益显著,且符合电网公司节能减排和降低购电费用的要求。

4.2.3 误差场景数目对帕累托前沿的影响分析

取10,50,100,200个误差场景时获得的帕累托前沿如图4所示。

表4列出了含50,100,200个误差场景的动态经济调度计算结果,均由算法3求得。通过表4不同误差场景数下的对比结果可知:①计算出的购电费用最小点相同,煤耗最小点相同;②当购电费用最小时,煤耗存在一定差异,当煤耗最小时,购电费用存在一定差异;③帕累托前沿上的其他点也存在一定的差异。造成这种差异的原因为:①为适应不同的误差场景,发电机出力不同造成的误差;②计算帕累托前沿时解耦本身造成的误差;③两端点的误差会影响到帕累托前沿上的其他点。

经检验可以看到,经过场景削减后取10个误差场景获得的动态经济调度结果对于剩余的190个误差场景均可行。对于取50,100,200个误差场景的情况,也有类似的结果。产生这种现象的原因在于,风电接入最大功率占最大负荷的百分比仅为3.68%,系统有足够的旋转备用应对风电功率的波动。

4.3案例2

此案例对应于风电接入最大功率占最大负荷的百分比为8.13%,经过场景削减同样分别生成10,50,100,200个误差场景的情况。在这种情况下,其计算规模、计算性能及误差场景数目对帕累托前沿的影响规律与案例1获得的结果类似。

但对应于10,50,100,200个误差场景的情况,其动态经济调度结果对于剩余的误差场景则并不完全可行。取10个误差场景获得的动态经济调度结果对于剩余的190个误差场景有4个不可行;取50个误差场景获得的动态经济调度结果对于剩余的150个误差场景有2个不可行;取100和200个误差场景获得的动态经济调度结果对于剩余的误差场景均可行。产生这种现象的原因在于,风电接入最大功率占最大负荷百分比较高,需要考虑如何选择误差场景数目。取100和200个误差场景时获得的帕累托前沿如图5所示。

为实现常规发电机出力在预测场景与误差场景之间的增减,需要快速调整发电机的出力。因此,场景越多计算结果越精确,调度决策对于系统中由于风电并网带来的不确定性的适应能力也越强。

5 结论

1)应用场景法和NBI法将多目标随机动态经济调度问题转化为一系列单目标非线性规划问题,进而用非线性原对偶内点法求解,可得具有对角加边结构的高维线性修正方程组,便于实施解耦计算。

2)运用异步块迭代法对具有对角加边结构的高维线性修正方程组进行解耦,从而大幅度降低了存储需求,增强了适应大系统、多场景的计算能力。

3)误差场景数目越多,计算结果越精确,优化方案对于系统中由于风电并网带来的不确定性的适应能力也会越强。

随机多目标 篇3

关键词：市场,广告预算,广告累积效应,随机优化模型

引言

企业每年都会投入数千万的资金用于推广和宣传自己的产品。 一个有效的优化模型可以帮助企业在合理的预算下获得更好的预期收益[1,2,3]。 针对多阶段多商品广告预算问题(以下简称MAB问题),文献[4]提出了一种确定性优化模型,并且在实际的案例中得到了有效的运用。 在确定性模型中,每一个未知的参数都被一个估量代替,问题的所有参数都被视作已知。 但是,在决定广告预算的过程中许多参数都是未知的。 这种不确定性也是不能忽视的。

本文针对MAB问题提出了一种随机优化模型,解决实际问题中参数的不确定性。

一、单一情形随机优化问题

这里首先引入单一情形的随机优化模型。 给出成本函数F(x,ζ),定义随机问题PS:

其中,FS(x):=E[F(x,ζ)]决策向量,可行域为,随机向量ξ的概率分布于集合中。这里可以看到,不确定参数只在成本函数中出现,并没有出现在可行域中。

二、MAB问题随机优化模型的建立

建立MAB模型的目的是使得广告带来的收益最大化, 这里我们需要区分商品不需要做广告就可以产生的基础销售量和由广告带来的销售量。 文献[5,6]中给出了广告行业里许多关键的概念。

文献[4]中提出了MAB问题的确定性优化模型,而解决这一问题更实际的方法是需要把以下随机因素考虑进去:广告的饱和水平,广告收益的递减率以及利润函数中的白噪声。 本文所介绍的随机优化模型,就是在确定性优化模型的基础上将这些因素考虑在内。

(一)广告的累积效应函数

广告的效应可以分为直接的效果和潜在的效果,而且这些效果都具有累积性,在进行广告投入的时候,一定要考虑往期广告的累积效应[7,8]。 广告累积效应函数将之前的广告投入和延期生效的广告效应考虑在内。 正如文献[10]中提到的,一个典型的方法是用Broadbent累积效应模型:

其中,是初始的累积效应值,yt是t阶段的累积效应值,δ是广告效应的保留率,gt是t阶段的广告投入。

(二)利润函数

本节在文献[4]中确定性利润函数的基础上引入随机MAB问题的利润函数P(g,ξ),定义如下:

其中:

对上述模型的分析如下:

误差 ε 的产生是由于模型中遗漏了某些变量和一些其他的随机扰动;ξ 是一个包含了问题中所有的随机参数的随机向量;αtjk表示t阶段产品j在媒介k上的广告饱和度,βtjk表示阶段产品j在媒介k上广告效应的递减速率,γtijk表示产品j对与产品i的交叉影响;等式(10)、(11)分别表示初始和最终的广告累积效应,等式(11)则是为了表达上的便利而引入的虚拟条件。

从统计学的角度看,P(g,ξ)等价于一个非线性回归问题,其中广告的投入向量g是自变量,利润P=P(g,ξ)是因变量。 因此,P是一个随机变量。 最常见的形式就是回归分析中估计独立随机变量的条件期望,即:E[P|g]=P(g,ξ)

因此,为了计算广告投入以产生最佳期望收益,我们就需要解决以下优化问题:

其中,定义了可行的广告投资的集合,同时注意到中并没有不确定参数。 那么这个问题就又回到了带有确定可行域的单一情形随机优化问题中去。

(三)目标函数

定义成本函数:

其中,ξω定义了随机向量 ξ 的各种现实情形[9]。 赋值每一种情形 ω 的权重为Wω, 定义有限支持向量集,从而有[11,12]。

(四)随机优化模型

定义MAB的随机优化模型:

其中,。注意到随机参数只出现在成本函数中,假定可行域为为线性的,则有:,表示广告预算。

显然,可以证明该问题是一个凸优化问题[10]。

三、结果比较

当是一组相互独立的随机变量时,

在文献[4]中,PS被确定性问题PEV近似代替,PEV定义为:

则有,

例给出成本函数:.则:

给出一定的参数值后运算结果如图:

我们可以看出,出于最优化的考虑,FEV(x)相比于FS(x)的近似效果并不理想。

结论

本文在文献[4]的基础上,将原问题中参数的不确定性加以考虑, 提出了随机优化模型解决MAB问题通过简单的计算比较, 随机优化模型的效果要优于确定性模型。

参考文献

[1]M.Garrat,R.Alcaraz,and S.Cohen.From brilliant to actionable:It takes technical brilliance and constant questioning to achieve the truly actionable in marketing ROI.In Esomar-Congress,2011

[2]P.Doyle and J.Saunders.Multiproduct advertising budgeting.Marketing Science,1990,9(20),pp97-113

[3]J.D.Little.Models and managers:The concept of a decision calculus.Management Science,2004,50(12),pp1841-1853

[4]C Beltran-Royo,H.Zhang,L.A.Blanco,and J.Almagro.Multistage multiproduct advertising budgeting.European Journal of Operational Research,2013,225(1),pp179-188

[5]P.J.Danaher and R.T.Rust.Determining the optimal return on investment for an advertising campaign.European Journal of Operational Research,1996,95,pp511-521.

[6]D.M.Hanssens,L.J.Parsons,and R.L Schultz.Market response models:Econometric and time series analysis.Kluwer academic publishers,second edition,2001

[7]S.Broadbent.One way TV advertisements work.Journal of the Market Research Society,1979,21(3),pp139-165

[8]D.M.Hanssens,L.J.Parsons,and R.L Schultz.Market response models:Econometric and time series analysis.Kluwer academic publishers,second edition,2001

[9]J.R.Birge and F.Louveaux.Introduction to Stochastic Programming.Springer,2nd edition,2011

随机多目标 篇4

到目前为止,多变量辨识方法比较多,比如最小二乘算法、随机梯度算法、辅助变量算法、极大似然算法等等,这么多算法究竟哪个计算量更小,收敛性能更好,实用性更强,谁优谁劣究竟如何评判,这是问题的关键。不过理论上很难从多方面判断算法的优劣和实用性,可通过具体的仿真例子为辨识算法的选择作理论上的指导。本文用随机梯度搜索原理或随机逼近原理,研究了遗忘梯度算法,子系统遗忘梯度算法和递阶遗忘梯度算法,并通过仿真比较研究它们的估计精度的高低和收敛性能的好坏,使算法能更好的指导实际的应用。

1 多变量系统辨识模型

一个状态空间描述的多变量系统,化为传递矩阵描述,经过参数化可以化为下列形式的差分方程[1],

(1)式中u(t)∈瓗r和y(t)∈瓗m分别是系统输入和输出向量,v(t)∈瓗m为零均值白噪声向量瓗1和瓗m×r分别是待辨识的系统参数和参数矩阵。

定义参数矩阵θ、参数向量为α、输入信息向量φ(t)和输出信息矩阵ψ(t)分别为

瓗

瓗ψ(t)∶=[y(t-1),y(t-2),…,y(t-n)]∈瓗m×n因此,可以得到随机系统辨识模型,

2 梯度辨识算法

辨识模型(2)式既包含一个参数向量α∈瓗n,又包含了一个参数矩阵θT∈瓗m×n0,使得传统的辨识方法,如递推最小二乘算法(RLS)和随机梯度算法(SG)等,不能直接辨识这模型的参数向量和参数矩阵[2]。一种直观的方案是,将该系统重新参数化,得到一个包含参数向量α与参数矩阵θ所有元的高维参数向量∈瓗n+mn0。另一种方案是,把辨识模型(2)式按照输出的数目分解为m个子系统,然后分别辨识每个子系统的参数向量。还有一种方案是,采用递阶辨识原理,直接辨识参数向量α∈瓗n与参数矩阵θ∈瓗n0×m,这种方法计算量最小在这三种辨识方案中,都可以采用最小二乘优化原理和随机梯度搜索原理(或随机逼近原理),来研究相应的辨识方法。限于篇幅,本文只介绍基于随机梯度的辨识方法。下面分别介绍这几种辨识方案和具体辨识算法,并比较各种方法的优缺点。

2.1 随机梯度算法

先引入一些符号。设Im是m×m单位阵:表示Kronecker积或直积,col[X]表示将矩阵X的列按次序排成的向量。

定义扩展的参数向量和扩展的输入输出信息矩阵分别为

瓗。于是辨识模型(2)式可以转化为熟悉的形式,

极小化准则函数可得带遗忘因子的随机梯度算法,简称遗忘梯度算法(FG),

其中(t)是的估计,λ为遗忘因子称为时变收敛因子。算法初始值可取为一个小实参数向量,如(0)=1(mnr+n)×1/p0,p0=106。遗忘梯度算法比SG算法有较快的收敛速度,而且具有克服数据饱和和跟踪时变参数的能力。关于算法的收敛性,可参考文献[3,4]证明。

2.2 子系统随机梯度算法

FG算法虽然比较直观,思路也简单,但是参数向量维数很高,且信息矩阵Υ(t)含有大量的零元素,这使得算法计算量很大。为克服这些缺点,下面研究分子系统的辨识算法,基本思想是将一个多输入多输出系统分为多个多输入单输出子系统,然后分别辨识每个分子系统的参数向量。

令表示向量y(t)的第j个元表示信息矩阵ψ(t)第j行表示信息矩阵v(t)第j个元表示参数矩阵θ第j列。辨识模型(2)式可以分解为m个子系统,

(5)式中第j个子系统的信息向量和参数向量分别为

定义和极小化m个准则函数

可以得到每个子系统的带遗忘因子随机梯度算法,简称子系统遗忘梯度算法(S-FG),

因为在每一个子系统都要估计α向量一次,所以α的估计可取m次估计的平均值,子系统随机梯度法避免了出现大量零元的信息矩阵,参数维数也由原来的n+mn0降低到n+n0,计算量大大减少,所以S-FG算法是优于遗忘梯度算法。

2.3 递阶随机梯度算法

尽管子系统随机梯度算法降低了参数维数,减少了计算量,但是由于每一个子系统都包含参数向量α,所以要(重复)估计m次α,增加了不必要的计算量。而下面要讨论的递阶随机梯度算法,既不包含大量零元,又不必(重复)估计参数向量α,计算量是最小的。

递阶随机梯度算法是基于递阶辨识原理提出的[1]。根据式(2)定义准则函数,

根据负梯度搜索原理,得到估计α和θ的带遗忘因子的递阶随机梯度辨识算法(FFHSG),简称递阶遗忘梯度算法(HFG),

这里α(t)和θ(t)分别是α和θ的估计。

3 仿真实验与算法比较

考虑2输入2输出一阶仿真系统,其中α=-0.8,θT=[2.0,1.0;1.0,2.0]。仿真时,{u1(t)}和{u2(t)}采用零均值单位方差不相关可测的持续激励信号序列,{v1(t)}和{v2(t)}采用零均值方差分别为σ2(1)=0.42和σ2(2)=0.52的白噪声序列,对应两个输出通道的信噪比分别为δns(1)=50.99%和δns(2)=51.51%。用本文三个算法估计这个系统参数,这里由于篇幅所限,只给出λ=0.98时参数估计相对误差误差随t变化曲线如图1所示。

从这个仿真例子的仿真结果,可以得出结论:这三个算法有一个相同点,即随着递推计算步骤的增加,参数估计误差(总的趋势)不断减小,即数据量越大,参数估计精度越高;当λ相同时,S-FG算法收敛最快,FG算法收敛最慢,HFG算法介于二者中间,而且平稳性也介于两者之间;对于很大的t(数据长度),三个算法最终的估计误差曲线相差不大,参见图1,但递阶遗忘梯度算法计算量最小。

4 结论

本文针对多变量系统研究了三个随机梯度型辨识方法,并对这三个算法的性能和特点进行了分析比较。

摘要：针对工业中广泛存在的多变量系统,研究了辨识这类系统的遗忘梯度辨识算法,分子系统遗忘梯度辨识算法和递阶遗忘梯度辨识算法,对这三种算法的计算量进行了比较分析,并给出了仿真例子。

关键词：最小二乘,随机梯度,递阶辨识,多变量系统

参考文献

[1]Feng Ding,Chen Tongwen.Hierarchical gradient-based identification of multivariable discrete-time systems.Automatica,2005;41(2):315—325

[2]方崇智,萧德云.过程辨识.北京:清华大学出版社,1988

[3]丁锋.随机梯度算法的收敛性分析.清华大学学报,1999;39(1):83—86

随机多目标 篇5

1 基于多期随机优化的个人财务计划模型

1.1 模型的相关描述

基于多期随机优化的个人财务计划包括养老金计划、债务计划、子女升学计划等等，如果有一项计划未完成，就需要受到惩罚，惩罚的程度是Pd，如果在规定的时间t内完成相关的计划，那么Dtn为0，如果在规定的时间t内没有完成相关的计划，那么Dtn为1。考虑到个人超额消费的问题，其效用是一种消费金额函数，那么模型的表达方式为：

在以上公式中，λ1，λ2均为权重，是根据个人偏好的差异，按照超额消费效用反应出消费的金额，那么U，期末财富效用可以使用期末财富数量来进行标示，即U (wTn)=wTn，那么此时就可以将上述的函数简化为：

以上的约束为：

在以上的公式中，Xhitn是在第t年资产i的投资额，smtn即超额消费，mbtn即第t年借入的资金量，Stn是第t年工资的收入，Ltn是第t年末未支付的部分，Ctn即核心消费，rtn是第t年初风险资产收益，gtn是第t年目标支付，rbtn是在第t年中拆入资金的成本情况，以上的参数是一种随机参数。

1.2 情景的生成

个人财务计划模型是建立在经济环境不确定的基础之上，模型核心问题就是处理未来不确定性的因素，未来不确定性因素可以由基本消费、工资水平、资产收益等变量来反应，因此，为了建立一个动态的模型，整个状态的构建必须能够以未来信息变化以及时间变化来反映，这就可以反映出整个情况的结构。图1是一种情景树，每个节点代表不同的状态，每一种状态都反映着阶段中随机变量的结果，每个路径都是一种情景，在t=0时，是初始状态，在初始状态的变量是已知的，在t=1时，会出现一些可能的状态，在情景树的节点之后会有很多的后续节点，可以对揭示的信息进行模拟。

一般情况下，我国居民的投资方向大多以国债、存款、股票、基金等为主，因此，本文考虑的资产种类也以存款、债券、股票为主，分析近年来资产物价的变动情况以及工资的上涨率，可以建立好相关的模型，并建立一种向量自回归的模型，使用该种模型就能够分析出情景树中不同节点的基本消费、资产收益以及工资的增长情况。

1.3 基于个人采取计划的多起随机优化模型遗传算法

个人的财务计划可以按照以下的程序进行：将基本的生活开支安排好，并规划出未偿债务的支出情况，如果资金不足，那么就难以满足生活的开支需要和设定的目标需求，那么此时就可以对外借款，如果资金情况难以满足生活开支需求以及设定目标需求，那么就可以将股票、存款、债券等投资项目进行出售或者以此来进行抵押贷款，如果在出售和抵押贷款完成后，但是资产依然不足，那么就可以将相关的支出进行延迟处理，按照个人财务计划的步骤，遗传算法的计算方式如下：

第一，选择好初始的生成种群解，设置好决策变量。

第二，计算出父代种群个体适应度函数F (X)，适应函数的设计为：

在上式中，M是一个正数。公式设计好之后，就可以选择操作，使用旋轮法来进行运算。

第三，进行交叉操作，并确定好交叉操作的父代，完成之后，会产生两个后代，再验证两个后代的实际可行性，如果后代可行，那么则可以使用后代来代替父代。

第四，进行变异操作，从i=1到种群大小m，重复以下的过程，在[0, 1]中随机产生r，如果r<变异概率，那么就选择个体vi为变异父代，M是一个足够大的正数，那么在空间RM中就会随机出现变异方向d，如果v+M检验不可行，那么就要将M设置成0到M之间的一个随机数，此时就会得到新的个体，再进行检验，知道结果可行为止。

第五，计算出种群的适应度函数，选择适宜的数字作为目标值，进行重复计算，直到结果可行为止。

2 仿真计算

某个家庭想要制定一项长期的采取计划，其初始投资存款为10万元，股票有20万元，国债有10万元，家庭年收入为8万元。在家庭负债方面，有房贷，房贷年偿还金额为2万元，偿还期限共计20年，孩子教育费用每年2万元，共计要支付10年，必备的生活开支为1.6万元，根据通货膨胀的情况，该家庭购买了养老保险，每年需要支付1.5万元，在60岁以后，每年可以得到2万元的养老保险，那么怎样才能够保证基本支出的同时满足家庭消费欲望，就需要根据其债务以及退休的计划，根据其资产的情况做好投资计划。

根据以上的公司求得各项工资、资产收益以及核心消费，按照遗传算法的步骤来求解，将初始种群m设定为50，净化代数T设定为100，求解的结果详见表1:

参考文献

[1]Mulvey J M, Vladimirou H.S tochast ic netw ork opti-mizat ionmodels for investment planning[J].A nnals of Op erati onsResearch, 1989, 20 (2) :187-217.

[2]津巴威廉T, 马尔维约翰M.全球资产与负债管理建模[M].顾娟, 张永乐译.北京:经济科学出版社, 2003.530-565, 566-586, 639-665.Ziemba W T, Mulvey J M.Worldw ide asset and l iabi li tymodel ing[M].Translat ed by Gu J, Zhang Y L.Beijing:Eco-nomic S cience Press, 2003.530-565, 566-586, 639-665.

[3]金秀, 冯英洁, 黄小原.基于多期随机优化的个人财务计划模型[J].东北大学学报, 2005 (09) .

随机多目标 篇6

矩量法(MOM)被看作是数值"精确解"，被广泛应用于电磁散射问题的分析中，但每次计算只能得到一个频率点的电流分布和雷达散射截面(RCS)，而在日渐发展的遥感技术和雷达目标识别中，需要随机分布目标的宽带RCS以产生一维距离像。与其它频域方法一样，为了获得目标的宽带RCS，应用MOM就必须在频带内的每个频率间隔点上逐点计算，当目标的RCS随频率变化剧烈时，必须以很小的频率间隔计算才能得到精确的频率响应，这就意味着在整个频带内矩阵方程求解次数的增加。为了克服这个缺点，E.H.Newman和G.J.Burke分别通过内插阻抗矩阵[1]和使用基于模型参数估计[2]的方法获得了宽带数据。一种类似的技术--渐近波形估计(AWE)技术被提出,并首先用于超大规模集成电路的时域分析[3]。近年来，AWE技术被逐渐应用到电磁问题的分析中[4,5,6,7,8]。

本文将AWE技术应用到MOM中，计算了二维随机分布的理想导体柱的宽带RCS。数值计算表明，本方法与传统MOM的逐点计算结果吻合良好，而计算效率却显著提高。

1 渐近波形估计技术

对于任意形状二维理想导体柱的电磁散射问题，其电场积分方程(EFIE)为：

式中为导体柱截面的周线，为导体表面上的电流密度，为入射电场，为第二类零阶Hankel函数，和分别表示场点和源点的位置矢量，k、分别为自由空间波数和波阻抗。

将理想导体柱的周界划分成N段，选择脉冲函数作为基函数，函数作为检验函数，应用MOM可将EFIE化成矩阵方程Z(k)I(k)=V(k)(2)

式中Z为阻抗矩阵，V为激励向量，它们的元素表达式分别为：

为入射平面波与轴间的夹角。

应用矩量法求解式(2)，只能得到一个频率点的导体表面电流。为了得到给定频带内的导体表面电流分布，就必须以一定的频率间隔重复求解式(2)。AWE技术通过将展开成关于的泰勒级数得到频带内的导体表面电流分布，即

展开系数的表达式为:

式中Z(i)表示Z(k)的i阶导数，V(n)表示V(k)的n阶导数。

为了扩大泰勒级数的收敛半径，可通过Padé逼近将展开成有理函数,即

确定，将其代入式(8)，便得到给定频带内任意频率点的表面电流分布，进而计算出散射场和宽带RCS。显然，在AWE方法中，只需进行一次矩阵求逆运算，便可得到给定频带内电流密度的分布，这正是AWE方法提高计算效率的原因所在。

2 数值计算与比较

为了验证应用AWE方法分析二维随机分布导体目标宽带电磁散射特性的有效性，本文对如下两种情况的宽带RCS进行了计算。

(1)36个正方形理想导体柱随机分布在边长为0.0103米的正方形平面内，每个正方形理想导体柱的边长为λ0/20(λ0为中心频率对应的波长)，根据蒙特-卡罗方法[9]生成的二维理想导体柱分布状态如图1所示，选择入射波为TM波，E zi=e-jkx，中心频率0f=3 5GHz,Padé逼近中的L=4,M=3,20～50GHz范围内的RCS频率响应如图2所示。

(2)16个正方形理想导体柱随机分布在边长为0.0102米的正方形平面内，每个正方形理想导体柱的边长为λ0/10，分布状态如图3所示，其他参数同(1)，20～50GHz范围内的RCS频率响应如图4所示。

图3 N=16,a=0/1 0,l=0.0 10 2 m图4 RCS与频率间的关系(参见右栏)

在上述两种情况下，每个正方形理想导体柱离散单元尺寸选择为λ0/40，从以上两个算例的RCS频率响应曲线来看，用Padé逼近的AWE技术得到的结果与MOM逐点计算的结果有较好的一致性，而计算效率则显著提高，表1为上述两种情况所用的计算时间。

表1计算时间及与传统矩量法的比较*(参见右栏)

3 结论

本文给出了一种分析二维随机分布导体目标宽带电磁散射特性的方法，文中对两种情况下的随机分布理想导体柱的宽带RCS进行了计算，结果表明，本文的计算结果与矩量法逐点计算的结果吻合良好，而计算效率显著提高。对于任意给定的频带，要想获得精确的宽带RCS，需要使用多个频率展开点的AWE技术。

参考文献

[1]Newman E H.Generation of wide-band data from the method of moments by interpolating the impedance matrix[J].IEEE Trans.Antennas Propagation,1988,36(12):1820-1824.

[2]Burke G J,Miller E K,Chakrabarti S,et al.Using model-based parameter estimation to increase the efficiency of computing electromagnetic transfer functions[J].IEEE Trans.Magn.,1989,25(4):2807-2809.

[3]Pillage L T,Rohrer R A.Asymptotic waveform evaluation for timing analysis[J].IEEE Trans.Computer-Aided Design,1990,9(4):352-366.

[4]Polstyanko S V,Dyczij-Edlinger R,Lee J F.Fast frequency sweep technique for the efficient analysis of dielectric waveguides[J].IEEE Trans.Microwave Theory Tech.,1997,45(7):1118-1126.

[5]Li M,Zhang Q J,Nakhla M.Finite difference solution of EM fields by asymptotic waveform techniques[J].IEE Proc.Part H:Microwaves,antennas and propagation,1996,143(6):512-520.

[6]Erdemli Y E,Gong J,Reddy C J,et al.Fast RCS pattern fill using AWE technique[J].IEEE Trans.Antennas and Propagation,1998,46(11):1752-1753.

[7]孙玉发,徐善驾.渐近波形估计技术在三维电磁散射问题快速分析中的应用[J].电子学报,2002(60)794-796.

[8]施长海,孙玉发.二维电大导体目标宽带雷达散射截面的快速计算[J].电波科学学报,2004(6):325-328.

随机多目标 篇7

闭环供应链是由正向和逆向供应链经整合后形成的一个封闭的供应链系统, 其目的是对物料进行循环利用, 减少污染排放和资源浪费, 同时改善供应链整体的绩效。因此, 闭环供应链除了包含传统供应链的特征, 还对环境保护具有重要意义, 近年来越来越受到政府和企业的重视, 也逐渐成为学界研究的热点[1,2,3]。

目前, 我国废弃物回收和处理已初具规模, 但仍存在许多亟待解决的问题, 如:回收责任主体不明确;具有正规资质的回收处理商在整个循环体系中作用不明显;缺乏有力的监管体系;市场不规范等等。

为解决这些问题, 必须借鉴其他发达国家和地区的先进经验, 建立符合中国国情的闭环供应链管理系统。为此, 我国学者从不同角度做了大量的研究工作。达庆利等综述了国内外的研究成果, 探讨了逆向物流系统结构、设计原则和定位问题, 并提出了多个有价值的研究方向[4]。徐家旺等就闭环供应链管理中的网络设计、供应链协调和不确定性等多个方面的问题进行了讨论[5]。

闭环供应链系统具有极强的复杂性和不确定性, 可供研究的问题和角度非常多。目前, 已有不少学者从系统学、运筹学和经济学等角度出发, 采用实证分析、优化理论和博弈分析等方法来研究闭环供应链管理中的各种决策问题, 包括闭环供应链管理体系设计[6,7]、回收网络定位和优化[8]、产品定价决策以及协调机制设计等[9,10]。

在前人的研究基础之上, 首先设计了图1所示的管理体系, 然后运用网络均衡理论和变分不等式方法详细研究闭环供应链系统中各利益相关者的决策行为和均衡条件, 随后建立了随机需求下多商品流闭环供应链均衡模型, 最后通过算例对模型进行了检验。

本文模型的基础理论是变分不等式和网络均衡理论, 它们都已被广泛应用于经济领域的均衡问题研究。 有关供应链网络均衡建模, 国外以Ann Nagurney[11]为首的团队做了许多的工作, 国内较早采用该理论的有滕春贤[12]、徐兵[13]、王志平[14]等学者, 模型的基本框架已经比较成熟, 未来的研究工作将主要集中在推广和应用上。

与已有的文献相比, 本文主要创新点体现在:

①闭环供应链结构完整。供应链往往具有复杂的结构, 而以往大多文献采用的结构都较为简单。文献[11]运用网络均衡理论研究了逆向供应链均衡, 但没有将正、逆向供应链整合起来。文献[15]、文献[16]采用同样方法研究了闭环供应链均衡, 但供应链中只含制造商和零售市场, 而且都没有研究补贴机制的影响。

②考虑了多商品流和随机市场。文献[12]研究了随机需求下的多商品供应链网络均衡问题, 但没有研究逆向供应链的情况。文献[17]、文献[18]采用超网络来描述闭环供应链复杂的结构, 并研究了补贴机制对供应链运作的影响, 然而文中都假设市场中只存在一种商品, 且需求固定。

1 系统描述和假设条件

设计如图1所示的多级闭环供应链系统:制造商处于供应链的中心位置, 负责对原材料进行加工。大部分原材料最终转化为合格商品, 剩余部分是废弃物料 (如边角料、不合格零部件、次品等) , 转卖给回收处理商进行资源再生。

第三方回收处理商聚集于再生资源中心。他们从市场中回收废弃商品, 从企业中回收废弃物料, 然后循环处理, 将再生材料卖给制造商, 废渣则被无害化处理。

管理部门位于供应链外部, 通过法规和政策调控供应链系统的运作。环保基金由管理部门组织建立, 主要用于回收处理费补贴, 其资金主要来源于制造商缴纳的环保税。

为便于下一步建模讨论, 作如下假设:

①再生原材料与初始原材料在使用性能上无差异, 原材料和商品按比率转换, 转换率为1;

②市场中存在多类别的商品, 商品可互相替代, 其差异体现在工艺和品牌上, 每类商品包含的原材料比率相等;

③制造商的生产能力有限, 用α表示原材料使用效率, 0≤α<1;

④回收处理商具备回收网络和处理技术, 循环处理能力有限, 用β表示资源再生比率, 0≤β<1;

⑤供应链中每个行为主体都处于非合作竞争状态, 目标都是追求自身的利益最大化;

⑥模型中出现的生产、销售和交易成本函数均为连续凸函数。

2 多商品闭环供应链网络均衡模型的建立

考虑一个由S个原材料供应商、I个制造商、J个零售商、N个回收处理商和市场组成的闭环供应链, 市场中存在L类商品。下面对各级决策者, 即原材料供应商s、制造商i、零售商j、回收处理商n以及市场的决策行为进行研究, 并推导基于变分不等式的均衡条件。其中, s∈{1, …, S}, i∈{1, …, I}, j∈{1, …, J}, n∈{1, …, N}。

2.1 原材料供应商行为分析及市场均衡模型

供应链中所有原料供应商生产同一种原材料。用qs表示供应商s的原材料产量, qsi表示供应商s卖给制造商i的原材料数量, 两者满足产销平衡条件$q{}_s={\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{si}$; QS是由所有qsi组成的向量;双方的交易价格为ρsi;供应商的生产成本与自身产量相关, 记为cs (qs) , 交易成本 (含运输成本) 则与交易量相关, 记为$\widehat{c}{}_{si}(q{}_{si})$。

此时, 可建立供应商s的最优化模型:

$\begin{array}{l}\max R{}_s={\displaystyle \sum _{i=1}^{{\rm I}}\rho }{}_{si}q{}_{si}-c{}_s(q{}_s)-{\displaystyle \sum _{i=1}^{{\rm I}}\widehat{c}}{}_s(q{}_{si})\\ \text{s}.\text{t}.q{}_{si}\ge 0‚i=1,\cdots ,{\rm I}\end{array}$

由假设条件⑤, 原材料供应市场均衡表现为供应商之间非合作竞争的Nash均衡。根据假设条件⑥, 生产成本函数和交易成本函数均为连续凸函数, 原材料供应市场的Nash均衡等价于如下变分不等式:

$\begin{array}{l}{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{s=1}^S\langle }}\frac{\partial c{}_s(q{}_{si}^{*})}{\partial q{}_{si}}+\frac{\partial \widehat{c}{}_s(q{}_{si}^{*})}{\partial q{}_{si}}-\rho {}_{si}^{*},q{}_{si}-q{}_{si}^{*}\rangle \ge 0,\\ \forall Q{}^S\in R{}_{+}^{S{\rm I}}(1)\end{array}$

其内在经济意义是:只有当制造商支付的单位价格等于供应商的边际成本时, 交易才会发生;否则, 双方交易量为零。

2.2 制造商行为分析及均衡模型

制造商都具备生产L类商品的能力。用qni表示制造商i从回收处理商n处采购的原材料数量, qwin表示i卖给n的废弃物料量;qli表示i生产的第l类商品的产量, qlij表示制造商i卖给零售商j的商品l的数量, 满足产销平衡条件qli=${\displaystyle \sum _{j=1}^{J}q}{}_{ij}^l$;QL1、QW、QN分别是由所有qlij、qwin、qni组成的向量;再用ρlij表示制造商i卖给零售商j的商品l的批发价格, ρin表示制造商i卖给回收处理商n的废弃物料的价格, ρni则表示制造商i从回收处理商n处购买原材料的价格。令αi为衡量制造商i的“绿色”生产能力和管理水平的指标; τl为制造商为单位商品l缴纳的环保税; 制造商与零售商、供应商、回收处理商的交易成本函数分别记作$\widehat{c}{}_{ij}^l(q{}_{ij}^l),\widehat{c}{}_{is}(q{}_{si}),\widehat{c}{}_{in}(q{}_{ni},q{}_{in}^w)$, 生产成本函数为cli (qli) 。

此时, 可建立制造商i的最优化模型:

$\begin{array}{l}\max R{}_i={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\rho }}{}_{ij}^{l}q{}_{ij}^l+{\displaystyle \sum _{n=1}^{{\rm N}}\rho }{}_{in}^{w}q{}_{in}^w-{\displaystyle \sum _{l=1}^{L}c}{}_i^l(q)\\ -{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\widehat{c}}}{}_{ij}^l(q{}_{ij}^l)-{\displaystyle \sum _{s=1}^S\widehat{c}}{}_{is}(q{}_{si})-{\displaystyle \sum _{n=1}^{{\rm N}}\widehat{c}}{}_{in}(q{}_{ni},q{}_{in}^w)\\ -{\displaystyle \sum _{s=1}^S\rho }{}_{si}q{}_{si}-{\displaystyle \sum _{n=1}^{{\rm N}}\rho }{}_{ni}q{}_{ni}-{\displaystyle \sum _{l=1}^L\tau }{}^{l}q{}_i^l\\ \text{s}.\text{t}.{\displaystyle \sum _{l=1}^{L}q}{}_i^l\le \alpha {}_i({\displaystyle \sum _{s=1}^{S}q}{}_{si}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni})(2)\\ {\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{in}^w\le (1-\alpha {}_i)({\displaystyle \sum _s^{S}q}{}_{si}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni})(3)\\ q{}_{si}\ge 0,q{}_{ni}\ge 0,q{}_{in}^w\ge 0,q{}_{ij}^l\ge 0,\forall s,i,j,n,l\end{array}$

约束 (2) 表示商品产量不多于最大可能总产量, 与之对应的拉格朗日乘子为λ1i, 向量λ1= (λ11, …, λ1i, …, λ1I) ;约束 (3) 说明卖出的废弃物料量不超过总废弃量, 对应的拉格朗日乘子为λ2i, 向量λ2= (λ21, …, λ2i, …, λ2I) 。

与2.1节同理, 在非合作竞争情形下, 制造商竞争达到均衡等价于如下变分不等式成立:

$\begin{array}{l}\begin{array}{l}{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\langle }}}\frac{\partial c{}_i^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\frac{\partial \widehat{c}{}_{ij}^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\tau {}^l+\lambda {}_{1i}^{*}-\rho {}_{ij}^{l*},q{}_{ij}^l-q{}_{ij}^{l*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{n=1}^{{\rm N}}\langle }}\frac{\partial \widehat{c}{}_{in}(q{}_{in}^{w*})}{\partial q{}_{in}^w}+\lambda {}_{2i}^{*}-\rho {}_{in}^{w*},q{}_{in}^w-q{}_{in}^{w*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{s=1}^S\langle }}\frac{\partial \widehat{c}{}_{is}(q{}_{si}^{*})}{\partial q{}_{si}}+\rho {}_{si}^{*}-\alpha {}_i\lambda {}_{1i}^{*}-(1-\alpha {}_i)\lambda {}_{2i}^{*},q{}_{si}-q{}_{si}^{*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{n=1}^{{\rm N}}\langle }}\frac{\partial \widehat{c}{}_{in}(q{}_{ni}^{*})}{\partial q{}_{ni}}+\rho {}_{ni}^{*}-\alpha {}_i\lambda {}_{1i}^{*}-(1-\alpha {}_i)\lambda {}_{2i}^{*},q{}_{ni}-q{}_{ni}^{*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}\langle }\alpha {}_i({\displaystyle \sum _{s=1}^{S}q}{}_{si}^{*}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni}^{*})-{\displaystyle \sum _{l=1}^{L}q}{}_i^{l*},\lambda {}_{1i}-\lambda {}_{1i}^{*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}\langle }(1-\alpha {}_i)({\displaystyle \sum _{s=1}^{S}q}{}_{si}^{*}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni}^{*})-{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{in}^{w*},\lambda {}_{2i}-\lambda {}_{2i}^{*}\rangle \ge 0,\end{array}\\ \forall (Q{}_1^L,Q{}^w,Q{}^S,Q{}^{{\rm N}},\lambda {}_1,\lambda {}_2)\in R{}_{+}^{{\rm I}J+{\rm I}{\rm N}+S{\rm I}+{\rm N}{\rm I}+2{\rm I}}(4)\end{array}$

其内在经济意义为:当商品批发价格等于边际成本时, 商品交易才会发生; 否则, 商品交易量为0; 当废弃物料的交易价格等于边际成本, 制造商才会卖出废弃物料。

2.3 零售商行为分析及均衡模型

零售商根据市场需求来决定商品价格和订货量等。用clj (qj) 表示零售商j的销售成本 (含商品的展示和库存成本) , 其中qj为j的商品总量;零售商j与制造商i的交易成本记为$\widehat{c}{}_{ji}^l(q{}_{ij}^l)$。

零售商面对的市场需求dlj (ρl1j) 为随机变量, 表示当零售商j对商品l制定的价格为ρl1j时, 市场对零售商j的反映。dlj (ρl1j) 的密度函数为flj (x, ρl1j) , 概率分布函数Flj (x, ρl1j) =F (dlj≤x) =∫x0flj (x, ρl1j) dx.

记$s{}_j^l={\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l$为零售商j从所有制造商处批发的商品l的总量。显然, 销售量不可能超过供应量和市场需求量, 即不大于min{slj, dlj}。

令Δl+j=max{0, slj-dlj}为零售商存货量大于销售量的那部分商品量, ηl+j表示商品剩余惩罚系数, ηl+j≥0; Δl-j=max{0, dlj-slj}为存货量小于销售量的那部分商品量, ηl-j表示缺货惩罚系数, ηl-j≥0。分别用el+j和el-j表示Δl+j和Δl-j的期望值, el+j=∫slj0 (slj-x) flj (x, ρl1j) dx, el-j=∫+∞slj (x-slj) flj (x, ρl1j) dx, 故总惩罚期望为:

$E(\eta {}_j^{l+}\Delta {}_j^{l+}+\eta {}_j^{l-}\Delta {}_j^{l-})=\eta {}_j^{l+}e{}_j^{l+}(s{}_j^l\rho {}_{1j}^l)+\eta {}_j^{l-}e{}_j^{l-}(s{}_j^l\rho {}_{1j}^l)$

此时, 可建立零售商j的最优化模型:

$\begin{array}{l}\max R{}_j={\displaystyle \sum _{l=1}^L\{}E(\rho {}_j^l\cdot \min \{s{}_j^l,d{}_j^l\})-E(\eta {}_j^{l+}\Delta {}_j^{l+}+\eta {}_j^{l-}\Delta {}_j^{l-})\\ -c{}_j^l(q{}_j)-{\displaystyle \sum _{i=1}^{{\rm I}}\widehat{c}}{}_{ji}^l(q{}_{ij}^l)-{\displaystyle \sum _{i=1}^{{\rm I}}\rho }{}_{ij}^{l}q{}_{ij}^l\}\\ \text{s}.\text{t}.q{}_{ij}^l\ge 0,\forall i,j,l\end{array}$

同理, 零售商竞争达到均衡等价于如下变分不等式成立:

$\begin{array}{l}\begin{array}{l}{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\langle }}}\eta {}_j^{l+}F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*})\\ -(\eta {}_j^{l-}+\rho {}_{1j}^{l*})(1-F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*}))\\ +\frac{\partial c{}_j^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\frac{\partial \widehat{c}{}_{ji}^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\rho {}_{ij}^{l*},q{}_{ij}^l-q{}_{ij}^{l*}\rangle \ge 0,\end{array}\\ \forall Q{}_1^L\in R{}_{+}^{{\rm I}JL}(5)\end{array}$

其内在经济意义为:当有市场交易时, 零售商j的销售价格以$1-F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*})$的概率等于以下四项之和, 分别是: 剩货成本 (概率为$F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*}))$、 缺货成本 (概率为$1-F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*}))$、边际成本 (包括销售和交易成本) 、零售商批发商品的价格。

2.4 顾客选择行为分析及市场均衡模型

顾客的选择行为能直接影响到企业的决策。 在正向供应链中, 顾客选择是否购买商品。在随机市场假设下, 由文[12]和文[19], 零售商和消费者之间发生交易的均衡条件是随机经济均衡条件, 可用下式表示:

$d{}_j^l(\rho {}_{1j}^{l*})\{\begin{array}{l}\le {\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^{l*}a.e.\rho {}_{1j}^{l*}=0\\ ={\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^{l*}a.e.\rho {}_{1j}^{l*}>0\end{array}(6)$

其中, a.e.表示几乎处处成立。

而在逆向供应链中, 顾客选择行为发生在废旧商品回收市场, 他们选择是否卖出自己的废旧商品。由于卖出废旧商品给顾客自身带来的收益非常有限, 而即使闲置也并不意味着损失, 所以是否卖出废旧商品很大程度上是一种个人心理选择行为。故假设, 只有当回收价格达到顾客的预期值时, 交易才会发生。即回收处理商n从市场中回收的商品l的数量ql2n取决于回收处理商n给出的回收价格ρl2n.用单调递增函数a (q) 刻画回收市场的顾客选择行为:

$a{}^l(q{}_{2n}^l)\{\begin{array}{l}\ge \rho {}_{2n}^{l*}‚q{}_{2n}^{l*}=0\\ =\rho {}_{2n}^{l*}‚q{}_{2n}^{l*}>0\end{array}(7)$

此外, 回收量不会超过市场的销售量, 即:

${\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{2n}^l\le {\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^{J}q}}{}_{ij}^l,\forall i,j,n,l(8)$

综合式 (6) 、式 (7) 、式 (8) , 零售和回收市场都达到均衡等价于如下变分不等式成立:

$\begin{array}{l}\begin{array}{l}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\langle }}{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l-d{}_j^l(\rho {}_{1j}^{l*}),\rho {}_{1j}^l-\rho {}_{1j}^{l*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{l=1}^L\langle }}a{}^l(q{}_{2n}^{l*})-\rho {}_{2n}^{l*}+\lambda {}_{3l}^{*},q{}_{2n}^l-q{}_{2n}^{l*}\rangle \\ +{\displaystyle \sum _{l=1}^L\langle }{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^{J}q}}{}_{ij}^{l*}-{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{2n}^{l*},\lambda {}_{3l}-\lambda {}_{3l}^{*}\rangle \ge 0,\end{array}\\ \forall (Q{}_2^L,\lambda {}_3,\rho )\in R{}_{+}^{{\rm N}L+L+JL}(9)\end{array}$

其中, QL2、ρ分别是由所有ql2n、ρl1j组成的向量; 向量λ3= (λ31, …, λ3l, …, λ3L) , λ3l是对应式 (8) 的拉格朗日乘子。

2.5 回收处理商行为分析及均衡模型

用βn表示回收处理商n的处理能力;无害化处理的单位成本为$\overline{c}$;环保基金给予废旧商品l的补贴为sl;废弃物处理成本与处理的总量相关, 记为cn (q) , 与制造商交易的成本记为$\widehat{c}{}_{ni}^l(q{}_{ni},q{}_{in}^w)$, 市场回收成本记为$\widehat{c}{}_n^l(q{}_{2n}^l)$。

此时, 可建立回收处理商n的最优化模型:

$\begin{array}{l}\max R{}_n={\displaystyle \sum _{i=1}^{{\rm I}}\rho }{}_{ni}q{}_{ni}+{\displaystyle \sum _{l=1}^{L}s}{}^{l}q{}_{2n}^l-{\displaystyle \sum _{i=1}^{{\rm I}}\rho }{}_{in}^{w}q{}_{in}^w-{\displaystyle \sum _{l=1}^L\rho }{}_{2n}^{l}q{}_{2n}^l\\ -c{}_n(q)-\overline{c}(1-\beta {}_n)({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{in}^w+{\displaystyle \sum _{l=1}^{L}q}{}_{2n}^l)\\ -{\displaystyle \sum _{i=1}^{{\rm I}}\widehat{c}}{}_{ni}(q{}_{ni,}q{}_{ni}^w)-{\displaystyle \sum _{l=1}^{L}c}{}_n^l(q{}_{2n}^l)\\ \text{s}.\text{t}.{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ni}\le \beta {}_n({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{in}^w+{\displaystyle \sum _{l=1}^{L}q}{}_{2n}^l)(10)\\ q{}_{ni}\ge 0‚q{}_{2n}^l\ge 0‚q{}_{in}^w\ge 0‚\forall i,n,l\end{array}$

式 (10) 说明卖掉的再生材料总量不超过再生处理所得到的总量, λ4n是与之对应的拉格朗日乘子, 向量λ4= (λ41, …, λ4l, …, λ4L) 。当条件⑥成立时, 回收处理商竞争达到均衡等价于如下变分不等式成立:

$\begin{array}{l}\begin{array}{l}{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{i=1}^{{\rm I}}\langle }}\frac{\partial \widehat{c}{}_{ni}(q{}_{ni}^{*})}{\partial q{}_{ni}}+\lambda {}_{4n}^{*}-\rho {}_{ni}^{*},q{}_{ni}-q{}_{ni}^{*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{i=1}^{{\rm I}}\langle }}\frac{\partial c{}_n(q{}_{in}^{w*})}{\partial q{}_{in}^w}+\frac{\partial \widehat{c}{}_{ni}(q{}_{in}^{w*})}{\partial q{}_{in}^w}+\overline{c}(1-\beta {}_n)+\rho {}_{in}^{w*}\\ -\beta {}_n\lambda {}_{4n}^{*},q{}_{in}^w-q{}_{in}^{w*}\rangle +{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{l=1}^L\langle }}\frac{\partial c{}_n(q{}_{2n}^{l*})}{\partial q{}_{2n}^l}+\rho {}_{2n}^{l*}\\ +\overline{c}(1-\beta {}_n)-s{}^{l*}-\beta {}_n\lambda {}_{4n}^{*},q{}_{2n}^l-q{}_{2n}^{l*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}\langle }\beta {}_n({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{in}^{w*}+{\displaystyle \sum _{l=1}^{L}q}{}_{2n}^{l*})-{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ni}^{*},\lambda {}_{4n}-\lambda {}_{4n}^{*}\rangle \ge 0,\end{array}\\ \forall (Q{}^{{\rm N}},Q{}^W,Q{}_2^L,\lambda {}_4)\in R{}_{+}^{{\rm N}{\rm I}+{\rm I}{\rm N}+{\rm N}L+{\rm N}}(11)\end{array}$

其内在经济意义是:只有当再生资源的价格等于加工成本和潜在机会成本之和时, 回收处理商才会卖出再生材料;否则, 交易量为0。

2.6 闭环供应链均衡模型

在市场机制调节下, 闭环供应链将达到均衡。此时, 原材料供应商、制造商、零售商、市场和回收处理商分别满足变分不等式 (1) 、 (4) 、 (5) 、 (9) 、 (11) 。

因此, 闭环供应链网络均衡状态表现为 (1) 、 (4) 、 (5) 、 (9) 、 (11) 之和, 即变分不等式 (12) 成立:

$\begin{array}{l}\begin{array}{l}{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{s=1}^S\langle }}\frac{\partial c{}_s^r(q{}_{si}^{r*})}{\partial q{}_{si}^r}+\frac{\partial \widehat{c}{}_s^r(q{}_{si}^{r*})}{\partial q{}_{si}^r}+\frac{\partial \widehat{c}{}_{is}(q{}_{si}^{*})}{\partial q{}_{si}}-\alpha {}_i\lambda {}_{1i}^{*}\\ -(1-\alpha {}_i)\lambda {}_{2i}^{*},q{}_{si}^r-q{}_{si}^{r*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\langle }}}\frac{\partial c{}_i^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\frac{\partial \widehat{c}{}_{ij}^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\frac{\partial c{}_j^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}+\frac{\partial \widehat{c}{}_{ji}^l(q{}_{ij}^{l*})}{\partial q{}_{ij}^l}\\ +\eta {}_j^{l+}F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*})-(\eta {}_j^{l-}+\rho {}_{1j}^{l*})(1-F{}_j^l({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l,\rho {}_{1j}^{l*}))\\ +\tau {}^l+\lambda {}_{1i}^{*},q{}_{ij}^l-q{}_{ij}^{l*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{i=1}^{{\rm I}}\langle }}\frac{\partial c{}_n(q{}_{in}^{w*})}{\partial q{}_{in}^w}+\frac{\partial \widehat{c}{}_{ni}(q{}_{in}^{w*})}{\partial q{}_{in}^w}+\frac{\partial \widehat{c}{}_{in}(q{}_{in}^{w*})}{\partial q{}_{in}^w}\\ +\overline{c}(1-\beta {}_n)+\lambda {}_{2i}^{*}-\beta {}_n\lambda {}_{4n}^{*},q{}_{in}^w-q{}_{in}^{w*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{n=1}^{{\rm N}}\langle }}\frac{\partial \widehat{c}{}_{in}(q{}_{ni}^{*})}{\partial q{}_{ni}}+\frac{\partial \widehat{c}{}_{ni}(q{}_{ni}^{*})}{\partial q{}_{ni}}+\lambda {}_{4n}^{*}-\alpha {}_i\lambda {}_{1i}^{*}\\ -(1-\alpha {}_i)\lambda {}_{2i}^{*},q{}_{ni}-q{}_{ni}^{*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}{\displaystyle \sum _{l=1}^L\langle }}\frac{\partial c{}_n(q{}_{2n}^{l*})}{\partial q{}_{2n}^l}+a{}^l(q{}_{2n}^{l*})+c(1-\beta {}_n)+\lambda {}_{3l}^{*}-s{}^{l*}\\ -\beta {}_n\lambda {}_{4n}^{*},q{}_{2n}^l-q{}_{2n}^{l*}\rangle \\ +{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{l=1}^L\langle }}{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ij}^l-d{}_j^l(\rho {}_{1j}^{l*}),\rho {}_{1j}^l-\rho {}_{1j}^{l*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}\langle }\alpha {}_i({\displaystyle \sum _{s=1}^{S}q}{}_{si}^{*}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni}^{*})-{\displaystyle \sum _{l=1}^{L}q}{}_i^{l*},\lambda {}_{1i}-\lambda {}_{1i}^{*}\rangle \\ +{\displaystyle \sum _{i=1}^{{\rm I}}\langle }(1-\alpha {}_i)({\displaystyle \sum _{s=1}^{S}q}{}_{si}^{*}+{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{ni}^{*})-{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{in}^{w*},\lambda {}_{2i}-\lambda {}_{2i}^{*}\rangle \\ +{\displaystyle \sum _{l=1}^L\langle }{\displaystyle \sum _{i=1}^{{\rm I}}{\displaystyle \sum _{j=1}^{J}q}}{}_{ij}^{l*}-{\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{2n}^{l*},\lambda {}_{3l}-\lambda {}_{3l}^{*}\rangle \\ +{\displaystyle \sum _{n=1}^{{\rm N}}\langle }\beta {}_n({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{in}^{w*}+{\displaystyle \sum _{l=1}^{L}q}{}_{2n}^{l*})-{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{ni}^{*},\lambda {}_{4n}-\lambda {}_{4n}^{*}\rangle \ge 0,\end{array}\\ \forall (Q{}^S,Q{}_1^L,Q{}^{{\rm N}},Q{}^W,Q{}_2^L,\rho ,\lambda {}_1,\lambda {}_2,\lambda {}_3,\lambda {}_4)\\ \in R{}_{+}^{S{\rm I}+{\rm I}J+{\rm N}{\rm I}+{\rm I}{\rm N}+{\rm N}L+JL+2{\rm I}+L+{\rm N}}(12)\end{array}$

为便于讨论, 将式 (12) 简写为〈G (X*) , X-X*〉≥0, X= (QS, QL1, QN, QW, QL2, ρ, λ1, λ2, λ3, λ4) 。分析式 (12) 可知, 由于各函数项都是连续函数, 因此只要可行集有界, 该均衡问题就有解。若各函数项也都是凸函数, 则式 (12) 有解且唯一。

3 算例及分析

考虑一个由1个原材料供应商、2个制造商、两个零售商、1个回收处理商组成的闭环供应链, 市场中有2种商品。为简便起见, 在不影响模型性质的前提下, 假设所有交易成本为0。

定义相关函数如下:

(1) 原料供应商的生产函数:

$c{}_s(q{}_s)=2.5({\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{si}){}^2+{\displaystyle \sum _{i=1}^{{\rm I}}q}{}_{si}+2$

(2) 制造商的生产函数:

$c{}_i^1=2.5({\displaystyle \sum _{j=1}^{J}q}{}_{ij}^1){}^2+{\displaystyle \sum _{j=1}^{J}q}{}_{ij}^1({\displaystyle \sum _{j=1}^{J}q}{}_{i^{\prime }j}^1+2)$,

$c{}_i^2=3({\displaystyle \sum _{j=1}^{J}q}{}_{ij}^1){}^2+{\displaystyle \sum _{j=1}^{J}q}{}_{ij}^1({\displaystyle \sum _{j=1}^{J}q}{}_{i^{\prime }j}^1+2)$,

i′∈{1, …, I}, i′≠i;

(3) 零售商的销售成本函数:

$\begin{array}{l}c{}_1^l=\frac{1}{2}({\displaystyle \sum _{i=1}^2{\displaystyle \sum _{l=1}^{2}q}}{}_{i1}^l){}^2\\ c{}_2^l=({\displaystyle \sum _{i=1}^2{\displaystyle \sum _{l=1}^{2}q}}{}_{i2}^l){}^2\end{array}$

再令零售商j的需求函数dlj (ρl1j) 在[0, bj/ρl1j]内服从均匀分布, 则有:

$f{}_j^l(x,\rho {}_{1j}^l)=\{\begin{array}{cc}\rho {}_{1j}^l/b{}_j,& 0\le x\le b{}_j/\rho {}_{1j}^l\\ 0,& \text{其}\text{他}\end{array}$

$F{}_j^l(x,\rho {}_{1j}^l)=\{\begin{array}{cc}x\rho {}_{1j}^l/b{}_j,& 0\le x<b{}_j/\rho {}_{1j}^l\\ 1,& b{}_j/\rho {}_{1j}^l\le x\end{array}$

dlj (ρl1j) =E (dlj) =bj/2ρl1j

(4) 市场回收偏好函数为:

$a{}^1=\frac{1}{2}({\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{2n}^1)+5$

$a{}^2=\frac{1}{2}({\displaystyle \sum _{n=1}^{{\rm N}}q}{}_{2n}^2)+3$

(5) 回收处理商的回收处理成本函数为:

$c{}_n=0.1({\displaystyle \sum _{i=1}^{2}q}{}_{in}^w+{\displaystyle \sum _{l=1}^{2}q}{}_{2n}^l)+3$

由于资源和能力都是有限的, 再根据以上函数的性质可知, 均衡问题 (12) 有唯一解。当设置参数为bj=10, ηl+j=ηl-j=1, c=0.2, r1=r2=5, α1=95%, α2=90%, θ1=θ2=80%, β=90%时, 求解闭环供应链均衡模型, 并针对两种补贴政策进行灵敏度分析。即统一补贴和差异补贴:统一补贴是指对不同类的商品给予相同补贴, 差异补贴是指对不同类商品给予不同补贴。

(1) 当提供统一补贴时, 发现补贴额与回收率同时递增。

而当补贴增加到大致3.5位置时, 回收率均超过60%, 达到4以后, 回收率均超过80% (见图2) , 实现了管理目标。

(2) 当提供差异补贴时, 发现不同类商品在回收市场存在竞争现象。

固定商品1的补贴, 逐渐增加对商品2的补贴, 得到图3。

上图表明, 在差异补贴政策下, 商品2的回收率随着补贴s2的增加而上升, 商品1的回收率则随之下降, 但幅度不大, 且始终保持在80%以上。这说明废旧商品的回收率也受到其他类商品补贴的影响, 因此它们在回收市场中也存在竞争现象。原因在于, 各类商品的市场回收偏好函数不同, 且回收处理成本也不同。

从供应链整体来看, 原材料总产量持续减少, 这是因为, 当补贴额增加后, 回收处理商在回收市场中就有能力提供更高的回收价格来激励顾客卖出废旧商品, 从而获得更多的再生原材料。均衡价格表明, 再生材料比原始材料更有价格优势, 进而造成供应商处的原材料总供应量减少。与此同时, 商品总供应量和生产量增加了, 这说明供应链整体运行良好。由于篇幅关系, 仅用图4反映上述情况。

总的来讲, 补贴能提高商品的回收率, 但要注意到商品在回收市场的竞争现象。因此在制定政策时, 不能仅仅考虑补贴的额度, 还要重视补贴在不同商品间的分配问题。

此外, 本算例假设β=0.9, 虽然目前技术能力还没有达到这个水平, 但是均衡解说明, 当资源再利用能力提高后, 可以明显减少自然资源消耗, 实现资源的充分循环利用, 有利于社会可持续发展。

4 结论

以具有多类商品流和随机市场的闭环供应链为背景, 通过建立网络均衡模型探讨了各决策者的行为与均衡条件, 分析了废旧商品回收处理的补贴机制设计问题, 得到了具有实际意义的结论。算例结果说明:补贴政策对闭环供应链系统的影响显著;异质商品间存在明显竞争;积极的回收政策可以提高资源循环利用率。

本文可为政府和企业提供决策参考, 但还存在许多不足之处。未来可开展的研究有:结合实际案例进行实证研究, 以检验算例得到的结论, 并完善理论模型; 采用多层规划研究上层为政府的最优政策设计问题; 供应链系统运行效率和政策效率评价问题。

摘要：针对包含供应商、制造商、零售商、顾客和回收处理商的闭环供应链, 在多商品流和随机需求假设下, 利用网络均衡理论, 分析了每类决策者的行为, 得到了基于变分不等式的闭环供应链均衡模型, 并给出了经济解释。对补贴政策进行了灵敏度分析, 实验表明:当补贴额增加时, 回收率也相应增加;异质商品在销售和回收过程中都存在竞争。