蒙特卡罗方法及应用

2024-08-16

蒙特卡罗方法及应用（共7篇）

蒙特卡罗方法及应用篇1

1.引言

近些年来, 随着电子计算机的迅速发展与广泛使用, 计算数学获得了新的发展, 增添了许多新内容, 蒙特卡罗方法就是其中一种。该方法的出现, 极大地促进了科学文明的进步, 开辟了人们研究问题的新途径。

在程序编制和计算机算法设计中常采用的算法一般是确定性算法, 也就是算法的每一个计算步骤都是确定的。但是在很多情况下, 当算法在执行过程中遇到一个选择时, 随机性选择经常比最优选择更节省时间且效率极高。我们把这种允许算法在执行过程中随机选择下一个计算步骤的算法叫做概率算法。概率算法可以很大程度地降低算法的复杂程度。

蒙特卡罗算法、拉斯维加斯算法和舍伍德算法都是一些常用的概率算法。这些算法的基本特征是:对要求问题的同一个例子, 用同一概率算法运算若干次, 所得到的结果可能完全不同。但是通过多次执行反复求解, 会使正确性和精确性达到令人满意的效果, 而失败或误差的概率则可以接近无限小[1]。

本文的主要研究内容:

(1) 蒙特卡罗方法的基本思想及其基本特点。

(2) 蒙特卡罗方法求解问题的基本过程。

(3) 蒙特卡罗方法的应用领域。

(4) 蒙特卡罗方法在计算曲线积分中的应用。

2.蒙特卡罗方法的基本思想及其基本特点

2.1蒙特卡罗方法的提出

二十世纪四十年代的二战中美国有个研制原子弹的计划———“曼哈顿计划”, 该计划的成员J.冯·诺伊曼和S.M.乌拉姆首先提出蒙特卡罗方法。著名数学家冯·诺伊曼采用世界有名的赌城———摩纳哥的Monte Carlo命名这一方法。而实际上, 在这之前, 蒙特卡罗方法就已经存在。早在1777年, 法国著名数学家布丰 (Georges Louis Leclere de Buffon, 1707-1788) 就提出用投针实验的方法计算圆周率π。这被认为是蒙特卡罗方法的起源。

2.2蒙特卡罗方法概述

蒙特卡罗方法又称为随机抽样技术、统计模拟法, 是一种随机模拟方法, 它是以概率和统计理论方法为基础的一种计算方法, 是使用随机数 (或更常见的伪随机数) 解决很多计算问题的方法。将所求解的问题同一定的概率模型联系起来, 用电子计算机实现统计模拟或抽样, 从而获得问题的近似解。为了象征性地表明这一方法的概率统计特征, 所以借用赌城蒙特卡罗命名。

2.3蒙特卡罗方法基本思想

当所要求解的问题是某种随机事件出现的概率, 或者是某个随机变量的期望值时, 通过某个“实验”的方法, 以这种事件出现的频率估计这个随机事件的概率, 或者得到的这个随机变量的某些数字特征, 并将它们作为问题的解。

3.蒙特卡罗方法求解问题的基本过程

蒙特卡罗方法求解问题的过程大致分为三个主要的步骤。

3.1描述或构造概率过程

对于那些本身就有随机性质的问题, 比如粒子的输运问题, 主要是要对这个概率的过程进行正确的描述和模拟;对于确定性的问题, 比如说计算数学上的定积分, 就一定要事先人为地构造一个概率过程, 而且这个概率过程的一些参量刚好是所要求的问题的解。也就是要把不具有随机性质的问题转化为有随机性质的问题。

3.2实现从已知概率分布抽样

构造了概率模型以后, 因为各种概率模型都可以看做是由各种各样的概率分布而组成的, 因此产生已知概率分布的随机变量 (或随机向量) , 就成为实现蒙特卡罗方法模拟实验基本的手段, 这是蒙特卡罗方法被称为随机抽样的原因。最基本、最简单、最重要的一个概率分布是 (0, 1) 上的均匀分布 (或称矩形分布) 。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样, 也就是一个具有此种分布的相互独立的随机变数序列。产生随机数的问题, 就是从这个分布的抽样问题。在计算机中可以用物理方法产生随机数, 但是价格比较昂贵, 且不能重复, 使用不方便。另一种方法是用数学的递推公式产生。这样产生的序列, 与真正的随机数序列不同, 所以称为伪随机数, 或伪随机数序列。不过, 经过多种统计检验表明, 它与真正的随机数序列或随机数具有很相似的性质, 所以可将其作为真正的随机数应用。由已知分布随机抽样有很多的方法, 不同于从 (0, 1) 上的均匀分布抽样的是, 此类方法都是凭借于随机序列从而得到实现, 也就是说, 其前提都是产生随机数。因此, 实现蒙特卡罗方法的基本工具是随机数。

3.3建立各种估计量

一般来说, 概率模型被构造出来并可以从中抽样以后, 一定要确定一个随机变量, 作为所求的问题的解, 我们把它叫做无偏估计。建立各种估计量, 也就是对模拟实验的结果进行观察和记录, 从而把问题的解求出来。

4.蒙特卡罗方法的应用领域

通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的问题。对于那些因为计算太过复杂而很难得到解析解或根本没有解析解的问题, 蒙特卡罗方法是一种有效的求解出数值解的方法。

除了在数学方面得到很好的应用外, 蒙特卡罗方法在其他很多领域也得到很好的应用, 比如说蒙特卡罗方法在金融学、工程学、宏观经济学、地质学、生物医学、计算物理学 (如粒子输运计算、空气动力学计算、量子热力学计算) 及计算机科学等广泛的领域都得到非常成功的应用。

下面探讨蒙特卡罗方法在数学领域中的应用。

5.蒙特卡罗方法在计算曲线积分中的应用

在物理的很多研究时, 要求平面或空间中的一条可求长度的曲线形物体的质量, 或者一个质点在平面或者空间中沿曲线L从A点运动到B点时力F所做的功, 这时就要用到曲线积分。

平面或空间曲中的曲线积分有两种类型: 第一型曲线积分和第二型曲线积分。下面着重探讨如何用蒙特卡罗法计算第一型曲线积分。

5.1第一型曲线积分的定义

设L为平面上可求长度的曲线段, f (x, y) 为定义在L上的函数。对曲线L做分割T, 它把L分成n个可求长度的小曲线段Li (i=1, 2, …n) , Li的弧长记为Δsi, 分割T的细度为, 在Li上任取一点 (ξi, ηi) (i=1, 2, …, n) , 若有极限,

且J的值与分割T与点 (ξi, ηi) 的取法无关, 则称此极限为f (x, y) 在L上的第一型曲线积分, 记作:

5.2第一型曲线积分的一般计算方法

设有光滑曲线

函数f (x, y) 为定义在L上的连续函数, 则

5.3第一型曲线积分的蒙特卡罗计算法

对于L上的曲线积分, 其中。令f (t) =, 则

在[α, β]上构造一个均匀分布的随机变量t, 则密度函数为

又令

则

(2) 式给我们传达了:F (t) 在t∈[α, β]上的数学期望就是 (1) 式的曲线积分[2]-[4]。

现在只要抽取概率密度为p (t) 的随机变量t的容量为N的随机样本ti (i=1, 2, …, n) , 由大数定理就可以得到:

5.4实例

设L是半椭圆

求第一型曲线积分。

首先, 由曲线积分方法得:

这样的积分很难直接求得结果, 这里考虑用蒙特卡罗法计算:

用计算机程序模拟10次得到随机表如下表所示:

模拟1000次得到随机表如下表所示:

模拟10000次得到随机表如下表所示:

模拟20000次得到随机表如下表所示:

模拟30000次得到随机表如下表所示:

由以上的模拟数据可以看出, 模拟次数越多, 得出的积分均值越趋于稳定。

5.5三维空间的第一型曲线积分

上面讨论的是平面上曲线段的曲线积分, 下面开始探讨蒙特卡罗法如何应用于三维空间中的曲线积分。对于三维空间中的光滑曲线段:

函数f (x, y, z) 为定义在L上的连续函数, 则

类比平面上的曲线积分, 很容易等到三维空间中, 曲线积分的蒙特卡罗算法:

6.结论和展望

与其他的计算方法相比较, 蒙特卡罗方法有下面几个优点:

(1) 问题的维数不影响该方法的收敛速度。但是很多的数值计算方法, 比如计算N重积分时, 达到相同的误差, 点数和维数的幂次成正比。

(2) 受问题的条件限制的影响不大。

(3) 程序的结构简单, 在计算机上实现蒙特卡罗计算时, 程序的结构简单清晰, 占用内存单位较少。

当然, 蒙特卡罗方法并不是完美的, 它也有其缺陷的地方:

(1) 伪随机数的均匀性影响随机变量取值从而影响结果。

(2) 收敛速度慢。蒙特卡罗模拟的收敛速度为O (n-1/2) , 一般不容易得到精确度较高的近似结果。对于维数比较少 (三维以下) 的问题, 不如其他方法好。

(3) 如误差是概率误差, 误差与点数 (抽样数) 的平方根成反比, 而其他数值方法的误差是确切的, 一般情况下, 误差与点数成反比, 因此, 对于维数不高的问题, 数值方法可以给出误差很小的结果, 而且计算量小。

蒙特卡罗方法与其他数值方法都有其本身的一定适用范围, 把二者结合起来, 取长补短, 发挥它们各自的长处, 目前在国外已受到普遍重视, 是蒙特卡罗方法发展的一个重要方向。

摘要：本文主要介绍蒙特卡罗方法的基本原理及思想特征, 详细讨论了蒙特卡罗方法在计算曲线积分等实际问题中的应用, 突出了蒙特卡罗方法的过程及优势。

关键词：蒙特卡罗方法,曲线积分,教学应用

参考文献

[1]王晓东.算法设计与分析[M].北京:清华人学出版社, 2008.

[2]黎锁平.运用蒙特卡罗方法求解随机性问题[J].甘肃工业大学学报, 2001, 27 (2) :95-97.

[3]姚贵平, 等.重积分近似计算方法的讨论[J].内蒙古农业大学学报, 2000, 21 (4) :98-101.

[4]宫野.计算多重积分的蒙特卡罗方法与数论网格法[J].大连理工大学学报, 2001, 41 (1) :20-23.

[5]尹增谦, 等.蒙特卡罗方法及其应用[J].物理与工程, 2002, 12 (3) :45-49.

[6]易大义, 陈道琦.数值分析引论[M].杭州:浙江大学出版社, 1998, 165-170.

[7]Sobol I M.Asymmetric convergence of approximations of the Monte Carlo method[J].Computational Mathematics and Mathematical Physics, 1993, 33:1391-1396.

蒙特卡罗方法及应用篇2

凝析气藏作为一种特殊类型的油气藏类型, 在世界气田开发中占有特殊重要的地位。塔里木盆地作为我国油气资源的重要产地, 凝析气分布广泛, 成因复杂, 以富、特富凝析气藏为主, 具有很高的经济价值。由于凝析气藏特殊的流体相态特征, 导致开发难度要比一般油气藏复杂、建设投资高、风险大。因此为了避免决策失误, 对于凝析气藏的项目开发应该做好开发方案的经济评价工作, 加强开发方案的风险分析。

经济评价风险分析方法主要有盈亏平衡分析、敏感性分析、蒙特卡洛模拟分析等。目前我国石油行业经济评价一般采用盈亏平衡和敏感性分析。传统的盈亏平衡分析和敏感性分析都是考虑单因素变动对评价指标影响, 是一种简化的分析方法。在现实生产中, 评价指标往往受到多个因素共同变化的影响, 所以传统的方法对于风险分析具有局限性和不合理性。相比于传统的方法, 蒙特卡洛模拟能综合考虑多种风险因素, 不仅能得到评价结果, 也能得到对应的发生概率, 风险结果更可靠、更贴近实际、更能体现经济评价的合理性。因此, 本文采用了基于蒙特卡洛模拟方法进行凝析气藏开发方案的经济评价风险分析, 并编程实现了该风险分析工具。结合凝析气藏油气同采的特点, 该方法考虑了对凝析气藏开发重要的建设投资、经营成本、凝析油和天然气产量、销售价格风险因素。通过模拟各风险因素的随机分布数值, 根据现金流量法计算经济评价指标值及其概率。经济评价指标采用税后财务内部收益率 (FIRR) 和净现值 (FNPV) 。最后利用量化的经济评价指标概率分布结果来判断开发方案的风险情况。

二、蒙特卡洛模拟风险分析

(一) 蒙特卡洛模拟分析原理

蒙特卡洛模拟分析又称统计模拟法、随机抽样技术。它是一种随机模拟方法, 是一类通过随机抽样和统计试验来求近似解的数值方法。蒙特卡洛模拟的基本思路:首先建立分析指标值的计算模型;分析与指标值相关的影响因素, 并分析建立影响因素的概率分布模型;依据模拟次数, 模拟产生相应批次的影响因素随机数, 带入分析指标值的计算模型, 计算出相应的指标值;最后对模拟试验结果加以分析统计, 计算评价指标值的统计特征量或概率分布。模拟次数越多, 概率分布越接近真实的分布。

(二) 蒙特卡洛模拟风险分析工具实现步骤

结合凝析气藏开发建设项目特点, 凝析气藏开发方案经济评价蒙特卡洛模拟风险分析具体实现步骤如下:

1. 确定经济评价指标模型。

油田开发建设项目经济评价指标主要包括:财务内部收益率 (FIRR) 和财务净现值 (FNPV) , 利用现金流量法计算。财务内部收益率是指能使项目计算期内净现金流量现值累计等于零时的折现率, 计算公式如下所示。

式中:CI—现金流入量

CO—现金流出量

(CI-CO) t—第t期的净现金流量

n—项目评价周期

财务净现值 (FNPV) 是指按设定的折现率, 一般采用基准收益率 (ic) 计算的项目期内净现金流量的现值之和, 计算公式如下所示。

现金流量法计算的关键内容是要计算项目评价周期内各年现金流入和现金流出值。针对油气开发建设项目, 现金流入包括:营业收入、增值税销项税、补贴收入、回收油气资产余值以及回收流动资金;现金流出包括:利用探井评价井投资、建设投资、流动资金、经营成本、增值税进项税、营业税金及附加、增值税、维持运营投资以及弃置费。

2. 分析影响经济评价指标的风险因素。

由于凝析气藏油气同采的特点, 该类油气藏的开发建设经济评价不同于一般油气项目。凝析气藏是介于油藏和纯气藏之间的复杂类型的特殊油气藏。凝析气田开发, 具有投资大、经营成本高、风险大的特点, 并且更为突出的特点是油气同采。油气同采的特点决定了油和气产量及其销售价格都会对项目最终效益产生影响, 因此油和气产量及价格都是重要的不确定性因素。并且不同油气比的项目, 油、气对项目的影响程度是不同的。结合凝析气藏开发项目特点和实际情况, 我们分析影响FIRR和FNPV的主要因素包括:建设投资、石油产量、石油销售价格、天然气销售产量、天然气价格以及经营成本六个因素。

建设投资:包括开发井投资和地面工程投资。对于油气开发建设项目, 建设投资非常大, 投资估算存在不确定性。建设投资作为现金流出的一个重要项目, 同时也会影响其他项目。受建设投资影响的项目包括:现金流入的回收油气资产余值、现金流出的增值税和弃置费。

石油产量:凝析气田开采伴随着产出凝析油。由于受到反凝析影响, 凝析油产量具有不确定性。对于凝析油产量高的气藏, 凝析油产量对项目的效益影响很大。石油产量直接影响现金流入中的营业收入, 同时影响开发方案的操作成本估算, 从而影响现金流出的经营成本。间接导致现金流入的增值税销项税, 现金流出的增值税、营业税金及附加等项目变化。

石油销售价格:石油销售价格受国际油价影响, 波动较大, 具有不确定性, 往往是影响项目效益最敏感的因素之一。石油价格直接影响营业收入和吨油特别收益金的计算。并间接影响其他项目计算。

天然气产量:凝析气藏开发主要商品是天然气, 由于受到油藏地质条件、开发方案等影响, 天然气实际产量往往与设计产量存在一定的偏差。因此在利用设计的方案产量进行评价时需要考虑天然气产量的不确定性。天然气产量与石油产量类似, 直接影响着营业收入和操作成本。

天然气销售价格:目前我国天然气销售价格是由政府定价决定, 对于未来年份价格估计需要考虑到不确定性。天然气价格直接影响营业收入。

经营成本:油气经营成本受到生产工艺、产量、措施、市场、管理等各方面因素影响, 不确定性很大。经营成本的变化直接影响流动资金和回收流动资金。

3. 建立风险因素随机模型。

结合开发项目实际情况和历史开发数据, 考虑各风险因素的随机模型。风险变量的概率分布模型主要包括均匀分布、三角分布和正太分布等。随机模型选取方法:当风险变量只能获得一个范围值时, 可采用均匀分布公式来描述;当风险变量除取得范围值外, 还知道最可能值, 则用三角分布公式来描述, 三角分布是采用较多的;当风险变量获得少量的随机值, 则根据多数风险变量具有正态分布或对数正态分布的特征, 可模拟为正太或对数正太分布公式来描述。

4. 模拟计算。

首先根据开发方案的基础参数计算出项目基本的FIRR和FNPV。然后根据风险因素随机函数产生随机数, 利用随机数与开发方案基础数据计算评价指标, 并根据设定的模拟次数 (N次) 完成计算, 得到N个FIRR和FNPV。

5. 分析统计。

对得到的N个FIRR和FNPV进行数理统计, 分析计算出评价指标的均值、标准差、最大值、最小值、经济风险值等数据, 绘制概率分布图。设定界限值, 计算大于或小于界限值概率。

三、方法应用分析

塔里木盆地某开发建设项目, 评价周期为20年, 其中建设周期为2年, 运营周期为18年。资本金比例为55%, 流动资金贷款比例为70%, 经营成本资金率为20%。原油商品率为96%, 天然气商品率为94%。石油销售价格为3577元/吨, 天然气销售价格为1060元/千方。税后基准收益率为12%, 还款方式采用等额本金, 还款期为8年。建设投资分为两年投入, 总额为58亿元, 其中开发井投资为46亿, 地面工程投资为12亿。维持运营投资48亿元。生产经营参数根据开发方案计算给出评价周期内每年的产量、产液量、开井数等参数。其中评价周期内石油总产量为898万吨, 天然气产量为16亿立方米。

为了分析计算, 本文所有的影响因素均采用三角分布模拟。经过专家估算, 风险因素的模拟参数为:建设投资模拟建设期总投资, 最小值为56亿元, 最大值为60亿元, 最可能值为58亿元;石油价格模拟生产期第一年销售价格, 最小值为3400元/吨, 最大值为3650元/吨, 最可能值为3577元/吨;石油产量模拟生产期平均产量, 数据为生产期年均产量, 最小值为43万吨/年, 最大值为48万吨/年, 最可能值为45万吨/年;天然气价格模拟生产期第一年销售价格, 最小值为800元/千方, 最大值为1200元/千方, 最可能值为1060元/千方;天然气产量模拟生产期平均产量, 数据为生产期年均产量, 最小值为0.6亿方/年, 最大值为0.9亿方/年, 最可能值为0.78亿方/年。模拟次数为10000次。通过Excel宏编程, 实现了基于蒙特卡洛模拟的凝析气藏开发方案经济评价风险分析计算工具。利用该工具进行模拟, 模拟结果如图1和表1所示。

计算得出, 该项目方案的税后投资财务内部收益率为14.32%, 投资财务净现值为4.97亿元。通过模拟计算, 净现值小于4.97亿元的概率为43.68%, 均值为53514万元, 虽然净现值最小值为负值, 但是通过计算小于0的概率只有1.8%, 经济风险=27188/53514=0.51>0.5, 说明项目经济风险较低;内部收益率超过14.32%的概率是57.1%, 大于税后基准收益率12%的概率为98.2%, 最大值为18.43%, 最小值为10.38%, 平均内部收益率为14.55%高于税后基准收益率, 经济风险=1.31/14.55=0.09<0.12, 项目经济效益较低。由此可以得出, 该项目的净现值大于0、收益率大于基准收益率的可能性很大, 因此该项目盈利能力可能性很大, 项目经济风险相对较小, 但是该项目的经济效益较低, 评估需谨慎。

四、结论

油气开发项目存在很大风险和不确定性, 因此油气开发项目经济评价必须要做好风险分析工作。凝析气藏由于其油气同采特点, 经济评价风险分析考虑的风险因素与一般油气藏不同, 需充分考虑凝析气藏石油和天然气产量及价格对方案风险的影响作用。本文采用了蒙特卡洛模拟方法进行凝析气藏开发方案经济评价风险分析, 通过以上应用分析表明, 蒙特卡洛模拟方法更加适合应用于油气开发项目风险分析, 它实现了多风险因素同时变化对评价指标的影响, 通过对评价指标概率计算, 提高了决策可靠性。

参考文献

[1]郭平, 李士伦, 杜志敏, 孙雷等.凝析气藏开发技术现状及问题[J].新疆石油地质, 2002 (3)

[2]孙志道等.凝析气藏早期开发气藏工程研究[M].北京:石油工业出版社, 2003

[3]孙龙德, 宋文杰, 江同文.塔里木盆地凝析气田开发[M].北京:石油工业出版社, 2003

[4]朱伟, 刘益超.石油建设项目投资风险评价方法探讨[J].天然气工业, 2000 (6)

[5]雷利安, 赵玉萍.应用蒙特卡洛法进行油气开发方案风险评价[J].石油天然气学报, 2006 (2)

[6]张明泉, 钟雄., 蒙特卡洛模拟在油田开发经济评价风险中的应用[J].西南石油大学学报, 2012 (4)

控制系统蒙特卡洛方法应用分析篇3

蒙特卡洛方法是利用随机数进行统计试验,求得统计特征值(如均值、概率等)作为待解问题的数值解,所作的统计试验称为蒙特卡洛仿真或者蒙特卡洛模拟。该方法是众多研究不确定性传递的方法之一,目的是要确定不确定性(包括系统参数变化、建模误差、外界扰动等)如何影响系统性能,其常用于复杂的非线性模型以及具有不确定参数的模型。

蒙特卡洛方法可以解决各种类型的问题,一般来说,该方法处理的问题可以分为两类:第一类是确定性的数学问题。包括计算多重积分、求解逆矩阵、线性代数方程组、积分方程和微分算子的特征值等;第二类是随机性问题,对于这类问题,有时可表示为多重积分或某些函数方程,并进而考虑用随机抽样方法求解。然而,一般情况下都不采用这种间接模拟方法,而是采用直接模拟方法,即根据实际物理情况的概率法则,首先建立一个与求所解有关的概率模型,使所求的解为所建立模型的概率分布或数学期望,然后对该模型进行随机抽样观察,即产生随机变量,最后用其算术平均值作为所求解的近似估计值。

当前,蒙特卡洛仿真及分析方法广泛应用于控制系统设计及分析中,包括制导与控制系统设计缺陷的分析、轨道及特定飞行器参数的识别[1];飞机鲁棒性分析中实时稳定裕度的测量[2];鲁棒自动着陆控制律的优化设计[3,4];电传飞控系统的飞行品质研究[5]等方面。本文对控制系统蒙特卡洛仿真的工作流程及关键问题进行分析,并给出一种控制系统的蒙特卡洛仿真方案。

1 蒙特卡洛仿真工作流程

Matlab环境下的蒙特卡洛仿真工作流程如图1所示。

蒙特卡洛工作流程为:(1)首先确定系统仿真参数,仿真次数可以通过理论分析方法或试验方法确定;(2)根据不确定性模型生成全部仿真所需的随机偏差和初始状态;(3)以预设的不确定偏差及初始状态进行仿真调用,并记录所需的结果及状态数据;(4)仿真全部结束后,对记录的数据进行分析,根据分析结果对仿真模型、控制系统设计进行相应的改进。在蒙特卡洛仿真过程中,需要重点考虑以下问题。

1.1 不确定性建模

大部分实际系统中总是包含有不确定性,蒙特卡洛仿真中需要采用随机数进行统计试验,所以需要建立系统不确定输入的概率模型,即确定不确定输入的分布类型。对于部分不确定性,根据经验或理论分析,可以知道其概率分布形式,但大部分情况概率分布是未知的。对于这类不确定性,较为普遍的方法是采用不确定性参数的标称值均值的正态分布作为不确定性的概率模型。但有些情况下,不确定性明显不服从正态分布,这就需要根据观测样本进行不确定性建模。假定实际系统不确定性的观测数据为互相独立的随机变量,不确定性建模步骤为:

(1)选择概率分布类型。选择概率分布的依据是随机变量的知识以及占有的数据。很多情况下,可以根据先验知识以及对产生该随机变量过程的认识决定或者拒绝哪种分布。另外,数据对概率选择具有决定性作用。一般有了足够的数据就可以选择出较好的概率分布,至少可以作出经验分布或者半经验分布,而且可以用数据检验所选择概率分布的正确性。对于统计量连续分布的情况,可以采用点估计法、直方图法和概率图法确定随机变量的概率分布类型。

(2)参数估计。初步确定概率分布类型后,需要确定出该分布的参数,常用的参数估计方法为矩估计法和极大似然估计法。

(3)拟合优度检验。选择参数分布类型,并估计具体参数后,应当进行分布的拟合优度检验,即检验所选择的具体分布是否符合给定的数据分布。拟合优度检验是一种统计假设检验,用于检验观测样本是否服从于概率分布函数的理论分布。常用的拟合检验方法包括检验法和Kolmogorov-Smirnov检验法。

1.2 仿真调度机制

蒙特卡洛仿真通常需要将一个仿真任务执行成千上万次,要耗费大量的时间,使得系统设计周期过长,因此并行计算是进行蒙特卡洛仿真的一个必然选择。Matlab的分布计算工具箱提供了利用计算机集群进行并行计算的功能。

分布计算工具箱通过MATLAB分布式计算引擎可支持开发能在多台计算机上运行的分布式MATLAB计算应用程序,无需改换原有的MATLAB开发环境。用户可先在MATLAB中完成计算程序的开发与调试,然后利用该工具箱将算法分解为多个相互独立或相关的任务。MATLAB分布式计算引擎在多个远程MATLAB进程中同步调度和处理这些任务,相比单一的MATLAB进程执行,运行时间大大缩短。

MATLAB分布计算工具,仿真系统由以下部分组成:(1)客户端:定义项目,将项目分解成多个任务,是用于进行MATLAB编程的计算机;(2)项目管理器:负责项目管理和任务分发,用户可以通过项目管理器将项目提交到项目队列以及将任务分发到工作机;(3)工作机:工作机接收项目管理器分配的任务,并进行计算。仿真系统各部分之间的交互如图2所示。

Matlab分布计算工具应用分为分布式执行以及并行计算,在应用分布计算工具进行蒙特卡洛仿真时,蒙特卡洛仿真作为用户提交的项目,可以自然地分为各个单次仿真任务,各单次仿真相互独立,因此属于分布式执行。

1.3 蒙特卡洛分析方法

对控制系统进行蒙特卡洛仿真本身不是目的,蒙特卡洛方法中重要的是应对蒙特卡洛仿真的结果进行充分有效的分析,使其能更好地服务于控制系统的设计及性能优化。对控制系统可以开展以下几方面的蒙特卡洛仿真分析。

(1)概率统计。首先确定任务成功的准则,然后通过蒙特卡洛仿真可以得到任务成功概率,能够直接对控制系统的性能进行评价。

(2)仿真模型验证。针对仿真模型中的某些错误和缺陷,单次仿真通常难以发现。因此对控制系统仿真模型,加入随机扰动及随机初始状态,进行蒙特卡洛仿真。在进行系统验证时,为充分暴露仿真系统的缺陷,设定的随机扰动的幅值可大于正常扰动幅值。通过对仿真统计结果和失败案例的概率分析,找出仿真系统的设计缺陷。

(3)控制系统灵敏度分析。采用灵敏度分析技术,确定特定输出与随机输入变量之间的灵敏度。采用最小二乘估计回归建立线性模型,使用多重测量的相关系数评估灵敏度。采用线性回归方法搜索仿真输入空间,找出对仿真输出影响最大的变量。

(4)定量的系统鲁棒性分析。采用蒙特卡洛方法计算系统的随机鲁棒性,得到定量的稳定性描述。

(5)控制器参数优化。在参数最优化计算中,当目标函数为非线性、存在多峰值或维数较高时,往往受到限制。这时应用蒙特卡洛方法比较有效,包括随机模拟算法,序贯随机模拟算法及随机搜索法等[6]。

(6)系统稳定性及飞行品质分析。对于电传飞控系统,通常考核它的动态特性和飞行品质,即输入驾驶员杆指令,考察时间响应、短周期阻尼、CAP、延迟时间等品质指标。通过蒙特卡洛仿真,可获得阶跃响应时间常数、稳态误差、周期阻尼、CAP、延迟时间等指标的散布点,计算各品质指标的数学期望、方差,可获得飞行品质的概率估计值。

2 控制系统蒙特卡洛仿真方案

本文建立一种基于MATLAB分布计算工具的蒙特卡洛仿真系统,主要由客户端以及分布式计算机网络构成,仿真系统逻辑结构如图3所示。

2.1 客户端

客户端需安装并启用分布式计算引擎服务,并安装MATLAB分布计算工具箱,用于创建分布式仿真作业。用户在客户端计算机的Matlab中构建Matlab仿真模型,创建蒙特卡洛仿真作业,并将作业提交给项目管理器,由分布式计算网络进行并行的仿真计算。客户端工作流程如图4所示。

最终的系统仿真结果会返回到客户端,因此蒙特卡洛分析工作也在客户端计算机进行。

2.2 分布式计算网络

分布式计算网络中由仿真管理节点和若干仿真节点组成。该网络中的计算机都需要安装并启用MATLAB分布式计算引擎服务,管理节点运行作业管理器进程,仿真节点运行工作机进程。仿真管理节点负责仿真调度,它接收客户端提交的作业,将作业所包含的仿真任务自动分配到各个仿真节点进行计算,收集各工作机的计算结果。所有任务完成后将仿真结果返回到客户端。分布计算网络环境下的蒙特卡洛仿真整体工作流程如图5所示。

3 结语

可知蒙特卡洛方法本质上是一种基于参数随机抽样的数值计算方法,通过这种方法得出的结果可信度基于两个条件:(1)应建立可靠的数学模型,并将此数学模型与实际数据进行检验,尽可能在不使模型过于复杂的情况下令模型符合实际数据;(2)对模型中的随机量分布有明确的了解,进而针对不同的分布采取合理的随机抽样方法。此外,如果条件许可,应该对计算得出的统计结果与实际对照。蒙特卡洛方法由于能够验证复杂系统的环境变化、子系统、部件及全部参数同时变化下的系统特性,可以进行参数的遍历性试验,因此具有较高的可信度,目前在国内外航空航天领域获得广泛应用,尤其在新型号飞行器控制系统的研制过程中,成为一种有效的、必备的检测和验证工具。

摘要：鉴于蒙特卡洛方法在控制系统设计及分析中的广泛应用,对蒙特卡洛仿真的工作流程及关键技术问题进行分析,并给出一种基于MTTLAB分布式计算引擎的控制系统蒙特卡洛仿真方案。

关键词：控制系统,蒙特卡洛,分布式计算,项目管理器

参考文献

[1]PEGGY S,WILLIAMS.A monte carlo dispersion analysis of the X-33simulation software[Z].NASA Dryden Flight Research Center,AIAA,2001.

[2]JOHN T,BOSWORTH,SUSAN J STACHOWIAK.Real-time stability margin measurements for X-38robustness analysis[Z].NASA Dryden Flight Research Center,2005.

[3]GERTJAN LOOYE,HANS-DIETER JOOS,DEHLIA WILLEMSEN.Application of an optimization-based design process for robust autoland control laws[Z].AIAA,2001.

[4]GHN,LOOYE,SANCHO,et al.Design of robust autoland control laws using mu-Synthesis[C].AIAA Guidance,Navigation,and Control Conference and Exhibit,2002.

[5]杨云,张平.电传飞控系统的蒙特卡洛分析与设计[J].航空学报,2008,(29):18-23.

蒙特卡罗方法及应用篇4

随机潮流是当前电力系统稳态分析和规划决策的重要工具,与传统确定性潮流相比,它能够考虑各种不确定因素,从而更全面深刻地揭示系统的运行特性及薄弱环节[1,2]。随着分布式电源、电动汽车的大规模接入,储能技术的发展及电力市场化改革的稳步推进,电网将呈现出较大的随机性和高维特征。因此,加快建立适应新环境的电网分析和评估方法具有重要意义[3,4,5]。

随机潮流自20世纪70年代由Borkowska[6]提出以来,国内外学者做了大量研究并提出了多种计算方法,大致可归纳为三类:模拟法、近似法和解析法。以蒙特卡洛仿真(MCS)为代表的模拟法[7]通过简单随机采样(simple random sampling,SRS)得到输入变量的样本,再对每个采样点进行确定性潮流计算,最后统计获得各状态量的分布情况。整个过程原理简单且适用性广,在保证样本规模足够大的情况下能够获得很高的精度,不足之处是仿真次数多,耗时长,一般用于评估其他方法的优劣。近似法,如点估计法[8,9]以及解析法中的快速傅里叶变换[10]、半不变量法[11]等都是利用随机变量的数字特征并结合数值方法来求解状态量的分布,有效地减少了计算规模和时间,但当考虑输入变量相关性时,公式会非常复杂且高阶矩的误差较大。文献[12]创建了一种新的序列运算理论,简化了卷积运算,但理论架构的特殊性制约了其大规模地推广应用。文献[13]中的无迹变换可直接通过求解输出变量的均值和方差得到其概率分布,避开了复杂的非线性变换。然而由于该方法以高斯分布为变换基础,具有一定的局限性。因而近年来,一些学者又开始以MCS为基础,试图通过优化抽样或简化模型等手段来解决耗时的难题,比如文献[14]采用拉丁超立方采样技术,显著地提高了采样效率。文献[15]提出扩展准蒙特卡洛方法,弥补了拉丁超立方采样不能保证序列低差异性的缺陷。而线性MCS的引入可以避免反复的潮流迭代。MCS的计算精度除了受采样个数的限制,还与输入变量的模型密切相关。目前,大多数文献均假设各输入量服从正态分布,局限性较大。文献[16]用双峰分布和双Weibull分布分别拟合负荷及风速模型,适用于一些特殊场合,但缺乏通用性。文献[17]建立了负荷的K均值模型,用聚类的方法处理不同时段的负荷数据,模型隐含地计入了节点负荷间的相关性,但聚类过程及显著性水平检验需要大量时间。

本文提出一种基于混合高斯模型(GMM)的改进MCS随机潮流计算方法,对于有测量数据的输入变量,通过GMM直接建立其精确的概率模型,而在缺失大量数据的情况下,根据输入变量特点选取常用分布进行拟合。在MCS基础上引入均匀设计抽样(UDS)得到计及相关性的输入变量样本,并利用多重线性化潮流方程保存精度、简化计算。

1 输入变量的随机模型

负荷的多样性及分布式电源的接入使得输入变量(节点注入功率)的概率密度函数(probability density function,PDF)难以用单一的标准分布准确拟合(如图1所示)。尤其在中低压网络,节点负荷的聚集程度下降,采用正态模型会有较大偏差。

GMM理论上能够平滑任何类型的输入变量分布,利用测量数据和期望最大化(expectation maximization,EM)算法可以方便地确定参数。目前,GMM作为研究热点广泛应用于机器学习、图像分割、模式识别等领域,本文主要用于模型的建立。

1.1 GMM

GMM是由若干个高斯(正态)分布线性叠加而成,对于一维随机变量X,其PDF定义为[18]:

式中:为第i个成分的分布;ωi,μi和σi分别为混合高斯函数第i个成分的比例、均值和标准差。其中,比例参数ωi必须满足归一性条件,即。

GMM所含的高斯分布个数n直接影响模型的精确程度,n越大,模型的偏差越小,但待求参数和耗时增大。图1给出了用GMM(n=4)近似实际的负荷样本分布。

1.2 参数估计与EM算法

GMM建模的关键在于高斯成分个数n及每个成分所含参数ωi,μi和σi的确定。为简化分析,假设n为定值,一般参数估计问题可直接利用最大似然函数法,如GMM的最大对数似然函数为[19]:

式中:NL为数据点xi(i=1,2,…,NL)的个数。

要找到一组合适的参数ωk,μk和σk(k=1,2,…,n),使式(3)取最大。然而由于数据点xi所属成分类别未知,无法对式(3)直接求导解方程来求得最大值,因此这里引入EM算法[20],其过程主要分为以下两步。

步骤1(期望):估计数据由每个高斯成分生成的概率,即每个成分被选中后的期望。对于数据xi,它由第k个成分生成的概率为:

步骤2(最大化):假设式(4)求出的γ(i,k)就是正确的“数据xi由成分k生成的概率”,根据式(5)至式(7)用最大似然法得到对应的模型参数。

式中:。

重复迭代上面两步,直到似然函数式(3)的值收敛为止。为了保证取到全局最优解,可先用K均值聚类法[21]得到粗略结果作为参数初值,再用EM算法进行细致迭代。

高斯成分个数n的确定可应用赤池信息量准则(AIC)[22],它是衡量模型拟合优良性的一种标准,可表示为:

式中:p为参数个数;L为似然函数。ηAIC的值越小,表明所取的模型越好,它可以权衡模型复杂度与数据拟合准确性之间的矛盾。

综上,GMM所有参数求解的过程包括:①设高斯成分个数n取2或3;②根据随机变量X的测量数据,采用EM算法得到ωi,μi和σi(i=1,2,…,n)的值;③由式(8)求得GMM的AIC值η(n)AIC;④返回②,n=n+1;⑤重复以上过程直到|η(n+1)AIC-η(n)AIC|<ξ。对于无测量数据或测量数据存在缺失的随机变量模型的建立方法,可参见附录A。

2 具有相关性的输入变量样本的产生

2.1 UDS技术

UDS是由方开泰和王元教授于1978年提出的一种新的抽样方法[23],它主要基于空间填充的思想,考虑试验点在试验范围内均匀散布。相比SRS,其优势有:①覆盖同样大小的样本空间,UDS的采样规模更小;②UDS的稳健性更好。设有m个服从[0,1]区间均匀分布的随机变量Y1,Y2,…,Ym,N为样本数,则利用UDS产生样本矩阵Ym×N的步骤如下。

步骤1:设Λ为正整数集合的一个无穷子集,n∈Λ。生成一个正整数向量H1×m=[h1,h2,…,hm],其中h1=1,1<hj<n且对任意i≠j,hi≠hj。

步骤2:从多项分布中选取m个独立同分布样本,构成向量B1×m=[b1,b2,…,bm]。

步骤3:通过SRS生成样本矩阵Em×N=[e1,e2,…,ej,…,em]T(ej=[ej1,ej2,…,ejN]),其中,随机变量e服从[-0.5,0.5]上的均匀分布。

步骤4:样本矩阵Ym×N=[y1,y2,…,yi,…,ym]T(yi=[yi1,yi2,…,yiN])中的元素yij为:

式中:{·}表示保留小数部分;i=1,2,…,m;j=1,2,…,N。

以[0,1]上均匀分布的两个随机变量y1和y2为例,令N=5,分别用UDS和SRS产生y1和y2的样本各三次,样本点分布如图2所示。可以明显看出,UDS产生的样本点在空间散布的更均匀,当采样点数量相同时,UDS比SRS覆盖的样本空间更大,采样效率更高。

2.2 边缘变换和相关性处理

产生服从给定分布的输入变量样本是MCS随机潮流计算的基础,这里借助边缘变换:设输入变量w的累积分布函数(cumulative distribution function,CDF)记作F(w),可以证明[24],u=F(w)服从[0,1]区间上的均匀分布。对于变量y∶U(0,1),其CDF为v=(y),也是[0,1]上的均匀分布。因此,输入变量w可由式(10)的变换求得,即边缘变换。附录B图B1给出了边缘变换的图形解释,它能实现随机变量从均匀分布到任意分布的变换。

式中:Φ为标准正态分布的CDF。

输入变量w1和w2(有测量数据)之间的相关性可用相关系数表示[25]:

式中:分别为w1和w2的样本均值。

而对无测量数据的输入变量,其相关系数则作为输入参数直接给定。这样便形成一个n1×n1阶的相关系数矩阵Σ=(ρij),n1为输入变量个数。

为了使由边缘变换得到的输入变量样本满足式(11)的相关性,引入中间变量v′,v′服从高斯分布,其作用是将变量的相关性变换转移到高斯域内进行。

由UDS产生的样本矩阵Yn×N=[y1,y2,…,yi,…yn]T,其中yi~U(0,1)。先由式(12)得到对应的高斯变量vi:

由于随机变量之间的相关系数在边缘变换过程中会发生改变,需要对相关系数矩阵Σ做相应修正,对于从高斯分布到均匀分布的变换,Σ的非对角元素应修正为[26]:

修正后的相关系数矩阵记为Σ*=(ρ*ij)。通过对Σ*进行Cholesky分解可以得到具有相关性的高斯样本矩阵Vcorr[27]:

式中:V=[v1,v2,…,vi,…,vn]T(vi=[vi1,vi2,…,viN]);L为上三角阵,Σ*=LLT。

将Vcorr中的样本代入式(15)(式(12)的逆向过程),得到具有相关性的均匀分布样本:

最后再做一次边缘变换,即可产生具有相关性的输入变量wi的样本:

式中:为输入变量wi的CDF。

对于GMM拟合的输入变量,上述过程可以简化。设GMM第k个高斯成分的参数为(ωk,μk,σk),ωk可理解为第k个成分被选中的概率。在由式(14)得到具有相关性的标准高斯样本vicorr=[vi1corr,vi2corr,…,viNcorr]后,根据比例系数ω大小,随机选取每个样本点所对应的高斯成分,假设样本点vijcorr(j=1,2,…,N)对应第k个高斯成分,直接应用式(17)可得输入变量wi的样本wicorr=[wi1corr,wi2corr,…,wiNcorr],从而避免了反函数求解的困难。

实际样本的相关系数矩阵Σ′与Σ之间存在一定偏差,这是由于边缘变换是非线性的,变换过程中会产生截断误差。文献[28]给出了特定分布下系数的修正方法,通常情况下偏差都很小,可忽略不计。

3 多重线性蒙特卡洛仿真

由于潮流方程的非线性,在对每个输入变量的样本点做确定性潮流计算时,需要反复迭代,这占了MCS超过95%的耗时。文献[29]采用交流线性化模型,将潮流方程在基准运行点处进行泰勒展开并忽略2阶以上的高次项,可得:

式中:W为节点注入功率;X和Z分别为状态变量和支路潮流;X0,Z0,W0分别为基准运行点处的状态变量、支路潮流和节点注入功率;S0和T0为灵敏度矩阵,S0=J0-1,T0=G0J0-1,其中J0为雅可比矩阵,。

但当输入变量的变化范围较大时,采用单点线性化模型会引起较大的截断误差,为了在提高MCS计算速度的同时不至于损失太多精度,下面引入一种多重线性化的方法[30]。

设Ptot为系统总的有功功率,有

式中:NN为PQ节点个数;Pi为节点i的负荷有功功率;NG为PV节点个数;PGj为接入节点j的发电机输出的有功功率。

由于负荷功率和新能源出力均为随机变量,可知Ptot也是随机变量,其PDF根据式(20)得到。附录C图C1大致给出了Ptot的PDF曲线,将其等间距地划分为RT个区域,RT一般取6~8。在区域R1至RT内选取相应的基准运行点,并将潮流方程在各点处分别线性化,得到

式中:Si为第i个区域对应的灵敏度矩阵;Xi0为第i个区域的状态向量:f(·)为节点功率方程;为第i个区域的基准功率向量。

的求解步骤如下:①由第1节和第2节的建模和抽样方法产生具有相关性的输入随机变量样本S;②根据式(19),计算各样本的Ptot并判断其所对应的区域;③将属于同一区域的样本取平均值,即得到各个区域的基准功率向量。

综合以上分析,采用本文所提基于GMM的多重线性MCS方法求解电力系统随机潮流的流程如图3所示。

4 算例分析

依托IEEE 30小型节点系统对所提方法的性能进行评估,算例的结构及线路参数见文献[31]。在节点25,29处分别接入两个10 MW的小型风电场(无监测数据),假设风速均服从尺度参数为8.09、形状参数为2.17的Weibull分布,且风电机组采用恒功率因数控制方式,输出特性可表示为文献[32]中的分段函数。

节点14处还接有一容量为2 MW的太阳能光伏电站,其无功出力为0,输出有功近似满足附录A式(A2)的β分布,其中形状参数α=β取0.9,光伏系统在最强光照条件下的输出功率PTmax为1.5 MW。

表1给出了用GMM拟合的负荷模型参数,一般高斯成分个数n取2或3即可达到AIC的收敛要求。表中只列出了有功功率的参数,假定负荷无功功率与有功功率之间的功率因数角恒定,有

式中:μP和σP分别为有功功率高斯成分的均值和标准差;φ为功率因数角。

除表1中节点4,7,18,21外,其余节点的负荷功率均服从高斯分布,均值即为文献[31]的负荷数据,标准差取均值的10%。考虑节点负荷之间的相关性,可将系统大致分为两个区域,区域1包含节点16至20,相关系数ρ1=0.6,区域2(节点15、节点23至26)的相关系数ρ2=0.8,两个区域间存在弱相关性(ρ1,2=0.2)。

采用本文方法对上述算例进行仿真计算,N=1 000。为了比较负荷的GMM及高斯模型对随机潮流结果的影响,利用式(23)和式(24),将表1中的节点负荷模型替换为均值和方差与原GMM相等的高斯分布:

高斯模型和GMM两种模型下节点26电压幅值及支路19-20无功功率的PDF如图4所示。

由图4可以看出,虽然两者的波动范围基本一致,但图线之间存在一定偏差,卡方检验[33]可用于拟合优度的量化分析,卡方值越小说明函数的拟合效果越好,利用MATLAB中的chi2gof函数得到GMM的卡方值分别为128.63和103.06,小于相应的高斯分布的卡方值478.07和584.52。进一步求出越限指标,则GMM对应的电压、无功功率的越限概率分别为0.76%和1.42%,而高斯模型的结果为0.42%和2.62%,从而导致低估或高估系统运行的风险。而且建模产生的误差会代入后面的潮流计算中,在误差传递作用下放大(或缩小),具体情况视方法本身而定。

输入随机变量样本的选取会直接影响MCS的计算精度。本文2.1节介绍了UDS技术及采样特点,为了评估其有效性,假设以N=20 000的MCS所得随机潮流结果为准确值,用μa和σa分别表示输出变量的准确均值和均方差,μs和σs分别表示仿真得到的输出变量的均值和均方差,则输出变量均值和均方差的相对误差计算公式可表示为:

式中:下标r表示输出变量的类型(包括电压幅值、相角,支路有功功率、无功功率)。

图5比较了不同采样规模下UDS和SRS的平均电压相对误差曲线。图中:εμmean,U为电压均值的平均相对误差;εσmean,U为电压均方差的平均相对误差。可以看出,样本容量相同时,无论是精度还是误差稳定性,UDS均优于传统的SRS,并且采样规模越小,优势越明显。

另外,对潮流方程的简化处理也会不可避免地引入截断误差。以节点25和支路10-22为例,三种潮流模型(线性、非线性、多重线性)得到的电压幅值和支路有功的CDF如图6所示(图中“非线性”即为常规MCS)。

由图6可知,采用本文的多重线性潮流模型可以获得较高的精度,其图线与简化前非线性潮流方程的结果几乎完全吻合,而单点线性模型的曲线在两端附近的误差较大,这是由于截断误差与输入随机变量样本点和基准点的偏离程度成正相关。

引入方差和的根均值(ARMS)来衡量输出变量概率分布的计算精度[34]:

式中:ξr为ARMS指标;Cr,P,i和Cr,M,i分别为所提方法和常规MCS得到的输出变量CDF上第i个点的值。

表2列出了系统部分输出变量在不同潮流模型下的ARMS值。ARMS值越大,说明所提方法与常规MCS所得输出变量的概率分布偏差越大。可以看出,多重线性化方法所得ARMS值大多小于1%,相较于单点线性化的结果,精度上有了质的提高。

在主频为2.63GHz、运行内存为2GB的intel i3计算机上,比较所提方法和MCS的计算耗时如表3所示。

当样本规模较小(N=50)时,两种方法的耗时主要由样本的生成时间组成,由于所提方法在输入变量模型和采样技术上较MCS复杂,计算时间略微长些。而N=500,1 000时,样本规模逐渐成为MCS耗时的决定因素,随着N的增加,MCS的计算时间呈线性增长,而所提方法耗时的增加并不明显且远小于相同采样规模下的MCS。

IEEE 118节点系统测试结果详见附录D。

5 结语

针对当前随机潮流方法无法同时兼顾计算精度和耗时、建模手段过于单一等不足,本文提出一种计及输入变量相关性的基于GMM的改进MCS随机潮流方法,该方法具有以下特点:①对于已知测量数据的输入变量,不管其分布类型及复杂程度,都能用GMM精确拟合,通用性较好;②利用UDS技术提高采样效率,结合边缘变换和Cholesky分解的方法可以方便地产生任意关联的输入变量样本;③建立多重线性潮流模型,在大幅提升计算速度的同时,最大限度地减小因潮流方程线性化引起的截断误差。

算例测试结果说明,在相同精度条件下,本文方法的计算效率相较于传统MCS有了显著提高,尤其适用于实际大型输配电系统在考虑新能源接入、负荷波动等不确定因素下的运行分析、安全评估及规划调度,具有较好的工程应用价值。

蒙特卡罗方法及应用篇5

关键词：三维公差综合,拟蒙特卡罗,旋量参数,广义逆矩阵,Halton序列

0 引言

公差设计不仅有助于降低产品生产成本,还影响产品质量以及装配成功率[1]。公差设计所涉及的因素较多,如何合理分配公差即公差综合是一个复杂的问题。目前有关公差综合的研究主要针对线性或二维尺寸公差的优化设计,对于计算机辅助三维公差综合的研究还不多[2,3]。

在三维CAD系统中,确定零件配合特征与装配功能要求的关系并进行公差综合方法的研究,对于辅助设计人员进行公差设计、提高产品质量有着重要意义。近些年来,一些研究人员在三维公差表示的基础上,把蒙特卡罗方法用于三维公差综合。文献[4]针对产品设计阶段预定义的装配功能要求,利用Jacobian矩阵建立与装配功能要求相关的3D公差链,并以此为基础采用蒙特卡罗法进行三维公差综合仿真研究。王太勇等[5]探讨了采用蒙特卡罗法进行三维尺寸及公差设计问题的实现方法。在特征自由度三维公差表示模型的基础上,许本胜等[6]采用蒙特卡罗法进行了三维公差综合仿真研究。然而,在进行三维公差综合仿真时,以上研究侧重于阐述蒙特卡罗法进行计算机仿真的流程,在具体的应用中简单地采用伪随机数进行仿真计算,对于蒙特卡罗方法中随机数的产生机制考虑不足,导致在大样本重复运算时会出现结果不稳定、收敛速度慢的问题,不能很好地发挥蒙特卡罗方法的应用效果。

拟蒙特卡罗方法是一种采用拟随机序列的蒙特卡罗方法,相比采用伪随机数的蒙特卡罗方法,拟蒙特卡罗方法不但具有更快的收敛速度,在仿真分析中所获得的结果也更为稳定和可靠。本文针对计算机辅助三维公差综合中仿真计算收敛速度慢以及结果不稳定的问题,在获得三维公差综合仿真等式的基础上,研究基于拟蒙特卡罗法的三维公差综合仿真方法,并通过实例进行具体的分析和阐述。

1 拟蒙特卡罗方法

拟蒙特卡罗方法的基本思想为:首先建立与问题相关的概率模型或随机过程,然后通过对概率模型或随机过程的观察或抽样试验来计算所求随机参数的统计特征,最后给出所求解的近似值。与蒙特卡罗方法不同的是,在对概率模型或随机过程进行随机抽样的过程中,拟蒙特卡罗方法采用的是拟随机序列。相较伪随机序列,拟随机序列的点列分布更加均匀,在仿真模拟时可以加快收敛速度并提高计算结果的稳定性和可靠性。

图1、图2所示分别为二维伪随机序列和二维拟随机序列。

图1、图2中的二维随机序列横纵坐标分别为区间[0,1]上的随机数。图2中的拟随机序列为二维Halton序列,与图1中的伪随机序列相比,它在函数空间的分布更加均匀。随机序列对应的点列在函数域上分布越均匀,其偏差相应越小,仿真计算时的收敛速度就越快,计算准确度也就越高。拟蒙特卡罗方法的计算误差可用Koksma-Hlawka不等式来确定[7]:

$\begin{array}{l} | \frac{1}{Ν} \sum_{k = 1}^{Ν} f (x_{k}) - \int_{Ι^{d}} f (u) d u | \leq \\ V (f) D_{Ν}^{*} (x_{1}, x_{2}, \dots, x_{Ν}) \end{array}$

其中,Id表示d维空间[0,1),随机点列x1,x2,…,xN∈Id;d、N分别为序列元素的维数和个数;V(f)是函数f在空间中Hardy-Krause意义下的有界变差函数;D*N(x1,x2,…,xN)定义如下:

D*N(x1,x2,…,xN)=DN(B;

$\begin{array}{l} Ρ) = \\ \sup_{B \in Ι^{d}} | \frac{A (B; Ρ)}{Ν} - λ_{d} (B) | \end{array}$

A(B; $Ρ) = \sum_{n = 1}^{Ν} χ_{B} (x_{n})$

其中,P是空间Id中的包含点集x1,x2,…,xN;对于任意属于空间Id的子集B,λd是d维勒贝格测度;χB是B的特征函数,A(B;P)是计算满足xn∈B(1≤n≤N)的点的个数的函数。

综上可得前N阶拟蒙特卡罗序列的误差收敛阶为O(N-1(lgN)d),而蒙特卡罗序列的误差收敛阶为O(N-1/2),因而拟随机序列的收敛速度要高于伪随机序列的收敛速度[7]。

2 基于拟蒙特卡罗的三维公差综合仿真

2.1 旋量参数和旋量矩阵

传统的尺寸与极限配合的公差表示方法已逐渐向新的产品几何技术规范(geometrical product specifications,GPS)过渡,出现了一些新的公差表示方法,如特征自由度[8]、旋量参数等[9,10]。这些公差表示方法本质上有着共同之处,即在三维CAD系统中,零件可看作是由点、线、面等基本几何特征构成的,将公差看作是构成零件几何特征微观上的三维几何变动量。为了规范和方便对公差问题的分析和描述,ISO/TC213延伸了特征的概念,定义了名义特征、名义派生特征、实际特征等特征类型,如图3所示。

从图3可看出,测量得到的实际特征为拟合特征或拟合派生特征,由于存在加工误差,其与名义特征或名义派生特征之间不可避免存在偏差,因而,从测量角度看,名义特征或名义派生特征在允许偏差范围(即公差)内是变动的,变动量可通过旋量参数来表示。对于选定的参考坐标系OiXYZ,旋量参数包括:名义特征或名义派生特征以原点Oi为中心的沿X、Y和Z轴平动方向的变动ui、vi、wi和绕X、Y和Z轴转动方向的变动αi、βi和γi。三维公差为名义特征或名义派生特征在三维空间允许的变动范围,相应地可通过旋量参数所构成的旋量矩阵来表示,记为Ti=[Diθi],其中Di=[uiviwi],θi=[αiβiγi]。

需说明的是,如果名义特征或名义派生特征沿某个方向的变动不产生新的特征,则对应的旋量参数记为零。如图3中,圆柱轴线名义派生特征沿轴线方向平动和绕轴线转动不产生新的特征,因而当选定轴线为参考坐标系OiXYZ的Z轴时,其三维公差对应的旋量矩阵为Ti=[Diθi],其中Di=[uivi 0],θi=[αiβi 0]。

2.2 三维公差综合仿真等式

产品是由若干零件装配而成的,对于规定的装配功能要求,与之相关的各零件配合面形成闭合的装配特征链,如图4所示。

图4中,产品的装配功能要求可视为装配特征链中首末配合特征之间的相对变动量,记为T0=[D0θ0],其中D0=[u0v0w0],θ0=[α0β0γ0]。装配特征链中第i个配合特征的变动量Ti(i=1,2,…,n,n为装配特征链配合特征个数,且n>1)对装配功能要求的累积影响量记为T0_i=[D0_iθ0_i],D0_i=[u0_iv0_iw0_i],θ0_i=[α0_iβ0_iγ0_i]。T0、Ti间的关系可以利用机器人学中刚体的坐标转换方法建立,为了减少坐标转换的计算工作量,将参考坐标系与装配基准坐标系各坐标轴取为同向。原点Oi在装配基准坐标系中的位置记为(Xi, Yi, Zi),利用机器人学中刚体坐标变换的知识,有

TT0_i=J0_iTTi (1)

$J_{0_i} = [\begin{matrix} 1 & 0 & 0 & 0 & Ζ_{i} & - Y_{i} \\ 0 & 1 & 0 & - Ζ_{i} & 0 & X_{i} \\ 0 & 0 & 1 & Y_{i} & - X_{i} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}] (2)$

从而:

$Τ_{0}^{Τ} = \sum_{i = 1}^{n} Τ_{0_i}^{Τ} = J [Τ_{1} Τ_{2} \dots Τ_{n}]^{Τ} (3)$

系数矩阵J称为雅可比矩阵,有

J=[J0_1J0_2 … J0_i … J0_n] (4)

由于装配特征链中配合特征数目n>1,因而J为6×6n阶长方阵。式(3)中,将Ti(i=1,2,…,n)视为未知量,对于给定的T0,式(3)有解但解不唯一,为获得唯一解,取解的最小范数作为约束条件,其现实意义为:各零件配合特征理想公差视为零,获得解与理想值之间具有范数意义下的最小偏差。记x=[T1T2 … Tn]T,y=TT0,由广义逆矩阵的知识,存在唯一的矩阵J+,满足:

$\min_{J x = y} ∥ x ∥ = ∥ J^{+} y ∥$

矩阵J+称为长方阵J的广义逆矩阵,J+可通过下式进行计算[11]:

J+=JT(JJT)-1 (5)

综合式(1)～式(5)可得三维公差综合模型如下:

[T1T2 … Tn]T=J+TT0 (6)

式(6)即为三维公差综合仿真等式。

2.3 三维公差综合仿真

拟蒙特卡罗方法的关键在于拟随机序列的生成,目前已经提出了若干拟随机序列,如Faure序列、Sobol序列以及Halton序列等。本文采用Halton序列进行仿真计算。记s维Halton序列为X=(x1,x2,…,xs),通过选定某个十进制整数m和以下步骤来生成s维Halton序列。

(1)选择s个基(base)b1,b2,…,bs,这s个基可以使用前s个素数(2,3,5,7,…)。

(2)对每个基bi,将m转换成bi进制:

$m = \sum_{k = 1}^{n_{i}} a_{k} b_{i}^{k - 1} i = 1, 2, \dots, s (7)$

其中,ni为将m转换成bi进制后的位数,ak为将m转换成bi进制后第k位上的数字。转换后得到一个随机向量(m(b1),m(b2),…,m(bi),…,m(bs)),其中m(bi)为m对应的bi进制整数。

(3)将m(bi)按位反序排列(如120反序排列为021),然后在前面加小数点得到新的值并将其转换成十进制数,记为xi,有

$x_{i} = \sum_{k = 1}^{n_{i}} a_{k} b_{i}^{- k} i = 1, 2, \dots, s (8)$

式中,ak为小数点后的第k位数字。

(4)置m←m+1,重复步骤(2)、步骤(3)直至得到所需的随机向量个数。

由上述步骤生成Halton序列,进行适当数学变换可得到服从其他分布的随机序列,如正态分布转换方法为:将s维Halton序列中各随机变量xi相加,依据中心极限定理,当s充分大时,有

$\frac{1}{\sqrt{n / 12}} (\sum_{i = 1}^{s} x_{i} - \frac{n}{2}) ~ Ν (0, 1) (9)$

通过式(9)获取标准正态分布随机变量,且当s=12时可达到较好的精度,从而要得到服从正态分布的随机变量x～N(μ,σ2),作如下变换即可:

$x = (\sum_{i = 1}^{12} x_{i} - 6) σ + μ (10)$

综合上述三维公差综合模型并利用Halton序列生成特定分布随机数的方法,基于拟蒙特卡罗的模拟仿真思路为:首先,根据特定分布函数,通过Halton序列获取装配功能要求T0各旋量参数的随机数;然后,由式(6)计算得到装配特征链中各配合特征旋量参数的值,重复计算便得到各配合特征旋量参数的样本;最后,对所得样本进行统计分析,得到样本均值、样本方差等统计量。图5所示为基于拟蒙特卡罗的三维公差综合仿真流程。

3 实例分析

图6a所示齿轮泵齿轮的精度等级为5级,齿轮分度齿圈径向跳动公差为Fr=0.016mm,综合考虑补偿热变形、保证正常润滑等因素,要求啮合轮齿侧向间隙Fc最小值为0.03mm,最大值为0.1mm。

3.1 装配功能要求与配合特征旋量参数

轮齿啮合的接触形式为线接触(图6b),装配功能要求(functional requirement of assembly,FRA)即侧向间隙可看作是两轮齿名义啮合线A1、C1空间位置发生了相对变动。装配链中首末配合特征旋量矩阵T0=[0 v0w0α0β0γ0],各旋量参数与Fr、Fc的关系为

${\begin{cases} v_{0} = F_{c} 0.03 m m \leq v_{0} \leq 0.1 m m \\ w_{0} = F_{r} / 2 - 0.016 / 2 m m \leq w_{0} \leq 0.016 / 2 m m \\ α_{0} = F_{c} / 2 r 0.03 / r \leq α_{0} \leq 0.1 / r \\ β_{0} = F_{r} / 2 b - 0.008 / b \leq β_{0} \leq 0.008 / b \\ γ_{0} = F_{r} / b - 0.03 / b \leq γ_{0} \leq 0.1 / b \end{cases} (11)$

式中,r为齿轮分度圆半径;b为齿轮宽度。

图7为与装配功能要求相关的各零件配合特征及三维公差对应的旋量参数。

3.2 三维公差综合仿真

图7中,{Oi}(i=1,2,…,6)为装配特征链中各配合特征参考坐标系,{Oi}与装配基准坐标系{O0}各坐标轴方向相同,各原点坐标为:O1=O6=(0,0,0),O2=(-b,0,r),O3=(-b,0,r),O4=(-b,0,-r),O5=(-b,0,-r)(b=19mm,r=12.5mm)。利用式(2)、式(4)得到系数矩阵J为

$J = [\begin{matrix} 0 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 0 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 \\ 0 & 1 & 0 & - r & 0 & 0 & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & - r & 0 & 0 & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & 0 & 0 & - b \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & b & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}]$

由式(6),利用MATLAB软件计算得到齿轮泵各配合特征三维公差域内旋量参数与侧向间隙各变动量之间的关系。以分度圆特征C1[0 v6w6α6β6γ6]为例,通过计算得到以下等式:

$\begin{array}{l} v_{6} = \frac{r γ_{0} + v_{0}}{2 (11 + r^{2})} \\ w_{6} = \frac{b β_{0} + w_{0}}{2 (11 + b^{2})} \end{array}$

$\begin{array}{l} α_{6} = \frac{α_{0}}{2} \\ β_{6} = \frac{(11 + 2 b^{2}) β_{0} - 2 b w_{0}}{4 (11 + b^{2})} \\ γ_{6} = \frac{(11 + 2 r^{2}) γ_{0} + 2 r v_{0}}{4 (11 + r^{2})} \end{array}} (12)$

在进行三维公差综合仿真时,首先需对装配功能要求T0进行模拟采样。在批量生产的情况下,尤其当装配特征链中配合特征数目n较大时,可认为T0服从正态分布。这里假定Fc、Fr服从正态分布。对于Fc～N(μ,σ2),根据设计要求,取侧向间隙均值μ=(0.03mm+0.1mm)/2=0.065mm;方差σ用于描述实际结果与设计目标相偏离程度,可由6σ设计原则来确定[12],取σ=0.0117mm。同样可得到Fr～N(0,0.00272)。由式(10),利用12维Halton序列分别获取Fc、Fr的随机数,进而通过式(11)即可得到T0中各旋量参数v0、w0、α0、β0及γ0的随机数。

以旋量参数γ6为例,分别采用蒙特卡罗和拟蒙特卡罗方法,利用获得的旋量参数v0、w0、α0、β0及γ0的样本通过式(12)进行仿真计算,在不同容量样本和相同容量样本情况下进行两次实验,得到的仿真统计结果见表1和表2。

表1中,γ6(p)为采用蒙特卡罗方法计算的样本均值,γ6(q)为采用拟蒙特卡罗方法计算的样本均值。由表1可看出,在不同容量样本情况下进行多次仿真计算,采用拟蒙特卡罗方法得到的标准差S较小,说明采用拟蒙特卡罗方法获得的结果波动较小。另外,由表1也可看出,当增大样本容量时拟蒙特卡罗方法的收敛速度更快,尤其当样本容量较大时,计算结果更为稳定。

由表2可看出,在相同容量样本情况下进行多次仿真计算,采用拟蒙特卡罗方法得到的结果的波动较小,表明采用拟蒙特卡罗方法得到的结果更为稳定。

由上述分析可知,拟蒙特卡罗方法更为高效,获得的结果也更稳定和可靠,因而能更好地运用于三维公差综合模拟仿真。取样本容量M=5000,通过拟蒙特卡罗法进行公差综合仿真计算,以分度圆特征C1为例,其结果见表3。

表3中, $\bar{X}$ 一栏为几何特征C1各旋量参数的三维公差分配结果,S为样本标准差,Xmax、Xmin分别为最大值、最小值。各旋量参数最大值、最小值之差即为C1在参考坐标系{O1}中各方向的变动范围,差值越小,表明对相应方向上C1的几何变动量要求越严格;标准差对应样本分布集中程度,其值越小,样本分布越集中,相应几何变动量越难以加工保证,在设计阶段应予以重视。

4 结语

本文在利用旋量参数和旋量矩阵建立三维公差综合仿真等式的基础上,详细讨论了用拟蒙特卡罗法进行三维公差综合仿真的方法。通过将所阐述的方法应用于齿轮泵装配体的三维公差综合,得到泵体、齿轮轴及轮齿配合特征旋量参数统计量。由于三维公差设计问题的复杂性,如何将仿真计算得到的配合特征各旋量转化成具体的公差大小还有待进一步的研究。通过对实例结果的分析可以看出,相比蒙特卡罗方法,采用拟蒙特卡罗方法进行三维公差综合仿真时的收敛速度更快,获得的结果更为稳定和可靠,可以简单高效地应用于三维公差综合的仿真计算,从而验证了本文所阐述方法的有效性和实用性。

参考文献

[1]Hong Y S,Chang T C.A Comprehensive Review ofTolerancing Research[J].International Journal ofProduction Research,2002,40(11):2425-2459.

[2]王恒,宁汝新.计算机辅助公差设计与分析的研究现状及展望[J].航空制造技术,2006(3):73-75,102.Wang Heng,Ning Ruxin.Present Status and Prospect ofComputer Aided Tolerance Design and Analysis[J].Aero-nautical Manufacturing Technology,2006(3):73-75,102.

[3]孙晓燕,黄克正,张勇.计算机辅助三维公差设计的研究现状[J].工具技术,2007(41):15-19.Sun Xiaoyan,Huang Kezheng,Zhang Yong.PresentState and Review of Computer Aided 3D ToleranceDesign[J].Tool Engineering,2007(41):15-19.

[4]LaperrieáL,Kabore T.Monte Carlo Simulation of Toler-ance Synthesis Equations[J].International Journal ofProduction Research,2001,39(11):2395-2406.

[5]王太勇,熊越东,路世忠,等.蒙特卡罗仿真法在尺寸及公差设计中的应用[J].农业机械学报,2005,36(5):101-104.Wang Taiyong,Xiong Yuedong,Lu Shizhong,et al.Application of Monte Carlo Method to Dimensionand Tolerance Design[J].Transactions of the Chi-nese Society for Agricultural Machinery,2005,36(5):101-104.

[6]许本胜,王灿,景晖,等.基于特征自由度变动的公差综合建模与仿真[J].中国机械工程,2011,22(24):2981-2985.Xu Bensheng,Wang Can,Jing Hui,et al.ToleranceSynthesis Modeling and Simulation Based on DOF ofGeometric Variations of Features[J].China Me-chanical Engineering,2011,22(24):2981-2985.

[7]雷桂媛.关于蒙特卡洛及拟蒙特卡洛方法的若干研究[D].杭州:浙江大学,2003.

[8]黄美发,钟艳如.CAD系统中并行公差的建模方法[J].中国机械工程,2004,15(18):1623-1625.Huang Meifa,Zhong Yanru.A New Model for Con-current Tolerancing in CAD System[J].China Me-chanical Engineering,2004,15(18):1623-1625.

[9]胡洁,熊光楞,吴昭同.基于变动几何约束网络的公差设计研究[J].机械工程学报,2003,39(5):20-26.Hu Jie,Xiong Guangleng.Study on Tolerance De-sign Based on Variational Geometric ConstrainsNetwork[J].Chinese Journal of Mechanical Engi-neering,2003,39(5):20-26.

[10]张为民,陈灿,李鹏忠,等.基于雅可比旋量法的实际工况公差建模[J].计算机集成制造系统,2011,17(1):77-83.Zhang Weimin,Chen Can,Li Pengzhong,et al.Toler-ance Modeling in Actual Working Condition Based onJacobian-Torsor Theory[J].Computer IntegratedManufacturing Systems,2011,17(1):77-83.

[11]姚俊,张玉春.工程矩阵方法[M].2版.北京:国防工业出版社,2007.

蒙特卡罗方法及应用篇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

蒙特卡罗方法, 是一种通过对随机变量的统计试验, 随机仿真求解数学、物理、工程技术问题近似的数学方法。通常作决策分析时, 蒙特卡罗可以随机仿真各种变量间的动态关系, 解决某些具有不确定性的复杂问题, 被公认为是一种经济而有效的方法。

它可以通过以下步骤来实现:

(1) 根据实际问题, 构造模拟的数学模型。

(2) 确定描述研究对象行为的概率分布, 构造仿真模型仿真的前提是可以用概率分布来对受不确定性影响的参数加以描述, 在蒙特卡罗仿真中, 生成大量的项目假想情况来反映实际项目的特征。

(3) 根据确定的概率分布进行随机抽样, 即进行数字仿真。

(4) 根据随机数字仿真结果和项目的要求, 统计研究对象事件发生的频数, 并运用数理统计知识求取各种统计量。

1 投标报价风险分析中变量基本类型概率分布的确定

在工程项目经济评价中, 由于不能给出随机变量变化的范围而不能反映出未来的变化特性, 因此通常只能借助主观估计的数据和信息来分析风险变量的变化特性。为提高精度, 需要引用概率分布来描述风险变量的变化规律。

分布是一种基于对样本数据分析的概率分布, 它的概率分布完全由样本的均值和偏差来予以确定, 分布的价值在于其曲线的多样性, 通过改变参数a和b来改变曲线的类型, 使得曲线的类型丰富多变, 它可以用于表示向任一方向偏斜, 无偏斜, 或偏斜十分严重的数据。这就为我们构造不同风险变量的概率分布模型, 提供了便利的方法。

分布的概率密度函数为:

式中, P (x) 为密度函数;a, b为分布参数;β (a, b) 为函数。

通过现有的蒙特卡罗计算程序, 我们由己知的工程造价资料如分部分项工程费等, 计算出a, b的值从而确定出该成本因素的β分布曲线。如图1所示。

2 仿真误差和仿真次数的确定

在蒙特卡罗仿真模拟时, 模拟次数越多越正确, 误差也越小, 但实际上模拟次数过多不仅费用高, 整理计算结果费时费力, 故模拟次数过多也无必要。但模拟次数过少, 随机数的分布就不均匀, 影响模拟结果的可靠性。

理论上我们可以用数理统计的知识对预测值的期望的精确度进行估计, 由此确定模拟次数。假设模拟的随机变量记为ξ1, ξ2, ……, ξn, 在给定模拟次数N和置信度p时, 在N足够大时, 采用中心极限定理估计法估计误差。依据中心极限定理估计法估计, 以ξ軃代替以E (ξ) 的误差。若是ξ1, ξ2, ……, ξn, 总体ξ的容量为N的子样, 则ξ1, ξ2, ……, ξn, 不仅相互独立, 且与总体同分布。当E (ξ) =μ<∞, D (ξ) =σ2>0, 则由中心极限定理可知, 当N充分大时, 统计量似地服从N (0, 1) 分布, 则对于给定的置信度β有:

即:

取, 得到, 此为在置信度为β的误差表达式。

我们从以上结果可知, 仿真的误差取决于模拟变量的标准差σ和模拟次数N, 由于对于确定的模拟变量来说σ是一个确定的值, 所以只有通过增加模拟次数N才能达到减少误差的目的。

在实际操作中, 投标造价分析人员不必麻烦地确定置信度和计算误差, 只要对模拟结果的曲线进行分析, 当模拟所得的曲线越光滑, 模拟的精度越高。在可以接受的运行的速度下, 通过观察曲线的光滑程度, 基本可以确定是否需要提高模拟次数再次进行模拟。

3 结语

蒙特卡罗法是投标报价风险分析工作中的一种有效方法, 但其实验和计算过程较为繁琐。由于计算机技术的发展, 现代的蒙特卡罗方法, 已经不必亲自动手做实验, 而是借助计算机的高速运转能力, 使得原本费时费力的实验过程, 变得简单快速, 这为蒙特卡罗方法在现代项目管理中应用提供了技术基础。

目前, 发达国家已经把蒙特卡罗模拟方法列入项目管理的常规方法。有关计算机应用软件也己经有许多种产品, 如C语言、VB、VC、Matlab等, 在造价风险分析中常运用Matlab软件进行蒙特卡罗模拟计算。计算机软件与蒙特卡罗方法的结合使得投标报价风险分析更为高效、准确, 为造价人员进行风险分析提供了一条科学有效的途径。

摘要：介绍一种确定预测价格概率分布的技术——蒙特卡罗方法。投标企业在投标报价阶段可以用它来评价某一预测价格, 并由此确定降价幅度、是否低于成本价等投标报价决策中急需解决的问题。

关键词：蒙特卡罗法,投标报价,风险分析

参考文献

[1]杨小俊.工程量清单招标模式下的成本分析与风险防范.建筑, 2004 (10)

[2]许高峰.工程项目投标的优化报价模型.系统工程理论与实践, 2004 (2)

[3]石道元, 常茜惠.蒙特卡罗模拟技术在投资风险决策中的应用.中国管理信息化, 2007 (10)

【蒙特卡罗方法及应用】推荐阅读：

蒙特卡罗模拟08-16