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

摘要

本发明公开了一种基于全局灵敏度分析的天然气网络储气配置方法，主要步骤为：2)建立服从正态分布的气负荷概率模型F(x)，并抽样得到节点气负荷样本矩阵X。2)将负荷的每组样本代入天然气网络的稳态能流计算，得到输出变量yED，并计算输出变量yED关于节点气负荷的全局灵敏度指标Si。3)确定最大全局灵敏度指标式中，imax为节点气负荷变量编号；在负荷变量imax所在节点设置储气装置。基于稀疏多项式混沌展开方法，本发明可以快速建立输出变量关于输入不确定变量的表达形式，从而有效提高全局灵敏度分析效率。

权利要求

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

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

T

所述天然气网络参数包括气源参数、管道参数、节点气负荷A＝[x1,x2,···,xN] 、气负荷的秩相关系数和平衡节点压力初始值；

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

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

2.1)设定天然气网络中有1个平衡节点和N‑1个非平衡节点；平衡节点编号记为1；非平衡节点均接有气负荷；天然气网络节点所在位置集合记为 网络中与任一节点a直接相连的所有节点构成的集合为 设定抽样次数为NED；设定N维独立标准正态变量ξ＝[ξ1,Tξ2,···,ξN]；

2.2)建立服从正态分布的气负荷概率模型，其累积分布函数F(x)为：式中，x表示气负荷；X为气负荷某一实际取值；μx为气负荷的均值；气负荷的均值μx等于气负荷基准值；σx为气负荷的标准差；

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

2.3.1)建立节点气负荷A＝[x1,x2,···,xN]的秩相关系数矩阵CS，即：式中，ρsij为节点气负荷xi和节点气负荷xj的秩相关系数，i, 当i＝j时，ρsij＝1；

当i≠j时，0≤ρsij≤1；

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

2.3.3)将气负荷的秩相关系数转化为等效线性相关系数，转换公式如下：式中，ρij为在标准正态空间中节点气负荷xi和节点气负荷xj的等效线性相关系数；

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

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

CX＝LL；            (6)

2.3.5)基于独立标准正态变量样本矩阵Xξ和下三角矩阵L，建立线性相关的标准正态变量样本矩阵X'，即：X'＝LXξ；             (7)式中，矩阵X'中各元素的线性相关性满足矩阵CX；矩阵X'中各元素的等效秩相关性满足矩阵CS；

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

xim＝Fi (Φ(x'im))；        (8)‑1

式中：Fi (·)表示变量xi的累积分布函数的反函数；Ф(·)表示标准正态分布变量的累积分布函数；x'im为矩阵X'中的元素，表示第i个变量的第m个样本；

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

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

3.2)令某一管道两端节点编号为a和b，建立流经管道的气流量方程：式中，Fa‑b为由a流向b的天然气流量；pa和pb分别为节点a和b的压力；sa‑b为方向变量，当pa>pb，sa‑b＝1，否则sa‑b＝‑1；D为管道a‑管道b的直径；f为管道的摩擦系数；L为管道长度；SG为天然气的相对密度；a≠b；

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

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

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

基于公式(10)计算得到的节点a注入流量与实际值之间存在的差异，计算节点a注入流量的不平衡量△Fa：式中，Fa由节点a的气负荷决定，即Fa＝‑xa；

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

ΔF＝[ΔF1 ΔF2…ΔFN]            (12)

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

结合牛顿法，以节点压力p＝[p1,p2,···,pN]作为状态变量，建立天然气网络中关于状态变量的修正方程组△F：ΔF＝JΔp；                  (13)T

式中，△p＝[△p1,△p2,···,△pN]为状态变量p的修正量向量；J为雅克比矩阵，具体表示为：式中，表示偏微分；

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

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

Δp＝‑J ΔF；               (15)根据所得修正结果对节点压力进行更新，即：(k+1) (k) (k)

p ＝p +Δp  ；            (16)式中，k表示迭代次数；

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

根据以下收敛条件判断是否结束计算：I)当修正量△p满足max{|△p|}<λ，且迭代次数k≤Tmax时，迭代结束，并输出目标变量y；

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

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

ED

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

ED

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

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

式中，f(ξ)表示输出变量y的样本向量y 关于输入变量A的原始表达式；

为变量 的s阶正交多项式； 为多项式展开系数；a＝T T

[a0,a1,···,ap‑1] ；H(ξ)＝[1,H1(ξ),···,Hp(ξ)] ；多项式混沌展开形式(17)的最高展开阶数为p；

公式(17)的展开项数M如下所示：式中，符号！表示阶乘；

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

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

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

4.2.1)设置初始参数，包括输出变量y的样本向量y 、正交多项式矩阵H(ξ)、最大允许ED误差ε、初始迭代次数k＝0、展开系数a(k)＝0、残差r(k)＝y ‑Ha(k)和项数索引集T

4.2.2)确定出正交多项式矩阵H中与残差r(k)最相关的列C(k)＝r(k)H；确定行向量C(k)中数值最大的项cmax，元素cmax所在列数记为j(k)；

4.2.3)令k＝k+1，将j(k)添加至集合索引集Λ(k)＝Λ(k‑1)∪j(k‑1)；

ED

4.2.4)基于最小二乘法，求解展开系数a(k)＝arg mina||y ‑H(k)a||2；其中，H(k)为矩阵H中索引集Λ(k)对应列数的向量构成的矩阵；

ED

4.2.5)更新残差r(k)＝y ‑Ha(k)，判断残差r(k)是否满足||r(k)||2>ε；若是，则返回步骤

4.2.2)，若否，则计算结束，并输出优化结果a(k)；其中，||·||2表示二范数；

4.3)计算输出变量y的均值E[·]和方差D[·]，即：式中，i＝(i1,i2,···,iN)表示多阶下标；|i|＝i1+i2+···+iN；

4.4)结合Sobol’计算原理以及公式(17)，计算输出变量关于输入变量的全局灵敏度，即：

式中： 表示变量 的全局灵敏度指标；下标集合 表示含下标(i1,i2，···,is)中一个或多个元素； 表示公式(17)中任意变量的任意阶的多项式的平方；

集合 如下所示：

5)确定最大全局灵敏度指标 式中，imax为节点气负荷变量编号；在负荷变量imax所在节点设置储气装置。

说明书

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

技术领域

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

背景技术

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

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

发明内容

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

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

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

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

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

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

[0010] 2.1)设定天然气网络中有1个平衡节点和N‑1个非平衡节点。 平衡节点编号记为1。非平衡节点均接有气负荷。天然气网络节点所 在位置集合记为 网络中与任一节点a直接相连的所有节点构成 的集合记为 设定抽样次数为NED。设定N维独立标准正态变量 T
ξ＝[ξ1,ξ2,···,ξN]。

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

[0012]

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

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

[0015] 2.3.1)建立节点气负荷A＝[x1,x2,···,xN]T的秩相关系数矩阵CS， 即：

[0016]

[0017] 式中，ρsij为节点气负荷xi和节点气负荷xj的秩相关系数，i,  当i＝j时，ρsij＝1。当i≠j时，0≤ρsij≤1。

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

[0019]

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

[0021]

[0022] 式中，ρij为在标准正态空间中节点气负荷xi和节点气负荷xj的 等效线性相关系数。

[0023] 基于等效线性相关系数ρij，建立气负荷的等效线性相关系数矩 阵CX，即：

[0024]

[0025] 2.3.4)对矩阵CX进行Cholesky分解，得到下三角矩阵L：

[0026] CX＝LLT。    (6)

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

[0028] X'＝LXξ。    (7)

[0029] 式中，矩阵X'中各元素的线性相关性满足矩阵CX。矩阵X'中 各元素的等效秩相关性满足矩阵CS。

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

[0031] xim＝Fi‑1(Φ(x'im))。    (8)

[0032] 式中：Fi‑1(·)表示变量xi的累积分布函数的反函数。Ф(·)表示标 准正态分布变量的累积分布函数。x'im为矩阵X'中的元素，表示第i 个变量的第m个样本。

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

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

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

[0036] 3.2)令某一管道两端节点编号为a和 对于本 文所研究的低压配气网络(压力小于等于75mbar)，建立流经管道 的气流量方程：

[0037]

[0038] 式中，Fa‑b为由a流向b的天然气流量；pa和pb分别为节点a 和b的压力；sa‑b为方向变量，当pa>pb，sa‑b＝1，否则sa‑b＝‑1；D为 管道a‑b的直径；f为管道的摩擦系数；L为管道长度；SG为天然气 的相对密度，标准工况下为0.5548。

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

[0040]

[0041] 式中，Fa为节点a的注入流量；

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

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

[0044]

[0045] 式中，Fa由节点a的气负荷决定，即Fa＝‑xa。

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

[0047] ΔF＝[ΔF1 ΔF2 … ΔFN]T    (12)

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

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

[0050] ΔF＝JΔp；    (13)

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

[0052]

[0053] 式中，表示偏微分。

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

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

[0056] Δp＝‑J‑1ΔF；    (15)

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

[0058] p(k+1)＝p(k)+Δp(k)；    (16)

[0059] 式中，k表示迭代次数。

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

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

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

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

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

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

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

[0067]ED

[0068] 式中，f(ξ)表示输出变量样本y 关于输入变量X的原始表达式。 为变量 的s阶正交多项式。 为多项式展开系数。a＝
T T
[a0,a1,···,ap‑1] 。H(ξ)＝[1,H1(ξ),···, Hp(ξ)]。多项式混沌展开形式(17)的最高展开阶数为p。

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

[0070]

[0071] 式中，符号！表示阶乘。

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

[0073]

[0074] 式中：ξ(i)表示变量ξ的第i组样本。

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

[0076] 4.2.1)设置初始参数，主要包括输出变量向量yED、正交多项式 矩阵H(ξ)、最大允ED许误差ε、初始迭代次数k＝0、展开系数a(k)＝0、 残差r(k)＝y ‑Ha(k)和项数索引集[0077] 4.2.2)确定出正交多项式矩阵H中与残差r(k)最相关的列C(k)＝r(k)TH。确定行向量C(k)中数值最大的项cmax，元素cmax所在列数记 为j(k)。

[0078] 4.2.3)令k＝k+1，将j(k)添加至集合索引集Λ(k)＝Λ(k‑1)∪j(k‑1)。

[0079] 4.2.4)基于最小二乘法，求解展开系数a(k)＝arg mina||yED‑H(k)a||2。 其中，H(k)为矩阵H中索引集Λ(k)对应列数的向量构成的矩阵。

[0080] 4.2.5)更新残差r(k)＝yED‑Ha(k)，判断残差r(k)是否满足||r(k)||2>ε。 若是，则返回步骤4.2.2)，若否，则计算结束，并输出优化结果a(k)。 其中，||·||2表示二范数[0081] 4.3)计算输出变量y的均值E[·]和方差D[·]，即：

[0082]

[0083]

[0084] 式中，i＝(i1,i2,···,iN)表示多阶下标。|i|＝i1+i2+···+iN。

[0085] 4.4)结合Sobol’计算公式和公式(17)，计算输出变量关于输入 变量(ξi1,ξi2,···,ξis)的全局灵敏度，即：

[0086]

[0087] 式中： 表示变量 的全局灵敏度指标；下标 集合 表示含下标(i1,i2，···,is)中一个或多个元素； 表示公式(17)中任意变量
的任意阶的多项式的平方。

[0088] 集合 如下所示：

[0089]

[0090] 5)确定最大全局灵敏度指标 式中，imax为节点气负荷变量编号。在负荷变量imax所在节点设置储气装置。

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

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

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

附图说明

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

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

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

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

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

具体实施方式

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

[0100] 实施例1：

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

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

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

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

[0105] 概率抽样的主要步骤如下

[0106] 1)设定天然气网络中有1个平衡节点和N‑1个非平衡节点。平 衡节点编号记为1。非平衡节点均接有气负荷。天然气网络节点所在 位置集合记为 网络中与任一节点a直接相连的所有节点构成的 集合记为 设定抽样次数为NED。设定N维独立标准正态变量ξT
＝[ξ1, ξ2,···,ξN]。

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

[0108]

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

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

[0111] 3.1)建立节点气负荷A＝[x1,x2,···,xN]T的秩相关系数矩阵CS， 即：

[0112]

[0113] 式中，ρsij为节点气负荷xi和节点气负荷xj的秩相关系数，i,  当i＝j时，ρsij＝1。当i≠j时，0≤ρsij≤1。

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

[0115]

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

[0117]

[0118] 式中，ρij为在标准正态空间中节点气负荷xi和节点气负荷xj的 等效线性相关系数。

[0119] 基于等效线性相关系数ρij，建立气负荷的等效线性相关系数矩 阵CX，即：

[0120]

[0121] 3.4)对矩阵CX进行Cholesky分解，得到下三角矩阵L：

[0122] CX＝LLT。    (6)

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

[0124] X'＝LXξ。    (7)

[0125] 式中，矩阵X'中各元素的线性相关性满足矩阵CX。矩阵X'中 各元素的等效秩相关性满足矩阵CS。

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

[0127] xim＝Fi‑1(Φ(x'im))。    (8)

[0128] 式中：Fi‑1(·)表示变量xi的累积分布函数的反函数。Ф(·)表示标 准正态分布变量的累积分布函数。x'im为矩阵X'中的元素，表示第i 个变量的第m个样本。

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

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

[0131] 4.2)令某一管道两端节点编号为a和 对于本 文所研究的低压配气网络(压力小于等于75mbar)，建立流经管道 的气流量方程：

[0132]

[0133] 式中，Fa‑b为由a流向b的天然气流量；pa和pb分别为节点a 和b的压力；sa‑b为方向变量，当pa>pb，sa‑b＝1，否则sa‑b＝‑1；D为 管道a‑b的直径；f为管道的摩擦系数；L为管道长度；SG为天然气 的相对密度，标准工况下为0.5548。

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

[0135]

[0136] 式中，Fa为节点a的注入流量；

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

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

[0139]

[0140] 式中，Fa由节点a的气负荷决定，即Fa＝‑xa。xa为节点a的气 负荷。

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

[0142] ΔF＝[ΔF1 ΔF2 … ΔFN]T    (12)

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

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

[0145] ΔF＝JΔp；    (13)

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

[0147]

[0148] 式中，表示偏微分。

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

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

[0151] Δp＝‑J‑1ΔF；    (15)

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

[0153] p(k+1)＝p(k)+Δp(k)；    (16)

[0154] 式中，k表示迭代次数。

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

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

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

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

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

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

[0161]

[0162] 式中，f(ξ)表示输出变量样本yED关于输入变量X的原始表达式。 为变量 的s阶正交多项式。 为多项式展开系数。a＝
T T
[a0,a1,···,ap‑1] 。H(ξ)＝[1,H1(ξ),···, Hp(ξ)]。多项式混沌展开形式(17)的最高展开阶数为p。

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

[0164]

[0165] 式中，符号！表示阶乘。

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

[0167]

[0168] 式中：ξ(i)表示变量ξ的第i组样本。

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

[0170] 5.2.1)设置初始参数，主要包括输出变量向量yED、正交多项式 矩阵H(ξ)、最大允ED许误差ε、初始迭代次数k＝0、展开系数a(k)＝0、 残差r(k)＝y ‑Ha(k)和项数索引集[0171] 5.2.2)确定出正交多项式矩阵H中与残差r(k)最相关的列C(k)＝r(k)TH。确定行向量C(k)中数值最大的项cmax，元素cmax所在列数记 为j(k)。

[0172] 5.2.3)令k＝k+1，将j(k)添加至集合索引集Λ(k)＝Λ(k‑1)∪j(k‑1)。Λ(k‑1)为第k‑1次迭代的项数索引集，j(k‑1)为第k‑1次迭代的元素cmax所在 列数。

[0173] 5.2.4)基于最小二乘法，求解展开系数a(k)＝arg mina||yED‑H(k)a||2。 其中，H(k)为矩阵H中索引集Λ(k)对应列数的向量构成的矩阵。min 表示取最小值。

[0174] 5.2.5)更新残差r(k)＝yED‑Ha(k)，判断残差r(k)是否满足||r(k)||2>ε。 若是，则返回步骤5.2.2)，若否，则计算结束，并输出优化结果a(k)。 其中，||·||2表示二范数。

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

[0176]

[0177]

[0178] 式中，i＝(i1,i2,···,iN)表示多阶下标。|i|＝i1+i2+···+iN。

[0179] 5.4)结合Sobol’计算公式和公式(17)，计算输出变量关于输入  变量的全局灵敏度，即：

[0180]

[0181] 式中： 表示变量 的全局灵敏度指标；下标 集合表示含下标(i1,i2，···,is)中一个或多个元素； 表示公式(17)中任意变量
的任意阶的多项式的平方。

[0182] 集合 如下所示：

[0183]

[0184] 式中，α为集合 表示的元素。

[0185] 6)确定最大全局灵敏度指标 式中，imax为节点气负荷变量编号。在负荷变量imax所在节点设置储气装置。

[0186] 实施例2：

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

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

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

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

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

[0192] 5)确定最大全局灵敏度指标 式中，imax为节点气负荷变量编号；在负荷变量imax所在节点设置储气装置。

[0193] 实施例3：

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

[0195] 1)设定天然气网络中有1个平衡节点和N‑1个非平衡节点。平 衡节点编号记为1。非平衡节点均接有气负荷。天然气网络节点所在 位置集合记为 网络中与任一节点a直接相连的所有节点构成的 集合记为 设定抽样次数为NED。设定N维独立标准正态变量ξT
＝[ξ1, ξ2,···,ξN]。

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

[0197]

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

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

[0200] 3.1)建立节点气负荷A＝[x1,x2,···,xN]T的秩相关系数矩阵CS， 即：

[0201]

[0202] 式中，ρsij为节点气负荷xi和节点气负荷xj的秩相关系数，i,  当i＝j时，ρsij＝1。当i≠j时，0≤ρsij≤1。

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

[0204]

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

[0206]

[0207] 式中，ρij为在标准正态空间中节点气负荷xi和节点气负荷xj的 等效线性相关系数。

[0208] 基于等效线性相关系数ρij，建立气负荷的等效线性相关系数矩 阵CX，即：

[0209]

[0210] 3.4)对矩阵CX进行Cholesky分解，得到下三角矩阵L：

[0211] CX＝LLT。    (6)

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

[0213] X'＝LXξ。    (7)

[0214] 式中，矩阵X'中各元素的线性相关性满足矩阵CX。矩阵X'中 各元素的等效秩相关性满足矩阵CS。

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

[0216] xim＝Fi‑1(Φ(x'im))。    (8)

[0217] 式中：Fi‑1(·)表示变量xi的累积分布函数的反函数。Ф(·)表示标 准正态分布变量的累积分布函数。x'im为矩阵X'中的元素，表示第i 个变量的第m个样本。

[0218] 实施例4：

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

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

[0221] 2)令某一管道两端节点编号为a和 对于本文 所研究的低压配气网络(压力小于等于75mbar)，建立流经管道的 气流量方程：

[0222]

[0223] 式中，Fa‑b为由a流向b的天然气流量；pa和pb分别为节点a 和b的压力；sa‑b为方向变量，当pa>pb，sa‑b＝1，否则sa‑b＝‑1；D为 管道a‑b的直径；f为管道的摩擦系数；L为管道长度；SG为天然气 的相对密度，标准工况下为0.5548。

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

[0225]

[0226] 式中，Fa为节点a的注入流量；

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

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

[0229]

[0230] 式中，Fa由节点a的气负荷决定，即Fa＝‑xa。

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

[0232] ΔF＝[ΔF1 ΔF2 … ΔFN]T    (4)

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

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

[0235] ΔF＝JΔp；    (5)

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

[0237]

[0238] 式中，表示偏微分。

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

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

[0241] Δp＝‑J‑1ΔF；    (7)

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

[0243] p(k+1)＝p(k)+Δp(k)；    (8)

[0244] 式中，k表示迭代次数。

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

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

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

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

[0249] 实施例5：

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

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

[0252]ED

[0253] 式中，f(ξ)表示输出变量样本y 关于输入变量X的原始表达式。 为变量 的s阶正交多项式。 为多项式展开系数。a＝
T T
[a0,a1,···,ap‑1]。H(ξ)＝[1,H1(ξ),···, Hp(ξ)]。多项式混沌展开形式(公式(1))的最高展开阶数为p。

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

[0255]

[0256] 式中，符号！表示阶乘。

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

[0258]

[0259] 式中：ξ(i)表示变量ξ的第i组样本。

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

[0261] 2.1)设置初始参数，主要包括输出变量向量yED、正交多项式 矩阵H(ξ)、最大允许ED误差ε、初始迭代次数k＝0、展开系数a(k)＝0、 残差r(k)＝y ‑Ha(k)和项数索引集[0262] 2.2)确定出正交多项式矩阵H中与残差r(k)最相关的列C(k)＝r(k)TH。确定行向量C(k)中数值最大的项cmax，元素cmax所在列数记 为j(k)。

[0263] 2.3)令k＝k+1，将j(k)添加至集合索引集Λ(k)＝Λ(k‑1)∪j(k‑1)。

[0264] 2.4)基于最小二乘法，求解展开系数a(k)＝arg mina||yED‑H(k)a||2。 其中，H(k)为矩阵H中索引集Λ(k)对应列数的向量构成的矩阵。

[0265] 2.5)更新残差r(k)＝yED‑Ha(k)，判断残差r(k)是否满足||r(k)||2>ε。 若是，则返回步骤2.2)，若否，则计算结束，并输出优化结果a(k)。 其中，||·||2表示二范数。

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

[0267]

[0268]

[0269] 式中，i＝(i1,i2,···,iN)表示多阶下标。|i|＝i1+i2+···+iN。

[0270] 4)结合Sobol’计算公式和公式(1)，计算输出变量关于输入变 量的全局灵敏度，即：

[0271]

[0272] 式中： 表示变量 的全局灵敏度指标；下标 集合 表示含下标(i1,i2，···,is)中一个或多个元素； 表示公式(17)中任意变量
的任意阶的多项式的平方。

[0273] 集合 如下所示：

[0274]

[0275] 实施例5：

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

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

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

[0279] 2.1)输入基础数据

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

[0281] 2.2)参数初始化

[0282] 设天然气网络中有且仅有一个平衡节点，非平衡节点均接有气 负荷，总数目N＝10， 表示N个节点所在集合，网络中与任一节 点a直接相连的所有节点构成的集合记为[0283] 设在抽样步骤中，基于拉丁超立方的抽样次数NED＝1000，展开 系数求解步骤中，‑5
最大允许误差ε＝10 。

[0284] 设在采用牛顿法求解天然气系统稳态能流模型步骤中，最大迭 代次数为Tmax＝‑6100，收敛判据为λ＝10 。

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

[0286] 3.1)建立气负荷概率模型

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

[0288]

[0289] 3.2)气负荷概率模型抽样

[0290] 定义天然气网络中气负荷变量：X＝[x1,x2,···,x10]T，采用秩相关 系数描述气负荷之间的相关性，已知秩相关系数矩阵CS为：

[0291]

[0292] 定义独立标准正态变量ξ＝[ξ1,ξ2,···,ξ10]T，利用拉丁超立方抽样 方法产生关于变量ξ，规模为N×NED的概率模型样本矩阵Xξ：

[0293]

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

[0295]

[0296] 对矩阵CX进行Cholesky分解，求得下三角矩阵L：

[0297]

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

[0299]

[0300] 式中：X’中各元素的线性相关性满足矩阵CX，等效的秩相关性 满足矩阵CS。

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

[0302]

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

[0304]

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

[0306] 4.1)建立管道流量方程

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

[0308]

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

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

[0311]

[0312] 4.3)计算雅克比矩阵J

[0313] 雅克比矩阵J为：

[0314]

[0315] 4.4)计算修正量

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

[0317]

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

[0319]

[0320] 4.5)收敛性判断

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

[0322]

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

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

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

[0326]

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

[0328]

[0329] 展开项数M为：

[0330]

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

[0332]

[0333] 式中：ξ(i)表示变量ξ的第i组样本。

[0334] 5.2)稀疏展开系数求解

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

[0336]

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

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

[0339]

[0340]

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

[0342]

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

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

[0345] 8)实验效果

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

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

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