蒙特卡洛模拟方法

2024-09-20

蒙特卡洛模拟方法（共7篇）

蒙特卡洛模拟方法篇1

0 引言

科学技术的进步促使机床朝着高速和高精度方向发展，机床精度性能的设计与提高变得日益重要。国内外学者在机床精度设计领域开展了大量的研究，取得了一定的进展。文献[1,2,3]运用多体运动学理论建立了机床的误差传递模型；文献[4,5]分析了加载时工作零件的变形量，将其转变为雅可比旋量修正量，通过对雅可比旋量公差模型在实际工况下的修正，建立了基于雅可比旋量和实际工况的装配体公差数学模型；文献[6]给出了三维公差累积运动学模型的一般表示式，研究了该模型在公差优化设计中的应用。但上述研究都将机床零部件间的几何变量约束统一用6项公差变动要素（3项平动误差与3项转动误差）进行描述，而公差变动要素与零部件公差之间的关系并不明确，无法建立反映零部件公差与机床精度映射关系的数学模型，也就无法有效地对机床零部件公差进行合理分配，只能根据并不完善的公差模型来进行公差分配，可靠性低。因此，建立能用于工程实际且定义严格的公差模型，对机床的精度分析与分配研究具有重要意义。

公差建模的关键是可以对满足公差的公差变动要素作出正确的解释[7]，即公差变动要素如何在公差域中变动。文献[8,9,10]运用基于数学定义的公差建模方法，用变动不等式和约束不等式严谨地描述了公差变动要素与公差间的关系，建立了不同类型公差的数学模型，但由于约束不等式中变动要素存在变动顺序不确定性，难以得到公差变动要素实际变动区间与公差间的具体函数关系，因此该模型无法用于直接指导精度设计。

本文以平面尺寸公差与平面度公差为研究对象，在深入研究基于数学定义的公差分析方法的基础上，将蒙特卡洛模拟与响应面方法应用于公差建模，考虑公差原则及约束条件等因素对公差变动要素的影响，建立变动要素实际变动区间带宽与公差间的响应面模型，以解决现有精度设计方法实用性不强的问题，并通过典型实例验证该方法的有效性。

1 基于小位移旋量的公差数学表示

三维空间中，物体表面可抽象为点、线、面等基本要素，并有3个沿坐标轴方向的平动自由度和3个绕坐标轴的转动自由度。用d=(u,v,w）和θ=（α，β，δ）分别表示平动与转动自由度的微小变动矢量，则这两组矢量组成的合成矢量D=（θ，d）=（α，β，δ，u,v,w）称为小位移旋量（small displacement torsors,SDT），α、β、δ、u、v、w为SDT的旋量参数。

公差域表示公差实际表面脱离名义表面变动的范围或区域，而物体表面的点、线、面等基本要素的变动量均可以通过各自的SDT表示。据此，Bourdet在1996年首次将SDT引入到公差领域，提出了基于SDT的公差数学表示方法[11]，即用旋量参数的变动来描述公差实际表面在公差域中的变动。旋量参数也称为公差变动要素。表1列出了几种形位公差对应的公差域和SDT。

2 基于蒙特卡洛模拟与响应面方法的公差建模

公差分为尺寸公差、形状公差和位置公差。针对尺寸公差和位置公差，分析时可以假设变动后公差实际表面的形状保持不变，即平面变动后依旧为平面，直线变动后依旧为直线，但对于形状公差，这一假设并不成立[12]。为此，本文从两类公差中各取一例，详细阐述基于蒙特卡洛模拟与响应面方法的公差建模步骤，其他类型的公差建模可依此类推。

2.1 平面尺寸公差建模

2.1.1 求解公差变动不等式与约束不等式

平面尺寸公差规定的公差域形式如图1所示。图中，TU、TL分别为上下偏差，公差T=TU+TL;D为基本尺寸；局部坐标系的z轴与平面法向平行。对于矩形平面，坐标系原点在平面中心位置；对于不规则平面，可用一等效矩形替代，则原点在等效矩形的中心位置。z=0表示名义公差平面；z(x,y）表示实际变动平面。由表1可知，变动平面z(x,y）对应的SDT为（α，β，0,0,0,w)T。尺寸公差用数学表示如下：

其中，z(x,y)=dz+xβ+yα为实际变动平面方程，dz为实际变动平面中心点与参考坐标系原点在z方向上的位移；Δz(x,y）为变动平面z(x,y）上任意两点的z坐标值之差。由于变动后公差实际表面的形状保持不变，仍为矩形平面，极值情况必然发生在平面的4个顶点处，因此只需考虑变动平面顶点的z向坐标值是否在公差域内即可。

由图1可知，4个顶点的x、y向坐标分别为（a,b）、（-a,b）、（a,-b）、（-a,-b），由实际变动平面方程可知，4个顶点的坐标分别为

由式（1）可得

由式（2）、式（3）可得旋量参数的变动不等式[8]为

约束不等式[8]为

其中，x、y分别取4个顶点的x、y向坐标值。

变动不等式式（4）定义了旋量参数的理想变动区间。约束不等式的存在使旋量参数不可能同时取最大值，参数的实际变动区间比理想变动区间要小，如果按旋量参数的理想变动区间进行精度设计，则会提高零件加工精度的要求，增加制造成本。而旋量参数实际变动区间与公差间的具体关系尚不明确。

2.1.2 蒙特卡洛模拟法求解旋量参数实际变动区间

蒙特卡洛模拟法也称为随机模拟方法，它以概率论和数理统计为基础，通过对随机变量的统计实验和随机模拟来求解问题近似解。由于多种因素综合作用的影响，产品公差实际表面变动的形式在公差域内是无法预知的，具有随机性，其分布的规律随加工条件而定，因此，与公差对应的SDT旋量参数也是随机量[13]。由式（4）～式（6）可知，旋量参数的变动顺序不同，其最终的取值也会不同。本文根据变动不等式和约束不等式，运用蒙特卡洛模拟法，按给定的参数分布类型和变动顺序进行模拟试验，生成相应的参数随机数，模拟实际公差表面的变动，对满足约束条件的随机数加以保留，不满足约束条件的则剔除，当产生足够多满足条件的随机数样本后，对样本进行分析，求解出对应旋量参数的变动区间带宽，再根据变动区间带宽和均值得到旋量参数的实际变动区间。求解步骤如下：

(1）根据研究对象的具体情况，确定旋量参数的理想概率分布模型。本文假定旋量参数的分布函数均为正态分布：

式中，μ、σ分别为正态分布的均值和均方差。

(2）根据式（4）确定各参数理想分布的均值与方差。平面尺寸公差SDT中非零的旋量参数为α、β、w，由于正态分布曲线的有效分布范围为6σ（变量位于区间（-3σ，3σ）内的概率达99.73%），且分布的中心与变动区间的中心重合，因此由式（4）可知，各旋量参数理想分布的均值与方差分别为

(3）根据抽样规则对旋量参数进行抽样。3个非零旋量参数共有6种变动顺序，为了尽可能地模拟零件的实际加工工况，对每一种变动顺序都进行抽样。以变动顺序α→β→w为例，其抽样流程如图2所示。

图2中，k的初始值为1;K为要求抽取的合格样本数，为保证模拟精度，K>20 000;K1、K2为常量，根据具体研究对象而定。设置判别式k1<K1、k2≤K2的目的是为了防止前面的参数抽样值过大，导致后面的参数无论取何值都无法满足约束方程，进入死循环。由式（4）～式（6）可知，旋量参数实际变动区间带宽只与公差值有关，而与上下偏差的大小无关。因此，本文在抽样的过程中令TL=TU。

按6种变动顺序编程进行抽样后，α、β、w均得到容量为6 K的实际变动区间带宽样本。

(4）旋量参数实际分布函数的假设检验。由于约束不等式式（5）、式（6）的限制，旋量参数的实际分布类型并不一定与理想分布类型相同，故需对实际分布类型进行假设检验。本文采用χ3拟合检验法，以参数α为例，其基本思路为：在α分布未知时，先根据抽样获得的观测值对α的分布类型作出假设：

式中，F0（α）为假设的α分布函数。

将实数轴分为k个互不相交的区间（ai,ai+1](i=1,2，…，k），其中a1、ak+1可分别取-∞、+∞，区间的划分视具体情况而定。在假设H0下计算概率：

计算理论频数。由于α的样本容量为6 K，因此理论频数为6 Kpi。

计算样本观察值落入区间（ai,ai+1]中的个数fi(i=1,2，…，k),fi称为实际频数。

计算统计量

当样本容量6 K≥50时，统计量χ2近似服从自由度为k-r-1的χ2分布（r为确定F0（α）表达式中被估计参数的个数，如F0（α）为正态分布，则r=2）。

对于给定的显著性水平ε（本文取ε=0.05），确定临界值χε2(k-r-1），如果χ2>χε2(k-r-1），则拒绝假设H0，如果χ2≤χε2(k-r-1），则接受假设H0。

(5）采用极大似然估计法估计实际分布函数的参数。如果经第（4）步假设检验验证α也属于正态分布N（μ，σ2），则根据α的样本，采用极大似然估计法估计α分布函数的参数，得到μ与σ2的极大似然估计量为

式中，为α样本的均值。

由于变动不等式与约束不等式均具有对称性，故旋量参数实际分布的均值近似等于理想分布的均值。

(6）求解旋量参数的实际变动区间。根据参数的实际分布函数类型与均方差，查表得到参数的实际变动区间带宽：

其中，G表示相对分布系数，如果旋量参数为正态分布，则G等于1。旋量参数的实际变动区间为

2.1.3 运用响应面法建立旋量参数实际变动区间带宽与公差间的响应面函数

响应面方法是数学方法与统计方法相结合的一种建模方法，主要用于解决如何建立复杂系统输入（变量）与输出（响应）之间近似函数关系的问题。其基本思路是，依据试验设计原则选定设计空间中的试验点，用多组试验点及其对应的响应值拟合响应面曲面，从而构造具有明确表达形式的显式函数[14,15]。

通过前文分析可知，只要确定了旋量参数实际变动区间带宽即可确定参数的实际变动区间。对于任意的平面尺寸公差均可通过蒙特卡洛模拟法得到与之对应的旋量参数实际变动区间带宽的响应值，但两者间的函数关系并不明确，无法用严格的数学公式表达。为此，本文采用响应面方法，建立公差与旋量参数实际变动区间带宽间的响应面模型，进而确定公差与旋量参数实际变动区间之间的响应关系，为后续的精度分析与精度分配研究打下基础。建模的步骤如下：

(1）试验设计。根据具体的研究对象确定公差最大变动范围（Tmin,Tmax），采用正交试验法，在（Tmin,Tmax）内等间隔地取n个值Ti(i=1,2，…，n），并以Ti作为研究对象的公差，运用蒙特卡洛模拟法求解旋量参数α、β、w的实际变动区间带宽Fiα、Fiβ、Fiw(i=1,2，…，n）。试验完成后，每个旋量参数均得到n个实际变动区间带宽响应值。

(2）以公差T为设计变量，旋量参数实际变动区间带宽Fy(y=α，β，w）为待求性能函数，（Ti,Fiy)(i=1,2，…，n;y=α，β，w）为建模样本，运用不含交叉项的二次型多项式构造Fy与T之间的响应面函数

式中，ci(i=0,1,2）为待求系数。

采用最小二乘法求解C=[c0c1 c2]T的无偏估计：

(3）验证模型精度。响应面生成后，采用复相关系数R2对响应面模型进行验证，其表达式为

式中，^Fiy、Fiy分别为第i个试验点对应的响应面预测值与试验计算值；为试验计算值的均值。

复相关系数R2代表了响应面预测值与真值之间的差异程度，在0～1之间取值，值越大，表示差异度越小。

2.2 平面度公差建模

形位公差数学建模的目的是确定形位公差域边界的变动范围及公差变动要素的表示方法。由工程实践可知，对某几何要素给出形位公差的原因是：当其他的公差（如尺寸公差）存在，但仍不能满足对几何要素形位变动限制的要求时，需要给出公差值较小的形位公差来进一步限制。实际上，尺寸公差与形位公差是同一实际几何要素反映的两种不同概念的公差，是同一问题的两个方面。因此，在研究形位公差的建模时不能离开尺寸公差而单独处理，需要进一步研究形位公差与尺寸公差的关系，即要考虑公差原则。

由GB4249-1996可知，公差原则包含独立原则和相关要求，相关要求中又分为包容要求、最大实体要求和最小实体要求。本文以独立原则为例，对平面度公差进行建模。

2.2.1 求解平面度公差变动不等式与约束不等式

由平面度的定义可知，若平面度大小为TF，则平面上所有的点均必须位于距离为TF的两平行平面所形成的公差域内。用Bottom和Top分别表示两平行平面，Bottom面表示材料边内的平面，Top面表示材料外的平面，平面度公差规定的形式如图3所示。

图3中，+TSU、-TSL分别表示尺寸公差TS的上下偏差；TF表示平面度公差；平面A为基准平面，平面B为变动平面。设Bottom与Top平面对应的SDT分别为（αB，βB,0,0,0,wB)T、（αT，βT,0,0,0,wT)T，则Bottom面和Top面的变动方程分别为

由独立原则可知，平面B必须在尺寸公差域（ST域）内，即与平面A相距D-TSL～D+TSU的区域。此时，对平面度公差域（FT域）的唯一限制是在FT域内存在一个平面，它同时也在ST域。采用基于数学定义的公差分析方法求解得到Bottom面旋量参数的变动不等式为

约束不等式为

其中x、y分别取Bottom面4个顶点的x、y向坐标值。

由于Bottom面和Top面互相平行，故有

通过式（7）～式（10）确定了Bottom和Top面的变动范围，但变动平面在Bottom和Top面限制的区域中如何分布并不确定。为此，本文采用二元线性回归法求解变动平面旋量参数的值，其基本思路如下：已知Bottom和Top面的方程，在名义平面内设置m×n个取样点，使之成m×n网格分布，通过随机模拟法对每个取样点在Bottom和Top面限制的范围内取z值，这样就得到Q(Q=m×n）个平面度区域内的模拟点。根据Q个模拟点的值，采用二元线性回归算法计算得到变动平面旋量参数的值。旋量参数的表达式如下：

2.2.2 蒙特卡洛模拟法求解旋量参数实际变动区间

由于可以假设变动后公差实际表面的形状保持不变，故用蒙特卡洛模拟法求解尺寸公差旋量参数实际变动区间时，通过模拟抽样得到的旋量参数值即可作为公差的实际旋量参数。与尺寸公差不同之处在于，形状公差只确定了公差特征要素的分布区域，通过蒙特卡洛模拟法抽样只能确定公差特征要素分布区域的上下边界（即Bottom和Top面），而特征要素在区域内如何分布仍不确定，需要再采用二元线性回归法求解公差的实际旋量参数。求解步骤如下：

(1）根据研究对象的具体情况，确定旋量参数的理想概率分布模型。本文假定旋量参数的分布函数均为正态分布：

(2）根据Bottom面变动不等式式（7）确定各旋量参数理想分布的均值与方差。由式（7）可知，旋量参数αB、βB、wB理想分布的均值与方差分别为

(3）根据抽样规则对旋量参数进行抽样。Bottom面的3个非零旋量参数同样有6种变动顺序，需对每一种变动顺序都进行抽样。以变动顺序αB→βB→wB为例，其抽样流程如图4所示。

按6种变动顺序进行抽样后，α、β、w均得到容量为6 K的实际变动区间带宽样本。

第（4)～(6）步与2.1.2节所述过程相似，在此不再重复阐述。

2.2.3 运用响应面法建立旋量参数实际变动区间带宽与公差间的响应面函数

由式（7）可知，与尺寸公差响应面方法建模的不同之处在于，平面度公差建模过程中包含尺寸公差和平面度公差2个设计变量，需同时考虑这两者对变动区间带宽的影响，建模的步骤如下：

(1）试验设计。根据具体的研究对象确定尺寸公差TS和平面度公差TF的最大变动范围（TSmin,TSmax）、（TFmin、TFmax）。采用正交试验法，在（TSmin,TSmax）、（TFmin、TFmax）内分别等间隔取nS、nF个值，得到n(n=nS×nF）组（Ti,S,Ti,F)(i=1,2，…，n）值，以（Ti,S,Ti,F）作为研究对象的尺寸公差和形位公差，运用蒙特卡洛模拟法求解旋量参数α、β、w的实际变动区间带宽Fiα、Fiβ、Fiw(i=1,2，…，n）。试验完成后，每个旋量参数均得到n个实际变动区间带宽响应值。

(2）以公差TS、TF为设计变量，旋量参数实际变动区间带宽Fay(y=α，β，w）为待求性能函数，（Ti,S,Ti,F,Fiy)(i=1,2，…，n;y=α，β，w）为建模样本，运用不含交叉项的二次型多项式构造Fy与TS、TF之间的响应面函数：

采用最小二乘法求解C=[c0c1 c2 c3 c4]T的无偏估计

(3）验证模型精度。与2.1.3节所述过程相似，在此不再重复阐述。

3 实例分析

以平面度公差为例，如图3所示，设a=40mm、b=30mm,3个旋量参数α、β、w的理想概率均服从正态分布，TS、TF的最大变动范围分别为0.4～1mm、0.05～0.4mm。在最大变动范围内均等间隔取10个值，得到100个试验点Ti,S,Ti,F(i=1,2，…，100）。运用MATLAB软件，按2.2.2节所述的蒙特卡洛模拟法求解旋量参数实际变动区间带宽步骤编写计算程序，计算得到每个试验点对应的旋量参数实际变动区间带宽值Fiy(i=1,2，…，100;y=α，β，w）。将实际变动区间带宽与由变动不等式式（7）确定的理想变动区间带宽进行对比分析，如图5～图7所示。图中，残余带宽为理想变动区间带宽与实际变动区间带宽的差值。

令下降比例等于残余带宽除以理想变动区间带宽。下降比例均值为100个试验点对应的下降比例的均值。下降比例均值反映了将蒙特卡洛模拟法应用于精度设计时对提高精度设计技术经济性的效果。表2列出了3个旋量参数实际变动区间带宽下降比例均值。

由表2可知，在考虑约束条件及参数变动顺序后，运用蒙特卡洛模拟法求解得到的旋量参数实际变动区间带宽较理想变动区间带宽下降比例的均值分别为22.49%、22.85%、20.34%，说明采用蒙特卡洛分析方法可以显著提高精度设计的技术经济性，且更加符合零件生产实际。

以100组（Ti,S,Ti,F,Fyi)(i=1,2，…，100;y=α，β，w）值为建模样本，分别建立旋量参数实际变动区间带宽Fy(y=α，β，w）与公差TS、TF间不含交叉项的二次型多项式响面模型，模型表达函数如下：

用复相关系数R2验证3个响应面模型的精度，如表3所示。由表3可知，3个响应面模型的R2值均大于0.985，表明本文建立的公差响应面模型具有很高的精度。

4 结语

本文将蒙特卡洛模拟法良好的求解非确定性问题的能力与响应面方法优秀的建模性能结合起来，提出了一种基于蒙特卡洛模拟与响应面方法的公差建模方法。在对平面尺寸和形状公差域变动关系进行深入研究的基础上，采用蒙特卡洛模拟法求解公差与公差变动要素变动区间之间的响应关系，再运用响应面方法建立了两者之间的响应面数学模型。以平面度公差建模为例，验证了该方法可以显著提高精度设计的技术经济性。进一步将该方法应用到机床的精度设计当中，对提高机床精度设计方法的实用性，降低机床的生产成本具有重要意义。

蒙特卡洛模拟方法篇2

期权是最重要的金融衍生工具之一。合理定价是期权发挥功能的基础。对于欧式期权, 我们已经有经典的Black-Scholes公式[1], 由于美式期权允许期权持有人在期权到期日之前的任何时刻执行, 无法使用此公式为其定价, 所以美式期权的定价一直是一个极具挑战性的话题。美式期权一般不存在封闭解, 即使在某些特殊情况下存在封闭解, 也会因含有大量的奇异积分而在实际中很难直接应用, 因而寻求有效的数值方法就成共识。已有许多研究者对此进行了研究, 提出了许多有用的方法, 如:蒙特卡洛法[2,3]、二叉树法[4]、有限元法[5]、有限差分法等。本文在此基础上, 把蒙特卡洛模拟法与二叉树法、有限差分法二种数值方法进行比较, 得出它的优缺点;介绍它的基本原理, 并以单个标的资产美式看跌期权为例提供一个具体的算法实现步骤。

二、几种数值方法的分析比较

二叉树法和有限差分法采用逆向求解的方法均不适合处理具有多个标的资产的期权定价问题。这是因为当期权的标的资产不只一种时, 这两种方法会因为离散节点数量呈指数的急剧增长而变得不可行, 而且这两种方法也无法处理期权的收益依赖于标的资产价格历史信息的期权定价问题。与二叉树法和有限差分法相比, 蒙特卡洛模拟法具有两个明显的优势:一是比较灵活易于实现和改进[3];二是估计的标准误差及收敛速度与解决问题维数独立, 从而能够很好地解决基于多标的资产的高维衍生证券的定价问题[8]。但是它也存在明显的不足:一是由于它采用的是正向求解的方法, 无法计算出在每个时刻继续持有期权的期权收益, 从而无法比较在该时刻立即执行期权的收益与继续持有期权的期权收益, 进而无法决定是立即执行期权还是继续持有期权, 所以直到几年前, 人们还认为此法只适合具有固定执行时间的欧式期权定价, 直到近年来, 随着数理金融学的发展, 出现了一些运用蒙特卡洛方法模拟美式期权定价的算法, 其中影响最大的是由Longstaff和Schwartz提出的最小二乘蒙特卡洛法, 该方法已成为目前使用蒙特卡洛模拟美式期权定价的标准方法;二是模拟结果不太精确, 具有较大的波动性, 在路径数目较小时表现尤甚。为提高精度, 要么增加模拟路径, 要么减少波动方差, 如文献[3]对如何减少模拟方差进行了研究。

三、蒙特卡洛模拟法的算法实现步骤

蒙特卡洛模拟方法的基本原理是:在有限个离散的时间点上, 根据标的资产价格的模拟样本路径在每个时刻的截面数据, 利用最小二乘回归法求得继续持有期权的期望收益, 并将其与该时刻立即执行期权的收益相比较, 如果后者大于前者, 则立即执行期权, 否则, 继续持有期权。

下面以单一标的资产美式看跌期权定价为例, 说明蒙特卡洛法的算法实现步骤, 使我们更好地理解和使用本方法。

设股票价格为S= (St) 0≤t≤T, 以该股票为标的资产的无红利支付美式看跌期权的执行价格为K, 到期日为T, r为无风险利率, 设标的资产价格服从几何布朗运动其中μ为标的资产价格的期望收益率, σ为波动率, dω为标准的布朗运动。假设期权市场是风险中性的。

(一) 生产标的资产价格样本路径

将[0, T]分成N等分, 则每个小区间的长度为在风险中性市场下, 用r代替μ, 则根据Ito公式可将改写为dlnS=r-σ22 dt+σdω, 从而可得为服从均值为0, 标准差为1的标准正态分布的随机样本。

从而可得这样由此式可得标的资产价格的一条样本路径Lj (S0, S1, …, ST) , j∈{1, …, M}, M为模拟样本路径的数量, 经过M次模拟, 可以得到样本路径矩阵。

(二) 计算每条样本路径的最优执行时间和期权收益

对美式看跌期权, 在任意时刻i的执行价值为Xi=m ax (K-Si) 。由于美式期权可以提前执行, 这使得我们在决定最优执行时间时, 必须权衡该时刻立即执行期权的即时收益与继续持有该期权的期望收益, 则i时刻的期权价值为只有通过逆向求解的方法得到, 因为它依赖于下一步执行期权的决策。换句话说, 为了知道现在的最优策略, 我们必须首先知道未来的最优策略, 而该策略由间接地定义。这正是用蒙特卡洛模拟方法模拟美式期权定价的难点所在, 也是为什么过去认为蒙特卡洛方法不适合模拟美式期权定价的主要原因, 因此, 直接应用模拟是不可行的的, 为此我们需要近似上述期望。蒙特卡洛方法就是以多项式函数近似中的条件期望, 其中多项式各系数通过最小二乘回归方法得到[2,3], 然后, 从后往前倒推, 逐步得到各时刻的γi, 直到0时刻, 就得到0时刻的期权价值。

(三) 对每条样本路径的期权收益贴现并求均值

经过M次模拟后, 得到M条样本路径, 以及每条样本路径上最优执行时间的期权收益。由于每条样本路径的执行时间不同, 对期权收益的贴现因子也不同, 所以, 必须按相应的贴现因子贴现, 然后求均值即得模拟值。

四、蒙特卡洛模拟算法的收敛性和计算效率

蒙特卡洛算法的有效性包括收敛性和计算效率两个方面。它们主要取决于基函数数目, 离散时间点数目N, 样本路径数目M, 以及模拟随机数的产生方法。关于收敛性问题, Longs taff和Schw artz证明, 令样本路径数目M→∞, 只要基函数数目足够大, 则模拟期权价值收敛于真实期权价值的ε领域内 (ε为任意正数) 。蒙特卡洛算法的一个重要特点是它优异的计算效率, 尤其在多个标的资产美式期权定价模拟方面, 与处理多维美式期权定价模拟的随机网格方法相比, 该算法的计算效率要高出600多倍。

摘要：针对美式期权的定价问题, 把蒙特卡洛模拟法与二叉树法、有限差分法二种数值方法进行比较, 得出它的优缺点;介绍蒙特卡洛模拟法的基本思想, 并以单个标的资产美式看跌期权为例说明其具体算法步骤。

关键词：美式期权,蒙特卡洛法,期权定价

参考文献

[1]F Black, , M Scholes.The pricing of optionsand corporate liabilities[J].J Political Economy, 1973, 81:637-654.

[2]Longstaff, Schwartz E S.Valuing American optionsby simulation:a simple least squares approach[J].Review of Financial Studies, 2001, 14 (1) :113-148.

[3]郑承利美式期权的几种蒙特卡洛仿真定价方法比较[J]系统仿真学报2006, 18 (10) :2929-2935

[4]约翰赫尔著张陶伟译期权、期货及其衍生证券[M]华夏出版社北京2000

[5]Hull J Option futures and other derivative securities[M]2nd Edition Englewood Cliffs:Prentice Hall1993

[6]张铁美式期权定价问题的数值方法[J]应用数学学报200225 (1) :113-122

[7]李玉立金朝嵩美式看跌期权的差分格式[J]重庆建筑大学学报200426 (4) :100-114

蒙特卡洛模拟方法篇3

近年来，电力市场化改革[1]和以风电为代表的可再生能源接入[2]给电力系统运行和规划带来更多不确定性，随机潮流PLF(Probabilistic Load Flow）计及这些不确定因素的影响，可综合分析系统的运行风险和薄弱节点，因而受到广泛关注[3,4,5]。

PLF首先由Borkowaska于1974年提出[6],40年来，国内外学者对其进行了大量优化与改进，目前主要可分为解析法[7,8,9,10]和模拟法[11,12,13]两大类。解析法可分为卷积法[7]、半不变量[8,9]、点估计[10]。卷积法利用输入变量分布函数的卷积运算得到输出变量的概率分布，计算量大。半不变量通过线性化潮流方程建立的线性关系将卷积运算转化为半不变量的代数运算，计算效率高但准确性差，得到的变量概率分布在首尾可能出现小于零的部分。点估计是已知输入变量的若干阶矩信息求取输出变量矩特征的数学方法，计算效率高，但高阶矩误差较大。

除卷积法外，解析法计算效率高，但误差均较大。而模拟法可准确地对实际物理过程建模，且建模质量不受变量分布性质、系统非线性的影响，因而比其他PLF方法更准确，其中蒙特卡洛模拟MCS(Mont Carlo Simulation）应用最为广泛。基于简单随机抽样SRS(Simple Random Sampling）的MCS收敛慢，获得准确解的计算代价大，因此，如何在保证精度的同时提高计算效率一直是MCS的研究热点。研究表明，样本的差异性[14]是影响MCS计算精度的主要因素。文献[11,12,13]采用各种基于拉丁超立方采样LHS(Latin Hypercube Sampling）的方法改进MCS，但由于无法保证序列的低差异性，因此无法克服MCS收敛性的瓶颈。而以Sobol序列为代表的低差异序列LDSs(Low-Discrepancy Sequences)MCS在计算精度和效率上均优于LHS，已在电子电路设计中取得了应用[14]，但在随机潮流中的应用还没有见诸报道。

另一方面，为得到含相关性的样本序列，随机排序[15]、Cholesky分解[16]、遗传算法[17]等改进排序方法均在MCS中得到应用。随机排序在样本较小时效果差；Cholesky分解仅针对相关系数矩阵正定的情况遗传算法有迭代操作，结果较准确但计算耗时长。因此，本文提出采用奇异值分解的方法对样本进行排序，计算量与Cholesky分解相当，并可处理相关系数矩阵非正定的情况。

为提高MCS的计算效率并处理相关系数矩阵非正定的情况，本文提出一种基于Nataf变换和准蒙特卡洛模拟NQM(Nataf transformation based Quas Monte Carlo simulation method）的PLF方法。首先利用Sobol方法产生服从均匀分布的LDSs，而后利用基于奇异值分解的Nataf变换生成服从输入变量联合分布的序列，再经过若干次确定型潮流计算得到输出变量的概率分布特性。对IEEE 30节点系统和我国某实际大区域电网的仿真验证了本文方法的有效性。

1 基于奇异值分解的Nataf变换

Nataf变换是在已知输入变量边缘分布的情况下重构联合分布的数学方法[18]。现将其过程简述如下：对任意n维输入变量X=[X1，…，Xi，…，Xn]T设其相关系数矩阵为[ρXij]n×n，其中ρXij定义如式（1所示。

其中，σi和σj分别为Xi、Xj的标准差；cov(Xi,Xj）为变量Xi、Xj的协方差；ρXij为相关系数。随机向量X可用Nataf变换转化为标准正态分布向量Z=[Z1，…，Zi，…，Zn]T:

其中，Fi为Xi的累积分布函数；Φ为标准正态分布的累积分布函数。设变换以后Z的相关系数矩阵为[ρZ ij]n×n，由Nataf变换性质，[ρXij]n×n与[ρZij]n×n可根据式（3）相互转化：

其中，fXi Xj(Xi,Xj）为Xi、Xj的联合概率密度函数；μ、σ分别为相应变量的期望值与标准差；Zi Zj(Zi,Zj，ρZ ij）为相关系数为ρZij的二维标准正态分布概率密度函数。本文采用文献[19]中方法对式（3）进行求解以求得[ρZij]n×n。而由于[ρZ ij]n×n可能为非正定或非满秩阵，此时其Cholesky分解不存在，无法采用文献[16]的方法对样本进行排序。但一般相关系数矩阵都为对称阵，存在奇异值分解，本文通过证明定理1说明奇异值分解排序方法的过程。

定理1设K为n维独立标准正态分布向量，矩阵由向量Z的相关系数矩阵ρZ的奇异值分解生成：

其中，UρZ为酉矩阵；∑ρZ为由矩阵奇异值构成的对角矩阵，并按从大到小排列。则由式（5）定义的Z*相关系数矩阵等于ρZ:

证明如下。

对其期望值满足式（6):

蒙特卡洛模拟方法篇4

传统的边坡稳态评价方法是以安全系数作为评价指标, 计算结果的正确性取决于计算模型与实际的吻合程度, 缺点是将安全指标常数化, 忽略了岩土体物理力学参数的不确定性和离散性。相对而言, 基于概率论和数理统计的边坡可靠性分析方法更加严谨, 更加科学。这里运用蒙特卡洛模拟方法对雅砻江某边坡进行稳定性分析。

1 边坡可靠性蒙特卡洛法模拟

1.1 Monte Carlo模拟基本原理

蒙特卡洛 (Monte Carlo) 模拟方法是一种依据统计抽样理论, 利用计算机研究随机变量的数值方法。在目前的可靠度计算中, Monte Carlo模拟方法广泛应用于土木工程、机械工程、水电工程以及边坡工程等领域。对于一般的边坡稳定性问题, 根据岩土结构、破坏机理和受力状况, 可以建立如下的状态函数:

K=F (x1, …, xm) (1)

其中, x1, …, xm分别为重度、粘聚力、摩擦系数、孔隙水压力、荷载强度、降雨强度等随机变量。

它们具有一定的分布 (大多服从正态分布或对数正态分布) , 其统计值为已知。设状态函数为安全系数, 且随机地从各个随机变量xi (i=1, 2, …, m) 的母体中抽取一个具有相同分布变量x′1, …, x′m, 由式 (1) 求得一个安全系数的随机样本K′。如此重复, 直至达到预期精度的充分次数N, 就可得到N个相对独立安全系数样本值K1, …, KN。安全系数所表征的极限状态为K=1, 构造随机变量为:

设在N次试验中, 出现Ti=1 (1, 2, …, N) , 即K≤1的次数为M, 则斜坡的破坏概率为:

$Ρ_{f} = \frac{Μ}{Ν}$ (4)

当N足够大时, 由安全系数的统计样本值K1, …, KN可能较精确地拟合安全系数的概率分布F (k) , 并估计其分布参数。其均值和标准差分别为:

$μ_{Κ} = \frac{1}{Ν} \sum_{i = 1}^{Ν} Κ_{i}$ (5)

$σ_{Κ} = [\frac{1}{Ν - 1} \sum_{i = 1}^{Ν} (Κ_{i} - μ_{Κ})^{2}]^{\frac{1}{2}}$ (6)

由中心极限定理, 在N≥50的条件下, 破坏概率的积分表达式为:

$Ρ_{f} = \frac{1}{\sqrt{2 π} σ_{Κ}} \int_{0}^{1} e^{- \frac{(Κ - μ_{Κ})}{2 σ_{Κ}^{2}}} d Κ$ (7)

由式 (4) 和式 (7) 计算破坏概率可得到可靠性指标为:

β=Φ-1 (1-Pf) (8)

其中, Φ (β) 为标准正态分布函数。

$Φ (β) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{β} e^{- \frac{x^{2}}{2}} d x$ (9)

1.2 Monte Carlo法模拟的基本流程

Monte Carlo法模拟的基本流程如图1所示。

1.3 误差估计与模拟次数

由于随机抽样是概率为P的伯努里试验, 所以Pf的期望值为: $E (Ρ_{f}) = E (\frac{Μ}{Ν}) = p$ 。

Pf的方差为: $D (Ρ_{f}) = \frac{p (1 - p)}{Ν}$ ;

Pf的标准差为: $σ_{Ρ}_{f} = \sqrt{\frac{p (1 - p)}{Ν}} (10)$

在实际应用中p为未知, 可用计算Pf作为p的估计值, 则有:

$\overset{\land}{σ}_{Ρ_{f}} = \sqrt{\frac{p (1 - p)}{Ν}}$ 。

从式 (10) 可见, 随机抽样计算结果的标准差σPf与 $\sqrt{Ν}$ 成反比, 若需将其标准差降低一个数量级, 则试验次数必须增大百倍。

现确定随机抽样需要的试验次数。按中心极限定理, 当N充分大时, 有:

$\frac{Ρ_{f} - p}{σ_{Ρ_{f}}} \sim Ν (0, 1)$ 。

因此, 显著水平为α的置信区间为:

$Ρ {| \frac{Ρ_{f} - p}{σ_{Ρ_{f}}} | \leq t_{\frac{α}{2}}} = 1 - α$ 。

其中, $t_{\frac{α}{2}}$ 为N (0, 1) 的 $\frac{α}{2}$ 的百分位值, $t_{\frac{α}{2}}$ 满足下式:

$\frac{1}{\sqrt{2 π}} \int_{- \infty}^{t_{\frac{α}{2}}} e^{- \frac{t^{2}}{2}} d t = 1 - \frac{α}{2}$ 。

于是得到显著水平为α的置信区间为:

$| \frac{Ρ_{f} - p}{σ_{Ρ_{f}}} | \leq t_{\frac{α}{2}}$ 。

设p和Pf的绝对误差为ε, 则可得到:

$ε = | Ρ_{f} - p | \leq t_{\frac{α}{2}} \cdot σ_{Ρ_{f}} = t_{\frac{α}{2}} \sqrt{\frac{p (1 - p)}{Ν}} \approx t_{\frac{α}{2}} \sqrt{\frac{Ρ_{f} (1 - Ρ_{f})}{Ν}}$ 。

可得模拟次数为:

$Ν = \frac{t_{\frac{α}{2}}^{2} \cdot Ρ_{f} \cdot (1 - Ρ_{f})}{ε^{2}}$ (11)

在式 (11) 中, Pf较小, 模拟次数N较大, 随着N值的增大, 误差ε减小, 逐渐趋于收敛。对于给定误差ε和置信度1-α, 假定N取某一值, 由式 (11) 可确定相应的ε值, 如果ε≤ε0, 所取的N满足要求;如果ε>ε0, 可继续增加N, 直至满足ε<ε0, 可求出所给定精度的随机模拟解。

一些学者建议蒙特卡洛模拟次数N需要满足下式:

N≥100/Pf (12)

由于Pf是一个很小的数值, 这就要求N很大, 如Pf=0.1%, 则计算次数N将达到十万次, 这样大的模拟次数, 在进行计算机分析时, 将花费过多的时间。对于一般性的边坡工程, 取5 000～10 000, 已经完全可以满足工程评价的精度要求。

2 工程运用

2.1 滑坡体基本情况

雅砻江某滑坡体整体呈“U”字形, 后缘分布高程约2 420 m, 后缘有基岩出露, 形成高约4 m～7 m的陡壁, 基岩岩性为大理岩、变质砂岩, 前缘至江边, 江边无基岩出露, 上、下游侧均为冲沟。滑坡体沿江边长约485 m, 从江边到滑坡体后缘平面距离约660 m, 估算面积约28.8万m2, 高程在2 200 m以上, 平均估算厚度约21.67 m。

据钻孔平洞揭露, 物质组成主要由崩坡积物及中部滑坡堆积物、底部滑带等组成, 下伏基岩岩层产状一般为N15°～30°W NE∠45°～55°, 顺坡倾向。

2.2 蒙特卡洛模拟分析

本例并非采用GEOSLOPE软件原有蒙特卡洛模拟, 而是对原有岩土参数蒙特卡洛抽样10 000次, 并通过自编程序将抽样数据导入SLOPE软件逐次计算, 计算结果显示如图2, 图3所示。

通过Bishop方法计算得安全系数均值为1.249;通过Morgenstern-Price 法计算得安全系数均值为1.216。另外输入原有数据, 采用软件自带蒙特卡洛算法可得计算结果如图4, 图5所示:Bishop方法计算得安全系数均值为1.411;通过Morgenstern-Price法计算得安全系数均值为1.408。两种算法计算结果比较见表1。

3 结语

本文结合GEOSLOPE软件和抽样软件提出了一种新的边坡稳定性计算方法。利用抽样软件对力学参数进行蒙特卡洛抽样, 再通过自编软件将抽样数据导入GEOSLOPE软件进行计算和统计, 得到安全系数均值, 该值较原有软件计算结果小10%～15%, 计算结果更加精确, 可信度更高。

参考文献

[1]祝玉学.边坡可靠性分析[M].北京:冶金工业出版社, 1993.

[2]孟庆山, 雷学文.土工参数的统计方法及其工程应用[J].冶金科技大学学报, 1992 (1) :10.

[3]张天宝.土坡稳定分析与土工建筑物的边坡设计[M].成都:成都科技大学出版社, 1987.

[4]潘家铮.建筑物的抗滑稳定和滑坡分析[M].北京:中国水利电力出版社, 1980.

[5]Sarma, S.K..Stability Analysis of Embankments and Slopes[J].Geotechnique, 2001 (23) :51-52.

[6]崔政权, 李宁.边坡工程——理论与实践最新发展[M].北京:中国水利电力出版社, 1999.

蒙特卡洛模拟方法篇5

1 井间干扰概率评价模型及流程

为了衡量干扰程度大小,提出了干扰井数比例和干扰储量比例两个指标。其中干扰井数比例指评价区中发生井间干扰的气井个数占全部气井的比重,干扰储量比例指因井间干扰而重复开发的储量占评价区总储量的比重,是造成投入和产出不匹配的主要原因之一。

干扰井数比例=发生井间干扰的气井数/总气井数;

干扰储量比例=气井共同开采储量/评价区总储量。

1.1 两口气井井间干扰程度

为评估干扰程度大小,结合矿场经验做出如下假设和简化:(1)单口气井泄气面积随机性强,但大量气井统计的单井泄气面积符合一定概率分布;(2)单井动态储量与泄气面积之间线性正相关;(3)为提高储层纵向动用程度和单井产量,气井常采取分层压裂合层开采工艺,因而泄气面积简化为单井多层的综合平均泄气面积,不再分层测算;(4)均匀井网。

井间干扰评价的关键问题是确定泄气面积干扰程度以及泄气面积内气井动态储量干扰程度。设两口气井泄气半径为r和R,井距为L,则井间干扰存在以下两种可能性(图1)。第一种可能性,井网较稀疏,井距大于气井泄气半径之和(L>r+R),井间无干扰;第二种可能性,井网较密(R-r<L<R+r),井间存在干扰。此外,生产中一口气井泄气范围不可能覆盖临井,因此理论上不存在R>L+r的可能性。井间干扰面积可根据两口气井的泄气半径和井距确定[16]。

井控动态储量Re和泄气面积S具有一定的线性相关性:Re=ρ×S,其中两者斜率为井控动态储量丰度ρ。基于井控动态储量和泄气面积之间的统计规律可以得到干扰储量大小。井间干扰储量为:

式(1)中ρr和ρR分别为两口井的井控动态储量丰度,θr和θR分别为两个弧的角度。

1.2 井区井间干扰程度

针对特定评价区域(图2),井间干扰程度评价过程如下:首先,对所有气井Nt的井控动态储量求和的R1;其次,消除内部干扰井N干扰的井间干扰储量R2对井网加密效果的影响;再次,评价区域外储量R3在界定计算单元之外需要劈分出去,基于R3对井控储量进行修正;最后,根据评估井区总储量Rt测算储量控制程度、干扰井数比例和干扰储量比例大小。其中,储量控制程度

干扰井数比例

干扰储量比例

1.3 井间干扰概率模拟流程

气井泄气面积和井控动态储量不确定性大、随机性强,为了模拟加密后井间干扰程度的所有可能结果,引入了蒙特卡洛随机模拟方法。蒙特卡洛基本原理如下:通过对气井泄气面积和井控动态储量按其概率分布进行Num次随机抽样,并且对每一次抽样结果都按上述关系式进行一次井间干扰的运算处理,当模拟次数趋于无穷大时,这Num次运算的结果就构成了井间干扰的概率分布[16]。

井间干扰概率模拟流程如下:

(1)选择代表性区块,要求该区块各项参数随机分布规律与全气田相似,而且生产历史较长,气井泄气面积和井控动态储量基本稳定。

(2)结合地质认识和生产动态分析方法,分析确定代表性区块内各气井的泄气面积和井控动态储量大小。

(3)绘制气井泄气面积和井控动态储量概率分布图,绘制单井泄气面积与井控动态储量关系图。

(4)确定某一井网密度大小,根据单井泄气面积与井控动态储量概率分布规律,随机多次抽值,测算储量控制程度、干扰井数比例和干扰储量比例大小。

(5)变换不同井网密度,重复步骤(4)。

(6)整理随机模拟结果,统计分析不同井网密度时储量控制程度、干扰井数比例和干扰储量比例的随机概率分布和期望值,优化井网密度。

2 井间干扰概率规律分析—以苏×井区为例

苏里格气田苏×井区,自2002年先后开展了老井试采、加密钻井、干扰试验等多项试验研究工作。对全区82口气井进行分析预测,平均泄流面积0.27 km2,中值0.21 km2,远低于骨架井网600 m×1 200 m的井控面积0.72 km2,骨架井网不能充分动用地质储量,仍有加密潜力。

苏×加密区位于苏×井区中部,面积12 km2,天然气地质储量丰度1.348×108m3/km2,加密区气井泄气面积和井控动态储量分布规律与苏×井区基本一致(图4),符合统计规律,因此可以利用加密区有限单井样本预测整个区块的加密潜力。为此,模拟评价了井网密度从1口/km2不断加密到8口/km2时井间干扰程度变化规律,其中每一井网密度随机模拟500次。模拟结果如下:

2.1 储量控制程度概率

初始井网较稀疏时,加密能显著提高储量控制程度。随着井网逐渐加密,井间干扰增强,井控储量增量趋于平缓,井网密度达到4井/km2后趋于平稳(图5)。

2.2 干扰井数比例概率

干扰井数随井网加密而快速增加。1井/km2时,井间几乎无干扰,加密至2井/km2时干扰井数比例P50=24%左右,加密至4井/km2时干扰井数P50=92%,加密至8井/km2时干扰井数比例约为100%(图6)。

2.3 干扰储量比例概率

干扰储量随井网加密而逐渐增加,但相较干扰井数“慢热”。1井/km2时,几乎无干扰储量;加密至2井/km2时,干扰储量仍然较少,P50=0.6%;加密至4井/km2时干扰储量P50=17%,加密至8井/km2时干扰储量P50=94%(图7)。与干扰井数比例相比,井网加密后干扰储量比例增长速度更慢,同一井网密度干扰储量比例小于干扰井数比例,例如井网密度4井/km2时,干扰井数比例中值92%,干扰储量比例17%,此时虽然大部分气井发生井间干扰,但是彼此之间相互干扰储量规模有限。由于干扰储量更加直观反映井网控制储量和最终采出量,是决定井网经济收益的本质因素,因此合理井网优化应更加重视对干扰储量比例这一指标的考量。

2.4 动态分析与综合评价

苏×加密区评价面积12 km2,骨架井网11口井,加密20口气井后井网密度达到2.6口/km2。干扰概率模拟结果表明,井网加密至2.6井/km2时,干扰井数比例P50=64.2%,干扰储量比例P50=4.5%,井控储量P50=56%,区块采收率46%。

加密后,有7口加密井在投产前地层压力低于区块原始压力,表明至少14口井(45%)存在井间干扰。加密区气井动态分析表明,加密使得储量控制程度达到51%,采收率达到44%,动态评价结果与井间干扰概率模拟结果基本一致,表明井间干扰模拟方法具有一定的实用性。

3 结论

(1)低渗致密砂岩气藏开发具有较大的不确定性,井网加密风险大,基于地质统计原理和蒙特卡洛方法建立了储量控制程度概率、干扰井数比例概率和干扰储量比例概率模拟方法,并形成了井间干扰概率评价流程。

(2)井网加密后干扰井数比例和干扰储量比例均增加,但同一井网密度干扰储量比例小于干扰井数比例,合理井网取值应更加侧重于对干扰储量比例的考量。

蒙特卡洛模拟方法篇6

Black-Scholes模型的假设条件过于完美, Papapantoleon对Levy过程进行了概述[1], 指出一类无限可分、左极限右连续的随机过程, 包括了Brown运动、泊松过程、指数过程及其它无穷跳跃过程, 一般没有确切的密度函数解析式, 而特征函数形式可根据Laplace分解得到具体形式, 即φ (u) =E (eiuX) =e-tφX (u) , 特征指数φX (u) 可以表达为Levy-Khinchine公式:

$φ_{X} (u) = - i μ u + \frac{σ^{2}}{2} u^{2} + \int_{R} (1 - e^{i u x} + i u x 1_{| x | < 1}) k (d x) (1)$

(μ, σ, k (dx) ) 称为Levy过程特征三项, 分别表示线性漂移率、连续布朗运动扩散、纯跳跃测度, 这三项测度能全面概括随机过程的运动特征。有限的跳跃到达率为: $Ι = \int_{R} k (d x) = λ < \infty$ , 无穷到达率为:I=∞.

随机过程的一次、二次变差表达式为:

$\begin{array}{l} J_{1} = \int_{- \infty}^{+ \infty} | x | k (x) d x \\ J_{2} = \int_{- \infty}^{+ \infty} x^{2} k (x) d x (2) \end{array}$

Madan, Carr和Chang[2]建立的VG (C, G, M) 模型可以表示为正Gamma (C, M) 和负Gamma (C, G) 两个不同过程的代数和, 是一类无穷活动率过程, 存在二次变差。模型能够描绘无穷小跳跃的活动过程, Carr, Geman, Madan和Yor[3]对VG模型的进行了扩展, 在VG测度中加入了参数Y, 以此刻画不同活动水平下的随机过程。即:

$k_{C G Μ Y} (x) = {\begin{matrix} \frac{C e^{G x}}{| x |^{Y + 1}}, & x < 0 \\ \frac{C e^{- Μ x}}{x^{Y + 1}}, & x > 0 \end{matrix} (3)$

VG模型是Y=0的特殊情形。当Y<0时, 属于类似于复合泊松分布的有限活动率过程, 0<Y<1为无限活动率过程, 1<Y<2为无限活动率且不存在有限变差, Y>2为无限活动率不存在二次变差。CGMY模型能有较好的兼容效果, 能通过参数的设定得到不同活动率形式下的随机过程, 研究这个过程下的期权定价结果具有重要意义。测度的复杂性不仅给参数估计带来了一定的麻烦, 也使得计算机模拟变得非常不便。 Carr和Wu[4]提出了CGMY过程的模拟步骤, 便是基于Levy测度下的随机数生成规则, 通过差分N阶变换, 得到Kummer等式[5]的解, 属于合流超几何函数 (confluent hypergeometric function, CHF) 。Zhang[6]提出了两种模拟方案, 一是复合泊松过程模拟方法, 二是随机时变布朗运动 (Time-Changed Brownian Motion, TCBM) 。两者依赖CHF计算, 即Francesco[7]得出的解析式。Ji[8]对CGMY过程的模拟采用的Madan和Yor (2006) [9] 的方法。Zhang[6]和Ji[8]两种判断随机数淘汰的条件都是通过Kummer等式解的变形而来, 都包含CFH, 各有优势和不足, 前者极端值偏离, 后者尾部较重。

看涨期权而言, 通常以现货价值的上下浮动5%为界区分价外、价内和平价。由于CGMY测度复杂, 模拟中随机数的淘汰涉及CHF计算, 对参数的依赖程度高。Carr和Madan[10], Carr和Wu[11], Ji[8]快速傅立叶变换 (FFT) 数值方法便利而快捷, 但在虚值期权定价时失效。

基于上述现状, 本文试图解决虚值期权定价精度问题。CGMY模型特别是参数Y比较小的情形下, 探索提高模拟效率的方法, 引入复数合流超几何函数 (CCHF) , 改进计算效率, 克服模拟困难。

2 模型与参数

2.1 VG模型

基于两个相反Gamma过程的代数和建立的无穷纯跳跃VG模型, 以参数 (C, G, M) 表示。

其中Gamma (a, b) 分布的密度函数:

$f_{G} (x; a, b) = \frac{b^{a}}{Γ (a)} x^{a - 1} e^{- b x} (4)$

若G (a, b) 表示Gamma过程, 价格增量演变过程可以用等式描绘:

$Δ L_{Τ} = G_{Τ}^{1} (C, Μ) - G_{Τ}^{2} (C, G) (5)$

2.2 CGMY模型

作为无穷活动率下纯跳跃VG模型的拓展, 仅考虑纯跳跃来刻画资产价格演变过程, 那么它的特征三项中连续布朗运动扩散可以为0, 即 (w, 0, k (x) ) 。跳跃的活动到达率水平取决于Y的大小。Y在 (0, 1) 区间就能体现无穷小跳跃、无穷到达率水平。w为线性漂移率, 非随机趋势, 根据式 (3) Levy测度k (x) , 可以得到CGMY过程的跳跃频率和强度。

由于CGMY模型没有确切的密度函数解析式, 而特征函数表达式比较简洁, 通过Levy-Khinchine公式, 模型的特征函数为:

$φ_{C G Μ Y} (u) = e^{C Γ (- Y) [(Μ - i u)^{Y} + (G + i u)^{Y} - Μ^{Y} - G^{Y}]} (6)$

其中伽马函数:

$Γ (a) = \int_{0}^{\infty} x^{a - 1} e^{- x} d x (7)$

关于参数Θ下随机过程Xt的Laplace变换公式是:

$φ_{X_{t}} (u; Θ) = E (e^{i u X_{t}}) = e^{- t φ_{X} (u)} = l_{t}^{Θ} (φ_{X} (u)) (9)$

指数Laplace形式是自然底数e的指数形式, 往往称为特征指数、累积特征指数, 转化为资产价格对数收益率特别方便, 利率期限结果模型中用于计算零息债券的价格公式 (假设有价证券面值为一个单位) 即为E (e-tRt) =e-tR (λ) =lΘt (λ) , 是关于利率仿射指数变换的函数表达式。

另外, 令u=-i, 有φXt (-i;Θ) =E (eXt) =e-tφX (-i) , φX (-i) 是随机过程的指数线性漂移率、单位时间的预期对数收益率, 用于建立风险中性资产价格模型和鞅过程。

由CGMY模型的特征指数得到CGMY过程分布的累积函数:

$\begin{array}{l} C_{n} = \frac{1}{i^{n}} \frac{d^{n} \log (Φ (u))}{d u^{n}} |_{u = 0} \\ = C Γ (n - Y) [Μ^{Y - n} + (- 1)^{n} (G^{Y - n})] (10) \end{array}$

依据特征指数可以计算风险中性漂移率修正参数:

$\begin{array}{l} w = φ (- i) \\ = C Γ (- Y) [(Μ - 1)^{Y} + (G + 1)^{Y} - Μ^{Y} - G^{Y}] (11) \end{array}$

在非单位时间段t内, CGMY模型的参数变换关系:

$C_{t} = t C; G_{t} = G; Μ_{t} = Μ; Y_{t} = Y (12)$

为便于计算和模拟, 设立一个正的极小值ξ→0, 作为无穷小跳跃的门槛, 高于小跳门槛的以CGMY的Levy测度表示, 低于该门槛则以无穷小波动的连续布朗运动代替。-ξ<x<ξ称为小跳跃, |x|>ξ为CGMY的跳跃, 通过设定和排除趋近于正态分布的无穷小跳跃, 能更加方便地模拟随机数的分步以及进行数值计算。

2.3 特征函数的矩估计

有了矩条件表达形式, 参数的矩估计采用Y.Miyahara的方法[18]。

设z分布服从CGMY过程的增量分布, k阶原点矩表达式:

$\begin{array}{l} \hat{m}_{k} = E (z^{k}) \\ \hat{A} (u) = E (e^{i u z}) \\ \hat{A}^{(k)} (u) = d^{k} E (e^{i u z}) / d u^{k} \\ \hat{m}_{k} = \hat{A}^{(k)} (u) / i^{k} |_{u = 0} (13) \end{array}$

根据式 (7) 的矩条件: $\hat{C}_{n} = C Γ (n - Y) [Μ^{Y - n} + (- 1)^{n} (G^{Y - n})], Ν$ 取一个比较大的整数解四阶方程组, 参数结果表达式:

$\begin{array}{l} \hat{Y} = \frac{\ln (Re (\ln (\hat{A} (Ν))))}{\ln Ν} \\ \hat{A} (u) = E (e^{i u x}) \\ \hat{C} = \frac{Re (\ln (\hat{A} (Ν)))}{2 Ν^{\hat{Τ}} Γ (- \hat{Y}) \cos (\frac{σ^{2} \cdot \hat{Y}}{2})} (14) \end{array}$

代入下列方程组求解:

$\begin{array}{l} \hat{C} Γ (- \hat{Y}) (\hat{Y} - 1) \hat{Y} [G^{Y - 2} + Μ^{Y - 2}] = C_{2} - C_{1}^{2} \\ \hat{C} Γ (- \hat{Y}) \hat{Y} (\hat{Y} - 1) (\hat{Y} - 2) [G^{Y - 2} + Μ^{Y - 2}] \\ = C_{3} - 3 C_{1} C_{2} + 2 C_{1}^{3} (15) \end{array}$

方程 (15) 由于指数复杂, 普通计算很难有较好的准确性。可在Matlab中采用函数最小值法求参数。X表示参数集, f1 (X) 、f2 (X) 分别表示模型矩条件, 参数估计原理为

$X = a r g m i n (f_{1} (X) - Μ_{1})^{2} + (f_{2} (X) - Μ_{2})^{2} (16)$

通过上述低阶矩估计方法获取参数后, 进行期权的数值定价或蒙特卡罗模拟比较在CGMY模型下的实值期权与虚值期权之间的差异。

3 分布模拟原理

3.1 传统模拟计算技术

根据Madan和Yor[9]的模拟方法, 在随机泊松跳跃时间内先产生一个逆CGMY随机变量, 可称之为levy测度下的随机跳跃y, 再通过计算CHF (y) 作为判断条件进行拒绝-接受检验, 可以接受的随机数y作为随机时间增量生成漂移率为A的布朗运动, 步骤:

①计算以下参数:

A= (G-M) /2

B= (G+M) /2

$k = C \sqrt{π} / [(Γ ((Y + 1) / 2) \cdot 2^{Y / 2}]$

λ=2k/ (YεY/2)

d=kε1-Y/2/ (1-Y/2)

ε→0

②生成以λ为参数的指数分布a≤ΔT, 以此作为判断跳跃的时间, 反之退出模拟

③生成两个独立的均匀分布 :u1, u2～U (0, 1)

④生成列维逆正态分布, 即随机跳跃的幅度: $y = \frac{ε}{u_{1}^{\frac{2}{Y}}}$

⑤根据随机数生成接受-拒绝的判断条件 (以下称U法则) :

$u = e^{\frac{- A^{2} + B^{2}}{2}} \frac{Γ (\frac{Y + 1}{2})}{\sqrt{π}} U (\frac{Y}{2}, \frac{1}{2}, \frac{B^{2}}{2} y) > u_{2}$

⑥若满足上述判断, 以随机数作为跳跃时间幅度的累积:sum=sum+y

⑦把跳跃作为时变布朗运动的随机时间, 生成CGMY过程的增量: $Ζ = A \cdot s u m + ξ \sqrt{s u m}$ , 这样价格过程就可以表示成:

$S_{t} = S_{0} e^{(r - q) t + Ζ (t)} (17)$

U (a, b, x) 是一个合流超几何函数。

另外, Zhang[6]运用了类似的淘汰规则, 是基于Kummer方程的第二个解形式 (下简称H法则) :

sum=sum+yi1h (yi) >ui

$\begin{array}{l} h (y) = 2^{Y} \exp {\frac{(B^{2} - A^{2}) y}{2}} [F (\frac{Y}{2}, \frac{1}{2}, \frac{B^{2}}{2} y) \\ - \frac{\sqrt{B^{2} y} \cdot π \cdot Γ (\frac{Y + 1}{2})}{\sqrt{2} Γ (\frac{Y}{2}) Γ (\frac{3}{2})} F (\frac{Y + 1}{2}, \frac{3}{2}, \frac{B^{2}}{2} y)] (18) \\ F (a, b, x) = U (a, b, x) \end{array}$

3.2 复数合流超几何函数计算

特别地广义超几何函数表达式为:

$F (n, d; x) = \sum_{k = 0}^{\infty} \frac{\prod_{i = 1}^{j} \frac{Γ (n_{i} + k)}{Γ (n_{i})}}{\prod_{i = 1}^{m} \frac{Γ (d_{i} + k)}{Γ (d_{i})}} \cdot \frac{x^{k}}{k!} (19)$

Y越小产生的随机数y可能越大, 函数计算复杂度增强, 出现的极端值非常不稳定。考虑到计算过程受到参数影响, 对参数进行复变化可以收敛实数分布域, 本文采用复数变换合流超几何函数 (Kummer Complex函数) 代替传统CHF, 是Stepan[23]的一种广义傅立叶变换算法, 即:

$C (z; Θ) \equiv \int_{- \infty}^{\infty} e^{i z v} \cdot F (x; Θ) d v (20)$

z可以为一切复数, 通过变换, 可以简化相关计算。特征函数就是这类转换的一种特例。采用了傅立叶变换方式, 该函数可以扩展到复数领域参数计算, 反之也可以将虚数转换到实数空间, 是傅立叶逆转换, 表达式为:

$F (x; Θ) \equiv \frac{1}{2 π} \int_{i z - \infty}^{i z + \infty} e^{- i z k} \cdot C (z; Θ) d z (21)$

虽然需要通过正反两次变换, 但无论参数多大的实数督可以变成二元复数, 大大减少了CHF的计算时间和复杂度。

Stepan[23]指出通过函数计算对比F函数和U函数的值误差不超过10-6, 计算速度却快了很多。同样也解决了y非常大的时候广义超几何函数无结果的问题, 计算效率得到较大程度的提高。

为进行比较, 针对不同Y模拟生成1000个随机数, 并进行CGMY分布的模拟, 比较两种计算方法的平均时间。当Y越小时, 需要的时间越长, 随机数出现超极端值的概檬就越高, CHF函数的计算就越复杂, 可能导致死循环, 通过无效值比例可以看出模拟1000个随机数时, 就可能导致一次死循环。两种函数计算时间比较参见表1[10]引入转换因子, 将期权价格进行广义傅立叶变换, 期权价格的傅立叶形式为 (α为转换因子也称修正因子) :

$F (v) = φ (v) = \int_{- \infty}^{\infty} e^{i v k} e^{α k} C_{Τ} (k) d k (24)$

式中, k为执行价 (对数形式) , CT (k) 表示到期日的期权价格, 关于执行价序列的期权价格向量转换成关于v向量的傅立叶序列, 该序列有许多优点, 代入欧式期权定价公式并进行积分秩序变换, 密度函数积分成特征函数得到以下结果 (小写的价格和执行价均表示对数形式) :

F (v) =∫∞-∞eivkeαkCT (k) dk

=∫∞-∞eivkeαk∫∞ke-rt (es-ek) f (s) dsdk

通过推导得到了特征函数形式下的期权傅立叶形式的数值解:

$φ (v) = \frac{e^{- r t} Φ (v - (α + 1) i)}{α^{2} + α - v^{2} + i (2 α + 1) v} (25)$

只需要进行傅立叶逆变换就可以得到期权的实数解:

$C_{Τ} (k) = \frac{e^{- α k}}{2 π} \int_{- \infty}^{\infty} e^{- i v k} φ (v) d v (26)$

解析式代表了在特征函数下的价格封闭解形式, 连续函数给电脑计算却带来了不便, 而且现实中的期权执行价序列也是离散的多个有限点, 利用Trapeezoidal计算规则进行离散化可以得到:

$\int_{- \infty}^{\infty} e^{- i v k} φ (v) d v \approx \sum_{j = 0}^{Ν - 1} e^{- i v_{j} k} φ^{'} (v_{j}) η (27)$

特别注意到积分结果依赖于积分间隔η和积分上限a=Nη, 通过设立k、v和N的关系, 可以得到快速傅立叶形式:

$Y = \sum_{j = 0}^{Ν - 1} e^{- i \frac{2 π}{Ν} j} X (28)$

通过特征函数求期权价格的傅立叶形式, 并计算与执行价无关的序列X, 即 $h_{j} = φ^{'} (v_{j}) η e^{i \frac{Ν λ}{2} η j}$ , 再进行快速傅立叶逆变换, 是FRFT数值方法 (Fractional FFT) : $Y = \sum_{j = 0}^{Ν - 1} e^{- i \cdot 2 π j ε} X$ 的特殊情形, 令 $ε = \frac{1}{Ν}$ , 得到新的序列Y就是关于执行价的期权价值。通过CGMY的参数、不同的执行价、恒生指数、无风险利率、时间参数设定后, 历时不到2s, 数值计算结果如图1所示。

图1显示, 执行价在20000点以内, 实值期权 (in the money) 定价结果比较精确, 误差较小。以20000点为界的价外期权 (out of the money, 虚值期权) 深度虚值甚至出现负数, 不符合期权的实际价值, 需要进一步研究虚值范围内蒙特卡罗模拟定价效果, 假设K>S即为虚值期权。

4.3 蒙特卡罗模拟

通过第3节的模拟比较, 为了提高模拟效率, 一例以Kummer Complex函数计算合流超几何函数值。分别采用H和U淘汰法则来模拟生成CGMY分布如图2所示。

数据拟合统计检验采用MATLAB中t-test方法的统计检验见表3。

H0=1接受原假设, 通过t-test二维拟合检验。

现若进行期权定价只需要对C进行参数变换:

$C_{Τ} = C \cdot Τ = 254 C (29)$

风险中性修正的漂移率为:

CGMY过程的风险中性下模拟路径:

$S_{t} = S_{0} e^{(r - w) t + Ζ (t)} (31)$

欧式看涨期权蒙特卡罗模拟定价关系式:

$c (Κ) = E (\max (S_{t} - Κ, 0)) = \frac{1}{n} \sum_{i = 1}^{n} m a x ((S_{t} - Κ, 0)) (32)$

通过MATLAB编程实现模拟路径及期权价格计算, 设定路径数np=10000, 每条路径步数ns=100, 共需产生CGMY随机数1000000个, 历时小于30s. 而采用普通广义超几何函数进行计算, 因一周内产生死循环重新运行MATLAB数次, 故不做时间统计。通过模拟计算得到的期权价格参见图3。

如图3, 通过模拟比较, 无论哪种法则, 相对数值方法误差更小, 深度虚值领域也存在正的价值, 这比较符合期权理论和现实。相对数值结果, 模拟定价的效果要精确很多, 证明该方法的有效性。相比起数值方法, 这种模拟方法在虚值期权定价上更具先进性。同时, 两种法则中H法则下的模拟结果较为发散, 在恒生指数期权定价方面, Kummer函数显然比Francesco函数好。

5 结论

通过CGMY模型下数值和模拟两种常用定价技术的比较, 进一步提出改进该模型下蒙特卡罗模拟方法, 引入复数合流超几何函数计算, 以提高效率, 结果显示:①运用于恒生指数期权定价, 数值计算期权价格十分快捷和便利, CGMY模型下实值期权内定价精确有效, 但不适合虚值期权定价。②CGMY模型下蒙特卡罗模拟虚值期权定价比较有效, 深度虚值期权也不例外。相比起数值计算有很大的改进。误差较小, 速度较慢。③CGMY模拟中采用广义超几何函数或实数域内的合流超几何函数在Y<0.5的计算都是低效率的, 通过合流超几何函数的傅立叶变换使得计算变得快捷而且有效, 在Y越小时越明显, 采用复数域的合流超几何函数代替合流超几何函数计算可以提高计算速度。④Kummer微分方程的两个非线性相关解都属于随机数淘汰条件, 复数代替广义合流超几何函数, 在期权定价模拟采用Kummer方程第一类解, 相对第二类解稳定和收敛一些。

根据以上结论, 提出如下建议:①新函数计算速度不依赖于模型参数变化, 不因参数无穷小而影响计算速度, 极大提高随机数生成的效率, 从而避免随机数生成死循环。CGMY模型下, Y接近于零或者为负数的时候最为适用。②如果只从单向期权进行定价, 实值期权采用数值计算而虚值期权采用模拟的方法比较快捷和稳健。

蒙特卡洛模拟方法篇7

关键词：电力市场,蒙特卡洛,边际成本,强迫停运,负荷预测

0 引言

随着上世纪90年代开始的英格兰电力市场改革获得成功[1],世界范围内的电力工业正经历着市场化的改革,我国也在积极推进电力市场的改革进程[2]。目前,各国的电力市场大都采用的是分时竞价模式,该模式将一天24小时分为若干投标时段,每一时段对应着不同的电力价格。该价格由发电商报价和实时电力负荷确定。该市场模式的理论基础是基于系统短期发电边际成本(Short-run Marginal Cost of Power Generation)的实时电价理论[3]。该理论认为合理的实时电价应与系统短期发电边际成本接近。由此可见,分析系统短期发电边际成本是一项非常有意义的工作。对发电商来说,可以短期发电边际成本为基础进行策略报价[4],对市场监管者来说,短期发电边际成本是判断实时电价是否合理以及市场是否具有效率的重要依据。

电力市场中,市场参与者需要进行实时电价的短期或长期预测,而实时电价预测的重要依据之一就是系统短期发电边际成本。也就是说,市场参与者对未来某时间段的发电边际成本非常关注。而未来条件存在各种不确定性,如负荷需求预测的不确定性、发电机组随机强迫停运的不确定性、燃料价格的不确定性等。因此,考虑这些不确定性的电价预测将涉及短期发电边际成本(或边际报价)的概率学研究,不仅需要估计其期望和方差,而且需要确定其完全的概率分布。

目前报道的系统短期发电边际成本的概率学估计方法大多基于发电系统随机生产模拟技术。文献[5]指出,在实时定价中短期发电边际成本一般分为边际运行成本和边际缺电成本两部分。其中边际运行成本通常定义为负荷增加一单位后所需的发电运行成本(主要是燃料费用)的增量,边际缺电成本定义为负荷增加一单位后用户缺电成本的增量。文献[6]给出了一个估计发电边际成本概率分布的方法,但涉及负荷和边际发电成本概率分布的双变量累积量法表示,算法过于复杂且没有考虑边际缺电成本。文献[7-9]分别用发电系统随机生产模拟的累积量法和离散Markov状态空间法给出了发电边际成本均值和方差的近似解析表示。文献[10]利用微增技术,给出了随机生产模拟中系统边际成本均值的近似计算方法。文献[11]用发电随机生产模拟技术对边际缺电成本的概率学估计作了较为全面的研究。纵观这些研究,现有的方法大多基于发电系统随机生产模拟技术,该技术需要较详细的负荷资料以形成负荷持续曲线,而获取未来某一时间段的详细负荷资料有一定的困难;此外,现有的随机生产模拟技术只考虑机组强迫停运和负荷预测的不确定性,而没有计及燃料价格的不确定性。

本文考虑了电力生产过程中的各种随机因素,利用蒙特卡洛法对电能生产过程进行仿真。与传统的基于随机生产模拟技术的方法相比,该方法不需要详细的负荷资料,只需要进行简单的负荷预测,而且易于考虑电能生产过程中的各种随机因素,物理概念清晰,便于实现。

1 发电市场模型

假定某区域电力市场发电侧一共有n台机组竞价上网。市场中,交易中心以机组为出清单位,按统一的出清电价进行结算。市场某一投标时段的负荷需求为PL。以机组为研究对象,其成本函数为二次函数:

出力约束:

交易中心根据各机组的报价进行出清计算。如果各机组都按照发电成本进行报价(即按式(1)进行报价),此时的出清电价就是该时段的短期发电边际成本[12]。本文暂时只考虑系统负荷平衡,忽略机组的启停约束和电网的安全约束。交易中心代表所有电力用户向发电商购买电能,其目标是使电能采购费用f最小:

显然,这是一个求条件极值的问题,利用拉格朗日法可以求得该时段此时段短期发电边际成本和各机组出力:

需要注意的是,由式(6)决定的各机组出力Pi有可能不满足式(2)所给出的机组出力约束。此时,需要将越限机组的出力定为该机组的出力下限或上限,然后重新计算市场短期发电边际成本和各机组出力,直至所有机组都没有越限。

2 发电边际成本分析

虽然式(5)给出了某时段的短期发电边际成本,但以此为依据分析未来某时段短期发电边际成本时,该式中所有参数都是随机变量。参数ai,bi与燃料价格有关,是随机变量;PL为负荷预测的结果,也是随机变量。而且此时段并不是所有机组都能够投入运行,这些机组都存在随机强迫停运的问题。尤为重要的是,机组随机停运可能会导致市场中可用发电容量不足,系统中部分用户失去供电,此时就需要引入边际缺电成本这一概念进行分析,式(5)更是无法直接应用。考虑到这些随机因素后,系统短期发电边际成本就再也不是一个确定性的数值,无法直接由式(5)给出。固然可以从式(5)出发,推导短期发电边际成本的概率学表达式,但是这些算法非常复杂,且不便考虑边际缺电成本[6,7,8,9]。蒙特卡洛方法是进行随机模拟的有力工具[13]。为了分析系统短期边际发电成本,本文利用蒙特卡洛方法对电能生产过程进行仿真,模拟时可以很方便地考虑上文提到的各种随机因素。

2.1 随机因素分析

本文在模拟电能生产过程时主要考虑三种随机因素:机组的随机强迫停运,负荷的不确定性以及燃料价格的不确定性。

发电机组采用双状态模型,即正常运行状态和故障状态。用正常运行持续时间t1和检修时间t2描述。一般情况下工作时间和修复时间服从指数分布,即故障率λ和修复率µ是常数[14],得

式中:γ为服从均匀分布的随机数;tMTTF为平均工作时间;tMTTR为平均修复时间。

对未来某时段系统短期边际发电成本进行分析时,负荷数据由负荷预测获得。为了使模拟结果更精确,还需要在负荷预测的基础上分析负荷的概率分布情况。统计资料表明,负荷在预测值附近随机变动的概率分布属正态分布[15],且其标准差δ的百分值大体与系统容量的平方根成反比。例如,容量2 000 MW的系统,标准差δ约为3%;容量为20 000MW的系统,标准差δ约为1%。

与机组随机强迫停运与负荷的随机性一样,发电燃料价格的不确定性对短期边际发电成本的估计同样具有重要影响。燃料价格的随机波动会引起机组发电成本函数中参数的随机波动。因此,为了考虑燃料价格的不确定性,应把各机组成本函数式(1)中参数ai,bi描述为随机变量。显然,随机变量ai,bi的概率分布由发电燃料(主要是石油和煤炭)价格所服从的概率分布决定。而研究燃料价格所服从的概率分布是一项很有挑战性的课题,经济学家们也提出了相应的预测模型。本文不打算对此问题进行深入研究,此处假定所有机组采用同一种燃料,燃料价格考虑为一个随机变量,它有三个等概率的取值ρ1,ρ2,ρ3,均值为ρ2。

2.2 蒙特卡洛模拟

蒙特卡洛方法是进行随机模拟的有力工具,本文将该方法应用于发电模拟。为保证模拟结果有足够精度,需要进行足够次数的模拟,在模拟结束后求得系统的边际运行成本和边际缺电成本。计算边际缺电成本时需要知道系统失负荷概率和系统单位缺电损失。其中,系统失负荷概率可以由模拟结果获得,系统单位缺电损失可取为每度电的使用价值,它与社会的劳动生产率密切相关,可近似估计为系统的国民生产总值除以总发电量。这里假设为所有机组中最大单位运行成本的120%。算法流程如图1所示。

3 算例分析

为了验证算法的有效性,本文利用C++语言编了计算程序,并以某实际系统部分机组构成的区域电力市场为算例[12,13]进行了仿真计算。该区域电力市场由10台机组组成,总装机容量为1 980 MW,机组参数如表1所示。该表中di,ai,bi为机组成本函数式(1)中的相关参数。参数ai,bi的取值与燃料价格相关,为服从某种分布的随机变量,表1中给出的数值为变量ai,bi的期望EXai,EXbi。此数值可以预测燃料价格期望值为基础,通过分析各机组的耗量特性获取。表中所有机组的强迫停运率(Forced Outage Probability)一样,取为0.05;负荷服从正态分布,标准差δ为系统容量的3%;参数ai,bi有三个等概率的取值[0.9EXai,EXai,1.1EXai]和[0.9EXbi,EXbi,1.1EXbi]。

为了研究算法的收敛性,此处研究两种负荷水平下(负荷预测的期望分别为1 200 MW和1 600MW),系统短期发电边际成本期望与模拟次数之间的关系,结果如图2所示。

从图2可以看出本文所提出的仿真算法收敛性良好,两种负荷水平下,当模拟次数达到5 000次时都能可靠收敛,在以下的计算中,模拟次数均设为5 000次。

图3给出了不同负荷水平下的系统短期发电边际成本期望和边际缺电成本期望。从图3可以看出随着负荷的增大,系统短期发电边际成本期望逐步增大,其中边际缺电成本在总发电边际成本中所占的比例呈显著增加的趋势。此分析结果是符合客观事实的,在负荷水平较低的时候,系统失负荷概率不大,因此边际缺电成本在总发电边际成本中所占的比例比较小,而随着负荷的增大,系统是负荷概率逐渐增大,边际缺电成本所占比重呈显著增加趋势。

由于各种随机因素的存在,系统短期边际发电成本为随机变量。在进行电价预测时,仅仅掌握其期望是不够的,还要研究其概率分布函数。算法可以在对各次模拟结果进行统计的基础上求出系统短期边际发电成本的概率分布函数。图4给出的是负荷为1200 MW时系统短期边际发电成本的概率分布函数。

4 总结

【蒙特卡洛模拟方法】推荐阅读：

蒙特卡洛06-09

蒙特卡罗模拟08-16

卡洛斯-06-16

皮埃蒙特07-15