CN110717643A

CN110717643A - 一种基于全局灵敏度分析的天然气网络储气配置方法

Info

Publication number: CN110717643A
Application number: CN201910724123.2A
Authority: CN
Inventors: 赵霞; 胡潇云; 冯欣欣; 李欣怡
Original assignee: Chongqing University
Current assignee: Chongqing University
Priority date: 2019-08-07
Filing date: 2019-08-07
Publication date: 2020-01-21
Anticipated expiration: 2039-08-07
Also published as: CN110717643B

Abstract

本发明公开了一种基于全局灵敏度分析的天然气网络储气配置方法，主要步骤为：2)建立服从正态分布的气负荷概率模型F(x)，并抽样得到节点气负荷样本矩阵X。2)将负荷的每组样本代入天然气网络的稳态能流计算，得到输出变量y^ED，并计算输出变量y^ED关于节点气负荷的全局灵敏度指标Si。3)确定最大全局灵敏度指标

式中，i_max为节点气负荷变量编号；在负荷变量i_max所在节点设置储气装置。基于稀疏多项式混沌展开方法，本发明可以快速建立输出变量关于输入不确定变量的表达形式，从而有效提高全局灵敏度分析效率。

Description

一种基于全局灵敏度分析的天然气网络储气配置方法

技术领域

本发明涉及天然气网络稳态能流计算领域，具体是一种基于全局灵敏度分析的天然气网络储气配置方法。

背景技术

能流计算作为获取天然气状态分布的基本方法，是天然气网络运行和规划的基础。在实际的系统运行过程中，负荷波动等随机因素的变化通常会改变天然气网络的运行状态，甚至造成不可忽视的负面影响。通常情况下，运行人员会在天然气网络中配置储气以提高网络气压水平和运行稳定性，但出于经济性的考虑，在所有节点设置储气往往是不合理的。因此，定量评估系统中的不确定因素对于重要状态变量的影响，以识别影响系统状态的关键不确定因素，对于运行人员合理配置储气、提高天然气网络安全稳定运行具有重要的指导意义。

全局灵敏度分析是一种有效的不确定因素分析方法。它能够充分计及系统中不确定因素的概率分布特性和相关性，并量化它们对于输出变量波动的影响。目前Sobol’法是一种常用的全局灵敏度分析方法，传统的Sobol’法需要以大量确定性计算作为基础，计算效率较低。

发明内容

本发明的目的是解决现有技术中存在的问题。

为实现本发明目的而采用的技术方案是这样的，一种基于全局灵敏度分析的天然气网络储气配置方法，主要包括以下步骤：

1)获取天然气网络参数，并建立天然气网络的稳态能流模型。

所述天然气网络参数主要包括气源参数、管道参数、节点气负荷A＝[x₁,x₂,···,x_N]^T、气负荷的秩相关系数和平衡节点压力初始值。

2)基于天然气网络中节点气负荷的随机性和秩相关性，建立服从正态分布的气负荷概率模型F(x)，并基于该模型对所有气负荷进行抽样，得到节点气负荷样本矩阵X。

抽样得到节点气负荷样本矩阵X的主要步骤如下：

2.1)设定天然气网络中有1个平衡节点和N-1个非平衡节点。平衡节点编号记为1。非平衡节点均接有气负荷。天然气网络节点所在位置集合记为

网络中与任一节点a直接相连的所有节点构成的集合记为

设定抽样次数为N_ED。设定N维独立标准正态变量ξ＝[ξ₁,ξ₂,···,ξ_N]^T。

2.2)建立服从正态分布的气负荷概率模型，其累积概率函数F (x)为：

式中，x表示气负荷。X为气负荷某一实际取值。μ_x为气负荷的均值。气负荷的均值μ_x等于气负荷基准值。σ_x为气负荷的标准差。

2.3)基于拉丁超立方抽样法和高斯Copula函数，对气负荷A 进行抽样，主要步骤如下：

2.3.1)建立节点气负荷A＝[x₁,x₂,···,x_N]^T的秩相关系数矩阵C_S，即：

式中，ρ_sij为节点气负荷x_i和节点气负荷x_j的秩相关系数，i,

当i＝j时，ρ_sij＝1。当i≠j时，0≤ρ_sij≤1。

2.3.2)基于拉丁超立方抽样法，建立变量ξ的概率模型样本矩阵X_ξ，即：

2.3.3)将气负荷的秩相关系数转化为等效线性相关系数，转换公式如下：

式中，ρ_ij为在标准正态空间中节点气负荷x_i和节点气负荷x_j的等效线性相关系数。

基于等效线性相关系数ρ_ij，建立气负荷的等效线性相关系数矩阵C_X，即：

2.3.4)对矩阵C_X进行Cholesky分解，得到下三角矩阵L：

C_X＝LL^T。 (6)

2.3.5)基于独立标准正态变量样本矩阵X_ξ和下三角矩阵L，建立线性相关的标准正态变量样本矩阵X'，即：

X'＝LX_ξ。 (7)

式中，矩阵X'中各元素的线性相关性满足矩阵C_X。矩阵X'中各元素的等效秩相关性满足矩阵C_S。

2.3.6)将矩阵X'中的变量样本进行等概率转换，得到气负荷样本矩阵X。其中，矩阵X中任一元素x_im如下所示：

x_im＝F_i ^-1(Φ(x'_im))。 (8)

式中：F_i ^-1(·)表示变量x_i的累积分布函数的反函数。Ф(·)表示标准正态分布变量的累积分布函数。x'_im为矩阵X'中的元素，表示第i 个变量的第m个样本。

3)建立天然气稳态能流模型，将每一组气负荷样本作为已知量代入该模型并采用牛顿法求解，得到输出变量y的样本向量y^ED。

采用牛顿法求解天然气系统稳态能流模型，以得到y^ED的主要步骤如下：

3.1)选取矩阵X的第m列作为负荷变量的第m组样本作为天然气稳态能流模型的已知量，采用牛顿法求解天然气能流模型，得到输出变量向量

其中y_m表示气负荷的第m组样本对应的输出变量；输出变量y为任意节点压力或者管道流量；m＝1，2，…N_ED；采用牛顿法求解天然气网络能流模型的最大迭代次数为T_max，收敛判据为λ。

3.2)令某一管道两端节点编号为a和

对于本文所研究的低压配气网络(压力小于等于75mbar)，建立流经管道的气流量方程：

式中，F_a-b为由a流向b的天然气流量；p_a和p_b分别为节点a 和b的压力；s_a-b为方向变量，当p_a>p_b，s_a-b＝1，否则s_a-b＝-1；D为管道a-b的直径；f为管道的摩擦系数；L为管道长度；S_G为天然气的相对密度，标准工况下为0.5548。

仅考虑天然气网络中的管道元件，建立节点a的管道流量平衡方程：

式中，F_a为节点a的注入流量；

3.3)计算节点注入天然气流量的不平衡量

将天然气各节点压力的初值代入管道流量的计算公式，得到天然气网络中的管道流量以及各节点的注入流量。此时经过式(10) 得到的节点a注入流量与实际值之间存在差异，采用△F_a表示节点a 注入流量的不平衡量：

式中，F_a由节点a的气负荷决定，即F_a＝-x_a。

采用列向量△F统一表示所有节点注入流量的不平衡量：

ΔF＝[ΔF₁ ΔF₂ … ΔF_N]^T (12)

3.4)建立修正方程中的雅克比矩阵J

结合牛顿法可知，以节点压力p＝[p₁,p₂,···,p_N]^T作为状态变量，建立天然气网络中关于状态变量的修正方程组：

ΔF＝JΔp； (13)

式中，△p＝[△p₁,△p₂,···,Δp_N]^T为状态变量p的修正量向量；J 为雅克比矩阵。具体表示为：

式中，

表示偏微分。

3.5)计算管道节点压力修正量

基于式(12)和式(14)的计算结果，可以计算所有节点压力的修正量：

Δp＝-J^-1ΔF； (15)

根据所得修正结果可以对节点压力进行更新，即：

p^(k+1)＝p^(k)+Δp^(k)； (16)

式中，k表示迭代次数。

3.6)重复步骤3.4)和3.5)，在每次迭代过程中计算雅克比矩阵 J和节点压力修正量Δp。并根据以下收敛条件判断是否结束计算：

a.当修正量△p满足max{|Δp|}<λ，且迭代次数k≤T_max时，迭代结束，并输出目标变量y；

b.当max{|Δp|}<λ，且k＞T_max，停止迭代，输出“能流计算结果不收敛”，结束计算；

c.当max{|Δp|}≥λ，且k＞T_max；取k＝k+1，返回步骤3.4)继续迭代计算。

4)基于已有的负荷变量矩阵X和输出变量向量y^ED，采用稀疏多项式混沌展开方法计算输出变量y关于各节点气负荷的全局灵敏度指标S_i。

计算输出变量y^ED关于节点气负荷的全局灵敏度指标Si的主要步骤如下：

4.1)设置Hermite多项式作为负荷变量的最优正交多项式基函数，建立输出变量样本y^ED关于输入变量X的多项式混沌展开形式，即：

式中，f(ξ)表示输出变量样本y^ED关于输入变量X的原始表达式。

为变量

的s阶正交多项式。

为多项式展开系数。a＝[a₀,a₁,···,a_p-1]^T。H(ξ)＝[1,H₁(ξ),···, H_p(ξ)]^T。多项式混沌展开形式(17)的最高展开阶数为p。

公式(17)的展开项数M如下所示：

式中，符号！表示阶乘。

变量ξ中所有变量样本相应的各阶多项式矩阵如下所示：

式中：ξ⁽ⁱ⁾表示变量ξ的第i组样本。

4.2)计算公式(17)的稀疏展开系数，主要步骤如下：

4.2.1)设置初始参数，主要包括输出变量向量y^ED、正交多项式矩阵H(ξ)、最大允许误差ε、初始迭代次数k＝0、展开系数a(k)＝0、残差r_(k)＝y^ED-Ha_(k)和项数索引集

4.2.2)确定出正交多项式矩阵H中与残差r_(k)最相关的列C_(k)＝r_(k) ^TH。确定行向量C_(k)中数值最大的项c_max，元素c_max所在列数记为j_(k)。

4.2.3)令k＝k+1，将j_(k)添加至集合索引集Λ_(k)＝Λ_(k-1)∪j_(k-1)。

4.2.4)基于最小二乘法，求解展开系数a_(k)＝arg min_a||y^ED-H_(k)a||₂。其中，H_(k)为矩阵H中索引集Λ_(k)对应列数的向量构成的矩阵。

4.2.5)更新残差r_(k)＝y^ED-Ha(k)，判断残差r_(k)是否满足||r_(k)||₂>ε。若是，则返回步骤4.2.2)，若否，则计算结束，并输出优化结果a_(k)。其中，||·||₂表示二范数

4.3)计算输出变量y的均值E[·]和方差D[·]，即：

式中，i＝(i₁,i₂,···,i_N)表示多阶下标。|i|＝i₁+i₂+···+i_N。

4.4)结合Sobol’计算公式和公式(17)，计算输出变量关于输入变量(ξi₁,ξi₂,···,ξi_s)的全局灵敏度，即：

式中：

表示变量

的全局灵敏度指标；下标集合

表示含下标(i₁,i₂，···,i_s)中一个或多个元素；表示公式(17)中任意变量

的任意阶的多项式的平方。

集合

如下所示：

5)确定最大全局灵敏度指标

式中，i_max为节点气负荷变量编号。在负荷变量i_max所在节点设置储气装置。

本发明的技术效果是毋庸置疑的。基于稀疏多项式混沌展开方法，本发明可以快速建立输出变量关于输入不确定变量的表达形式，从而有效提高全局灵敏度分析效率。

综上所述，为合理量化天然气网络中不确定因素及其相关性对关键状态变量的影响，基于稀疏多项式混沌展开实现天然气网络的全局灵敏度分析，并以全局灵敏度指标为指导，选择储气配置节点，为提高系统安全运行状态提供支持。

本发明基于全局灵敏度分析实现对天然气网络的储气配置，考虑天然气网络的节点气负荷的不确定性，基于稀疏多项式混沌展开实现快速高效的全局灵敏度分析，以量化负荷不确定性对于输出变量的影响，并根据所得全局灵敏度指标指导储气的配置。

附图说明

图1为基于全局灵敏度分析的天然气网络储气配置方法流程图；

图2为本发明方法中基于正交匹配追踪计算稀疏多项式展开系数流程图；

图3为本发明方法所采用的天然气网络结构图；

图4为发明方法的流经管道7-8气流量关于所有气负荷的全局灵敏度仿真结果图；

图5为基于本发明方法选择最大灵敏度节点配置储气与选取管道首末节点配置储气的能流计算对比图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

参见图1，一种基于全局灵敏度分析的天然气网络储气配置方法，主要包括以下步骤：

概率抽样的主要步骤如下

1)设定天然气网络中有1个平衡节点和N-1个非平衡节点。平衡节点编号记为1。非平衡节点均接有气负荷。天然气网络节点所在位置集合记为

网络中与任一节点a直接相连的所有节点构成的集合记为

设定抽样次数为N_ED。设定N维独立标准正态变量ξ＝[ξ₁, ξ₂,···,ξ_N]^T。

2)建立服从正态分布的气负荷概率模型，其累积概率函数F(x) 为：

3)基于拉丁超立方抽样法和高斯Copula函数，对气负荷A进行抽样，主要步骤如下：

3.1)建立节点气负荷A＝[x₁,x₂,···,x_N]^T的秩相关系数矩阵C_S，即：

当i＝j时，ρ_sij＝1。当i≠j时，0≤ρ_sij≤1。

3.2)基于拉丁超立方抽样法，建立变量ξ的概率模型样本矩阵 X_ξ，即：

3.3)将气负荷的秩相关系数转化为等效线性相关系数，转换公式如下：

3.4)对矩阵C_X进行Cholesky分解，得到下三角矩阵L：

C_X＝LL^T。 (6)

3.5)基于独立标准正态变量样本矩阵X_ξ和下三角矩阵L，建立线性相关的标准正态变量样本矩阵X'，即：

X'＝LX_ξ。 (7)

3.6)将矩阵X'中的变量样本进行等概率转换，得到气负荷样本矩阵X。其中，矩阵X中任一元素x_im如下所示：

x_im＝F_i ^-1(Φ(x'_im))。 (8)

4)建立天然气稳态能流模型，将每一组气负荷样本作为已知量代入该模型并采用牛顿法求解，得到输出变量y的样本向量y^ED，主要步骤如下：

4.1)选取矩阵X的第m列作为负荷变量的第m组样本作为天然气稳态能流模型的已知量，采用牛顿法求解天然气能流模型，得到输出变量向量y^ED＝[y₁,y₂,···,yN_ED]^T；其中y_m表示气负荷的第m组样本对应的输出变量；输出变量y为任意节点压力或者管道流量； m＝1，2，…N_ED；采用牛顿法求解天然气网络能流模型的最大迭代次数为T_max，收敛判据为λ。

4.2)令某一管道两端节点编号为a和

式中，F_a为节点a的注入流量；

4.3)计算节点注入天然气流量的不平衡量

式中，F_a由节点a的气负荷决定，即F_a＝-x_a。x_a为节点a的气负荷。

采用列向量△F统一表示所有节点注入流量的不平衡量：

ΔF＝[ΔF₁ ΔF₂ … ΔF_N]^T (12)

4.4)建立修正方程中的雅克比矩阵J

ΔF＝JΔp； (13)

式中，△p＝[△p₁,△p₂,···,△p_N]^T为状态变量p的修正量向量；J 为雅克比矩阵。具体表示为：

式中，

表示偏微分。

4.5)计算管道节点压力修正量

Δp＝-J^-1ΔF； (15)

根据所得修正结果可以对节点压力进行更新，即：

p^(k+1)＝p^(k)+Δp^(k)； (16)

式中，k表示迭代次数。

4.6)重复步骤4.4)和4.5)，在每次迭代过程中计算雅克比矩阵 J和节点压力修正量Δp。并根据以下收敛条件判断是否结束计算：

a.当修正量△p满足max{|△p|}<λ，且迭代次数k≤T_max时，迭代结束，并输出目标变量y；λ为修正量阈值。

b.当max{|△p|}<λ，且k＞T_max，停止迭代，输出“能流计算结果不收敛”，结束计算；

c.当max{|△p|}≥λ，且k＞T_max；取k＝k+1，返回步骤4.4)继续迭代计算。T_max为迭代次数最大值。

5)基于已有的负荷变量矩阵X和输出变量向量y^ED，采用稀疏多项式混沌展开方法计算输出变量y关于各节点气负荷的全局灵敏度指标S_i，主要步骤如下：

5.1)设定Hermite多项式作为输出变量样本y^ED的最优正交多项式基函数，建立输出变量样本y^ED关于输入变量X的多项式混沌展开形式，即：

式中，f(ξ)表示输出变量样本y^ED关于输入变量X的原始表达式。为变量

的s阶正交多项式。

公式(17)的展开项数M如下所示：

式中，符号！表示阶乘。

变量ξ中所有变量样本相应的各阶多项式矩阵如下所示：

式中：ξ⁽ⁱ⁾表示变量ξ的第i组样本。

5.2)计算公式(17)的稀疏展开系数，主要步骤如下：

5.2.1)设置初始参数，主要包括输出变量向量y^ED、正交多项式矩阵H(ξ)、最大允许误差ε、初始迭代次数k＝0、展开系数a(k)＝0、残差r_(k)＝y^ED-Ha_(k)和项数索引集

5.2.2)确定出正交多项式矩阵H中与残差r_(k)最相关的列C_(k)＝r_(k) ^TH。确定行向量C_(k)中数值最大的项c_max，元素c_max所在列数记为j_(k)。

5.2.3)令k＝k+1，将j_(k)添加至集合索引集Λ_(k)＝Λ_(k-1)∪j_(k-1)。Λ_(k-1)为第k-1次迭代的项数索引集，j_(k-1)为第k-1次迭代的元素c_max所在列数。

5.2.4)基于最小二乘法，求解展开系数a_(k)＝arg min_a||y^ED-H_(k)a||₂。其中，H_(k)为矩阵H中索引集Λ_(k)对应列数的向量构成的矩阵。min 表示取最小值。

5.2.5)更新残差r_(k)＝y^ED-Ha(k)，判断残差r_(k)是否满足||r_(k)||₂>ε。若是，则返回步骤5.2.2)，若否，则计算结束，并输出优化结果a_(k)。其中，||·||₂表示二范数。

5.3)计算输出变量y的均值E[·]和方差D[·]，即：

5.4)结合Sobol’计算公式和公式(17)，计算输出变量关于输入变量

的全局灵敏度，即：

式中：

表示变量

的全局灵敏度指标；下标集合表示含下标(i₁,i₂，···,i_s)中一个或多个元素；

表示公式(17)中任意变量

的任意阶的多项式的平方。

集合

如下所示：

式中，α为集合

表示的元素。

6)确定最大全局灵敏度指标

实施例2：

一种基于全局灵敏度分析的天然气网络储气配置方法，主要包括以下步骤：

2)基于天然气网络中节点气负荷的随机性和秩相关性，建立服从正态分布的气负荷累积概率模型F(x)，并对气负荷A进行抽样，得到节点气负荷样本矩阵X。

5)确定最大全局灵敏度指标

式中，i_max为节点气负荷变量编号；在负荷变量i_max所在节点设置储气装置。

实施例3：

一种基于全局灵敏度分析的天然气网络储气配置方法，主要步骤见实施例2，其中，抽样得到节点气负荷样本矩阵X的主要步骤如下：

网络中与任一节点a直接相连的所有节点构成的集合记为

2)建立服从正态分布的气负荷累积概率模型F(x)，即：

3)基于拉丁超立方抽样法和高斯Copula函数，对气负荷概率模型F(x)进行抽样，主要步骤如下：

当i＝j时，ρ_sij＝1。当i≠j时，0≤ρ_sij≤1。

3.4)对矩阵C_X进行Cholesky分解，得到下三角矩阵L：

C_X＝LL^T。 (6)

X'＝LX_ξ。 (7)

x_im＝F_i ^-1(Φ(x'_im))。 (8)

实施例4：

一种基于全局灵敏度分析的天然气网络储气配置方法，主要步骤见实施例2，建立天然气稳态能流模型，将每一组气负荷样本作为已知量代入该模型并采用牛顿法求解的主要步骤如下：

1)选取矩阵X的第m列作为负荷变量的第m组样本作为天然气稳态能流模型的已知量，采用牛顿法求解天然气能流模型，得到输出变量向量其中y_m表示气负荷的第m组样本对应的输出变量；输出变量y为任意节点压力或者管道流量；m＝1， 2，…N_ED；采用牛顿法求解天然气网络能流模型的最大迭代次数为 T_max，收敛判据为λ。

2)令某一管道两端节点编号为a和

式中，F_a为节点a的注入流量；

3)计算节点注入天然气流量的不平衡量

将天然气各节点压力的初值代入管道流量的计算公式，得到天然气网络中的管道流量以及各节点的注入流量。此时经过式(2)得到的节点a注入流量与实际值之间存在差异，采用ΔF_a表示节点a 注入流量的不平衡量：

式中，F_a由节点a的气负荷决定，即F_a＝-x_a。

采用列向量ΔF统一表示所有节点注入流量的不平衡量：

ΔF＝[ΔF₁ ΔF₂ … ΔF_N]^T (4)

4)建立修正方程中的雅克比矩阵J

ΔF＝JΔp； (5)

式中，Δp＝[Δp₁,Δp₂,···,Δp_N]^T为状态变量p的修正量向量；J 为雅克比矩阵。具体表示为：

式中，

表示偏微分。

5)计算管道节点压力修正量

基于式(4)和式(6)的计算结果，可以计算所有节点压力的修正量：

Δp＝-J^-1ΔF； (7)

根据所得修正结果可以对节点压力进行更新，即：

p^(k+1)＝p^(k)+Δp^(k)； (8)

式中，k表示迭代次数。

6)重复步骤4)和5)，在每次迭代过程中计算雅克比矩阵J和节点压力修正量△p。并根据以下收敛条件判断是否结束计算：

a.当修正量△p满足max{|△p|}<λ，且迭代次数k≤T_max时，迭代结束，并输出目标变量y；

c.当max{|△p|}≥λ，且k＞T_max；取k＝k+1，返回步骤4.4)继续迭代计算。

实施例5：

一种基于全局灵敏度分析的天然气网络储气配置方法，主要步骤见实施例2，其中，计算输出变量y^ED关于节点气负荷的全局灵敏度指标S_i的主要步骤如下：

1)设定Hermite多项式作为输出变量样本y^ED的最优正交多项式基函数，建立输出变量样本y^ED关于输入变量X的多项式混沌展开形式，即：

为变量

的s阶正交多项式。

为多项式展开系数。a＝[a₀,a₁,···,a_p-1]^T。H(ξ)＝[1,H₁(ξ),···, H_p(ξ)]^T。多项式混沌展开形式(公式(1))的最高展开阶数为p。

公式(1)的展开项数M如下所示：

式中，符号！表示阶乘。

变量ξ中所有变量样本相应的各阶多项式矩阵如下所示：

式中：ξ⁽ⁱ⁾表示变量ξ的第i组样本。

2)计算公式(1)的稀疏展开系数，主要步骤如下：

2.1)设置初始参数，主要包括输出变量向量y^ED、正交多项式矩阵H(ξ)、最大允许误差ε、初始迭代次数k＝0、展开系数a(k)＝0、残差r_(k)＝y^ED-Ha_(k)和项数索引集

2.2)确定出正交多项式矩阵H中与残差r_(k)最相关的列C_(k)＝r_(k) ^TH。确定行向量C_(k)中数值最大的项c_max，元素c_max所在列数记为j_(k)。

2.3)令k＝k+1，将j_(k)添加至集合索引集Λ_(k)＝Λ_(k-1)∪j_(k-1)。

2.4)基于最小二乘法，求解展开系数a_(k)＝arg min_a||y^ED-H_(k)a||₂。其中，H_(k)为矩阵H中索引集Λ_(k)对应列数的向量构成的矩阵。

2.5)更新残差r_(k)＝y^ED-Ha(k)，判断残差r_(k)是否满足||r_(k)||₂>ε。若是，则返回步骤2.2)，若否，则计算结束，并输出优化结果a_(k)。其中，||·||₂表示二范数。

3)计算输出变量y的均值E[·]和方差D[·]，即：

4)结合Sobol’计算公式和公式(1)，计算输出变量关于输入变量

的全局灵敏度，即：

式中：表示变量

的全局灵敏度指标；下标集合

表示含下标(i₁,i₂，···,i_s)中一个或多个元素；

表示公式(17)中任意变量

的任意阶的多项式的平方。

集合

如下所示：

实施例5：

一种验证基于全局灵敏度分析的天然气网络储气配置方法的实验，主要步骤如下：

1)搭建具有11节点的天然气网络。

2)输入基础数据及初始化

2.1)输入基础数据

首先输入天然气网络的基本参数，包括：气源参数、管道参数、气负荷及其秩相关系数、平衡节点压力初始值。

2.2)参数初始化

设天然气网络中有且仅有一个平衡节点，非平衡节点均接有气负荷，总数目N＝10，表示N个节点所在集合，网络中与任一节点a直接相连的所有节点构成的集合记为

设在抽样步骤中，基于拉丁超立方的抽样次数N_ED＝1000，展开系数求解步骤中，最大允许误差ε＝10^-5。

设在采用牛顿法求解天然气系统稳态能流模型步骤中，最大迭代次数为T_max＝100，收敛判据为λ＝10^-6。

3)基于拉丁超立方和高斯Copula函数的概率变量抽样

3.1)建立气负荷概率模型

假设所有气负荷服从正态分布。以节点2气负荷为例，其负荷基准值为219m³/h，标准差取其大小的5％，其累积概率分布函数为：

3.2)气负荷概率模型抽样

定义天然气网络中气负荷变量：X＝[x₁,x₂,···,x₁₀]^T，采用秩相关系数描述气负荷之间的相关性，已知秩相关系数矩阵C_S为：

定义独立标准正态变量ξ＝[ξ₁,ξ₂,···,ξ₁₀]^T，利用拉丁超立方抽样方法产生关于变量ξ，规模为N×N_ED的概率模型样本矩阵X_ξ：

将气负荷的秩相关系数转化为等效的线性相关系数。根据推导，可以得到气负荷的等效线性相关系数矩阵C_X：

对矩阵C_X进行Cholesky分解，求得下三角矩阵L：

根据已有独立标准正态变量样本矩阵X_ξ和下三角矩阵L，得到线性相关的标准正态变量样本矩阵X’：

式中：X’中各元素的线性相关性满足矩阵C_X，等效的秩相关性满足矩阵C_S。

将矩阵X’中的变量样本进行等概率转换，便能得到气负荷样本矩阵X：

X的各列即表示气负荷的各组样本，以第一列X₁为例，给出气负荷第一组样本值：

4)采用牛顿法求解天然气稳态能流模型

4.1)建立管道流量方程

根据输入数据中已有管道压力的初始值，建立相应的管道流量方程，以管道1为例，得到管道流量F_1-2为：

4.2)计算节点注入天然气流量的不平衡量

根据节点负荷取值，获得节点流量不平衡量向量：

4.3)计算雅克比矩阵J

雅克比矩阵J为：

4.4)计算修正量

根据不平衡量向量△F和雅克比矩阵J，可以计算各状态变量(节点压力)的修正量，以第一次迭代后的结果为例，各节点压力修正量为：

根据修正结果更新节点压力，以第一次迭代后的结果为例，计算出第一次迭代后的节点压力：

4.5)收敛性判断

由计算结果可知，当迭代4次(k＝4)时，满足收敛条件，即 max{|△p|}＝6.6992×10^-9<λ。此时可以根据状态变量的数据求解任一待求输出变量，以流经管道7-8的流量(简写为F_7-8)为目标输出变量 y，得到y的样本y^ED：

5)构造稀疏多项式混沌展开表达式

5.1)独立标准正态变量样本的多项式矩阵构造

由于ξ中各变量相互独立且服从标准正态分布，取Hermite正交多项式作为各变量对应的最优正交多项式基函数。单变量Hermite 正交多项式的前三阶展开式分别为：

建立输出变量y关于输入变量的多项式混沌展开表达式。取最高展开阶数p＝2，可将展开形式可以改写为：

展开项数M为：

基于上述内容，结合建立X_ξ中各变量样本相应的各阶多项式，并以矩阵形式H(ξ)表示：

式中：ξ⁽ⁱ⁾表示变量ξ的第i组样本。

5.2)稀疏展开系数求解

采用正交匹配追踪方法可以快速定位并求解混沌多项式展开模型中的非零展开系数，最终得到稀疏多项式展开系数a：

6)基于稀疏多项式混沌展开的全局灵敏度计算

通过上述步骤，可以得到输出变量关于输入变量的稀疏多项式混沌展开方程，基于该方程，可以计算输出变量的均值和方差如下：

结合Sobol’计算公式和所得稀疏多项式混沌展开表达式，对于任意输入变量x_i，求解各输出变量关于它们的全局灵敏度指标：

7)基于全局灵敏度分析的储气配置

选择式(22)中最大的灵敏度指标，有S_imax＝0.5203，且有i_max＝8，这说明第八个输入不确定变量，即负荷x₈的不确定性对于输出变量的影响最大。实际上在天然气网络中，由于平衡节点编号为1，若将负荷变量x按照所接节点编号由小到大排列，可知x₈实际接入节点为9。综上所述，基于全局灵敏度指标可知，在节点9设置储气装置的效果理论上为最优。

8)实验效果

以附图3所示的11节点天然气网络为仿真对象，视F_7-8为关键输出变量，得到其关于所有气负荷的全局灵敏度指标。设计以下仿真算例以验证所提方法对储气配置具有指导作用的有效性：1.不配置任何储气；2.以全局灵敏度指标为参考，将储气配置在全局灵敏度指标最大的气负荷所在节点；3.不参考全局灵敏度指标，分别在管道首节点(节点7)配置储气；4.不参考全局灵敏度指标，分别在管道末节点(节点8)配置储气。分别在上述四种情况下考虑负荷的不确定性，并计算天然气网络的概率能流分布。假设储气提供的气流量保持为15m³/h不变。

F_7-8关于所有气负荷的全局灵敏度指标如附录图4所示，采用 F_i表示节点i所接气负荷；四种场景下F_7-8计算结果如附录图5所示。

由附录图4可见，当F_7-8是天然气网络中的关键状态变量时，节点9的气负荷，即p₉的全局灵敏度指标最大，这表明在所有不确定变量中，F_7-8受到p₉的影响最大。为提高系统运行安全稳定性，考虑在天然气网络中配置储气。基于全局灵敏度指标可知，在节点9 处设置储气最为有利。由附录图5可见，当储气点设置在节点9时， F_7-8的概率密度曲线向右偏移程度最大，并且曲线最为陡峭，表明储气的配置使得F_7-8均值增加，且其波动程度明显减小，这意味着在此处配置储气可以有效增强F_7-8的稳定性。相比之下，在节点7处设置储气虽然也能平抑F_7-8的波动性，但由于p₇的全局灵敏度指标小于p₉，因此在平抑效果上有所降低。而p₈的全局灵敏度指标在三者中最小，因此在节点8配置储气对于F_7-8的影响也是最小，观察附录图5可知，在节点8配置储气对于F_7-8的波动程度基本没有影响，与理论分析一致。上述仿真结果验证了所提方法的正确性和有效性。

Claims

1.一种基于全局灵敏度分析的天然气网络储气配置方法，其特征在于，主要包括以下步骤：

1)获取所述天然气网络参数，并建立天然气网络的稳态能流模型；

3)将负荷的每组样本作为已知量代入天然气网络的稳态能流模型，采用牛顿法求解所述稳态能流模型，得到输出变量y的样本向量y^ED；

4)基于气负荷样本矩阵X和输出变量y^ED，采用稀疏多项式混沌展开方法求解输出变量y关于各节点气负荷的全局灵敏度指标S_i；

5)确定最大全局灵敏度指标

2.根据权利要求1所述的一种基于全局灵敏度分析的天然气网络储气配置方法，其特征在于：所述天然气网络参数主要包括气源参数、管道参数、节点气负荷A＝[x₁,x₂,···,x_N]^T、气负荷的秩相关系数和平衡节点压力初始值。

3.根据权利要求1所述的一种基于全局灵敏度分析的天然气网络储气配置方法，其特征在于，抽样得到节点气负荷样本矩阵X的主要步骤如下：

1)设定天然气网络中有1个平衡节点和N-1个非平衡节点；平衡节点编号记为1；非平衡节点均接有气负荷；天然气网络节点所在位置集合记为网络中与任一节点a直接相连的所有节点构成的集合为

设定抽样次数为N_ED；设定N维独立标准正态变量ξ＝[ξ₁,ξ₂,···,ξ_N]^T；

2)建立服从正态分布的气负荷概率模型，其累积概率函数F(x)为：

式中，x表示气负荷；X为气负荷某一实际取值；μ_x为气负荷的均值；气负荷的均值μ_x等于气负荷基准值；σ_x为气负荷的标准差；

式中，ρ_sij为节点气负荷x_i和节点气负荷x_j的秩相关系数，

当i＝j时，ρ_sij＝1；当i≠j时，0≤ρ_sij≤1；

3.2)基于拉丁超立方抽样法，建立变量ξ的概率模型样本矩阵X_ξ，即：

式中，ρ_ij为在标准正态空间中节点气负荷x_i和节点气负荷x_j的等效线性相关系数；

3.4)对矩阵C_X进行Cholesky分解，得到下三角矩阵L：

C_X＝LL^T； (6)

X'＝LX_ξ； (7)

式中，矩阵X'中各元素的线性相关性满足矩阵C_X；矩阵X'中各元素的等效秩相关性满足矩阵C_S；

3.6)将矩阵X'中的变量样本进行等概率转换，得到气负荷样本矩阵X；其中，矩阵X中任一元素x_im如下所示：

x_im＝F_i ^-1(Φ(x'_im))； (8)

式中：F_i ^-1(·)表示变量x_i的累积分布函数的反函数；Ф(·)表示标准正态分布变量的累积分布函数；x'_im为矩阵X'中的元素，表示第i个变量的第m个样本。

4.根据权利要求1所述的一种基于全局灵敏度分析的天然气网络储气配置方法，其特征在于：建立天然气稳态能流模型，将负荷的每组样本作为已知量代入天然气网络的稳态能流模型，采用牛顿法求解所述稳态能流模型，得到输出变量y的样本向量y^ED的主要步骤如下：

1)选取矩阵X的第m列作为负荷变量的第m组样本作为天然气稳态能流模型的已知量，采用牛顿法求解天然气能流模型，得到输出变量向量y^ED＝[y₁,y₂,···,yN_ED]^T；其中y_m表示气负荷的第m组样本对应的输出变量；输出变量y为任意节点压力或者管道流量；m＝1，2，…N_ED；采用牛顿法求解天然气网络能流模型的最大迭代次数为T_max，收敛判据为λ；

2)令某一管道两端节点编号为a和b，建立流经管道的气流量方程：

式中，F_a-b为由a流向b的天然气流量；p_a和p_b分别为节点a和b的压力；s_a-b为方向变量，当p_a>p_b，s_a-b＝1，否则s_a-b＝-1；D为管道a-管道b的直径；f为管道的摩擦系数；L为管道长度；S_G为天然气的相对密度；a≠b；

式中，F_a为节点a的注入流量；

3)计算节点注入天然气流量的不平衡量；

将天然气各节点压力的初值代入管道流量的计算公式，得到天然气网络中的管道流量以及各节点的注入流量；

基于公式(10)计算得到的节点a注入流量与实际值之间存在的差异，计算节点a注入流量的不平衡量△F_a：

式中，F_a由节点a的气负荷决定，即F_a＝-x_a；

采用列向量△F统一表示所有节点注入流量的不平衡量，即：

ΔF＝[ΔF₁ ΔF₂ …ΔF_N]^T (12)

4)建立修正方程中的雅克比矩阵J

结合牛顿法，以节点压力p＝[p₁,p₂,···,p_N]^T作为状态变量，建立天然气网络中关于状态变量的修正方程组△F：

ΔF＝JΔp； (13)

式中，△p＝[△p₁,△p₂,···,△p_N]^T为状态变量p的修正量向量；J为雅克比矩阵。具体表示为：

式中，

表示偏微分。

5)计算管道节点压力修正量

基于式(12)和式(14)的计算结果，计算所有节点压力的修正量Δp：

Δp＝-J^-1ΔF； (15)

根据所得修正结果对节点压力进行更新，即：

p^(k+1)＝p^(k)+Δp^(k)； (16)

式中，k表示迭代次数。

6)重复步骤4)和5)，在每次迭代过程中计算雅克比矩阵J和节点压力修正量△p；

根据以下收敛条件判断是否结束计算：

I)当修正量△p满足max{|△p|}<λ，且迭代次数k≤T_max时，迭代结束，并输出目标变量y；

II)当max{|△p|}<λ，且k＞T_max，停止迭代，输出“能流计算结果不收敛”，结束计算；

III)当max{|△p|}≥λ，且k＞T_max；取k＝k+1，返回步骤4)继续迭代计算。

5.根据权利要求1所述的一种基于全局灵敏度分析的天然气网络储气配置方法，其特征在于：基于已有的负荷变量矩阵X和输出变量向量y^ED，采用稀疏多项式混沌展开方法计算输出变量y关于各节点气负荷的全局灵敏度指标S_i的主要步骤如下：

1)设置Hermite多项式作为负荷变量的最优正交多项式基函数，建立输出变量样本y^ED关于输入变量的多项式混沌展开形式，即：

式中，f(ξ)表示输出变量样本y^ED关于输入变量A的原始表达式；H_s(ξi₁,ξi₂,···,ξi_s)为变量(ξi₁,ξi₂,···,ξi_s)的s阶正交多项式；a₀,a_i,ai₁i₂,···,ai₁i₂···i_p为多项式展开系数；a＝[a₀,a₁,···,a_p-1]^T；H(ξ)＝[1,H₁(ξ),···,H_p(ξ)]^T；多项式混沌展开形式(17)的最高展开阶数为p；

公式(17)的展开项数M如下所示：

式中，符号！表示阶乘；

变量ξ中所有变量样本相应的各阶多项式矩阵如下所示：

式中，ξ⁽ⁱ⁾表示变量ξ的第i组样本；

2)计算公式(17)的稀疏展开系数，主要步骤如下：

2.2)确定出正交多项式矩阵H中与残差r_(k)最相关的列C_(k)＝r_(k) ^TH；确定行向量C_(k)中数值最大的项c_max，元素c_max所在列数记为j_(k)；

2.3)令k＝k+1，将j_(k)添加至集合索引集Λ_(k)＝Λ_(k-1)∪j_(k-1)；

2.4)基于最小二乘法，求解展开系数a_(k)＝arg min_a||y^ED-H_(k)a||₂；其中，H_(k)为矩阵H中索引集Λ_(k)对应列数的向量构成的矩阵。

2.5)更新残差r_(k)＝y^ED-Ha(k)，判断残差r_(k)是否满足||r_(k)||₂>ε；若是，则返回步骤2.2)，若否，则计算结束，并输出优化结果a_(k)；其中，||·||₂表示二范数

3)计算输出变量y的均值E[·]和方差D[·]，即：

式中，i＝(i₁,i₂,···,i_N)表示多阶下标；|i|＝i₁+i₂+···+i_N；

4)结合Sobol’计算原理以及公式(17)，计算输出变量关于输入变量(ξi₁,ξi₂,···,ξi_s)的全局灵敏度，即：

式中：Si₁i₂···i_s表示变量(ξi₁,ξi₂,···,ξi_s)的全局灵敏度指标；下标集合Ωi₁,i₂,···,i_s表示含下标(i₁,i₂，···,i_s)中一个或多个元素；

表示公式(17)中任意变量(ξi₁,ξi₂,···,ξi_s)的任意阶的多项式的平方；

集合Ωi₁,i₂,···,i_s如下所示：