全球能源互联网运行特性仿真的自适应变步长积分方法与流程

本发明涉及全球能源互联网宏观运行特性领域，具体涉及一种全球能源互联网运行特性仿真的自适应变步长积分方法。
背景技术：
：能源和环境问题是当今世界共同面临的重大挑战。长期过度开发使得传统化石能源趋于枯竭，能源地缘政治事件频发引发区域战争不断，碳排放超限造成气候变暖趋势加剧。大规模开发利用清洁、低碳能源，实现能源消费向低碳转型是关系人类可持续发展的重要基础。树立全球能源观，构建全球能源互联网，统筹全球能源资源开发、配置和利用，保障能源的安全、清洁、高效和可持续供应，已成为国际社会的共识。全球能源互联网是物理系统与信息系统的有机融合，同时也是能量流与信息流的统筹结合，其具有若干复杂特征。全球能源互联网具有复杂的动力学特性，其宏观运行特性受国际政治局势、经济发展、能源格局等各种因素的影响。如何量化分析国际政治、经济、能源间的关系，设计一种适用于全球能源互联网宏观运行特性仿真的自适应变步长积分方法，从而模拟和推演在不同国际政治、市场环境下全球能源发展格局及其变化情况，对于预测国际能源供需、优化全球能源配置及构建全球能源互联网过程中的宏观决策等均具有战略性指导意义。可以预见，全球能源互联网不仅仅涉及能源系统，而且涉及国际政治系统和经济系统，其影响因素众多，采用数学模型模拟和推演的过程中势必会存在诸多变量和微分方程，非线性程度极高，而且不同子系统间时间常数差异巨大。现有数学模型往往仅考虑其中某一个或部分子系统，无法将国际政治系统-经济系统-能源系统建立在统一框架下进行仿真；此外，现有数值积分方法也无法模拟秒级和年度级时间常数共存的情形，为了能够模拟秒级数据变化，势必要减小迭代步长，大大增加计算耗时。技术实现要素：本发明为了解决上述问题，提出了一种全球能源互联网运行特性仿真的自适应变步长积分方法，本发明提取能源、政治和经济系统中的主要影响因素，摒除次要影响因素；然后进行历史数据统计分析，建立金融市场、能源价格等变量关于全球政治事件的量化函数关系，获得国际政治系统-经济系统-能源系统间的数学模型，能够针对各类突发政治事件采用不同的初值条件，在保证计算精度的基础上有效减少计算耗时。为了使本领域技术人员能够更好的明白技术方案，首先进行名词解释如下：政治系统：主要是指会对全球一次能源和电力消费产生显著影响的政治因素集合，比如重大国际或区域战争、多国间政治博弈、政治突发事件、以及石油含量十分丰富的中东区域局势等，都是在本专利中需要考虑的政治系统的构成。经济系统：主要是指会对全球一次能源和电力消费产生显著影响的经济因素集合，比如、全球gdp总额、各产业gdp、国际金融市场、能源类股票价格变化、石油价格、电力价格、美元汇率以及各种突发事件如金融危机等，都是在本专利中需要考虑的经济系统的构成。能源系统：由于本专利主要考虑能源的宏观运行特性，因此本专利中能源系统主要包括一次能源消费量、电力消费量、石油供给量、主要产油区产量、线路网损率、资源开发成本、清洁能源消费量与装机量等，通过描述各因素间的因果关系来描述各系统的动态特性。为了实现上述目的，本发明采用如下技术方案：一种全球能源互联网运行特性仿真的自适应变步长积分方法，提取能源、政治和经济中的影响因素，根据各影响因素之间的因果关系，确定速率变量和状态变量，建立全球能源互联网宏观运行特性因果关系图，对不同类型的影响因素的历史数据进行标准化处理，依次选取每个变量，查找所有指向该变量的变量集合，依据历史数据确定变量间的代数方程，完成全球能源互联网系统变量间数学公式构建，将方程差分化处理，在突发事件发生前采用改进欧拉法进行迭代，在突发事件发生后采用四阶龙格库塔法进行迭代，以自适应步长进行仿真预测，预测推演出全球能源互联网的宏观运行特性。进一步的，确定政治系统、经济系统和能源系统三个子系统的重要内部变量以及子系统间的重要边界变量；作为一种优选方案：选取政治系统变量集合为：{政治突发事件、多国博弈、国际能源政策、中东政治局势}；选取经济系统变量集合为：{全球gdp总额、人均gdp、第一产业gdp、第二产业gdp、第三产业gdp、第一产业增长额、第一产业增长率、第二产业增长额、第二产业增长率、第三产业增长额、第三产业增长率、金融危机、石油价格、电力价格、美元汇率、资源开发成本}；选取能源系统变量集合为：{居民用电量、第一产业用电量、第二产业用电量、第三产业用电量、清洁能源消费量、新增清洁能源装机、电力消费量、网损率、石油消费量、能源消费总量、储产比、石油供给、opec产量、非opec产量}。进一步的，确定各个子系统内部变量之间的因果关系以及子系统间重要边界变量的因果关系，完成全球能源互联网宏观运行特性因果关系图，以确定的因果关系图为基础，依次分析每个因果关系间是否与要添加状态变量和速率变量，其余变量均为辅助变量，在因果关系图的基础上添加速率变量，形成模型流图。进一步的，如果因果关系间的分析结果为需要添加速率变量，则证明该因果关系满足微分关系，否则，说明该因果关系满足代数方程关系。进一步的，模型流图中的变量即为全球能源互联网宏观特性仿真中的全部变量集合，收集变量集合中所有变量的历史数据，保证所有数据的维数一致。进一步的，若因果关系满足微分关系，则依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}，如果该集合{s}中有一个变量为速率变量，则该变量x为状态变量，确定状态变量及速率变量之间的微分方程，然后依据历史数据确定系数取值。进一步的，若因果关系满足代数方程关系，依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}，如果该集合{s}中所有变量均为辅助变量，则该变量x也是辅助变量，然后依据历史数据确定变量间的代数方程，完成全球能源互联网系统变量间数学公式构建。确定仿真开始时刻、终止时刻和总仿真时间，并依据突发事件类型和初值条件选取仿真步长，仿真步长的值在突发事件发生前视稳态数据精度要求而定，在突发事件发生后，仿真步长视暂态数据精度和数据量要求而定，采用变步长策略。判断所建立的代数方程是否存在时间延迟，如果存在，根据变量种类自定义延迟时间，选取为仿真步长的整数倍。进一步的，若在仿真步长△t内能够获得m组自变量x的数据，则需要分别利用m组数据带入方程中求因变量或状态量值，然后计算平均值作为t+△t时刻因变量或状态量的值；若在仿真步长△t内未能够获得自变量x的数据，则需要采用插值法求取变量x在该时刻的值，然后带入微分方程方程或代数方程求取t+△t时刻状态量的值。与现有技术相比，本发明的有益效果为：本发明方法以全球能源互联网为载体，通过全球能源互联网宏观运行特性仿真将政治系统、经济系统和能源系统有机融合在一个仿真框架内，通过选取合适变量描述三个子系统，将重大政治和经济事件量化表达，能够有效描述三个子系统的动态特性，并且预测推演全球能源互联网的宏观运行特性。具体效果如下：1)各子系统变量选取典型恰当。政治、经济和能源系统联系复杂，变量众多，非线性程度高。本方法通过选取典型变量描述三个子系统，在保证仿真精度的同时能够有效降低模型维数，提高计算速度，工程适用性强。2)建立子系统内及子系统间的定量数学函数模型。通过进行历史数据统计和调研分析，建立金融市场、能源价格与需求等关于各类事件的量化函数，准确定量分析各类事件对于政治、经济和能源系统的显著影响。3)提出一种自适应变步长积分方法。在政治事件或战争发生前后，分别采用不同步长对微分方程和代数方程进行迭代，提高计算速度；根据各类变量时间尺度差异，采用均值法或插值法处理数据过剩或进行数据补充，方法简单准确，易于程序实现。附图说明构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。图1是本发明方法流程图；图2是本发明方法因果关系图；图3是本发明方法流图。具体实施方式：下面结合附图与实施例对本发明作进一步说明。应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属
技术领域：
的普通技术人员通常理解的相同含义。需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。正如
背景技术：
所介绍的，现有技术中存在现有数学模型往往仅考虑其中某一个或部分子系统，无法将国际政治系统-经济系统-能源系统建立在统一框架下进行仿真；同时现有数值积分方法也无法模拟秒级和年度级时间常数共存的情形，为了能够模拟秒级数据变化，势必要减小迭代步长，大大增加计算耗时的不足，为了解决如上的技术问题，本申请提出了适用于全球能源互联网宏观运行特性仿真的自适应变步长积分方法。解决全球能源互联网宏观运行特性仿真难题，提取能源、政治和经济系统中的主要影响因素，摒除次要影响因素；然后进行历史数据统计分析，建立金融市场、能源价格等变量关于全球政治事件的量化函数关系，获得国际政治系统-经济系统-能源系统间的数学公式；最后，提出一种适用于全球能源互联网宏观运行特性仿真的自适应变步长积分方法，能够精确求解国际政治系统-经济系统-能源系统数学模型，并能够针对各类突发政治事件采用不同的初值条件，在保证计算精度的基础上有效减少计算耗时。下面将按照步骤详细介绍：步骤1)：确定政治系统、经济系统和能源系统三个子系统的重要内部变量以及子系统间的重要边界变量。步骤2)：确定各个子系统内部变量之间的因果关系以及子系统间重要边界变量的因果关系，完成全球能源互联网宏观运行特性因果关系图。步骤3)：以上述确定的因果关系图为基础，依次分析每个因果关系间是否与要添加状态变量和速率变量，其余变量均为辅助变量。步骤4)：按照步骤3)分析结果，在因果关系图的基础上添加速率变量，形成模型流图。步骤5)：模型流图中的变量即为全球能源互联网宏观特性仿真中的全部变量集合。收集变量集合中所有变量的历史数据，注意所有数据的维数一致。如果在步骤3)中需要添加速率变量，则证明该因果关系能够写成微分方程，执行步骤7)；否则，说明该因果关系能够采用代数方程表示，执行步骤8)。步骤6)：对不同子系统间量纲不一致的数据做标准化处理，获得标准化数据。步骤7)：依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}。如果该集合{s}中有一个变量为速率变量，则该变量x为状态变量，确定状态变量及速率变量之间的微分方程，然后依据历史数据确定系数取值。步骤8)：依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}。如果该集合{s}中所有变量均为辅助变量，则该变量x也是辅助变量，然后依据历史数据确定变量间的代数方程，完成全球能源互联网系统变量间数学公式构建。步骤9)：确定步骤8)中所建立的代数方程是否存在时间延迟。如果存在，根据变量种类自定义延迟时间，一般选取为仿真步长的整数倍。步骤10)：确定仿真开始时刻、终止时刻和总仿真时间，并依据突发事件类型和初值条件选取合适的仿真步长。本专利采取如下取值策略：在突发事件发生前，各个社会经济变量将按照步骤7)和8)所求的公式平稳变动，不会发生数据突变，因此仿真步长一般选取为最小步长的10-100倍，视稳态数据精度要求而定；在突发事件发生后，模型中的变量将发生跃变或急剧变化，因此为了精细仿真该政治事件的影响，仿真步长一般选取为最小步长的1-10倍，视暂态数据精度和数据量要求而定。步骤11)：编写适用于全球能源互联网宏观运行特性计算的数值仿真程序，如利用matlab等，(当然，本领域技术人员也可以进行别的处理手段，只要保证其处理过程可行，能够得到结果即可)，这里主要说明对于各类方程的处理策略。对于微分方程而言，首先将方程差分化处理，然后本专利在突发事件发生前采用改进欧拉法进行迭代，提高计算速度；在突发事件发生后采用四阶龙格库塔法进行迭代，提高计算精度；对于代数方程而言，将按照是否存在时间延迟做相应处理。假设最小仿真步长为△t，根据方程类型具体说明如下：步骤11-1)：对于微分方程而言，首先将方程差分化处理。在突发事件发生前采用改进欧拉法进行迭代，仿真步长视稳态数据精度选取为最小步长的10-100倍。具体公式如下：步骤11-2)：在突发事件发生后微分方程采用变步长四阶龙格库塔法进行迭代，初始仿真步长视暂态数据精度选取为最小步长的1-10倍。公式如下：k1＝gi[t,xi(t)](3)其中，k1-k4是指在每个迭代步长内各个中间点的斜率值。可以看出，在每一步迭代步长中共计算了t、t+△t/2、t+2△t/3、t+△t四个中间点处的函数值，因此仍然属于四阶龙格库塔法。另外，所谓变步长的含义是指每一步迭代步长△t可以由下一步计算结果xi(t+△t)的误差大小来控制。误差计算公式为：当e(t+△t)>emax，下一步迭代步长缩减为上一步的一半；当e(t+△t)<emin，下一步迭代步长增加为上一步的两倍；否则，步长不变。其中emax、emin为预先设定的误差最大值与最小值。步骤11-3)：对于无延迟的代数方程yi(t)＝f[t,yj(t),x(t)]i＝1,2,…,h而言，需要利用t+△t时刻的代数量yj(t+△t)和状态量xj(t+△t)计算t+△t时刻的代数量值yi(t+△t)，见公式(8)。但是，yj(t+△t)和xj(t+△t)也通过预测获得，不可避免地会有预测误差。因此，为减小误差，通常采用前推回代法多次预测yj(t+△t)和yi(t+△t)，减小交接误差。yi(t+δt)＝f[t,yj(t+δt),x(t+δt)]i＝1,2,…,h(8)步骤11-4)：对于有延迟t的代数方程yi(t)＝f[t,yj(t)]i＝1,2,…,h而言，需要利用t+△t时刻的代数量yj(t+△t)计算t+△t时刻的代数量值yi(t+△t)。但是需要注意的是，若代数量yj(t)的值发生突变，则该突变量在延迟时间t后才能发挥作用，因此，当n△t<t(n是指代数量突变后所经历的迭代步长个数)时，仍然采用突变前该变量的历史数据做外推预测；当n△t>t后，采用突变后的值按照式(8)预测。步骤12)：由于各变量历史数据的时间尺度有差异，因此若自变量x的时间尺度小于仿真步长，也即在仿真步长△t内能够获得m组自变量x的数据，则需要分别利用m组数据带入方程中求因变量或状态量值，然后计算平均值作为t+△t时刻因变量或状态量的值；若自变量x的时间尺度大于仿真步长，也即在仿真步长△t内未能够获得自变量x的数据，则需要采用插值法求取变量x在该时刻的值，然后带入微分方程方程或代数方程求取t+△t时刻状态量的值。步骤13)：执行程序，输出计算结果并分析，结束程序。本申请的一种典型的实施方式中，如图1所示，具体包括：步骤1)：确定政治系统、经济系统和能源系统三个子系统的重要内部变量以及子系统间的重要边界变量。选取政治系统变量集合为：{政治突发事件、多国博弈、国际能源政策、中东政治局势}。选取经济系统变量集合为：{全球gdp总额、人均gdp、第一产业gdp、第二产业gdp、第三产业gdp、第一产业增长额、第一产业增长率、第二产业增长额、第二产业增长率、第三产业增长额、第三产业增长率、金融危机、石油价格、电力价格、美元汇率、资源开发成本}。选取能源系统变量集合为：{居民用电量、第一产业用电量、第二产业用电量、第三产业用电量、清洁能源消费量、新增清洁能源装机、电力消费量、网损率、石油消费量、能源消费总量、储产比、石油供给、opec产量、非opec产量}。步骤2)：确定各个子系统内部变量之间的因果关系以及子系统间重要边界变量的因果关系，完成全球能源互联网宏观运行特性因果关系图，见图2。步骤3)：以上述确定的因果关系图为基础，依次分析每个因果关系间是否与要添加状态变量和速率变量，其余变量均为辅助变量。如果需要添加速率变量，则证明该因果关系能够写成微分方程，执行步骤6)；否则，说明该因果关系能够采用代数方程表示，执行步骤7)。本实例中经济子系统中第一产业gdp、第二产业gdp和第三产业gdp三个变量可以作为状态变量，添加第一产业增长额、第二产业增长额、第三产业增长额作为速率变量，同时添加第一产业增长率、第一产业增长率、第一产业增长率作为辅助变量。步骤4)：按照步骤3)分析结果，在因果关系图的基础上添加速率变量，形成模型流图，见图3。步骤5)：模型流图中的变量即为全球能源互联网宏观特性仿真中的全部变量集合。收集变量集合中所有变量的历史数据，见表3，注意所有数据的维数一致。如果在步骤3)中需要添加速率变量，则证明该因果关系能够写成微分方程，执行步骤7)；否则，说明该因果关系能够采用代数方程表示，执行步骤8)。步骤6)：对不同子系统间量纲不一致的数据做标准化处理，获得标准化数据，见表4。步骤7)：依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}。如果该集合{s}中有一个变量为速率变量，则该变量x为状态变量，确定状态变量及速率变量之间的微分方程，然后依据历史数据确定系数取值。下面以第三产业gdp为例，说明其建立的微分方程，其中增长率r(t)不是常数，而是时间的函数：y＝y0+∫r(t)·ydt其中，y为第三产业gdp，y0为第三产业gdp初始值，r(t)为第三产业增长率。表1第三产业gdp增长率年份增长率年份增长率年份增长率年份增长率19962.4720010.9520068.22201110.7919977.1920024.72200712.8720122.8119980.99200312.2200810.0920133.1619994.49200412.02009-2.6120143.0620003.3720058.3920108.75步骤8)：依次选取每一个变量x，查找所有指向该变量的变量集合，记作{s}。如果该集合{s}中所有变量均为辅助变量，则该变量x也是辅助变量，然后依据历史数据确定变量间的代数方程，完成全球能源互联网系统变量间数学公式构建。下面以全球电力消费量为例，说明其建立的代数方程如下：y＝0.161x1+0.417x2+0.259x3+0.112x4-0.002x5+0.009x6+0.125x7-0.080其中，y为全球电力消费量，x1～x7分别表示第一产业用电量、第二产业用电量、第三产业用电量、居民总用电量、全球总gdp、电力价格和全球石油总消费量。步骤9)：确定步骤8)中所建立的代数方程是否存在时间延迟。如果存在，根据变量种类自定义延迟时间，一般选取为仿真步长的整数倍。步骤10)：确定仿真开始时刻、终止时刻和总仿真时间，并依据突发事件类型和初值条件选取合适的仿真步长。本专利采取如下取值策略：在突发事件发生前，各个社会经济变量将按照步骤7)和8)所求的公式平稳变动，不会发生数据突变，因此仿真步长一般选取为最小步长的10-100倍，视稳态数据精度要求而定；在突发事件发生后，模型中的变量将发生跃变或急剧变化，因此为了精细仿真该政治事件的影响，仿真步长一般选取为最小步长的1-10倍，视暂态数据精度和数据量要求而定。步骤11)：编写适用于全球能源互联网宏观运行特性计算的数值仿真程序，这里主要说明对于各类方程的处理策略。对于微分方程而言，首先将方程差分化处理，然后本专利在突发事件发生前采用改进欧拉法进行迭代，提高计算速度；在突发事件发生后采用四阶龙格库塔法进行迭代，提高计算精度；对于代数方程而言，将按照是否存在时间延迟做相应处理。假设最小仿真步长为△t，根据方程类型具体说明如下：步骤11-1)：对于微分方程而言，首先将方程差分化处理。在突发事件发生前采用改进欧拉法进行迭代，仿真步长视稳态数据精度选取为最小步长的10-100倍。具体公式如下：步骤11-2)：在突发事件发生后微分方程采用变步长四阶龙格库塔法进行迭代，初始仿真步长视暂态数据精度选取为最小步长的1-10倍。公式如下：k1＝gi[t,xi(t)](3)其中，k1-k4是指在每个迭代步长内各个中间点的斜率值。可以看出，在每一步迭代步长中共计算了t、t+△t/2、t+2△t/3、t+△t四个中间点处的函数值，因此仍然属于四阶龙格库塔法。另外，所谓变步长的含义是指每一步迭代步长△t可以由下一步计算结果xi(t+△t)的误差大小来控制。误差计算公式为：当e(t+△t)>emax，下一步迭代步长缩减为上一步的一半；当e(t+△t)<emin，下一步迭代步长增加为上一步的两倍；否则，步长不变。其中emax、emin为预先设定的误差最大值与最小值。步骤11-3)：对于无延迟的代数方程yi(t)＝f[t,yj(t),x(t)]i＝1,2,…,h而言，需要利用t+△t时刻的代数量yj(t+△t)和状态量xj(t+△t)计算t+△t时刻的代数量值yi(t+△t)，见公式(8)。但是，yj(t+△t)和xj(t+△t)也通过预测获得，不可避免地会有预测误差。因此，为减小误差，通常采用前推回代法多次预测yj(t+△t)和yi(t+△t)，减小交接误差。yi(t+δt)＝f[t,yj(t+δt),x(t+δt)]i＝1,2,…,h(8)步骤11-4)：对于有延迟t的代数方程yi(t)＝f[t,yj(t)]i＝1,2,…,h而言，需要利用t+△t时刻的代数量yj(t+△t)计算t+△t时刻的代数量值yi(t+△t)。但是需要注意的是，若代数量yj(t)的值发生突变，则该突变量在延迟时间t后才能发挥作用，因此，当n△t<t时，仍然采用突变前该变量的历史数据做外推预测；当n△t>t后，采用突变后的值按照式(8)预测。步骤12)：由于各变量历史数据的时间尺度有差异，因此若自变量x的时间尺度小于仿真步长，也即在仿真步长△t内能够获得m组自变量x的数据，则需要分别利用m组数据带入方程中求因变量或状态量值，然后计算平均值作为t+△t时刻因变量或状态量的值；若自变量x的时间尺度大于仿真步长，也即在仿真步长△t内未能够获得自变量x的数据，则需要采用插值法求取变量x在该时刻的值，然后带入微分方程方程或代数方程求取t+△t时刻状态量的值。例如，如果仿真步长为季度(0.25年)，也就是说，对于历史数据中以年为单位的数据而言数据量是不足的。就需要按照插值法求取该变量在0.25、0.5、0.75处的值，然后才能进行后续计算。下面以第三产业增长率为例，进行说明：1996-1998年的第三产业增长率分别为2.47％、7.19％、0.99％。下面利用插值法计算1997年和1998年各个季度的第三产业增长率。表21997-1998年各个季度第三产业增长率1997年第三产业增长率1998年第三产业增长率第一季度2.47第一季度7.19第二季度4.04第二季度5.12第三季度5.62第三季度3.05第四季度7.19第四季度0.99步骤13)：执行程序，输出计算结果并分析，结束程序。表3各变量历史数据表4各变量标准化后的历史数据以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。当前第1页12