CN104638644B - 一种含风电场电力系统动态随机最优潮流获取方法 - Google Patents

Abstract

本发明公开了一种含风电场电力系统动态随机最优潮流获取方法，本发明用于解决风速和负荷的随机性和时空相关性影响下的电力系统动态最优潮流。本发明首先建立风速的动态概率模型，分析风速的时间和空间相关性。然后，采用基于原对偶解耦内点法的确定性动态最优潮流计算得到最优调度方案。接着，在此调度方案下，基于半不变量法求解计及相关性的动态随机潮流，从而得到状态变量的概率分布，并据此调整机会约束的上下界。最后，迭代计算解得一组满足所有机会约束的最优调度方案。本发明能够有效处理输入变量随机性影响下的电力系统最优潮流问题，具有结果准确、实现方便的优点，所得结果对调度人员具有一定指导意义。

Description

一种含风电场电力系统动态随机最优潮流获取方法

技术领域

本发明属于电力系统运行分析和控制技术领域，特别涉及一种含风电场电力系统动态随机最优潮流获取方法。

背景技术

最优潮流(optimal power flow,OPF)是电力系统规划和运行的重要工具，考虑输入变量随机性的OPF模型主要分为2类：概率最优潮流(probabilistic optimal powerflow,POPF)模型和随机最优潮流(stochastic optimal power flow,SOPF)模型。前者是根据已知随机变量的概率分布，获得若干状态变量和控制变量的概率分布信息，属于随机性分析问题；而后者是建立满足若干机会约束条件的最优化模型，在优化过程中考虑了随机因素的影响，最终得到一组满足一定机会约束的最优解，属于随机规划问题。实际电力系统是一个动态变化的系统，因而有学者提出了动态最优潮流(dynamic optimal power flow,DOPF)模型，用于解决调度周期内的优化问题。然而传统DOPF模型忽略了输入变量的随机性，在新能源不断接入电网的形势下，亟需建立计及输入变量随机性的DOPF模型，并考虑其求解方法。

发明内容

发明目的：本发明的目的在于针对现有技术的不足，提供了一种考虑了输入变量随机性的含风电场电力系统动态随机最优潮流获取方法。

技术方案：技术方案：本发明提供一种含风电场电力系统动态随机最优潮流获取方法，包括以下步骤：

步骤1：建立风速的动态概率模型，分析风速的时间和空间相关性；

步骤2：建立基于机会约束规划的含风电场电力系统动态随机最优潮流(dynamicstochastic optimal power flow,下文中简称DSOPF)模型；

步骤3：不考虑输入变量的随机性，将输入变量的期望代入，应用原对偶解耦内点法求解确定性动态最优潮流(dynamic optimal power flow,下文中简称DOPF)，得到一组最优调度方案；

步骤4：在这组最优调度方案下，考虑输入变量的随机性和时空相关性，利用基于半不变量法的解析法求解动态随机潮流(dynamic probability power flow,下文中简称DPPF)，得到状态变量的概率分布；

步骤5：判断状态变量是否满足机会约束限制；如果满足机会约束限制，则停止计算，输出结果；否则，调整机会约束的上下界，并转至步骤3，进行迭代计算，直至找到一组满足所有机会约束限制的调度方案。其中，输出结果包括发电费用、各发电机的有功和无功出力曲线、各节点电压和支路潮流的数字特征和概率分布等

进一步，步骤1包括以下步骤：

步骤101：根据公式x(t)＝X(t)+ε(t)获得实际风速序列；式中，x(t)为实际风速序列，X(t)为t时刻的风速预测序列，ε(t)为时刻t的风速预测误差序列；由于实际风速既具有一定的随机性，也具有一定的规律性，所以实际风速序列可以看作风速预测序列和预测误差序列的叠加；

步骤102：自回归滑动平均(auto regression moving average，下文简称ARMA)模型是研究时间序列的重要方法，建立自回归滑动平均模型，(p，q)阶风速时间序列的自回归滑动平均模型如下：

式中，p和q分别表示AR模型和MA模型的阶数,x_t为时刻t的风速序列值，和θ_q分别为自回归(下文中简称AR)模型和滑动平均(下文中简称MA)模型的参数；a_t为时刻t的高斯白噪声序列；

步骤103：ARMA模型的扩充Yule-Walker方程下式所示，求解Yule-Walker方程得到ARMA模型对应时间序列的自相关函数ρ_t；

式中，G_q为ARMA序列的格林函数，σ_X和σ_a分别为风速序列和高斯白噪声序列的标准差；

步骤104：研究T个时段内风速序列的相关性，根据步骤103中的公式求得风速序列的时间相关系数矩阵如下：

步骤105：由于地理位置相近的风电场之间的风速在空间上存在较强的相关性，所以研究K座风电场间风速的相关性，则这K座风电场在时刻t风速的空间相关系数矩阵如下：

式中，t＝1,2,…,T，ρ_KK,t为第K和K座风电场在第t个时刻的风速的相关系数。

进一步，步骤2中DSOPF模型为：

将调度周期分为T个时段，每个时段持续时间为τ，且每个时段内的负荷和风速保持不变；DSOPF模型可表示为如下：

(1)目标函数

以发电费用最小作为目标函数，包括常规发电机组和风电场的发电费用：

$\min f=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}[\underset{i\∈S_G}{\mathrm{\Σ}}\mathrm{\τ}(a_i{\stackrel{\‾}{P}}_{\mathrm{Gi}\left(t\right)}^2+b_i{\stackrel{\‾}{P}}_{\mathrm{Gi}\left(t\right)}+c_i)+\underset{j\∈S_W}{\mathrm{\Σ}}\mathrm{\τ}{d}_{j}E\left({\tilde{P}}_{\mathrm{Wj}\left(t\right)}\right)]$

式中，上标“—”表示确定性变量，上标“～”表示随机变量，S_G和S_W分别为常规发电机组和风电场的节点集合，和场在t时刻的有功出力，E(·)表示求随机变量的期望，a_i,b_i,c_i为第i个常规发电机组的成本系数，d_j为第j个风电场的成本系数；

(2)静态约束

考虑到负荷和风速具有随机性，所得的状态变量(包括节点电压和支路潮流等)均是随机变量，因而采用机会约束描述状态变量的不等式约束。静态约束在所有时段均需要满足，为书写方便，以下变量均省略下标t(t＝1,2,…,T)。DSOPF的静态约束包括：

${\tilde{V}}_m\underset{n\∈S_B}{\mathrm{\Σ}}{\tilde{V}}_n(G_{\mathrm{mn}}\cos {\widetilde{\mathrm{\δ}}}_{\mathrm{mn}}+B_{\mathrm{mn}}\sin {\widetilde{\mathrm{\δ}}}_{\mathrm{mn}})-{\stackrel{\‾}{P}}_{\mathrm{Gm}}-{\tilde{P}}_{\mathrm{Wm}}+{\tilde{P}}_{\mathrm{Dm}}=0(m\∈S_B)$

$P_{\mathrm{Gi},\min }\≤{\stackrel{\‾}{P}}_{\mathrm{Gi}}\≤P_{\mathrm{Gi},\max }(i\∈S_G)$

$V_{i,\min }\≤{\stackrel{\‾}{V}}_i\≤V_{i,\max }(i\∈S_G)$

$\mathrm{Pr}\{V_{i^{\′},\min }\≤{\tilde{V}}_{i^{\′}}\≤V_{i^{\′},\max }\}\≥p_{V_{i^{\′}}}(i^{\′}\∈S_B[-\left]S_G\right)$

$\mathrm{Pr}\{Q_{\mathrm{Gi},\min }\≤{\tilde{Q}}_{\mathrm{Gi}}\≤Q_{\mathrm{Gi},\max }\}\≥p_{Q_{\mathrm{Gi}}}(i\∈S_G)$

$\mathrm{Pr}\{S_{{\mathrm{ll}}^{\′},\min }\≤{\tilde{S}}_{{\mathrm{ll}}^{\′}}\≤S_{{\mathrm{ll}}^{\′},\max }\}\≥p_{S_{{\mathrm{ll}}^{\′}}}({\mathrm{ll}}^{\′}\∈S_L)$

式中，和为节点m的电压幅值和相角，G_mn和B_mn分别是系统节点导纳矩阵中第m行第n列元素的实部和虚部，和分别是节点m的有功负荷和无功负荷，和分别为第i个常规发电机组和风电场的无功出力，为线路ll'的视在功率，下标“max”和“min”分别表示模型中对应物理量的上下界，Pr{·}表示不等式约束成立的概率值，为对应物理量的预设置信水平值，S_B为系统节点集合，S_L为系统线路集合

(3)动态约束

动态约束主要考虑发电机爬坡速率约束和供购电合同约束等：

$-R_{\mathrm{ampGi}}\≤{\stackrel{\‾}{P}}_{\mathrm{Gi}\left(t\right)}-{\stackrel{\‾}{P}}_{\mathrm{Gi}(t-1)}\≤R_{\mathrm{ampGi}}(i\∈S_G)$

$C_{i,\min }\≤\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\stackrel{\‾}{P}}_{\mathrm{Gi}\left(t\right)}\mathrm{\τ}\≤C_{i,\max }(i\∈S_G)$

式中R_ampGi为第i个常规发电机组在相邻时段间的最大爬坡功率，C_i,min和C_i,max分别为第i个常规发电机组在调度周期由合同确定的供电量上下限。

进一步，步骤3包括以下步骤：

步骤301：在不等式约束中引入松弛变量，将其转化为等式约束；

步骤302：引入障碍函数对松弛变量进行约束，构造拉到格朗日函数；

步骤303：得到最优解所满足的KKT条件；

步骤304：利用解耦方法求解KKT方程，得到最优解。

进一步，步骤4包括以下步骤：

步骤401：采用线性化交流潮流模型，以注入功率期望作为输入量对每个时段进行确定性潮流计算，并在基准运行点处对潮流方程进行泰勒展开，忽略2次以上的高次项，得到：

$\left\{\begin{array}{c}{X}_{\left(t\right)}=X_{0\left(t\right)}+S_{0\left(t\right)}\mathrm{\Δ}{W}_{\left(t\right)}\\ Z_{\left(t\right)}=Z_{0\left(t\right)}+T_{0\left(t\right)}\mathrm{\Δ}{W}_{\left(t\right)}\end{array}\right.$

式中，W_(t)为时段t的节点注入功率，X_(t)为节点状态变量，Z_(t)为支路状态变量，下标“0”表示基准运行点，S_0(t)和T_0(t)为灵敏度矩阵，J₀为基准点处的雅克比矩阵，

步骤402：将T个时段的输入变量和状态变量按照时间顺序组成一个新的矩阵，得到下式：

$\left\{\begin{array}{c}\mathrm{\Δ}{X}_1=S_1\mathrm{\Δ}{W}_1=S_1\mathrm{diag}\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{\mathrm{KT}}\end{array}\right\}\mathrm{\Δ}{W_1}^{*}\\ \mathrm{\Δ}{Z}_1=T_1\mathrm{\Δ}{W}_1=T_1\mathrm{diag}\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{\mathrm{KT}}\end{array}\right\}\mathrm{\Δ}{W_1}^{*}\end{array}\right.$

其中： $\left\{\begin{array}{c}\mathrm{\Δ}{X}_1={\left[\begin{array}{cccc}\mathrm{\Δ}{X}_{\left(1\right)}^T& \mathrm{\Δ}{X}_{\left(2\right)}^T& ...& \mathrm{\Δ}{X}_{\left(T\right)}^T\end{array}\right]}^T\\ \mathrm{\Δ}{Z}_1={\left[\begin{array}{cccc}\mathrm{\Δ}{Z}_{\left(1\right)}^T& \mathrm{\Δ}{Z}_{\left(2\right)}^T& ...& \mathrm{\Δ}{Z}_{\left(T\right)}^T\end{array}\right]}^T\\ \mathrm{\Δ}{W}_1={\left[\begin{array}{cccc}\mathrm{\Δ}{W}_{\left(1\right)}^T& \mathrm{\Δ}{W}_{\left(2\right)}^T& ...& \mathrm{\Δ}{W}_{\left(T\right)}^T\end{array}\right]}^T\\ S_1=\mathrm{diag}\left\{\begin{array}{cccc}{S}_{0\left(1\right)}& S_{0\left(2\right)}& ...& S_{0\left(T\right)}\end{array}\right\}\\ T_1=\mathrm{diag}\left\{\begin{array}{cccc}{T}_{0\left(1\right)}& T_{0\left(2\right)}& ...& T_{0\left(T\right)}\end{array}\right\}\end{array}\right.$

式中，ΔX_(t)、ΔZ_(t)和ΔW_(t)分别表示在时刻t的节点状态变量、支路状态变量和节点注入功率的扰动，σ_KT为ΔW₁中第KT个元素对应输入变量的标准差，diag{σ₁ σ₂ … σ_KT}表示对角线元素为σ₁,σ₂,…σ_KT，其余元素为0的方阵，ΔW₁ ^*为ΔW₁标准化后得到的输入变量组成的向量；

步骤403：考虑到ΔW₁ ^*的相关系数矩阵为C_KT，此相关系数矩阵包含时间和空间两部分相关性。将C_KT进行Cholesky分解得到下三角矩阵B_KT：

$C_{\mathrm{KT}}=B_{\mathrm{KT}}{B}_{\mathrm{KT}}^T$

步骤404：将ΔW₁ ^*中各随机变量表示为互不相关的随机变量ΔY的线性组合：

ΔW₁ ^*＝B_KTΔY

步骤405：进一步得到：

$\left\{\begin{array}{c}\mathrm{\Δ}{X}_1=S_1\mathrm{diag}\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{\mathrm{KT}}\end{array}\right\}{B}_{\mathrm{KT}}\mathrm{\ΔY}=S_{11}\mathrm{\ΔY}\\ \mathrm{\Δ}{Z}_1=T_1\mathrm{diag}\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{\mathrm{KT}}\end{array}\right\}{B}_{\mathrm{KT}}\mathrm{\ΔY}=T_{11}\mathrm{\ΔY}\end{array}\right.$

对于服从正态分布和离散分布的输入量可以采用常规数值方法求取ΔY的半不变量，对于服从其他分布函数或分布函数未知的输入变量需要采用蒙特卡罗抽样的方法求取半不变量，由此求得ΔY的各阶半不变量，其中，S₁diag{σ₁ σ₂ … σ_KT}B_KT＝S₁₁，T₁diag{σ₁ σ₂… σ_KT}B_KT＝T₁₁；

步骤406：由于ΔY中各变量互不相关，输出状态变量的各阶半不变量可由下式得到：

$\left\{\begin{array}{c}\mathrm{\Δ}{X}_1^{\left(k\right)}=S_{11}^{\left(k\right)}\mathrm{\Δ}{Y}^{\left(k\right)}\\ \mathrm{\Δ}{Z}_1^{\left(k\right)}=T_{11}^{\left(k\right)}\mathrm{\Δ}{Y}^{\left(k\right)}\end{array}\right.$

式中，和ΔY^(k)分别表示ΔX₁、ΔZ₁和ΔY中各随机变量的k阶半不变量所组成的矩阵，所述k取值为7，和分别为矩阵S₁₁和T₁₁中元素的k次幂所组成的矩阵；

步骤407：采用Gram-Charlier级数或Cornish-Fisher级数展开方法求得输出状态变量的概率分布。

进一步，步骤5包括以下步骤：

步骤501：机会约束可用统一的形式表示为：

$\mathrm{Pr}\{x_{\min }\≤\tilde{x}\≤x_{\max }\}\≥p_{\mathrm{rob}}$

其中，表示节点电压和支路功率组成的随机变量，的概率分布函数为F(x)，可以将上式转化为：

F(x_max)-F(x_min)≥p_rob

步骤502：判断步骤501中的公式是否成立，如果成立，则满足机会约束限制；否则，不满足机会约束。当不满足机会约束时，有以下两种情况：

情况1： $\left\{\begin{array}{c}F\left(x_{\min }\right)\≈0\\ F\left(x_{\max }\right)\≤p_{\mathrm{rob}}\end{array}\right.$

情况2： $\left\{\begin{array}{c}F\left(x_{\min }\right)\≥1-p_{\mathrm{rob}}\\ F\left(x_{\max }\right)\≈1\end{array}\right.$

步骤503：针对步骤502中2种不满足机会约束的情况，按照下式调整机会约束的上下界：

情况1： $\left\{\begin{array}{c}{x}_{\min 2}=x_{\min }\\ x_{\max 2}=x_{\max }\·\max \{1-\frac{F^{-1}\left(p_{\mathrm{rob}}\right)-x_{\max }}{x_{\max }},1-\mathrm{\α}\}\end{array}\right.$

情况2： $\left\{\begin{array}{c}{x}_{\min 2}=x_{\min }\·\min \{1+\frac{x_{\min }-F^{-1}\left(p_{\mathrm{rob}}\right)}{x_{\min }},1+\mathrm{\α}\}\\ x_{\max 2}=x_{\max }\end{array}\right.$

式中，x_max2和x_min2分别为调整后的状态变量上下界，F^-1(·)为的逆概率分布函数，α为调整系数，α取为0.1。

工作原理：本发明首先建立风速的动态概率模型，分析风速的时间和空间相关性。然后，采用基于原对偶解耦内点法的确定性动态最优潮流计算得到最优调度方案。接着，在此调度方案下，基于半不变量法求解计及相关性的动态随机潮流，从而得到状态变量的概率分布，并据此调整机会约束的上下界。最后，迭代计算解得一组满足所有机会约束的最优调度方案。

有益效果：本发明具有如下优点和技术效果：

(1)建立了一种基于机会约束规划的含风电场电力系统动态随机最优潮流模型，有效的解决风速和负荷的随机性和相关性影响下的电力系统动态最优潮流的获取问题；

(2)考虑风速的时空相关性，所得结果更接近真实结果；

(3)采用基于半不变量和原对偶解耦内点法的启发式算法进行求解，具有良好的收敛性，有效缩短了获得最优潮流的时间，从而可以快速将系统调至最优状体，所得结果具有工程实用性。

附图说明

图1为DSOPF的流程图；

图2为改进IEEE14节点系统结构图；

图3为风速序列自相关函数；

图4为风速预测曲线和日负荷曲线；

图5为节点电压幅值期望；

图6为节点电压幅值标准差；

图7为发电机有功出力曲线；

图8为发电机无功出力曲线；

图9为4种方案下节点电压幅值概率密度曲线。

具体实施方式

以下结合附图和实例对本发明的实施作进一步说明，但本发明的实施和包含不限于此。

如图1所示，以IEEE14节点系统为例采用本发明提供的方法计算DSOPF。如图2所示，IEEE14节点系统共有5台常规发电机组、11个负荷节点、17条线路支路和3条变压器支路，在节点10、11、13和14处分别接入额定功率50MW的风电场。将调度周期分为24个时段，每个时段持续时间为1h。以发电机有功出力上限的10％表示相邻时段允许的最大爬坡功率，发电机节点电压上下界设定为1.10和0.90，其余节点电压上下界设定为1.05和0.95，机会约束的置信水平均取为95％。各风电场均采用恒功率因数控制，功率因数均为0.98。

ARMA模型如下式所示：

x_t＝1.5047x_t-1-0.6635x_t-2+0.1150x_t-3+a_t-0.8263a_t-1+0.2250a_t-2

式中，a_t∈N(0,0.447 423²)，N(·)表示独立正态分布。

如图3所示，求得风速序列的自相关函数。如图4所示，4座风电场的风速预测曲线，假设风速序列服从正态过程，取风速误差的标准差为风速预测值的20％，同一时段4座风电场间风速的相关系数矩阵取为：

$C_K=\left[\begin{array}{cccc}1& 0.7& 0.7& 0.7\\ 0.7& 1& 0.7& 0.7\\ 0.7& 0.7& 1& 0.7\\ 0.7& 0.7& 0.7& 1\end{array}\right]$

风电场的有功出力-风速模型由下式表示：

式中，k₁＝P_r/(v_r-v_ci)；k₂＝-k₁v_ci；P_r为风电场的额定功率，取值为50MW；v_ci、v_r和v_co分别为风机的切入风速、额定风速和切出风速，取值分别为3m/s，14m/s和25m/s。

在原有负荷的基础上按比例增加200MW的有功负荷，以平衡风电场的有功出力，日负荷曲线如图4所示。假设负荷预测误差服从正态过程，算例中给定负荷数据作为负荷的预测值，以负荷预测值的5％作为负荷预测误差的标准差。仅考虑负荷预测的时间相关性，不考虑负荷的空间相关性。

所得各节点电压幅值的期望和标准差如图5和6所示，各常规发电机组的有功和无功出力曲线如图7和8所示。

从图5可得，系统除发电机节点外所有节点电压幅值期望均处于0.95～1.05之间，且与越限边界之间存在一定裕度，因而通过DSOPF计算使系统在保持良好的电压水平的基础上取得较好的经济性，同时能够经受负荷和风速预测的误差。

由图6可以看出，系统节点电压幅值的标准差在时段12～20内的较大，由于在这些时段内，负荷和风速预测值较大，导致负荷和风速的波动性较大，因而系统节点电压幅值具有较强的波动性。

同时，由图7和图8可以看出，常规发电机组在时段12～18内的出力较低，由于在此时段内，风电场的出力较高，而风电相比常规发电机组具有更强的经济性，因而在优先使用风力的前提下，常规发电机组的出力反而减小。

将本发明提供的方法与现有技术提供的方法的试验数据进行比较，讨论相关性对系统运行特性的影响：

方法1：计及相关性的DSOPF；即为本发明提供的方法的一种实施例；

方法2：不计及相关性的DSOPF；即为本发明提供的方法的另一种实施例；

方法3：在确定性DOPF计算结果基础上进行计及相关性的DPPF；即为现有技术中常用的方法

方法4：在确定性DOPF计算结果基础上进行不计及相关性的DPPF。即为现有技术提供的方法。

方法3和方法4实际上是将输入变量的期望代入DOPF模型，进行确定性计算得到优化调度方法，并在此调度方案下进行DPPF计算，而不考虑机会约束的限制。

4种方法下节点11在时段9的节点电压幅值的概率密度曲线如图9所示。方法3和方法4中节点电压幅值越界的概率均为50.0％，而方法1和方法2中节点电压越界的概率分别只有4.8％和4.9％，由此进一步佐证了本文所提SDOPF的有效性和实用性。由图9可以看出，方案1和3得到的节点电压幅值相比方案2和4具有更强的波动性，这是由于输入变量的相关性使得输入变量同时取得较大或较小取值的概率增加，进而使得节点电压处于低压和高压的概率增加，因而节点电压具有较强的波动性。波动性较强说明结果更加精确。

综上，本发明提供的方法获得的结果精确性更高。以上仿真结果验证了本发明所提的方法的有效性和实用性。该方法能在输入变量具有随机性和时空相关性的情况下，准确得到系统最优调度方案，且该调度方案能够经受和适应输入变量的波动，具有一定的工程实用价值。

Claims (3)

1.一种含风电场电力系统动态随机最优潮流获取方法，其特征在于：包括以下步骤：

步骤1：建立风速的动态概率模型，分析风速的时间和空间相关性；

步骤2：建立基于机会约束规划的含风电场电力系统动态随机最优潮流模型；

步骤3：不考虑输入变量的随机性，代入输入变量的期望，应用原对偶解耦内点法求解确定性动态最优潮流，得到一组最优调度方案；

步骤4：在这组最优调度方案下，考虑输入变量的随机性和时空相关性，利用基于半不变量法的解析法求解动态随机潮流，得到状态变量的概率分布；

步骤5：判断状态变量是否满足机会约束限制；如果满足机会约束限制，则停止计算，输出结果；否则，调整机会约束的上下界，并转至步骤3，进行迭代计算，直至找到一组满足所有机会约束限制的调度方案；

所述步骤1包括以下步骤：

步骤101：根据公式x(t)＝X(t)+ε(t)获得实际风速序列；式中，x(t)为实际风速序列，X(t)为t时刻的风速预测序列，ε(t)为时刻t的风速预测误差序列；

步骤102：建立自回归滑动平均模型，(p，q)阶风速时间序列的自回归滑动平均模型如下：

式中，p和q分别表示自回归模型和滑动平均模型的阶数,x_t为时刻t的风速序列值，和θ_q分别为自回归模型和滑动平均模型的参数；a_t为时刻t的高斯白噪声序列；

步骤103：自回归滑动平均模型的扩充Yule-Walker方程如下式所示，求解Yule-Walker方程得到自回归滑动平均模型对应时间序列的自相关函数ρ_t；

式中，G_q为自回归滑动平均序列的格林函数，σ_X和σ_a分别为风速序列和高斯白噪声序列的标准差；

步骤104：研究T个时段内风速序列的相关性；根据步骤103中的公式求得风速序列的时间相关系数矩阵如下：

步骤105：研究K座风电场间风速的相关性，则这K座风电场在时刻t风速的空间相关系数矩阵如下：

式中，t＝1,2,…,T，ρ_KK,t为第K和K座风电场在第t个时刻的风速的相关系数。

2.根据权利要求1所述的含风电场电力系统动态随机最优潮流获取方法，其特征在于：所述步骤3包括以下步骤：采用数学规划模型对确定性动态最优潮流的物理模型进行描述

步骤301：在不等式约束中引入松弛变量，将其转化为等式约束；

步骤302：引入障碍函数对松弛变量进行约束，构造拉到格朗日函数；

步骤303：得到最优解所满足的KKT条件；

步骤304：利用解耦方法求解KKT方程，得到最优解。

3.根据权利要求1所述的含风电场电力系统动态随机最优潮流获取方法，其特征在于：所述步骤4包括以下步骤：

步骤401：采用线性化交流潮流模型，以注入功率期望作为输入量对每个时段进行确定性潮流计算，并在基准运行点处对潮流方程进行泰勒展开，忽略2次以上的高次项，得到：

$\left\{\begin{array}{c}{X}_{\left(t\right)}=X_{0\left(t\right)}+S_{0\left(t\right)}\mathrm{\Δ}{W}_{\left(t\right)}\\ Z_{\left(t\right)}=Z_{0\left(t\right)}+T_{0\left(t\right)}{\mathrm{\ΔW}}_{\left(t\right)}\end{array}\right.$

式中，W_(t)为时段t的节点注入功率，X_(t)为节点状态变量，Z_(t)为支路状态变量，下标“0”表示基准运行点，S_0(t)和T_0(t)为灵敏度矩阵，J₀为基准点处的雅克比矩阵，

步骤402：将T个时段的输入变量和状态变量按照时间顺序组成一个新的矩阵，得到下式：

$\left\{\begin{array}{c}{\mathrm{\ΔX}}_i=S_1{\mathrm{\ΔW}}_1=S_{1}diag\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{KT}\end{array}\right\}{\mathrm{\ΔW}}_1^{*}\\ {\mathrm{\ΔZ}}_1=T_1{\mathrm{\ΔW}}_1=T_{1}diag\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{KT}\end{array}\right\}{\mathrm{\ΔW}}_1^{*}\end{array}\right.$

其中：

式中，ΔX_(t)、ΔZ_(t)和ΔW_(t)分别表示在时刻t的节点状态变量、支路状态变量和节点注入功率的扰动，σ_KT为ΔW₁中第KT个元素对应输入变量的标准差，diag{σ₁ σ₂ … σ_KT}表示对角线元素为σ₁,σ₂,…σ_KT，其余元素为0的方阵，ΔW₁ ^*为ΔW₁标准化后得到的输入变量组成的向量；

步骤403：将ΔW₁ ^*的相关系数矩阵C_KT进行Cholesky分解得到三角矩阵B_KT：

$C_{KT}=B_{KT}{B}_{KT}^T;$

步骤404：将ΔW₁ ^*中各随机变量表示为互不相关的随机变量ΔY的线性组合：

ΔW₁ ^*＝B_KTΔY；

步骤405：进一步得到：

$\left\{\begin{array}{c}{\mathrm{\ΔX}}_1=S_{1}diag\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{KT}\end{array}\right\}{B}_{KT}\mathrm{\Δ}Y=S_{11}\mathrm{\Δ}Y\\ {\mathrm{\ΔZ}}_1=T_{1}diag\left\{\begin{array}{cccc}{\mathrm{\σ}}_1& {\mathrm{\σ}}_2& ...& {\mathrm{\σ}}_{KT}\end{array}\right\}{B}_{KT}\mathrm{\Δ}Y=T_{11}\mathrm{\Δ}Y\end{array}\right.$

对于服从正态分布和离散分布的输入量可以采用常规数值方法求取ΔY的半不变量，对于服从其他分布函数或分布函数未知的输入变量采用蒙特卡罗抽样的方法求取半不变量，由此求得ΔY的各阶半不变量,其中，S₁diag{σ₁ σ₂ … σ_KT}B_KT＝S₁₁，T₁diag{σ₁ σ₂ …σ_KT}B_KT＝T₁₁；

步骤406：由于ΔY中各变量互不相关，输出状态变量的各阶半不变量由下式得到：

$\left\{\begin{array}{c}{\mathrm{\ΔX}}_1^{\left(k\right)}=S_{11}^{\left(k\right)}{\mathrm{\ΔY}}^{\left(k\right)}\\ {\mathrm{\ΔZ}}_1^{\left(k\right)}=T_{11}^{\left(k\right)}{\mathrm{\ΔY}}^{\left(k\right)}\end{array}\right.$

式中，和ΔY^(k)分别表示ΔX₁、ΔZ₁和ΔY中各随机变量的k阶半不变量所组成的矩阵，所述k取值为7，和分别为矩阵S₁₁和T₁₁中元素的k次幂所组成的矩阵；

步骤407：采用Gram-Charlier级数或Cornish-Fisher级数展开方法求得输出状态变量的概率分布。