一种基于Copula函数的时空相关风光序列模拟方法与流程

本发明属于风能发电和光伏发电技术领域，尤其是一种考虑时空相关性的风光序列的模拟方法。

背景技术：

发展低碳经济，实现经济和社会的可持续发展，已经成为人类社会的普遍共识，可再生新能源开发浪潮席卷全球。近年来，我国风电和光伏装机大幅增加，随着风能和太阳能的大规模利用，同一地区将会有多个风电场和光伏电站同时接入系统。由于新能源基地的大规模集中建设，风电场和光伏电站之间的距离较近，风速和光辐照度序列之间的时空相关性较强，这将会对场站的规划和运行产生一定的影响，而且随着风电场和光伏电站装机容量的不断增大和厂址密集程度的不断增加，这种影响会更加明显。因此，十分有必要研究如何准确地刻画多个风电场和光伏电站之间的时空相关性模型。

目前，基于风光相关性对电力系统规划和运行的影响的研究较多，研究内容包括电力系统可靠性评估、概率最优潮流以及电力市场等，都对电力系统的规划和运行具有指导意义。尽管取得了一系列成果，但大多数研究者仅针对风光序列单一的时间相关性或者空间相关性进行模拟生成技术的研究，而没有同时考虑风光序列的时空相关性。另外，很多研究者仅针对风光之间的时空相关性进行分析说明并得出相应的结论而未涉及序列的模拟生成技术。

因此，有必要提出一种基于copula的考虑时空相关性的风光序列的模拟方法，对于大规模集中建设的风电场和光伏电站的规划选址和装机容量的选取等方面具有重要意义。

发明目的

针对现有时空相关风光序列模拟技术研究的不足，本发明的目的在于提出一种基于copula函数的时空相关风光序列模拟方法，该方法考虑了多维风速和光辐照度序列之间的时空相关性，可以准确地刻画多个风电场和光伏电站风光序列之间的空间特性以及每个风电场和光伏电站风光序列自身的时间特性，模拟产生更符合实际的时空相关风光序列，为风电场和光伏电站的规划和运行提供指导，具有较强的通用性和工程实用性。

为了达到上述目的，本发明采用如下的技术方案：

一种基于copula函数的时空相关风光序列模拟方法，包括以下步骤：

1)通过风电场和光伏电站获取测量得到的n个风电场的历史风速数据和m个光伏电站的历史光辐照度数据，并根据实际需要确定需模拟的时间长度t；

2)采用核密度估计法来估计风电场的历史风速序列和光伏电站的历史光辐照度序列的边缘密度函数，从风电场和光伏电站历史风光序列中采样得随机样本{x1,x2,…,xn}，在任意点x处的核密度估计函数为：

其中，k(·)为核函数，h为窗宽，n为样本容量，xi为样本风速或光辐照度；

核密度估计法中最核心的两个问题是核函数的选取和窗宽的选取，根据实证，窗宽取为：

其中为样本标准差，n为样本容量；

当窗宽确定后，核函数的选取对核密度估计影响不大，选取高斯核函数，即

由核密度估计函数积分得多维历史风速序列和历史光辐照度序列的边缘分布函数；

3)基于风光空间相关性，根据n个风电场的历史风速数据和m个光伏电站的历史光辐照度数据利用极大似然估计法分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计；

根据copula理论，来自n个风电场和m个光伏电站的风光空间相关性变量{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}的联合概率密度函数为：

其中，为相应的不同类型(n+m)维阿基米德copula函数的概率密度函数，c(f1(x1,t),…,fn+m(xn+m,t))为相应的不同类型(n+m)维阿基米德copula函数的概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率密度函数，{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度；

基于风光空间相关性的对数似然函数为：

其中，θ为(n+m)维阿基米德copula函数的相关参数，f(x1,t,…,xn+m,t)为来自n个风电场和m个光伏电站的风光空间相关性变量{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}的联合概率密度函数，{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度，n为样本容量；

根据极大似然估计法原理，(n+m)维阿基米德copula函数的相关参数θ取使基于风光空间相关性的对数似然函数最大时对应的参数：

4)依据理论copula函数和经验copula函数之间欧式距离最短的原则，选取最符合n个风电场历史风速数据和m个光伏电站历史光辐照度数据之间风光空间相关性的(n+m)维copula函数c；

(n+m)维经验copula函数的表达式为：

其中，i(·)为指示函数，当括号内条件满足时为1，反之为0；{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度，n为样本容量；为(n+m)维样本的顺序统计量；

(n+m)维理论copula函数和经验copula函数之间的欧式距离为：

其中，c(·)为不同类型的(n+m)维理论copula函数，ce(·)为(n+m)维经验copula函数，n为样本容量；

根据copula理论，选取欧式距离最小的(n+m)维理论copula函数作为描述各个风电场和光伏电站历史风速序列和历史光辐照度序列之间风光空间相关性的copula函数；

5)基于风光空间相关性，利用蒙特卡洛模拟产生长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据{(z1,t,z2,t,…,zn,t,zn+1,t,…,zn+m,t),t＝1,2,…,t}；

借助风速和光辐照度的条件分布函数模拟产生风光空间相关性数据；由copula函数的sklar定理，风速或光辐照度的条件概率分布表示某一事件的概率，数值在0到1之间，根据蒙特卡洛模拟法原理，认为其与模拟产生的服从均匀(0,1)分布的随机数等同，两者建立等式关系，从而产生满足风光空间相关性的(n+m)维风光空间相关性数据{z1,z2,…,zn,zn+1,…,zn+m}：

a.首先产生(n+m)个服从均匀(0,1)分布的随机数{u1,u2,…,un,un+1,…,un+m}；

b.根据风速或光辐照度的条件概率分布特性：

其中，c0＝0，ci＝φ[f1(z1)]+…+φ[fi(zi)]，φ(·)为相应的阿基米德copula函数的生成元，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·|·)为风速或光辐照度的条件概率分布函数，{u1,u2,…,un,un+1,…,un+m}为服从均匀(0,1)分布的随机数，{z1,z2,…,zn,zn+1,…,zn+m}为所需求解的满足风光空间相关性的(n+m)维风光空间相关性数据；

已知随机数{u1,u2,…,un,un+1,…,un+m}，利用递归法求解方程组得{z1,z2,…,zn,zn+1,…,zn+m}；

c.重复以上两步t次，得到长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据{(z1,t,z2,t,…,zn,t,zn+1,t,…,zn+m,t),t＝1,2,…,t}；

6)基于风光时间相关性，根据各个风电场相邻时刻的历史风速数据和各个光伏电站相邻时刻的历史光辐照度数据利用极大似然估计法分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计；

根据copula理论，来自第i个风电场或光伏电站的风光时间相关性变量{(xi,t,xi,t+1),t＝1,2,…,ni}的联合概率密度函数为：

f(xi,t,xi,t+1)＝ci(fi(xi,t),fi(xi,t+1))fi(xi,t)fi(xi,t+1)

其中，为相应的不同类型二维阿基米德copula函数的概率密度函数，ci(fi(xi,t),fi(xi,t+1))为相应的不同类型二维阿基米德copula函数的概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率密度函数，{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度；

基于风光时间相关性的对数似然函数为：

其中，θi为二维阿基米德copula函数的相关参数，f(xi,t,xi,t+1)为来自第i个风电场或光伏电站的风光时间相关性变量{(xi,t,xi,t+1),t＝1,2,…,ni}的联合概率密度函数，{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度，ni为第i个风电场或光伏电站的样本容量；

根据极大似然估计法原理，二维阿基米德copula函数的相关参数θi取使基于风光时间相关性的对数似然函数最大时的参数：

7)依据理论copula函数和经验copula函数之间欧式距离最短的原则，分别选取最符合n个风电场历史风速序列和m个光伏电站历史光辐照度序列风光时间相关性的(n+m)个copula函数；

二维经验copula函数的表达式为：

其中，i(·)为指示函数，当括号内条件满足时为1，反之为0；{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度，ni为第i个风电场或光伏电站的样本容量；{xi(k),xi(j)}为二维样本的顺序统计量；

二维理论copula函数和经验copula函数之间的欧式距离为：

其中，ci(·)为不同类型的二维理论copula函数，cei(·)为二维经验copula函数，ni为第i个风电场或光伏电站的样本容量；

根据copula理论，选取欧式距离最小的二维理论copula函数作为描述每个风电场历史风速序列和每个光伏电站历史光辐照度序列风光时间相关性的copula函数；

8)基于风光时间相关性，借助一阶马尔科夫过程的特性以及风速和光辐照度的条件概率分布模拟产生(n+m)条长度为t的马尔科夫时间序列，即得到n个风电场和m个光伏电站时间长度为t的风速和光辐照度模拟数据；

由一阶马尔科夫过程可知，序列某时刻的状态只与最近时刻的状态有关，而与之前时刻的状态均无关；由copula函数的sklar定理可知，风速或光辐照度的条件概率分布表示某一事件的概率，数值在0到1之间，根据蒙特卡洛模拟法原理，认为其与模拟产生的服从均匀(0,1)分布的随机数等同，两者建立等式关系，由此产生满足风光时间相关性的t维马尔科夫时间序列{xi,1,xi,2,…,xi,t}：

根据一阶马尔科夫过程和风速或光辐照度的条件概率分布特性：

其中，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·|·)为风速或光辐照度的条件概率分布函数，ci(fi(xi,t),fi(xi,t-1))为相应的二维阿基米德copula函数的概率分布函数，{(zi,1,zi,2,…,zi,t),i＝1,2,…,n+m}为步骤5)所求得的长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据，{(xi,1,xi,2,…,xi,t),i＝1,2,…,n+m}为所需求解的(n+m)个t维马尔科夫时间序列；

已知由步骤5)求得的{(zi,1,zi,2,…,zi,t),i＝1,2,…,n+m}，利用递归法求解方程组得{xi,1,xi,2,…,xi,t}，最终求得(n+m)个t维马尔科夫时间序列{(xi,1,xi,2,…,xi,t),i＝1,2,…,n+m}，即得到需要的n个风电场和m个光伏电站时间长度为t的风速和光辐照度模拟数据。

本发明建立了基于copula函数的时空相关风光序列的蒙特卡洛模拟框架，考虑了当同一地区有多个风电场和光伏电站接入系统时多维风光序列之间的时空相关性，较为准确地刻画了这一地区实际风电场和光伏电站的出力模型，解决了时空相关风光序列的模拟问题。能应用于大规模集中建设的多个风电场和光伏电站的出力模拟，为风电场和光伏电站的规划选址和装机容量的选取提供指导。

附图说明

图1为基于copula函数的时空相关风光序列的蒙特卡洛模拟框架。

图2为历史数据和模拟数据的自相关系数图，其中，图2a)为美国nrelsolarradiationresearchlaboratorybaselinemeasurementsystem(bms)的风速自相关系数图，图2b)为美国nrelsolarradiationresearchlaboratorybaselinemeasurementsystem(bms)的光伏出力自相关系数图，图2c)为美国nrelnationalwindtechnologycenterm2tower(m2)的风速自相关系数图，图2d)为美国nrelnationalwindtechnologycenterm2tower(m2)的光伏出力自相关系数图。

图3为历史数据和模拟数据的四维经验分布函数值散点图。

具体实施方式

下面结合附图对本发明的具体实施方式及工作原理做进一步详细说明。

参见图1，一种基于copula函数的时空相关风光序列模拟方法，包括风光空间相关性数据的模拟生成和风光时间相关性数据的模拟生成两个部分，具体分为以下步骤：

1)通过风电场和光伏电站获取测量得到的n个风电场的历史风速数据和m个光伏电站的历史光辐照度数据，并根据实际需要确定需模拟的时间长度t。

2)分别对n个风电场历史风速序列和m个光伏电站历史光辐照度序列的边缘密度函数进行估计。具体方法如下：

本发明采用核密度估计法对风电场的历史风速序列和光伏电站的历史光辐照度序列的边缘密度函数进行估计，从风电场和光伏电站历史风光序列中采样得随机样本{x1,x2,…,xn}，在任意点x处的核密度估计函数为：

其中，k(·)为核函数，h为窗宽，n为样本容量，xi为样本风速或光辐照度。

核密度估计法中最核心的两个问题是核函数的选取和窗宽的选取，根据实证，窗宽取为：

其中为样本标准差，n为样本容量。

当窗宽确定后，核函数的选取对核密度估计影响不大，选取高斯核函数，即

由核密度估计函数积分得多维历史风速序列和历史光辐照度序列的边缘分布函数。

3)基于风光空间相关性，根据n个风电场的历史风速数据和m个光伏电站的历史光辐照度数据分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计。具体方法如下：

本发明采用极大似然估计法分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计，根据copula理论，来自n个风电场和m个光伏电站的风光空间相关性变量{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}的联合概率密度函数为：

其中，为相应的不同类型(n+m)维阿基米德copula函数的概率密度函数，c(f1(x1,t),…,fn+m(xn+m,t))为相应的不同类型(n+m)维阿基米德copula函数的概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率密度函数，{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度。

基于风光空间相关性的对数似然函数为：

其中，θ为(n+m)维阿基米德copula函数的相关参数，f(x1,t,…,xn+m,t)为来自n个风电场和m个光伏电站的风光空间相关性变量{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}的联合概率密度函数，{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度，n为样本容量。

根据极大似然估计法原理，(n+m)维阿基米德copula函数的相关参数θ取使基于风光空间相关性的对数似然函数最大时对应的参数：

4)选取最符合n个风电场历史风速数据和m个光伏电站历史光辐照度数据之间风光空间相关性的(n+m)维copula函数c。具体方法如下：

本发明通过计算各个理论copula函数和经验copula函数之间的欧式距离来确定最佳的copula函数，距离越小，拟合越好。

(n+m)维经验copula函数的表达式为：

其中，i(·)为指示函数，当括号内条件满足时为1，反之为0；{(x1,t,x2,t,…,xn+m,t),t＝1,2,…,n}为(n+m)维样本风速和光辐照度，n为样本容量；为(n+m)维样本的顺序统计量。

(n+m)维理论copula函数和经验copula函数之间的欧式距离为：

其中，c(·)为不同类型的(n+m)维理论copula函数，ce(·)为(n+m)维经验copula函数，n为样本容量。

根据copula理论，选取欧式距离最小的(n+m)维理论copula函数作为描述各个风电场和光伏电站历史风速序列和历史光辐照度序列之间风光空间相关性的copula函数。

5)基于风光空间相关性，利用蒙特卡洛模拟产生长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据{(z1,t,z2,t,…,zn,t,zn+1,t,…,zn+m,t),t＝1,2,…,t}。具体方法如下：

本发明借助风速和光辐照度的条件分布函数模拟产生风光空间相关性数据。由copula函数的sklar定理，风速或光辐照度的条件概率分布表示某一事件的概率，数值在0到1之间，根据蒙特卡洛模拟法原理，认为其与模拟产生的服从均匀(0,1)分布的随机数等同，两者建立等式关系，从而产生满足风光空间相关性的(n+m)维风光空间相关性数据{z1,z2,…,zn,zn+1,…,zn+m}：

a.首先产生(n+m)个服从均匀(0,1)分布的随机数{u1,u2,…,un,un+1,…,un+m}。

b.根据风速或光辐照度的条件概率分布特性：

其中，c0＝0，ci＝φ[f1(z1)]+…+φ[fi(zi)]，φ(·)为相应的阿基米德copula函数的生成元，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·|·)为风速或光辐照度的条件概率分布函数，{u1,u2,…,un,un+1,…,un+m}为服从均匀(0,1)分布的随机数，{z1,z2,…,zn,zn+1,…,zn+m}为所需求解的满足风光空间相关性的(n+m)维风光空间相关性数据。

已知随机数{u1,u2,…,un,un+1,…,un+m}，利用递归法求解方程组得{z1,z2,…,zn,zn+1,…,zn+m}。

c.重复以上两步t次，得到长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据{(z1,t,z2,t,…,zn,t,zn+1,t,…,zn+m,t),t＝1,2,…,t}。

6)基于风光时间相关性，根据各个风电场相邻时刻的历史风速数据和各个光伏电站相邻时刻的历史光辐照度数据分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计。具体方法如下：

本发明采用极大似然估计法分别对几类阿基米德copula函数即gumbelcopula函数、claytoncopula函数和frankcopula函数的相关参数进行估计，根据copula理论，来自第i个风电场或光伏电站的风光时间相关性变量{(xi,t,xi,t+1),t＝1,2,…,ni}的联合概率密度函数为：

f(xi,t,xi,t+1)＝ci(fi(xi,t),fi(xi,t+1))fi(xi,t)fi(xi,t+1)

其中，为相应的不同类型二维阿基米德copula函数的概率密度函数，ci(fi(xi,t),fi(xi,t+1))为相应的不同类型二维阿基米德copula函数的概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·)为历史风速序列或历史光辐照度序列的边缘概率密度函数，{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度。

基于风光时间相关性的对数似然函数为：

其中，θi为二维阿基米德copula函数的相关参数，f(xi,t,xi,t+1)为来自第i个风电场或光伏电站的风光时间相关性变量{(xi,t,xi,t+1),t＝1,2,…,ni}的联合概率密度函数，{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度，ni为第i个风电场或光伏电站的样本容量。

根据极大似然估计法原理，二维阿基米德copula函数的相关参数θi取使基于风光时间相关性的对数似然函数最大时的参数：

7)分别选取最符合n个风电场历史风速序列和m个光伏电站历史光辐照度序列风光时间相关性的(n+m)个copula函数。具体方法如下：

本发明通过计算各个理论copula函数和经验copula函数之间的欧式距离来确定最佳的copula函数，距离越小，拟合越好。

二维经验copula函数的表达式为：

其中，i(·)为指示函数，当括号内条件满足时为1，反之为0；{(xi,t,xi,t+1),t＝1,2,…,ni}为二维样本风速或光辐照度，ni为第i个风电场或光伏电站的样本容量；{xi(k),xi(j)}为二维样本的顺序统计量。

二维理论copula函数和经验copula函数之间的欧式距离为：

其中，ci(·)为不同类型的二维理论copula函数，cei(·)为二维经验copula函数，ni为第i个风电场或光伏电站的样本容量。

根据copula理论，选取欧式距离最小的二维理论copula函数作为描述每个风电场历史风速序列和每个光伏电站历史光辐照度序列风光时间相关性的copula函数。

8)基于风光时间相关性，模拟产生(n+m)条长度为t的马尔科夫时间序列，即得到n个风电场和m个光伏电站时间长度为t的风速和光辐照度模拟数据。具体方法如下：

本发明借助一阶马尔科夫过程的特性以及风速和光辐照度的条件概率分布模拟产生(n+m)条长度为t的马尔科夫时间序列。由一阶马尔科夫过程可知，序列某时刻的状态只与最近时刻的状态有关，而与之前时刻的状态均无关；由copula函数的sklar定理可知，风速或光辐照度的条件概率分布表示某一事件的概率，数值在0到1之间，根据蒙特卡洛模拟法原理，认为其与模拟产生的服从均匀(0,1)分布的随机数等同，两者建立等式关系。由此产生满足风光时间相关性的t维马尔科夫时间序列{xi,1,xi,2,…,xi,t}：

根据一阶马尔科夫过程和风速或光辐照度的条件概率分布特性：

其中，f(·)为历史风速序列或历史光辐照度序列的边缘概率分布函数，f(·|·)为风速或光辐照度的条件概率分布函数，ci(fi(xi,t),fi(xi,t-1))为相应的二维阿基米德copula函数的概率分布函数，{(zi,1,zi,2,…,zi,t),i＝1,2,…,n+m}为步骤5)所求得的长度为t的服从(n+m)维copula函数c的(n+m)维风光空间相关性数据，{(xi,1,xi,2,…,xi,t),i＝1,2,…,n+m}为所需求解的(n+m)个t维马尔科夫时间序列。

已知由步骤5)求得的{(zi,1,zi,2,…,zi,t),i＝1,2,…,n+m}，利用递归法求解方程组得{xi,1,xi,2,…,xi,t}，最终求得(n+m)个t维马尔科夫时间序列{(xi,1,xi,2,…,xi,t),i＝1,2,…,n+m}，即得到需要的n个风电场和m个光伏电站时间长度为t的风速和光辐照度模拟数据。

参见图2中图2a)、图2b)、图2c)和图2d)，基于美国能源部国家新能源实验室历史风光数据进行的案例分析，案例数据为美国nrelsolarradiationresearchlaboratorybaselinemeasurementsystem(bms)和nrelnationalwindtechnologycenterm2tower(m2)的四维历史风速和光辐照度序列。案例分析显示，利用本专利提出的模拟法产生的四个风速和光辐照度序列的自相关性与历史数据的自相关性基本一致，而利用传统的双参数韦布尔分布和贝塔分布方法产生的模拟数据的自相关性很弱，基本维持在零处稳定，与历史数据的自相关性相差甚远，从而验证了本专利提出的模拟方法产生的模拟数据更符合实际新能源场站历史风光数据的时间特性。

参见图3，历史数据和利用本专利提出的模拟法产生的模拟数据的四维经验分布函数值散点图比历史数据和利用传统方法产生的模拟数据的四维经验分布函数值散点图更接近于对角线，说明利用本专利提出的模拟法产生的模拟数据与历史数据具有更相似的空间特性，从而验证了本专利提出的模拟方法产生的模拟数据更符合实际新能源场站历史风光数据的空间特性。综上所述，本专利提出的模拟方法能准确刻画历史风光数据的时空特性，更符合实际新能源场站的模型，具有较强的通用性和工程实用性。