序贯蒙特卡洛模拟

序贯蒙特卡洛模拟（精选7篇）

序贯蒙特卡洛模拟 篇1

0 引言

风力发电的高比例接入对电力系统的安全经济运行产生巨大影响,风电所具有的时序性和相关性特征对电力系统随机生产模拟的要求更为严苛。随机生产模拟可以考虑各种不确定性因素,如电力负荷的随机波动、发电机组随机停运等情况,更深刻地描述了电力系统的生产过程[1]。目前电力系统随机生产模拟已有多种较为成熟的方法,主要分为解析法和模拟法两大类。解析法主要包括离散卷积积分法、等效电量法、累积量法等,在对模拟过程的时序性没有严苛要求的电力系统分析中因其计算迅速而占据主流地位。模拟法则主要以蒙特卡洛抽样法[2,3,4,5]为代表,进一步可细分为序贯、伪序贯和非序贯3种。序贯蒙特卡洛模拟法保留了系统运行的时序性,考虑了负荷变化与发电机组出力的时序变化特性,可以得到系统运行的详细有效信息,故更适合含风电场等可再生能源的电力系统的随机生产模拟。

含风电电力系统随机生产模拟目前存在的问题是:以多状态机组[6,7]为代表的模型无法严格考虑风电出力的时序性,而这种时序性是风电最显著、本质的特征,考虑准确与否对结果有至关重要的影响;以离散卷积为代表的模型在考虑机组运行动态约束方面有困难,所提供的生产模拟统计信息有限;在以蒙特卡洛为代表的模型中,对系统运行约束的细致模拟和其计算时间等方面有待改进。

针对以上问题,本文以电力系统接纳风电能力评估[8]为背景,提出基于序贯蒙特卡洛仿真的电力系统随机生产模拟滚动试探算法,讨论了负荷及风电出力的不确定性处理、生产费用计算、调峰评价指标等问题。对我国东北某省级电网实际系统进行了仿真分析,比较了风电并网前后对生产模拟指标及火电机组运行的影响,评估了强化系统约束、不同计算时间粒度对调峰不足系数、弃风率的影响,详细分析了全年及日平均各时段弃风差异及原因,验证了所提方法的合理性和实用性。

1 系统不确定因素模拟

1.1 负荷、风电不确定性处理

为准确地考虑负荷及风电的不确定性,使计算结果更加符合实际,引入随机偏差概念,即在负荷、风电原始值基础上叠加随机偏差。

因此,任意时刻t的实际负荷可表示如下:

其中,Lrt为t时刻的实际负荷;Lft为t时刻的原始负荷;δlt为负荷随机偏差,该偏差为服从均值为0、方差为σl2的正态分布的随机变量[9]。根据文献[10]的研究,δlt的标准差由下式计算:

其中,k一般取值为1。

风电实际出力由下式表示:

其中,Wrt为t时刻风电实际出力;Wft为t时刻风电出力期望值;δwt为t时刻风电随机偏差[11],服从如式(4)所示的方差的正态分布。

其中,WI为风电场总装机容量。

1.2 常规发电机组的运行状态模拟

常规发电机组用两状态模型描述。采用状态持续时间抽样法模拟每个元件状态转移和循环过程,认为正常运行持续时间τ1和故障修复时间τ2服从指数分布,计算公式如下所示。

其中,γ1、γ2为均匀分布的随机数;TTMTTF为平均工作时间;TTMTTR为平均修复时间。

2 电力系统随机生产模拟滚动试探模型

2.1 火电机组带负荷顺序

供热期:将供热机组分为2个容量段,第一分段为强迫出力,用必开可用供热机组的第一分段作为基荷修正负荷曲线;第二分段及其他火电机组按最小比耗量由小到大依次带负荷,在满足调峰充足的前提下保证经济性。当风电出力超过常规电源调节出力时,采取弃风措施,风电消纳受限[12],即调峰容量不足,记录该调峰不足容量值和时段;若某时刻所有可用发电机组全部满额出力仍不能满足负荷要求,即发电容量不足,记录该发电不足容量值和时段。

非供热期:除供热机组无需分段外,其他安排方式同上。

2.2 随机生产模拟过程

随机生产模拟过程的具体步骤如下。

步骤1根据风电出力波动选择计算时间粒度;系统火电机组按最小比耗量排序,根据元件随机故障/修复模型确定一年中元件的状态持续曲线;负荷、风电出力由1.1节确定。

步骤2用必开可用供热机组第一分段、联络线输送功率修正负荷,将每时段风电出力进行试探性全额消纳,即每时段首先将风电作为负的负荷修正负荷曲线,用变步长迭代算法安排可调节水电。

步骤3选择滚动的最小时间长度(文中设定为系统中最小的平均无故障时间(TMTTF),在一个滚动时长内假设所有机组正常运行,采用计及机组启停限制的多时段优先顺序法[13]进行机组组合计算。

步骤4根据序贯蒙特卡洛模拟的机组状态持续曲线查找此滚动时长内最早出现故障的机组和时刻,将该滚动时长内此时刻之前的时段的常规机组按前述规则带负荷。

步骤5以此时刻为下一个滚动时长起点,记录该时刻前一时刻各元件的运行状态及持续时间,作为下一滚动时长的初始状态。依此类推,进行下一个滚动时长的模拟,直至一年所有时段计算完毕。

步骤6统计年调峰不足系数、弃风率、动态成本、燃料成本、环境成本、可靠性等指标[14]。

步骤7重复上述模拟过程,直至满足收敛判据或达到给定模拟年数。

2.3 收敛条件

进行多年模拟直至某个收敛最慢的统计指标收敛,收敛判据如下式所示。

其中,X为某一收敛最慢的统计指标;σ(X)为X的标准差;E(X)为X的均值;N为模拟年数;ε为给定值,一般取0.03~0.08。

2.4 随机生产模拟流程

随机生产模拟流程如图1所示。流程中num取值0、1、2,根据一年时序的气候状况变化,num=0为供热期,num=1为非供热期,num=2为供热期。

3 评价指标计算

3.1 计及动态费用的生产成本分析

对于系统中的n台常规机组,按流程可以逐台计算一年的平均停机次数Sai,进而可以计算每台机组的动态费用。

3.1.1 动态费用计算

常规机组的动态费用Dc可按下式计算:

其中,Tci为第i台机组的停机费用;Qrci为第i台机组的热启动费用;Qlci为第i台机组的冷启动费用;Xi为第i台机组热启动次数,Yi为第i台机组冷启动次数,Xi+Yi=Sai。

假定常规机组启动方式为冷启动,则停机时间soff满足:

其中,Soff为设定的冷启动停机时间下限。

3.1.2 风力发电成本计算

风力发电成本包括建设投资回收费用及运行和维护费用。计算公式如下:

其中,Wc为风力发电总成本;I为投资回收成本;O为风力发电机组的运行费用;M为风力发电机组的维护费用;为年平均接纳功率。

3.1.3 总生产成本计算

全面计及燃料成本、环境成本、风电成本、动态发电成本的系统总生产成本为:

其中,Fc为燃料成本,Ec为环境成本。

其中,Pi,t为t时段第i台机组出力;ai、bi、ci为成本系数;T为一年内计算总时段数;m为环境污染物的种类;Ce为第e项环境污染物的单位环境价值,单位为元/(k W·h);PG为火电机组发电总量[15]。

风电的可避免费用为:

其中,Cf0为不含风电系统的总生产成本;Cf1为含风电系统的总生产成本。

3.2 调峰不足系数PPSCIC、电量不足期望EENS、电力不足概率PLOLP

一年模拟结束后,计算调峰不足系数PPSCIC:

其中,Nh为全年模拟调峰不足时段数。

其中,Ny为仿真总年数;Ct为t时段切负荷量;St取0或1,其中St=1表示出现失负荷。

4 算例分析

在基于序贯蒙特卡洛模拟法的电力系统随机生产模拟滚动试探模型基础上,通过编写Java程序实现模拟过程。计算环境为台式电脑,系统配置如下:CPU为Intel Pentium G2030双核,内存为3.43 G,操作系统为Windows 7(32 bit)。

4.1 算例系统

以中国东北某省级电网为例进行模拟评估。该省的常规电源总装机容量为17 446 MW,有52台抽气机组、13台凝气机组、2台水电机组。风电平均装机容量为4 793 MW,占总装机容量的21.6%。负荷最大值为11116 MW。每年4月15日到10月15日为非供热期,其余为供热期。供热机组供热期出力上限为额定出力乘以(1-Cv),出力下限为额定出力乘以Cb,Cv、Cb的取值详见表1,其他机组全年出力上下限不变。系统供热期所有可用供热机组全开。系统备用容量选取为网内最大一台火电机组容量,为600 MW,各机组爬坡速率按每分钟额定出力的1%计算。不同计算时间粒度下原始负荷和风电出力输入采用实际调度系统2014年的历史输出数据。

4.2 指标的收敛性分析

文中包含EENS、PPSCIC等统计指标。在评估流程中选择EENS统计指标满足式(6),即符合收敛性要求作为整个评估流程的结束标志。之所以选择EENS而未选择其他指标(以小时粒度为例),是因为经过多次算例测试发现EENS的收敛速度慢于其他指标。表2给出了算例接入风电时,分别模拟50 a、100 a、200 a的部分统计指标取值,以考察各种指标的收敛速度。图2给出了相应指标的收敛曲线。由图2可见,EENS指标收敛最慢。

4.3 场景分析

对以下4个场景进行比较分析(收敛判据取0.08)。

场景1:时间粒度为1 h,不加风电,考虑机组启停特性的系统。

场景2:时间粒度为1 h,加入风电实际出力,考虑机组启停特性的系统。

场景3:时间粒度为1 h,加入风电实际出力,不考虑机组启停特性的系统。

场景4:时间粒度为15 min,加入风电实际出力,考虑机组启停特性的系统。

4.3.1 不同场景下各项指标对比

(1)风电加入前后指标分析。

由场景1与场景2结果对比可以看出,风电加入前后对系统可靠性、成本及火电机组运行的影响。表3给出了2种场景的计算结果,得出如下结论。

a.风电场接入后,燃料成本减少20%,环境成本减少10%。

b.风电场接入后,动态成本增加了42%。

c.风电场的年平均可避免费用为5 200万元,风电并网后较并网前的总成本节省了2.7%。可见风电接入后,增加了动态成本,节省了燃料成本和环境成本,风电年平均可避免费用虽然不多,但减小总成本的同时,带来了可观的环境效益。

d.风电场接入后,平均每台火电机组利用小时数由风电并网前的3 765 h下降到3 502 h,降低了263 h。火电机组平均每小时运行台数由37台下降到34台,相当于风电平均每小时可替代3台火电机组运行。

(2)启停约束、不同计算时间粒度的影响分析。

表4给出了场景2、3、4的计算结果。

比较场景2和3发现:同一系统,启停约束加入后,调峰不足系数增加了10.49%,调峰情况更为严峻,且从图3可知供热期的平均调峰不足小时数大于非供热期,可见,供热期强迫出力加剧了调峰压力;弃风率增大了15%,机组开停机约束加强使得风电接纳减少,因此,在只考虑最大化风电接纳的情况下,一般不考虑机组启停约束。

比较场景2和4可得:时间粒度减小,弃风率减小、燃料成本降低、动态成本增加,可以更加细致真实地再现系统运行轨迹。虽计算时间有所延长,但能满足工程计算要求。因此,需要科学选择计算时间粒度,达到计算效率和时间的最佳结合。

4.3.2 场景2情况逐层分析

(1)年时间尺度弃风、开机台数分析。

表5给出了供热期与非供热期的风电接纳情况,可以看出供热期风电可用总量是非供热期的1.45倍,而风电利用小时数却少于非供热期,弃风率占到总弃风率的98.14%。这是由此系统的性质决定的。

图4给出了该省原始负荷减去联络线输送及供热强迫出力后的等效负荷以及风电可出力对比图,图5给出了该系统一年8 760 h风电可出力及风电接纳量数据。从图4和图5可以看出,该省供热期强迫出力大,使得供热期等效负荷小,间接缩小了供热期的风电出力空间,但该省风电出力在供热期相对充足,因此必然导致供热期弃风较多;而非供热期因风电接纳空间大、机组投入相对灵活且风力发电少,使得该时期弃风率较小,这与表4的分析结果相呼应。

图6给出了全年每小时各类型机组运行台数,发现供热期供热机组的开机台数变化不大,这是为了满足该省的供热需求,其运行台数的微小波动是由强迫停运造成的;非供热期2种类型机组的运行情况相对更加灵活。

(2)日时间尺度弃风分析。

图7、图8将统计年的指标纵向处理为全年日平均各时刻风电可出力值、等效负荷和弃风量,由图7可知风电高发期在22:00左右,而图8统计得到的等效负荷曲线变化趋势高峰期在10:00—12:00,第2个高峰期出现在16:00—18:00,第3个高峰出现在20:00—21:00;弃风量的高峰出现在23:00左右,其中夜间的弃风量明显大于白天,这与夜间负荷小且风电出力高有关。

5 结论

本文提出了综合考虑元件停运、启停限制、备用容量、爬坡速率等系统约束的基于序贯蒙特卡洛的电力系统随机生产模拟滚动试探算法,考虑了可再生能源风电的时序性和相关性,模拟过程细化且更加符合实际,具有灵活、简洁、提供的技术经济指标丰富等特点,在新能源电力系统技术经济评价中有独特的优势。根据需要,模拟时间粒度可以更小(如10 min或5 min),以期考虑运行中系统和元件更加细微的变化,提供更加详细的参考信息。

摘要：针对风电并入电网系统运行状况的变化,提出了基于序贯蒙特卡洛仿真的随机生产模拟滚动试探算法并应用于实际工程。该方法计及元件强迫停运、机组启停限制、运行经济性、系统供热期与非供热期、风电出力波动性、时序性等因素,将多时段优先顺序法融入计算体系,算法除提供传统经济性和可靠性指标外,新增了调峰不足系数、弃风率等调峰评价指标。对我国东北某省级电网实际数据进行计算分析,比较了风电并网前后对技术经济指标及火电机组运行的影响,重点评估了强化系统约束、不同计算时间粒度对调峰评价指标的影响,详细分析了每日各时段的全年平均的弃风差异及原因,验证了所提方法的合理性和工程应用价值。

关键词：序贯蒙特卡洛,随机生产模拟,风电,接纳能力,多时段优先顺序法,调峰不足系数,弃风分析

序贯蒙特卡洛模拟 篇2

关键词：风力发电,调峰裕度,序贯蒙特卡洛方法,调峰不足概率

0 引言

进入21世纪以来,风电蓬勃发展,装机容量一直保持高速增长,但风电的随机间歇性及预测精度低等特性,使风电并网面临诸多挑战。其中,风电并网使系统调峰问题凸显[1,2],各地均出现了不同程度的风电限电现象[3]。目前,国内外学者针对风电并网系统调峰问题作了大量研究。文献[4]考虑了风电预测对调峰的影响,设计了大规模风电运行的调度模式。文献[5]以调峰为控制目标,建立了风电场3层有功控制中的能量管理层控制模型。文献[6]结合吉林电网的实际情况,提出了计及调峰因素的日前发电计划协调优化方法。文献[7]基于确定性方法,从系统调峰角度评估了电网接纳风电的能力。文献[8-9]基于确定性方法,分别针对西北电网和京津唐电网,分析了风电并网后系统调峰特性和调峰能力计算方法。文献[10]采用概率性方法,考察了大规模风电出力对系统峰谷差的影响机理,分析了在各种外送模式及系统调峰协调方式下,大规模风电基地对系统调峰的影响规律,但未考虑风电预测和常规机组故障停运对系统调峰的影响。已有的研究多角度考察了风电并网后的系统调峰问题,取得了较为丰富的成果。但针对长期规划阶段中风电并网后的系统调峰裕度研究较少。

长期以来,一般采用实用的确定性方法研究和分析电力系统调峰裕度,采用调峰容量或调峰容量占总装机容量的比例来评价电力系统调峰裕度。由于风电并网后,系统纳入了大量不确定性因素,使用确定性方法评价电力系统调峰裕度很难满足含大规模风电的电源规划需求。

本文以计及风电随机波动和常规机组随机停运等不确定性因素为导向,基于序贯蒙特卡洛方法,提出了适用于大规模风电并网的系统调峰裕度评估方法,并在IEEE-RTS系统中对所述方法进行了有效性验证。

1 风电并网系统调峰裕度评估原理

1.1 风电并网后系统调峰需求的变化

系统负荷与人们生产生活密切相关,有较强的规律性。调峰需求如图1所示。

图1中,净负荷曲线是将风电作为负的负荷,在负荷曲线上剔除风电出力得到的曲线。电力系统日负荷曲线描述了一日24h的负荷变化情况。风电接入系统前,由于日前负荷预测较为准确,系统的日调峰需求表现为负荷的峰谷差,即Pmn。风电接入系统后,由于风电的随机间歇性使系统的日调峰需求具有较大的不确定性。根据风电出力的预测精度并结合图1,以下分3种情况讨论风电接入后其对系统调峰的影响。

1)如果日前能够准确预测风电出力,系统调峰需求则表现为净负荷峰谷差,即Pmn1,系统日最大负荷为净负荷峰荷。由于风电具有多时间尺度波动特性,风电出力曲线与负荷曲线叠加会出现正调峰、负调峰、过调峰[10]等多种情况。风电正调峰使净负荷峰谷差变小,可以有效减少该日常规机组调峰容量需求;风电负调峰与正调峰相反,会增大系统调峰需求。可见,即使实现次日风电出力准确预测,系统依然要面对风电负调峰带来的更大调峰需求。

2)实际中,短期风电出力预测较为准确,在风电发展规模适中且小于负荷峰谷差时(未出现风电过调峰现象),认为距离日前机组安排较近的次日负荷谷荷时段风电预测准确,距离日前机组安排较远的峰荷时段无法预测。文献[6]中吉林省电网针对风电并网进行日前机组安排时采用的就是该种风电预测处理方法。此时系统调峰需求表现为原始负荷峰荷与净负荷谷荷之差,即Pmn2,系统日最大负荷为原始负荷峰荷。此时不管风电出力曲线如何与负荷曲线叠加,系统调峰需求均比风电接入前有所增大,增加幅度取决于风电谷荷时段的出力大小。

3)考察不进行风电出力预测或风电出力预测完全不准时的系统调峰需求变化。进行日前机组安排时,既要考虑可能在系统日负荷达到峰值时的风电零出力情况,也要考虑可能在日负荷到达谷值时的风电满发情况。系统日调峰需求显著增加,变为日负荷峰谷差Pmn与风电接入容量PW之和,即Pmn3,系统日最大负荷为原始负荷峰荷。

通过对以上3种情况的分析可知,风电并网后的调峰需求与风电出力预测精度密切相关。对于规划阶段的系统调峰裕度应结合实际及未来电网中风电的预测精度进行考虑。

1.2 调峰机组

电力系统中常规发电机组的种类一般分为:水电机组、火电机组和核电机组。

水电机组按水库的调节性能可分为无调节、日调节、周调节、年调节和多年调节水电机组。除无调节水电机组外,由于启停迅速,有调节水电机组是理想的调峰机组。

火电机组一般分为常规火电机组、供热机组和燃气轮机机组。供热机组在热电联产时经济性高,一般担任基荷运行,不适合作为调峰机组。燃气轮机机组由于启动快、可以频繁启动的特点,适于作为调峰机组。常规火电机组运行在额定出力或接近额定出力时的经济性较好,并且常规火电机组启停周期长,不宜频繁启动。在风电大规模并网后系统调峰压力很大的情况下,一般以牺牲常规火电机组的部分经济性为代价,以热备用的形式参与调峰。文献[2]指出风电场群出力变化率在每分钟0~1.5%之间的概率约为99%;目前国内外超临界燃煤火电机组的调峰深度可达50%以上,出力调整速率约为每分钟3%~5%。所以本文认为常规火电机组出力调节速率能够跟上风电功率的变化速率,即系统中常规火电机组热备用容量可以满足等容量风电出力变化的调峰需求。

核电厂按其技术和经济特性要求,适于承担系统基本负荷,不适合作为调峰机组运行。

1.3 调峰容量比

火电机组运行时,其正常发电出力范围受最小技术出力和最大技术出力的约束,其最小技术出力一般达到额定出力的70%,少数机组可达60%或50%。即一台100 MW的机组,在正常运行时的最小发电容量一般为70 MW,其可灵活变动的出力范围为70~100 MW,体现了机组的调峰深度。这一灵活变动的容量数值与机组额定容量之比称为调峰容量比[11]。发电机组的调峰容量比RG为:

式中:PN为发电机组的额定发电出力;PGmax为发电机组的最大技术出力,多数火电机组的最大技术出力与额定出力相等;PGmin为发电机组的最小技术出力。

将调峰容量比概念泛化,可认为能够通过启停调峰的机组调峰容量比为1,不具备调峰能力,适于承担基荷运行的机组调峰容量比为0。即有调节水电机组、燃气轮机机组的调峰容量比为1,无调节水电机组、处于热电联产期供热机组和核电机组的调峰容量比为0。

系统综合调峰容量比RS为:

式中:PNi为第i台发电机组的额定发电出力;RGi为第i台机组的调峰容量比;n为系统的机组总台数。

系统综合调峰容量比可以在一定程度上反映系统的调峰裕度。电网的调峰能力是指正常运行的机组出力与运行机组的最小技术出力之差是否满足系统的负荷峰谷差要求[7]。调峰能否成功,取决于在满足系统日最大负荷时,该日投入的机组调峰容量是否满足该日的系统峰谷差。系统有效调峰容量比RES为:

式中:M为当日正常运行的机组集合,即当日系统峰荷时投运的常规机组集合。

1.4 调峰裕度评估实现步骤

本文采用序贯蒙特卡洛模拟方法,方法中没有考虑检修因素。

1)输入用于生成风速随机序列的模型和参数,用于生成各常规机组状态转移时间序列的模型和参数,常规机组的调峰容量比,以及负荷时序曲线等原始数据。

2)生成风速随机序列,进而得到系统风电出力时序曲线。生成各常规机组状态转移时间序列。由于调峰考察周期为1d,忽略常规机组的日内状态变化,将每日常规机组可用与否的状态取作常规机组状态转移时间序列中每日的第1个小时时标所对应的状态。

3)基于1.1节,选择风电预测精度处理方法,根据风电时序出力曲线及负荷时序曲线,计算每日的调峰需求Pmnx,角标x可取作1,2或3,分别对应1.1节中3种风电预测精度处理方法,同时计算每日最大负荷Pmax。

4)按常规机组调峰容量比从大到小排序。累加前m台可用机组容量,使总容量PM刚好满足该日系统最大负荷Pmax,并通过式(3)计算PM对应的系统有效调峰比RES,该系统有效调峰容量比即为该日最大系统有效调峰容量比。如果全部处于可用状态机组容量且小于该日系统最大负荷,则认为该日为一个失负荷日,直接进行下一日模拟。

5)判断模拟日调峰是否充裕,判断式为:

如果式(4)成立,则认为模拟日调峰充裕,否则调峰不足。

6)判断全年365日模拟是否结束。如果结束,计算年失负荷日概率,否则继续该年模拟。第i年失负荷日概率PDLOLPi为:

式中:Nd为全年模拟中失负荷日的数目。

7)进行多年模拟,直至年失负荷日概率收敛,收敛判据[12]为:

式中:N为模拟年数;E(·)和σ(·)分别用于求PDLOLPi的均值和标准差。

8)计算调峰不足概率。通过对评估方法的描述可知,判断调峰不足是建立在不失负荷的基础上,其概率可以用条件概率表示,以弱化失负荷日对调峰裕度指标的影响。调峰不足时的概率PPSCIP为:

式中:A表示系统不失负荷;B表示系统调峰不足;AB表示在系统不失负荷的基础上调峰不足;NPSCI为模拟过程中调峰不足次数。

将式(8)代入式(7)可得:

2 调峰裕度评估计算流程

调峰裕度评估的计算流程如图2所示。

3 调峰裕度评估模型

调峰裕度评估模型基本结构如图3所示,该模型不考虑输电网络约束。

3.1 常规机组状态转移模型

常规发电机组采用2状态模型,即正常运行状态和故障停运状态。对建立虚拟系统状态转移循环过程,本文采用最通用的状态持续时间抽样法[12]。一般情况下,认为正常运行持续时间和故障修复时间服从指数分布。正常运行持续时间τ1和故障修复时间τ2分别为:

式中:γ1和γ2为均匀随机数;TMTTF为故障前平均持续工作时间;TMTTR为修复前平均持续故障时间。

3.2 风速模型

模拟风速时间序列的方法较多[13,14,15,16],本文采用广泛应用于风速预测的自回归滑动平均(ARMA)模型[13]。该模型的一般表达式为:

式中:yt为在t时刻序列上的值;φ1,φ2,…,φn为自回归参数;θ1,θ2,…,θm为滑动平均参数;αt是一个均值为0、方差为σα2的正态白噪声过程,即αt∈N(0,σα2)。

时刻t对应的风速vt为:

式中:μ为平均风速;σ为风速标准差。

3.3 风力发电机出力模型

风力发电机出力与风速成非线性函数关系,且通常由风力发电机运行参数来描述。常用参数为切入风速、切出风速和额定风速,其分段函数为:

式中:vci,vco,vr分别为切入风速、切出风速和额定风速;Pr为额定功率;a,b,c为系数,获得方法可参见文献[17]。

4 算例分析

本文采用IEEE-RTS 32机系统[18],装机容量为3 405 MW,最大机组容量为400 MW,最小机组容量为12 MW,年最大负荷为2 850 MW。由于测试系统未给出机组调峰深度,本文根据文献[18]描述的机组类型给出各机组调峰容量比,如表1所示。

风电场由单机容量为1.5 MW的风电机组组成;风电机组采用GE公司1.5MW风电机组模型;切入风速、额定风速、切出风速分别为3 m/s,14m/s,25m/s;采用单一风电场风速数据,即假设风电场相关度为1;风速采用ARMA(3,2)模型[14];模型自回归参数φ1,φ2,φ3分别为1.709 0,0.908 7,0.094 8;滑动平均参数θ1和θ2分别为1.092 9和-0.289 2;正态白噪声αt的标准差为0.474 762;平均风速为6.9m/s,风速标准差为3.64 m/s。算例负荷时序曲线采用文献[18]中的年负荷时序曲线,算例中未考虑负荷的不确定性。

考察风电预测准确、风电出力只在谷荷时段预测准确和不进行风电预测3种情况(即1.1节中的3种情况),分别计算不同风电并网容量下的系统调峰不足概率指标,计算结果如图4所示。

由图4可以看出:随着风电并网容量的不断增大,系统调峰不足概率逐渐增大,如果不进行风电预测(情况3),随着风电接入规模的不断增大,调峰不足概率指标会迅速上升,在此种情形中若要保持每日谷荷时段均留有与风电装机容量等同的风电出力空间,会显著增加调峰需求;与情况3相比,风电预测部分准确(情况2)和风电预测准确(情况1)时,调峰不足概率指标增长缓慢;在同一风电并网容量下,风电预测部分准确时的调峰不足概率高于风电预测准确时的情况;风电预测部分准确时的调峰不足概率指标随风电并网容量增大的速度大于风电预测准确时的情况。

算例中阐释了风电功率预测和风电装机容量对系统调峰裕度的影响,积极主动提高风电预测精度将有效缓解调峰压力。影响系统调峰裕度的因素较多,电力系统负荷特性和系统拥有的灵活可调节容量也是影响系统调峰裕度的关键因素。进行有效的负荷管理,减小负荷本身的峰谷差,以及增加可灵活调节的电源,同样可以缓解系统的调峰压力,增大系统的调峰裕度。负荷特性和电源固有的灵活性已在本文建立的调峰裕度评估模型中进行了描述,如果系统负荷特性和电源结构发生变化,只需相应地更改负荷时序曲线和机组输入数据。

实际评估中选择何种风电预测精度描述,可以根据实际情况及调峰需求等指标是否易于求解等因素进行取舍。对于不同电源方案的横向比较,可选择任一种预测精度进行描述。

风速时序模型参数估计和多风电场风速相关性模拟的准确度对调峰充裕度有效评估十分重要。算例中采用单一风电场风速数据给出了风速时序模型参数。在实际应用时需要利用实测风数据来估计风速时序模型参数,并对系统内多风电场的风速相关性进行模拟。

5 结语

序贯蒙特卡洛模拟 篇3

国内外关于低压电力线信道传输的研究主要集中在两个方面,即基于低压电力线的高速OFDM数据通信的研究和基于低压电力线的室内家具智能控制系统的研究。由于低压电力线上负载的多样性和时变性,以及复杂的网络结构,要建立全面可靠的通用低压电力线载波信道传输模型就显得极为不易。基于前人的研究成果,在传输线理论的基础上,结合电路理论的二端口级联理论,推导出带负载电力线输出信号的数学表达式;并通过蒙特卡洛随机模拟的方法,模拟负载的随机接入和撤除对信号传输的影响,研究结果为探明低压电力线信号传输特性奠定一定的理论和实验基础,也可为目前兴起的路灯智能单灯控制节能系统的路由设置提供借鉴。

1低压电力线无负载模型

在没有负载时,低压电力线理论上是一条阻抗均匀分布的传输线,即可以认为它的参数是均匀分布的。对于均匀传输导线,其几何空间是一维的,因此传输线上的电压和电流只是时间t和距离x的函数。如图1所示,如果将x方向无限长的传输线看成无限多Δx长度传输线的级联,而每一段Δx长度的传输线又利用LC网络等效,那么x方向无限长的传输线就可以用无限多的级联网络表示[1]。

设均匀传输线的单位长度来回导线的电阻为R0,电感为L0;单位长度来回导线间的电导为G0,电容为C0。对于任意短的传输线Δx,可以用集中参数来描述。结合图1,由基尔霍夫电流、电压定律得到:

$\begin{array}{l}-i+(i+\text{Δ}i)+(u+\text{Δ}u)G{}_0\text{Δ}x+\\ C{}_0\text{Δ}x\frac{\partial (u+\text{Δ}u)}{\partial t}=0(1)\\ -u+iR{}_0\text{Δ}x+L{}_0\text{Δ}x\frac{\partial i}{\partial t}+(u+\text{Δ}u)=0(2)\end{array}$

基尔霍夫电流定律描述了电路中各支路电流之间的关系,基尔霍夫电压定律描述了电路中各支路电压之间的关系,它们都与电路元件的性质无关,而只取决于电路的连接方式。由此可得到均匀传输线的偏微分方程组:

$\begin{array}{l}-\frac{\partial i}{\partial x}=G{}_{0}u+C{}_0\frac{\partial u}{\partial t}(3)\\ -\frac{\partial u}{\partial x}=R{}_{0}i+L{}_0\frac{\partial i}{\partial t}(4)\end{array}$

在正弦稳态下,由于均匀传输线的参数R0,L0,G0,C0都是常数,所以其是线性的。因此,均匀传输线上的电流及线间电压是与电源同频率的时间变量t的正弦函数,而其幅值和相角则为空间变量x的函数[1]。对于这样的时间、空间函数,可以用含有空间变量x的向量来表示,即:

$\begin{array}{l}i(x,t)=\text{Ι}\text{m}[\sqrt{2}\dot{{\rm I}}(x)\text{e}{}^{\text{j}\omega t}](5)\\ u(x,t)=\text{Ι}\text{m}[\sqrt{2}\dot{U}(x)\text{e}{}^{\text{j}\omega t}](6)\end{array}$

将式(5)和式(6)分别代入式(3)和式(4)中,两边分别对x求导化简,可得:

$\begin{array}{l}\frac{-\text{d}{}^2\dot{{\rm I}}(x)}{\text{d}x{}^2}=Y{}_{\text{A}}\frac{\text{d}\dot{U}(x)}{\text{d}x}(7)\\ \frac{-\text{d}{}^2\dot{U}(x)}{\text{d}x{}^2}={\rm Z}{}_{\text{A}}\frac{\text{d}\dot{{\rm I}}(x)}{\text{d}x}(8)\end{array}$

式中:YA=G0+jωC0;ZA=R0+jωL0。令$\gamma =\sqrt{{\rm Z}{}_{\text{A}}Y{}_{\text{A}}}$,于是:

$\begin{array}{l}\gamma =\sqrt{(G{}_0+\text{j}\omega C{}_0)(R{}_0+\text{j}\omega L{}_0)}=\alpha +\text{j}\beta (9)\\ {\rm Z}{}_0=\frac{{\rm Z}{}_{\text{A}}}{\gamma }=\sqrt{\frac{R{}_0+\text{j}\omega L{}_0}{G{}_0+\text{j}\omega C{}_0}}(10)\end{array}$

式中:γ称为传播常数,它描述单位长度的均匀传输线中电流行波及电压行波的衰减和相位变化;α为衰减系数,描述单位长度的均匀传输线中电流行波及电压行波的衰减的变化;β为相位偏移系数,描述单位长度的均匀传输线中电流行波及电压行波的相位变化[2]。由式(9),可解得:

$\begin{array}{l}\alpha =\sqrt{\frac{1}{2}(R{}_{0}G{}_0-\omega {}^{2}L{}_{0}C{}_0)+\frac{1}{2}\sqrt{(R{}_0^2+\omega {}^{2}L{}_0^2)(G{}_0^2+\omega {}^{2}C{}_0^2)}}(11)\\ \beta =\sqrt{\frac{1}{2}(\omega {}^{2}L{}_{0}C{}_0-R{}_{0}G{}_0)+\frac{1}{2}\sqrt{(R{}_0^2+\omega {}^{2}L{}_0^2)(G{}_0^2+\omega {}^{2}C{}_0^2)}}(12)\end{array}$

从式(11),式(12)可以看出,当ω=0时,$\alpha =\sqrt{R{}_{0}G{}_0}‚\beta =0$;当ω=∞时,ωL0≫R0,ωC0≫G0,这样得到:

$\begin{array}{l}\alpha =\frac{R{}_0}{2}\sqrt{\frac{C{}_0}{L{}_0}}+\frac{G{}_0}{2}\sqrt{\frac{L{}_0}{C{}_0}}(13)\\ \beta =\omega {}^2\sqrt{L{}_{0}C{}_0}(14)\end{array}$

由此可见,对于低频信号,电力线对信号的影响主要在幅度上,而相位上变化很小,表现为电阻形式。随着频率的增加,信号幅度主要受电阻和接入电容、电感的影响,与信号本身的频率没有关系;而相位除了受电容、电感的影响外,还与信号自身的频率有关系,且成平方比的关系。

以空间任意一点的坐标为起点,假设起点的电压、电流相量分别为$\stackrel{\cdot }{{\displaystyle V{}_1}}$和$\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}‚$在距离电力线l处,可求得电力线的均匀传输线的正向传输方程为:

$\begin{array}{l}\stackrel{\cdot }{{\displaystyle V{}_2}}(l)=\stackrel{\cdot }{{\displaystyle V{}_1}}\frac{(\text{e}{}^{\gamma l}+\text{e}{}^{-\gamma l})}{2}-\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}{\rm Z}{}_0\frac{(\text{e}{}^{\gamma l}-\text{e}{}^{-\gamma l})}{2}\\ =\stackrel{\cdot }{{\displaystyle V{}_1}}\text{c}\text{h}\gamma l-\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}{\rm Z}{}_0\text{s}\text{h}\gamma l(15)\\ \stackrel{\cdot }{{\displaystyle {\rm I}{}_2}}(l)=\stackrel{\cdot }{{\displaystyle -V{}_1}}\frac{(\text{e}{}^{\gamma l}-\text{e}{}^{-\gamma l})}{2{\rm Z}{}_0}+\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}\frac{(\text{e}{}^{\gamma l}+\text{e}{}^{-\gamma l})}{2}\\ =\frac{\stackrel{\cdot }{{\displaystyle -V{}_1}}\text{s}\text{h}\gamma l}{{\rm Z}{}_0}+\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}\text{c}\text{h}\gamma l(16)\end{array}$

通过分析低压电力线无负载模型的过程及结果可以看到,电力线本身的特性会受到输入信号频率ω的影响。频率越高,影响就越明显,这可认为是电力线产生的集肤效应。

2低压电力线带负载模型

总结前人研究的成果,影响低压电力线传输特性的因素主要如下:

(1) 人为干扰因素,如电网上负载随机的接入、撤出,电机的停运、启动,家用电器的开、关,功率因数补偿电容的接入、撤出等;

(2) 自然干扰因素,如风吹日晒、雷电等引起的干扰;

(3) 低压电力线自带的50 Hz工频正弦信号的干扰等[3,4]。

其中,人为干扰因素对低压电力线传输特性影响最大[5]。理论上低压电力线是均匀平衡的传输线,实际上同一部分低压电力线在不同时间段、不同家庭用户、不同设备的接入和撤除时都会造成阻抗不匹配,所以信号会出现反射、驻波等复杂现象。这些现象的组合,使得信号的衰减随距离的变化关系变得非常复杂,甚至会出现近距离点的衰减比远距离点还大的现象。 结合实际情况,给出如图2所示的低压电力线带负载模型。在负载模型中,假设:

(1) 低压电力线包括主干线和分支线,分支线上任意连接负载;

(2) 电压信号可以在主干线和分支线上自由传输;

(3) 电力线上接入负载的数量是随机变化的。

根据电路二端口网络理论,电力线接入的负载可以看作一种级联的形式。如图3所示,均匀电力线终端接任意负载Z(i)L时,从主干线l(i)T处看进去的输入端阻抗等效到主干网络的阻抗为:

${\rm Z}{}_{\text{i}\text{n}}^{(i)}=\frac{\stackrel{\cdot }{{\displaystyle V{}_1}}(i)}{\stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}(i)}={\rm Z}{}_0\frac{{\rm Z}{}_{\text{L}}^{(i)}+{\rm Z}{}_0\tanh (\gamma l{}_{\text{B}}^{(i)})}{{\rm Z}{}_0+{\rm Z}{}_{\text{L}}^{(i)}\tanh (\gamma l{}_{\text{B}}^{(i)})}(17)$

所有的分支线路都可以用类似的方法等效到主干线上。如图4所示,当完成所有分支线等效到主干线上之后,主干线就形成了一个级联的二端口网络。如图5所示,引入二端口网络正向传输T参数模型,以一个端口的输入电压V${}_1^{(i)}$和输入电流I${}_1^{(i)}$为自变量,输出电压V${}_2^{(i)}$和输出电流I${}_2^{(i)}$为因变量,可以得到T参数的正向传输矩阵:

$\begin{array}{l}A{}^{(i)}=\frac{V{}_1^{(i)}}{V{}_2^{(i)}}|{}_{{\rm I}{}_2^{(i)}}=0B{}^{(i)}=\frac{{\rm I}{}_1^{(i)}}{V{}_2^{(i)}}|{}_{{\rm I}{}_2^{(i)}}=0\\ C{}^{(i)}=\frac{V{}_1^{(i)}}{{\rm I}{}_2^{(i)}}|{}_{V{}_2^{(i)}}=0D{}^{(i)}=\frac{{\rm I}{}_1^{(i)}}{{\rm I}{}_2^{(i)}}|{}_{V{}_2^{(i)}}=0\end{array}$

令${\rm N}{}_{\text{L}}^{(i)}=\frac{{\rm Z}{}_{\text{i}\text{n}}^{(i)}-{\rm Z}{}_0}{{\rm Z}{}_{\text{i}\text{n}}^{(i)}+{\rm Z}{}_0}$,其中NL成为反射系数,l(i)T是主干网上各个分支的距离。于是:

$\begin{array}{l}A{}^{(i)}=\text{e}{}^{\gamma l{}_{\text{Τ}}^{(i)}}\frac{1+{\rm N}{}_{\text{L}}^{(i)}\text{e}{}^{-2\gamma l{}_{\text{Τ}}^{(i)}}}{1+{\rm N}{}_{\text{L}}^{(i)}}(18)\\ B{}^{(i)}=\frac{{\rm Z}{}_0}{2}\text{e}{}^{\gamma l{}_{\text{Τ}}^{(i)}}(1-\text{e}{}^{-2\gamma l{}_{\text{Τ}}^{(i)}})(19)\\ C{}^{(i)}=\frac{\text{e}{}^{\gamma l{}_{\text{Τ}}^{(i)}}}{{\rm Z}{}_0}\frac{1+{\rm N}{}_{\text{L}}^{(i)}\text{e}{}^{-2\gamma l{}_{\text{Τ}}^{(i)}}}{1+{\rm N}{}_{\text{L}}^{(i)}}\frac{{\rm Z}{}_0+{\rm Z}{}_{\text{i}\text{n}}^{(i)}\tanh (\gamma l{}_{\text{Τ}}^{(i)})}{{\rm Z}{}_{\text{i}\text{n}}^{(i)}+{\rm Z}{}_0\tanh (\gamma l{}_{\text{Τ}}^{(i)})}(20)\\ D{}^{(i)}=\frac{\text{e}{}^{\gamma l{}_{\text{Τ}}^{(i)}}}{2}(1-\text{e}{}^{-2\gamma l{}_{\text{Τ}}^{(i)}})\frac{{\rm Z}{}_0+{\rm Z}{}_{\text{i}\text{n}}^{(i)}\tanh (\gamma l{}_{\text{Τ}}^{(i)})}{{\rm Z}{}_{\text{i}\text{n}}^{(i)}+{\rm Z}{}_0\tanh (\gamma l{}_{\text{Τ}}^{(i)})}(21)\end{array}$

所以,根据二端口网络的级联理论,总的T参数矩阵为:

$\mathit{{\rm T}}=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]={\displaystyle \prod _{i=1}^{{\rm N}}{\rm T}}{}_i={\rm T}{}_1{\rm T}{}_2{\rm T}{}_3\cdots {\rm T}{}_{{\rm N}}(22)$

那么电力线二端口网络的正向传输参数方程为:

$\left[\begin{array}{l}\stackrel{\cdot }{{\displaystyle V{}_1}}\\ \stackrel{\cdot }{{\displaystyle {\rm I}{}_1}}\end{array}\right]=\left[\begin{array}{cc}A& C\\ B& D\end{array}\right]\left[\begin{array}{l}\stackrel{\cdot }{{\displaystyle V{}_2}}\\ \stackrel{\cdot }{{\displaystyle {\rm I}{}_2}}\end{array}\right]=\mathit{{\rm T}}\left[\begin{array}{l}\stackrel{\cdot }{{\displaystyle V{}_2}}\\ \stackrel{\cdot }{{\displaystyle {\rm I}{}_2}}\end{array}\right](23)$

从建立起来的二端口网络模型可以看到,模型涉及到的参数有电力线主干线长度、分支线长度、传播系数、输入阻抗和接入负载的特性等。

3基于蒙特卡洛方法的模拟仿真

蒙特卡洛法方法又称随机模拟方法或统计实验方法,是一种和一般数值计算方法有本质区别的计算方法。它是通过不断产生随机数序列来模拟物理过程或其他随机过程。基本原理是,利用随机数进行实验,得到统计特征值作为问题的数值解。它模拟的是一个随机过程。低压电力线负载的接入、撤出正是一个随机过程。所以,利用蒙特卡洛方法模拟仿真低压电力线上负载的随机接入、撤出对信号传输的影响理论上是一种可行的方法。

3.1 蒙特卡洛仿真模型

根据文章前面两部分对电力线模型的分析,结合实际情况以及前人的研究成果,在进行蒙特卡洛模拟之前,做出如下假设:

(1) 输入信号为正弦信号Asin ωt,其幅度A恒定,角度ω恒定,信号随测试时间点变化;

(2) 电力线的传播系数[2]γ=0.007 7+0.003 5j为恒值;电力线的波阻抗特性Z0为恒值;主干线l(i)T长度变化范围:15～50 m;支线l(i)B长度变化范围:10～40 m;

(3) 随机接入或撤出的负载,无论是纯阻性、感性或容性负载,其被接入或撤出的概率同等。

蒙特卡洛方法模拟负载对信号传输影响流程图如图6所示,第一步是设置测试时刻,模拟不同时刻负载变化对信号传输的影响;第二步是产生该时刻的传输信号,这里采用的传输信号是正弦信号Asin ωt;第三步是产生随机接入或撤出设备的数量,在这里设备数量的随机变化是通过随机变化电力线上设备的总量来反映设备的接入或撤出;第四步,对于每一个负载,随机选择阻抗类型和阻抗值,随机选择其主干线长度和支线长度,计算该负载的传输矩阵;第五步,完成所有负载总传输矩阵的计算,输出此时的信号电压、电流;第六步,完成所有时刻的测试。

3.2 蒙特卡洛仿真结果与分析

为了尽最大可能描述负载对电力线传输信号的影响,这里从这几个角度进行测试:一是在同样的输入信号频率下,测试不同阻抗特性且阻抗值不同的负载对信号传输的影响;二是在同样的输入信号频率下,测试随机变化电力线负载数量对信号传输特性的影响;三是测试不同输入信号频率下,信号通过电力线随机变化的负载数量和不同阻抗值时的传输特性。

(1) 输入信号频率f=0.1,电力线上随机变化设备的数量范围为1～10。仿真结果如图7和图8所示。

图7描述的是纯电阻负载,其值随测试时刻t在设定范围内随机变化,阻抗与输入信号的频率无关;图8描述的是感性负载,其值随测试时刻t在设定范围内随机变化,特性阻抗与输入信号的频率有关jwL。对比图7和图8发现:图7仍保留着正弦信号的轮廓,而图8已经完全变成脉冲信号。感性设备对电路的影响要明显大于纯电阻对电路的影响,容性设备也是如此,这与前人的研究成果相吻合的。

(2) 输入信号频率f=0.1,电力线上随机变化设备的数量范围为1～10和1～50,电力线上的设备随机选择纯电阻设备、感抗设备和容抗设备。

相比图10,图9的图像质量明显要更好。图9的信号在一定程度上反映输入信号的轮廓,而图10则不能。当电力线上负载数量在小范围内变化时,采用一些特殊设备接收是可以实现通信的;当电力线上负载数量变化很大时,输入信号则基本完全消失,电力线信道已不能用做通信。所以,采用低压电力线作通信信道,设置中继是必须的。根据实际路况,让负载数量变化只在小范围进行,这样就可以在一定程度上保证信号的顺利传输。这在抄表行业上就是所谓固定的路由距离设置,而在智能路灯控制系统中,虽然可以采用自动路由协议,但也可以设置一个默认路由,这样可以节省自动路由的路径搜索时间,有利于提高应用效率。

(3) 输入信号频率f=20 Hz和f=100 kHz,电力线上随机变化设备的数量范围为1～10,电力线上的设备随机选择纯电阻设备、感抗设备和容抗设备。仿真结果如图11和图12所示。

根据前人研究的结果,低压电力线最佳通信频段在f=100 kHz左右,所以在这里分别测试了输入信号为:f=20 Hz和f=100 kHz。对比图11和图12,从通信的角度来看,图12的信号更可能还原出输入信号。图12的接收信号变化规律更遵循输入信号的变化规律,相比图11具有更加明显的单峰尖锐特性,这些因素都更有利于还原输入信号。在低压窄带通信中,采用扩频通信,把输入信号淹没在噪声之中,就是因为f=100 kHz左右时,信号更容易被接收和还原。

4结论

实验表明,低压电力线信道在不同阻抗特性的负载下,对信号传输的影响有很大的差别。纯电阻信道比纯电容或纯电感信道对信号的影响较小;随机接入或者撤出不同数量的设备,对信道的特性影响也有不同,数量变化范围大的,信道特性更加恶劣;同时,不同频率的传输信号,在信道中传输受到的影响也不一致,这提供了一个低压电力线传输频率的“频率窗口“。这表明采用二端口网络理论建立起的模型,是可行的;同时,采用的蒙特卡洛法来模拟也是行之有效的。

参考文献

[1]陈崇源.电路理论:端口网络与均匀传输线[M].武汉:华中科技大学出版社,1997.

[2]谢东垒,陈显圣,陈滟涛,等.均匀传输线传播常数的测试[J].汕头大学学报:自然科学版,2008,23(4):251-254.

[3]要趁红.基于扩频通信技术的远程抄表系统[D].西安:西安建筑科技大学,2008.

[4]张允霞.基于载波通讯的智能路灯控制系统的研究[D].哈尔滨:哈尔滨工业大学,2008.

[5]ZHAI Ming-yue.Transmission characteristics of low-voltage distribution networks in China under the smartgrids environment[J].IEEE Transactions on Power Deliv-ery,2011,26(1):173-179.

[6]KORKI M,HOSSEINZADEH N,VU H L,et al.A chan-nel model for power line communication in the smart grid[C]//Procee-dings of 2011IEEE/PES Power Systems Con-ference and Exposition.[S.l.]:IEEE,2011:1-7.

[7]ANATORY E J,THEETHAYI N,THOTTAPPILLI R,et al.Broadband power-line communication channel model:comparison between theory and experiments[C]//Proceed-ings of 2008IEEE International Symposium on Power LineCommunications and Its Applications.[S.l.]:IEEE,2008:322-324.

序贯蒙特卡洛模拟 篇4

一、蒙特卡洛法解决结构可靠性的基本原理

1.蒙特卡洛法基本理论及在结构可靠性中的应用范围

蒙特卡洛法是基于Bernoulli (贝努里) 大数定理与中心极限定理。

贝努里大数定理假设是N次独立重复试验中时间A发生的次数, P是每次试验中A发生的概率, 则:

在概率统计中, 事件A发生概率稳定于事件A在一次试验中发生的概率是指频率与P有较大偏差是小概率事件, 当试验次数足够多时, 可以用频率的近似值代替P这种稳定也成为依概率稳定。

而中心极限理论则是假设随机变量序列相互独立, 服从统一分布, 且有期望与方差:

基于上诉概率统计理论, 蒙特卡罗法是一种特殊的技术, 在无需任何物理实验情况下便可以产生大量的数据结果。对于结构可靠性分析, 只需通过建立结构构件分析所适用的概率统计模型, 并使用这种分布模型产生大量的数学数据对结构进行分析。模拟结构可靠性的基本流程可用与以下三种情况:

(1) 它可用于解决不可能或是极端复杂的封闭解析解。

(2) 它可以用于解决复杂的解析解如果一些简单的假设成立的话。

(3) 它可用于检查其他其他技术解答的答案。

2.结构可靠性指数定义

在开始结构可靠性分析前, 必须对“安全”与“可靠性”进行定义:

如果R指代抗力, 而Q指代荷载效应, 代表结构失效的可能性, 数学表达式为:

即当抗力小于荷载效应时, 结构被定义为失效的。

当进行结构可靠性分析时候, 可将所有变量转化为一种标准形式, 对于基本的变量R与Q, 这种标准形式可转化为:

通常称为简化变量, 通过公式 (6) 、 (7) , 极限状态方程可以表示为:

如果将结构可靠性指数定义为变异系数的反函数, 将公式 (8) 转化为二维线性模型, ZR定义为y轴, ZQ定义为x轴, 那么可靠性指数可以被理解为到简化变量函数直线的最短距离, 称之为β, 通过几何图形, 由以下方程计算出可靠性指数β:

虽然这种计算结构可靠性方法较为常用, 当然, 也存在其局限性, 当极限函数方程中, 变量间没有任何联系, 也有其他的方法可用于计算可靠性指数, 例如Hasofer-Lind所提出的可靠性指数的计算方法便是其中之一。

3.使用蒙特卡洛法计算结构构件可靠性指数

(1) 可靠性指数蒙特卡洛模拟计算流程

采用蒙特卡洛法计算结构可靠性指数步骤如下:

(1) 决定适用概率统计模型;

(2) 根据结构可靠性分析精度决定实验模拟次数;

(3) 建立对随机变量的抽样方法, 产生随机数据;

(4) 计算极限状态函数;

(5) 统计分析在失效的可能性;

(2) 可靠性指数蒙特卡洛模拟计算案例

计算使用美国ACI规范计算公式及材料, 研究案例考虑以下极限方程涉及到R (抵抗力) , D (恒荷载) , L (活荷载) 方程为:

恒荷载 (D) :平均值与标准值的比值λD, 变差系数Vd

活荷载 (L) :平均值与标准值的比值λL, 变差系数Vl

抗力 (R) :平均值与标准值的比值λR, 变差系数Vr

为各变量选择相符的概率模型, 在此案例中恒载与活载被定义为正态分布模型, 抗力被定义为对数分布模型。在概率模型定义完成后, 获取所需的分布概率模型数据, 各变量的平均数与方差如下:

在获得相关数据后, 需按照蒙特卡洛法选择模拟次数, 由于不同随机变量取值依据的随机数序列是独立分别产生, 各随机数序列的随机数在产生顺序上是随机的, 因此, 不同随机数序列的相同顺序上的随机数组合也是随机的。这每一组随机数组合对应的随机变量取值的组合模拟了一个状态, 称为一个蒙特卡洛模拟事件。

在概率统计中, 数据量越大, 获得的精确度越高, 但同时根据不同实验选择合适的模拟次数才能最有效地进行相关的计算, 在传统概率统计学中, 为事件发生的概率, 假设期望事件发生的概率低于X%, 将变异系数控制在Y%以内, 那么所需模拟的次数则为:

为变异系数, 由以下公式进行计算:

通常, 当工程处于正常使用极限状态时, 可靠度大约在0~1.5之间, 当工程结构承载力为极限状态时, 可靠度在2.7~4.2之间, 此案例中, 假定期望事件发生的概率为15%, 变异系数控制在10%以内, 因此计算模拟次数为566次, 实际模拟次数选择600次。

使用软件随机产生600组数据, 要求数值位于0和1之间。由于正态分布在结构可靠性分析中具有重要的作用, 通常考虑标准正态分布, 从而采用以下公式生成一组正态分布随机数:

即为累积分布函数的反函数。接下来为各变量产生随机模拟数据:

来计算极限状态函数:

使用软件对数据进行整理, 并计算事件发生的可能性:

将D、L、R、g在图表中进行整理, 与标准正态标量制成二维曲线, 得到结构可靠性指数以及失效概率。

二、结语

序贯蒙特卡洛模拟 篇5

随着近几年经济的迅猛发展, 各类大型工业投资项目每年不管规模上还是数量上都在一步步实现突破, 给市场经济带来了巨大的活力。然而, 一般企业由于缺乏评估专业技术人员, 还没有形成完善的项目技术经济评价体系对各工业投资项目风险的预测, 往往仅靠主观判断或工作经验, 导致投资项目的经济效益与预估值之间总存在一定的差距。若要工业投资项目的投资效益达到预期的指标, 科学、客观、全面的经济评价是尤为重要的。项目经济评价是根据国民经济和社会发展战略以及行业、地区发展规划的要求, 在做好产品市场需求预测及厂址选择、工艺技术选择等工程技术研究的基础上, 计算项目的效益和费用, 同时从多方案中比选出最佳方案, 保证项目在财务上具备可行性, 在国民经济上具备合理性, 同时实现社会效益的显著性[1]。在实际工程应用中, 我们往往仅注重盈亏平衡分析和敏感性分析, 但二者均没有考虑到不确定因素变动的随机性, 因此在应用上还存在着较大的盲目性。概率分析可以弥补这一点, 其所包含的解析法和蒙特卡洛法, 由于大多是针对投资项目实施过程中某个具体的问题进行分析, 所以在不同的方面各有其优势和劣势, 适用于不同的情况。蒙特卡洛法是一种以概率统计为指导的数值计算方法, 它利用数值模拟来解决随机变量在实际工作中所产生的问题, 本文就蒙特卡洛法及其如今在工业项目中的应用进行说明。

1 相关现状研究与评述

纵观近年来我国在工业领域的发展, 不难看出随着经济稳定、快速的增长, 项目的规模、投资和复杂度都已上升到了一个新高度, 工业发展已步入了高峰期。早在建国初期, 我国是从前苏联学习评价理念、沿用评价模式。评价指标由于当时经济水平不高使用的还是成本、劳动生产率等经济指标。随着改革开放相关政策的实施, 同时虚心向国外学习先进的评价理念, 我国的项目经济评价水平开始逐渐提升。目前, 我国项目经济评价有了坚实的理论基础, 在实际的运用中也有了较为全面的方法指导, 但一旦出现问题便无法应对, 不能制定完备的应急计划, 从而也无法采取风险防范措施。可以说, 我国目前工程在技术上虽取得了进步, 但在管理水平上仍有巨大的挑战, 仍缺乏一种可以对项目状态准确反映的判断系统和一套精确的量化指标。

国外的项目经济评价起源于资本主义初期, 早期的经济评价仅重视财务评价, 忽视了国民经济评价。随着经济的发展, 财务评价和国民经济评价渐渐成为彼此的补充[2]。资本主义时期, 经济评价仅从私人资本的角度评估项目, 力求项目的利润最大化。美国预算体系于20 世纪60 年代出台, 这一预算体系不仅完善了国民经济评价方法, 同时实现了经济评价的推广和应用。如今, 加州大学的Mulholand则是以项目管理过程中的信息数据经逻辑判断后建立出一个工程费用全周期动态分析模型;加拿大的Brenda博士采用蒙特卡洛模拟建立出的模型来解决工程进程中的状态分析。如上文所述, 国内外项目经济评价都拥有其漫长的发展史, 随着经济和社会的发展, 其理论和应用也在不断地转变和完善。但仍然不够具体、不够全面, 还存在着项目评定指标不完善、判定方法中的主观成分等问题有待解决。

2风险分析

我们通常将项目实际值与预期值之间常存在着的差异称为风险性。因此, 投资项目风险是指由项目实施过程中的不确定因素所导致投资项目实施后偏离预期结果而造成损失的可能性[3]。

投资项目风险分析是在项目实施之前, 依具体标准对项目实施过程中所可能存在风险的发生概率和影响程度进行分析, 并针对不同的风险状况制定出不同的回避措施或应对方法, 尽可能降低项目的不确定性。在前文中已提及, 盈亏平衡分析和敏感性分析均没有考虑到不确定因素变动的随机性, 对不确定因素发生概率的分析称为概率分析, 概率分析方法包括解析法和蒙特卡洛模拟法两种。

3蒙特卡洛法

蒙特卡洛法 (Monte Carlo Method) 得名于蒙特卡洛赌场, 因为与赌场中很多随机抽样类游戏很相似, 由20 世纪40 年代由洛斯阿拉莫斯国家实验室“曼哈顿计划”的成员John von Neumann、Stanis law Ulam和尼古拉斯率先提出[4]。当时正处于第二次世界大战, 主要是为核武器的研发和制作而工作。如今蒙特卡洛法因可以实现对结果变量概率分布的估算, 所以常被用于投资项目风险分析。蒙特卡洛法建立在概率论和数理统计的基础上, 若要求解某事件的发生概率, 或者是某随机变量的期望值时, 采用大量的随机抽样试验来模拟风险的发生情况, 用风险出现的频率来估计该事件的发生概率或其它量化指标, 进而可以实现对随机事件的风险分析。

3.1 蒙特卡洛法原理

根据概率的定义可知, 通过多次的试验, 在实验结果中某事件发生的频率可以用来估计概率的分布情况。蒙特卡洛法的工作原理正是基于该思路建立的。先是对项目相关的随机变量进行反复多次的随机抽样, 将得到的抽样组一一作为输入值代入到功能函数式中, 所输出的结果概率是一个近似值, 模拟试验的次数越高, 失效概率和可靠度指标就会越精确。

其中独立的随机变量xi ( i = 1 , 2 , 3 , … , n ) , 其对应的概率密度函数分别为f (x1) , f ( x2) , …, f ( xn) 功能函数式为表达式:

首先通过随机抽样得到n组随机数x1, x2, … , xn的抽样值, 然后带入功能函数y = g ( x1, x2, … , xn) 中, 如果随机数对应的功能函数值y ≤ 0 的有L组, 那么当n无穷大时, 即可得到y的结构失效概率及可靠指标。

综上所述, 不难看出此方法巧妙地避开了结构可靠度分析中的计算部分, 大大减轻了其难度。不用考虑该函数是否非线性、随机变量分布是否非正态。因此, 用数学模型来表示被试验的目标变量的表达式:

式中, xi ( i = 1 , 2 , 3 , …, n ) 是n个相互独立的随机变量;y是n个变量的函数, 即分析项目风险效果的评价指标。

该模型将重要的风险变量综合起来, 且各风险变量都有其概率分布。通过各变量的概率分布产生随机数并分别得到抽样值, 然后将各变量的抽样值带入评价模型进行函数值的计算。通过多次重复上述试验, 风险的总体状况得以呈现[5]。

3.2 蒙特卡洛法步骤

1) 编制风险列表。采用风险识别将所存在的敏感风险列入到标准化风险列表中。

2) 采用专家调查法, 根据专家的知识和经验, 对项目中各风险因素进行评价, 确定其发生概率和估计其对结果的影响。

3) 量化专家调查中的主观成分, 同时在风险因素组合安排时体现出来。

4) 采用大量的模拟试验可以得出项目风险的概率分布曲线, 从图中读出该投资项目总体风险情况, 为日后针对每一个环节的问题采取怎样的应对措施提供了有力的依据。

3.3 蒙特卡洛法优势

1) 蒙特卡洛法概念清晰, 计算简单, 可以解决很多复杂问题。

2) 概率收敛和问题维数无关。

3) 只要抽样的次数足够多, 模拟所得结果的精度就能得以保证, 因此该种方法往往被用作其他方法的检验结果使用。

4) 数值模拟的误差比较容易确定, 从而确定模拟的次数和精度。但是, 为了保证计算结果的精度, 蒙特卡洛法需要进行大量的实验, 工作量很大, 这也是该方法久久得不到广泛推广应用的原因。随着计算机技术的发展和数值模拟方法的改进, 该问题将会很快得以解决。

4蒙特卡洛法在工业中的应用

现代工业项目随着复杂度不断地提升, 可变因素和不确定因素也在随之增长。因此, 在进行项目的经济评价中, 一贯的评价方法满足不了社会现状的需要, 经济评价工作渐渐失去了意义。专业技术研究人员日益将注意力放在概率分析上, 概率分析方法在对项目不确定因素的分析和研究工作中渐渐走入了主流的队伍。但在工业投资项目中, 不确定因素往往不止一个, 而是有其复杂联系关系的多个因素。在解决这类问题时, 才体现出蒙特卡洛法其独特的优势。本文就我国工业项目中用蒙特卡洛模拟分析的部分文献和资料进行了粗浅的总结:

谢贤平、柴建设 (1995) 对矿产资源综合利用开发的风险进行分析, 结合蒙特卡洛模拟法实现计算机程序的编制, 并进行仿真, 得出其风险估计结果。杨作安 (2001) 通过结合矿山项目的实际生产状况, 利用计算机编程技术将蒙特卡洛模拟与净现值进行结合, 对矿山项目的经济动态实现跟踪。吴爱祥、张卫峰 (2000) 利用Mathematica并结合某金矿的实际情况完成了计算机程序的编写, 这种方法及其程序对其他矿山也具有一定的参考价值。刘晓明, 李洪岩 (2008) 将净现值作为指标, 运用蒙特卡洛模拟法完成了某市基础设施项目收益风险评估模型的构建[6]。卢识峰、刘纳兵 (2007) 结合某房地产的投资实例, 对蒙特卡洛模拟分析方法的步骤进行了详细说明。吴立寰 (2004) 运用蒙特卡洛模拟法对某市公路工程的备选方案进行模拟, 科学地为投资决策提供依据。在其他行业的应用, 连春光 (2009) 运用蒙特卡洛模拟分析法完成了对造船项目投资风险框架的构建, 在Matlab的辅助下, 研究和分析了原始数据的变量因素概率分布, 以净现值为输出指标, 得到了其仿真模型。林君晓、姜鹏飞等 (2006) 在天津某污水处理厂的污水处理实际应用中, 结合蒙特卡洛模拟给出的可行性评价, 为投资者提供了科学有效的决策依据。

5总结与展望

目前国内工业投资项目往往将不确定因素对项目的影响视为固定值, 该方法存在一定的局限性, 不仅不能准确地反映出经济评价指标的分布及其所带来的影响, 还影响了项目的正确决策。本文通过对目前工业投资项目经济评价现状、风险分析再到蒙特卡洛法的具体说明, 将蒙特卡洛法在经济评价中的独特优势凸现出来, 该方法解决了国内工业项目经济评价现状中所存在的关键问题。蒙特卡洛法通过对投资项目的风险及其所产生的影响进行分析估计和量化计算, 为项目投资和决策提供更全面、更科学、更合理的依据。展望日后的经济评价工作, 进一步建立不确定因素的概率分析模型和开发应用程序将成为研究工作的重点。

参考文献

[1]丁勇.电网电源开发的区域选择与厂址优化[D].南京:东南大学, 2005.

[2]蒋洪霞.德利公司垃圾渗沥液处理生产项目财务评价研究[D].沈阳:沈阳工业大学, 2006.

[3]姚斌.基于最大熵原理的投资项目风险分析[D].天津:天津大学, 2003.

[4]朱陆陆.蒙特卡洛方法及应用[D].武汉:华中师范大学, 2014.

[5]姚慧丽, 连春光, 白婧贤.基于Monte Carlo模拟的造船项目投资风险分析框架构建[D].南京:南京航空航天大学, 2012.

序贯蒙特卡洛模拟 篇6

电力系统风险评估中所采用的蒙特卡洛模拟法主要分为非序贯蒙特卡洛模拟法和序贯蒙特卡洛模拟法[1,2]。

1.1 非序贯蒙特卡洛模拟法

非序贯蒙特卡洛模拟法也叫状态抽样法[3], 在电力系统风险评估中被广泛的使用。假定每一个元件的状态有两种:工作、失效, 并且元件与元件之间的失效状态是两个独立事件, 那么每一元件可用一个在[0, 1]区间的均匀分布来模拟, 如果让Si代表元件i的状态, Qi代表其失效概率, 对于元件i来说在[0, 1]区间均匀分布的随机数Ri, 使

在抽样中选定一个系统状态后, 首先进行系统分析以判断其状态, 如果此时系统状态为失效, 那么对该状态的风险指标函数进行估计。如果增加抽样的次数使抽样的数量足够的多, 那么就可以将系统状态s的抽样频率作为其概率的无偏估计。在此基础上, 对系统的每一个状态的概率都进行无偏估计, 即可使用公式 (3) 至 (9) 计算系统失效频率、系统失效平均持续时间、系统失效概率以及系统其它风险指标。

由于在一个时间点上, 系统只能处于一个状态, 所以系统状态之间是互斥的, 可以得到系统的累计失效率为所有失效状态概率的加和:

在公式中 (5) 中, G表示状态的集合。累计失效率为:

公式 (6) 中, frs表示系统由状态r向状态s的转移频率。系统的累计失效频率的如下近似表达式:

在得到系统的累计失效频率之后, 系统处于失效状态的平均持续时间如下式所示:

如果使用公式 (7) 计算Ff, 那么平均失效持续时间在这里是一个近似值。

每个系统失效状态的任何其他风险指标函数都可通过系统风险分析技术得到, 例如负荷削减C (s) , 所有系统失效状态指标函数的数学期望公式如下所示:

1.2序贯蒙特卡洛模拟法

序贯蒙特卡洛法是按照时序, 在一个时间跨度上进行的模拟[4,5]。状态持续时间抽样法对元件状态持续时间的概率分布进行抽样, 该方法的实现步骤如下所示: (1) 设定初始状态, 一般情况下都假定所有的元件在开始时处于运行状态。 (2) 然后对每一元件在当前状态的持续时间进行抽样。设定状态持续时间的概率分布, 不同状态的设定不同的时间概率分布。公式 (10) 表示的是指数分布:

式中:Ff、Df、Pf表示的是频率、平均持续时间和系统失效概率;Ddk表示系统第k个停运状态的持续时间;Duj表示系统第j个运行状态的持续时间;Mdn表示模拟时间内系统失效状态出现的次数;Mup表示模拟时间内系统运行状态出现次数。

算法流程图如图4所示。

2 以海南省电网为例进行电力系统风险评估

2.1 海南省电网现状概述

根据2015年海南电网规划数据, 电网总装机容量将达到6850MW, 包含联网两回线路合计120MW容量, 各类型电源占比情况如图5所示。海南省电网最大负荷为5200MW, 全社会总用电量为290亿千瓦时。220k V线路形成了环绕海南日字形的双环网结构, 而各分区也形成了115k V供配电网架。

2.2 海南省电网风险评估结果及分析

2.2.1 系统风险评估结果

四种典型运行方式下系统风险评估结果如表1和表2所示。

根据上述评估结果可知:四种典型运行方式中只有夏大运行方式下电网存在过负荷风险;四种典型运行方式下电网均存在电压越限的风险, 其中冬大运行方式下风险最大, 冬小运行方式最小, 且电压越限的概率均接近于1。

2.2.2 系统脆弱性分析。

(1) 过负荷风险脆弱性分析。夏大运行方式下电网过负荷风险及概率分区分布情况如图6所示。由图6可知, 夏大电网系统过负荷风险主要分布于东方地区, 而过负荷概率主要分布于东方、澄迈地区, 从而澄迈地区夏大线路过负荷的风险值并不大, 但是概率值较高。 (2) 电压越限风险脆弱性分析。海南电网冬季运行方式下电压越限情况分区分布如图7所示。由图7可知, 冬大运行电压越限风险主要分布在三亚地区, 而相比于冬大方式, 冬小方式三亚地区电压越限风险所占比例增大。

3 结束语

文章以非序贯蒙特卡洛并施以加快算法对海南省电网进行了风险评估计算, 并对计算结果做了较为详细的分析, 验证了该算法技术在电力系统风险评估的有效性和可实践性。

摘要：文章首先概述了电力系统风险评估的概念及基本理论, 然后详细介绍了风险评估中最为常用的蒙特卡洛法, 主要分为序贯蒙特卡洛法和非序贯蒙特卡洛法。文章采用海南电网数据进行风险评估分析, 风险评估结果验证了电力系统风险分析在电力系统中的有效性、可实践性。

关键词：风险评估,模拟,有效性,可实践性

参考文献

[1]Fang Z.Computer Simulation and Monte-Carlo Method[M].Beijing:Publishing House of Beijing Industrial Institute, 1988.

[2]Rubinstein R Y.Simulation and Monte Carlo Method[M].New York:Wiley, 1981.

序贯蒙特卡洛模拟 篇7

期权是最重要的金融衍生工具之一。合理定价是期权发挥功能的基础。对于欧式期权, 我们已经有经典的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