CN105046369A - 一种基于能源中心的电气混联系统建模和优化调度方法 - Google Patents

Abstract

本发明公布了一种基于能源中心的电气混联系统建模和优化调度方法，本发明首先建立了电力网络、天然气网络和能源中心模型，电力网络和天然气网络通过能源中心耦合形成电气混联系统。然后以总的能源成本为目标函数、考虑了各种约束条件建立电气混联系统优化调度数学模型。用原-对偶内点法求解，求解过程中先后引入松弛变量和障碍参数将模型变为只有等式约束的模型，再引入拉格朗日乘子得到拉格朗日函数，之后用牛顿法求解其KKT条件形成的非线性方程组。构造的算例仿真结果表明本发明对电气混联系统的优化效果优于单独优化效果。

Description

一种基于能源中心的电气混联系统建模和优化调度方法

技术领域

本发明涉及一种基于能源中心的电气混联系统建模和优化调度方法，属于多能源综合控制利用领域。

技术背景

化石能源的广泛使用导致环境问题日益突出以及化石能源逐渐枯竭，传统的经济和社会发展模式难以持续，亟需建立更加高效、环保和可持续型的能源利用模式，能源中心(energyhub,EH)的提出使建立新的能源利用模式成为可能。能源中心被定义为由能源转化设备和储能设备构成、能够实现多种能源相互转化和存储的虚拟实体。能源中心作为多种能源形式的中间媒介，促使多种能源形式和能源系统的联系更加紧密，并考虑了经济性影响。能源中心输入端连接多种能源的供应，输出端连接多种能源的消费，不同能源形式在其间相互转化，这就为多种能源形式的优化调度提供了可能性。

相比于其他一次能源，天然气对环境的影响较小、储量丰富且易于存储；燃气轮机和联合循环电厂的日渐发展，使得天然气发电投资更少、效率更高、建设周期缩短且调控更加灵活方便；此外，“页岩气革命”的不断深化，将造成天然气价格显著下降。可以预见，天然气发电的比例将逐步提高，天然气也会凭借其诸多优势在未来能源利用模式中占据重要地位。天然气网络与电力网络的联系越来越密切，可以说未来的能源利用模式将是天然气网络与电力网络高度耦合的产物。

前人对电力网络和天然气网络的独立网络都有足够的研究，如最优潮流(optimalpowerflow,OPF)等经典问题，但少有考虑电力网络与天然气网络的综合优化调度问题。近年来，由于天然气在能源利用中的地位越来越突出，电力网络与天然气网络的综合研究引起了国内外学者的广泛关注。电气混联系统的研究主要分为协调规划和协调运行两个阶段，问题主要集中在不确定因素的处理、电负荷与电负荷的相互转化与影响、监管部门的协调规划与管理以及混联系统的动态行为与可靠性等方面。

电气混联系统优化调度本质上是非线性规划问题。本发明以总的能源成本作为目标函数，考虑了电力网络中有功、无功功率平衡约束，平衡节点相角平衡约束，发电机有功、无功出力约束，电压约束；天然气网络中流量平衡约束，压缩机状态量平衡约束，气源点输出约束和节点压力约束；能源中心的输入输出平衡约束和调度因子约束。本发明采用原-对偶内点法(primal-dualinteriorpointmethod，PDIPM)求解基于能源中心的电气混联系统优化调度问题，因为PDIPM有着收敛性好，计算效率高，鲁棒性强，初值选取不敏感，没有识别约束集困难等优点。最后对构造的算例进行编程仿真，结果分析验证了本发明方法的优势。

发明内容

发明目的：本发明针对现有技术所需解决的技术问题提供一种基于能源中心的新型能源供应形式即电气混联系统的建模和优化调度。

技术方案：本发明为实现上述目的，采用如下技术方案：

本发明为一种基于能源中心概念电气混联系统建模和优化调度方法，其特征在于所述方法依次按以下步骤实现：

1)获取电力网络的参数信息，包括：输电线路的首端节点和末端节点编号，支路π型等效电路的电阻、电抗，对地并联电导、电纳，变压器变比和阻抗，各节点负荷以及发电机输出有功、无功约束，各节点电压约束；

对于各电力网络中节点i：

$P_i=e_i\underset{j=1}{\overset{nb}{\Σ}}(G_{ij}{e}_j-B_{ij}{f}_j)-f_i\underset{j=1}{\overset{nb}{\Σ}}(G_{ij}{f}_j+B_{ij}{e}_j);$

$Q_i=f_i\underset{j=1}{\overset{nb}{\Σ}}(G_{ij}{e}_j-B_{ij}{f}_j)+e_i\underset{j=1}{\overset{nb}{\Σ}}(G_{ij}{f}_j+B_{ij}{e}_j);$

其中：P_i，Q_i分别为电力网络中节点i有功功率和无功功率；e_i，f_i分别为节点i电压向量的实部和虚部；G_ij，B_ij分别为节点导纳矩阵第i行第j列元素的实部和虚部。

2)获取天然气网络的参数信息，包括：输气管道的首端节点和末端节点编号，管道的长度、内径和传输效率等物理特性，压缩机的状态量，各节点气负荷以及气源点天然气输出约束，各节点压力约束；

引入包括管道流量方程、压缩机流量消耗方程及流量平衡方程；

理想条件下，管道k从节点i到节点j的流量值可以用下面的方程表示：

$f_k=f_{kij}=S_{ij}\×6.4774\frac{T_0}{{\mathrm{\π}}_0}\sqrt{S_{ij}\frac{({\mathrm{\π}}_i^2-{\mathrm{\π}}_j^2)D_k^5}{F_k{\mathrm{GL}}_k{T}_{ka}{Z}_a}};$

其中：f_kij为管道流量值； $S_{ij}=\left\{\begin{array}{cc}+1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j\&GreaterEqual;0\\ -1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j<0\end{array}\right.;$ F_k为管道摩擦系数；D_k为管道内径；G为气体重力系数；L_k为管道长度；π_i为节点i压力值；π_j为节点j压力值；π₀为标准压力值；T₀为标准温度值；T_ka为平均气体温度；Z_a为平均气体压缩系数；

对于高压网络的全紊流形式，流量方程可以进一步简化为：

$f_k=f_{kij}=M_k{S}_{ij}\sqrt{S_{ij}({\mathrm{\&pi;}}_i^2-{\mathrm{\&pi;}}_j^2)};$

其中： $M_k=\mathrm{\&epsiv;}\frac{18.73T_0{D}_k^{8/3}}{{\mathrm{\&pi;}}_0\sqrt{{\mathrm{GL}}_k{T}_{ka}{Z}_a}},$ ε为管道效率；

理想气体条件下，压缩机的能量消耗方程可以表示为：

$H_k=H_{kij}=B_k{f}_{Ck}\&lsqb;{\left(\frac{{\mathrm{\&pi;}}_j}{{\mathrm{\&pi;}}_i}\right)}^{Z_{ki}\left(\frac{\mathrm{\&alpha;}-1}{\mathrm{\&alpha;}}\right)}-1\&rsqb;;$

其中：f_Ck为通过压缩机的气体流量；π_i为气体注入压缩机压力；π_j为气体输出压缩机压力；Z_ki为压缩机入口的气体压缩系数；T_ki为压缩机输出温度；α为热量系数；η_k为压缩机效率；

转化为消耗的流量值：

τ_k＝α_Tk+β_TkH_kij+γ_TkH_k ² _ij；

每一个节点的流量平衡方程可以用下面的矩阵形式表示：

(A+U)f+w-Tτ＝0；

其中：A、U和T分别为支路、压缩机和压缩机消耗与节点有关的方向矩阵，f为支路流量值向量；w为各节点的气体注入向量；τ为各压缩机消耗流量值向量；

3)获取能源中心的参数信息，包括：各能源中心输入输出与电力网络及天然气网络的连接情况，各能源中心负荷以及各能源中心调度因子约束；

多输入多输出能源中心模型可以用矩阵方程描述：

其中：P和L分别为输入输出向量；矩阵C被称为耦合矩阵；

列出能源中心输入输出平衡方程：

$\left[\begin{array}{c}{L}_e\\ L_h\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\&eta;}}_T& {\mathrm{\&nu;\&eta;}}_{GTe}\\ 0& {\mathrm{\&nu;\&eta;}}_{GTh}+(1-\mathrm{\&nu;}){\mathrm{\&eta;}}_F\end{array}\right]\left[\begin{array}{c}{P}_e\\ P_g\end{array}\right];$

其中：η_T表示变压器效率；η_GTe表示燃气轮机气转电的效率；η_GTh表示燃气轮机气转热的效率；η_F表示燃气锅炉的效率；

4)以总的能源成本作为目标函数f(x)，选取控制变量u和状态变量根据各种约束条件列写等式约束和不等式约束，建立电气混联系统非线性规划数学模型并用PDIPM求解；

建立电气混联系统优化调度模型：

minf(x)；

s.t.h(x)＝0；

g_min≤g(x)≤g_max；

其中：h表示等式约束；g表示不等式约束；g_max，g_min分别表示不等式约束的上、下限；

5)引入广义松弛变量l和u，将模型中不等式约束转化为等式约束；

引入广义松弛变量l＝[l₁,…,l_r]^T,u＝[u₁,…,u_r]^T将广义不等式转化为等式：

g(x)+u-g_max＝0；

g(x)-l-g_min＝0；

应满足：u>0,l>0；

6)为了保证目标函数在可行域内取得最小值，对目标函数引入对数障碍参数μ，将模型转化为只含有等式约束的模型；

引入松弛变量和障碍参数将模型转化成只有等式约束的模型：

$\min .f\left(x\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}log\left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}log\left(u_j\right);$

s.t.h(x)＝0；

g(x)+u-g_max＝0；

g(x)-l-g_min＝0；

其中：扰动因子(障碍参数)μ＞0；

引入拉格朗日乘子得到拉格朗日函数：

$\begin{array}{c}Lag=f\left(x\right)-y^{T}h\left(x\right)-z^T\&lsqb;g\left(x\right)-l-g_{\min }\&rsqb;\\ -w^T\&lsqb;g\left(x\right)+u-g_{\max }\&rsqb;-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(u_j\right)\end{array};$

7)用拉格朗日乘子法求解上述等式约束模型，引入拉格朗日乘子z、w和y，得到拉格朗日函数，根据拉格朗日极值存在必要条件(KKT条件)列写非线性方程组；

根据拉格朗日极值存在条件(KKT条件)得到一组非线性方程组，用牛顿-拉弗逊法求解，其一阶修正方程为：

$\left[\begin{array}{cccccc}I& L^{-1}Z& 0& 0& 0& 0\\ 0& I& 0& 0& -{\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& I& U^{-1}W& 0& 0\\ 0& 0& 0& I& {\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& 0& 0& H^{\&prime;}& {\&dtri;}_x^{T}h\left(x\right)\\ 0& 0& 0& 0& {\&dtri;}_x^{T}h\left(x\right)& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}z\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}w\\ \mathrm{\&Delta;}u\\ \mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\end{array}\right]=\left[\begin{array}{c}-L^{-1}{L}_l^{\mathrm{\&mu;}}\\ L_z\\ -U^{-1}{L}_u^{\mathrm{\&mu;}}\\ -L_w\\ L^{\&prime;}\\ -L_y\end{array}\right];$

其中： $H^{\&prime;}=-\&lsqb;{\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w)\&rsqb;-{\&dtri;}_{x}g\left(x\right)(L^{-1}Z-U^{-1}W){\&dtri;}_x^{T}g\left(x\right);$ $L_x^{\&prime;}=L_x+{\&dtri;}_{x}g\left(x\right)\&lsqb;L^{-1}(L_l^{\mathrm{\&mu;}}+{\mathrm{ZL}}_z)+U^{-1}(L_u^{\mathrm{\&mu;}}-{\mathrm{WL}}_w)\&rsqb;;$ L_x，L_y，L_z，L_w，为KKT方程上次迭代的残差；和分别为h(x)和g(x)的海森矩阵；

8)定义对偶间隙Gap＝l^Tz-u^Tw和扰动因子其中σ＝[0,1]为中心参数；

9)用牛顿法求解上述非线性方程组，计算雅可比矩阵、海森矩阵和常数项等，求解修正方程组，得到各原变量和对偶变量的修正量，乘上步长对变量进行修正，直到对偶间隙Gap小于收敛精度(ε＝10^-6)，否则不收敛；

求解修正方程得到修正量，计算原变量和对偶变量步长：

$\left.\begin{array}{c}{\mathrm{\&alpha;}}_p=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-l_i}{{\mathrm{\&Delta;l}}_i},{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i},{\mathrm{\&Delta;u}}_i<0\right\},1\right\}\\ {\mathrm{\&alpha;}}_d=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-z_i}{{\mathrm{\&Delta;z}}_i},{\mathrm{\&Delta;z}}_i<0;\frac{-w_i}{{\mathrm{\&Delta;w}}_i},{\mathrm{\&Delta;w}}_i>0\right\},1\right\}\end{array}\right\}\left(i=1,2,\mathrm{...},r\right);$

按下式修正变量：

$\left[\begin{array}{c}{x}^{(k+1)}\\ l^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ l^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_p\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}u\end{array}\right];$

$\left[\begin{array}{c}{y}^{(k+1)}\\ z^{(k+1)}\\ w^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{y}^{\left(k\right)}\\ z^{\left(k\right)}\\ w^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_d\left[\begin{array}{c}\mathrm{\&Delta;}y\\ \mathrm{\&Delta;}z\\ \mathrm{\&Delta;}w\end{array}\right];$

10)输出优化调度结果，包括：电力网络发电机有功、无功出力和各节点压力，天然气网络气源点天然气输出和各节点压力值，能源中心各调度因子。

有益效果：本发明相对于现有技术而言：本发明首先建立了电力网络、天然气网络和能源中心模型，电力网络和天然气网络通过能源中心耦合形成电气混联系统。然后以总的能源成本为目标函数、考虑了各种约束条件建立电气混联系统优化调度数学模型。用原-对偶内点法求解，求解过程中先后引入松弛变量和障碍参数将模型变为只有等式约束的模型，再引入拉格朗日乘子得到拉格朗日函数，之后用牛顿法求解其KKT条件形成的非线性方程组。构造的算例仿真结果表明本发明对电气混联系统的优化效果优于单独优化效果。

附图说明

图1为本发明方法流程图；

图2为能源中心模型示意图；

图3为能源中心示例示意图；

图4为算例结果分析的不同能源中心输入结果对比图。

具体实施方案

下面结合附图对发明的技术流程进行详细说明：

电力网络模型实施例

电力网络的拓扑结构可以用导纳矩阵Y描述。电力网络的状态变量有节点电压、节点功率及支路功率等，这里用电压的直角坐标形式给出节点有功功率和无功功率：

$P_i=e_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{e}_j-B_{ij}{f}_j)-f_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{f}_j+B_{ij}{e}_j);$

$Q_i=f_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{e}_j-B_{ij}{f}_j)+e_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{f}_j+B_{ij}{e}_j);$

其中：P_i，Q_i分别为电力网络中节点i有功功率和无功功率；e_i，f_i分别为节点i电压向量的实部和虚部；G_ij，B_ij分别为节点导纳矩阵第i行第j列元素的实部和虚部。

天然气网络模型实施例

典型的天然气网络包括一个或多个气源点(天然气生产点和存储点)，一个或多个负荷(发电厂、其他的网络或存储点)，管道，压缩机和其他设备。用于天然气网络建模的基本要素有三个：1)管道；2)压缩机；3)注入点(节点或母线)。天然气模型两个最主要的部分是管道流量模型和压缩机模型。

天然气管道流量方程描述了天然气流量值与管道两端压力值以及相关的气体性质、管道特性和运行状态之间的关系。形成流量方程需要考虑很多因素，没有一个流量方程适用于所有的情况。理想条件下，管道k从节点i到节点j的流量值可以用下面的方程表示：

$f_k=f_{kij}=S_{ij}\&times;6.4774\frac{T_0}{{\mathrm{\&pi;}}_0}\sqrt{S_{ij}\frac{\left({\mathrm{\&pi;}}_i^2-{\mathrm{\&pi;}}_j^2\right)D_k^5}{F_k{\mathrm{GL}}_k{T}_{ka}{Z}_a}};$

其中：f_kij为管道流量值； $S_{ij}=\left\{\begin{array}{cc}+1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j\&GreaterEqual;0\\ -1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j<0\end{array}\right.;$ F_k为管道摩擦系数；D_k为管道内径；G为气体重力系数；L_k为管道长度；π_i为节点i压力值；π_j为节点j压力值；π₀为标准压力值；T₀为标准温度值；T_ka为平均气体温度；Z_a为平均气体压缩系数。

对于高压网络的全紊流形式，流量方程可以进一步简化为：

$f_k=f_{kij}=M_k{S}_{ij}\sqrt{S_{ij}({\mathrm{\&pi;}}_i^2-{\mathrm{\&pi;}}_j^2)};$

其中： $M_k=\mathrm{\&epsiv;}\frac{18.73T_0{D}_k^{8/3}}{{\mathrm{\&pi;}}_0\sqrt{{\mathrm{GL}}_k{T}_{ka}{Z}_a}},$ ε为管道效率。

压缩机安装在气体网络中为了传输气体和补偿由于摩擦阻力等造成的下游压力损失。压缩机的运行需要消耗大量的能量，在大型天然气网络中，驱动主要压缩机最经济的能量源是通过压缩机的天然气。对于压缩机的网络分析关键的方程是能量消耗，它是流过压缩机的气体流量和注入、流出压缩机的气体压力的函数。在理想气体条件下，加上经验的量化，绝热压缩机的能量消耗方程可以表示为：

$H_k=H_{kij}=B_k{f}_{Ck}\&lsqb;{\left(\frac{{\mathrm{\&pi;}}_j}{{\mathrm{\&pi;}}_i}\right)}^{Z_{ki}\left(\frac{\mathrm{\&alpha;}-1}{\mathrm{\&alpha;}}\right)}-1\&rsqb;;$

其中：f_Ck为通过压缩机的气体流量；π_i为气体注入压缩机压力；π_j为气体输出压缩机压力；Z_ki为压缩机入口的气体压缩系数；T_ki为压缩机输出温度；α为热量系数；η_k为压缩机效率。

转化为消耗的流量值：

${\mathrm{\&tau;}}_k={\mathrm{\&alpha;}}_{Tk}+{\mathrm{\&beta;}}_{Tk}{H}_{kij}+{\mathrm{\&gamma;}}_{Tk}{H}_{kij}^2;$

根据基尔霍夫第一定律：所有流入与流出某网络节点的流量代数和为零。故每一个节点的流量平衡方程可以用下面的矩阵形式表示：

(A+U)f+w-Tτ＝0；

其中：A、U和T分别为支路、压缩机和压缩机消耗与节点有关的方向矩阵，f为支路流量值向量；w为各节点的气体注入向量；τ为各压缩机消耗流量值向量。

能源中心模型实施例

能源中心被定义为由能源转化设备和储能设备构成、能够实现多种能源相互转化和存储的虚拟实体。能源中心之间通过能源传输设备相互连接，构成多能源网络系统。在能源中心内部，能源可能被负荷消耗，或者被转化为其他形式。

能源中心有多种能源形式α，β，…∈ε，每一种能源形式可以是能源中心的输入或输出。这些输入能量(输出能量)被定义为P_α,P_β,…，P_ω(L_α,L_β,…L_ω)。如图2所示。对于多输入多输出的能源中心可以用下面的矩阵方程描述：

其中：P和L分别为输入输出向量；矩阵C被称为耦合矩阵。耦合矩阵从数学上描绘了能量从能源中心输入到输出的分配。耦合矩阵中的元素为耦合因子。每一个耦合因子表示了一个特定的输入到一个特定的输出。

图3给出了一种能源中心的示例。能源中心输入由电力和天然气提供，能源在能源中心内部通过变压器、燃气轮机和燃气锅炉转化，输出供应电负荷和热负荷。其中ν为调度因子。据此可以列出能源中心输入输出平衡方程：

L_e＝η_TP_e+νη_GTeP_g；

L_h＝νη_GThP_g+(1-ν)η_FP_g；

写成矩阵形式：

$\left[\begin{array}{c}{L}_e\\ L_h\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\&eta;}}_T& {\mathrm{\&nu;\&eta;}}_{GTe}\\ 0& {\mathrm{\&nu;\&eta;}}_{GTh}+(1-\mathrm{\&nu;}){\mathrm{\&eta;}}_F\end{array}\right]\left[\begin{array}{c}{P}_e\\ P_g\end{array}\right];$

其中：η_T表示变压器效率；η_GTe表示燃气轮机气转电的效率；η_GTh表示燃气轮机气转热的效率；η_F表示燃气锅炉的效率。

电气混联系统优化调度模型实施例

根据上述模型，电力网络与天然气网络通过能源中心耦合，形成电气混联系统。以总的能源成本为目标函数，考虑电力网络、天然气网络及能源中心的各种约束，建立电气混联系统优化调度模型：

minf(x)；

s.t.h(x)＝0；

g_min≤g(x)≤g_max；

其中：h表示等式约束；g表示不等式约束；g_max，g_min分别表示不等式约束的上、下限。

进行选取变量，控制变量：

1)电力网络控制变量：发电机有功出力P_G；发电机无功出力Q_R；

2)天然气网络控制变量：气源点天然气供应量N_G；

3)能源中心控制变量：调度因子ν。

则本发明选取的控制变量为：

u＝[P_G1,...,P_Gng,Q_G1,...,Q_Gng,N_G1,…N_GnGS,ν₁,…,ν_nEH]^T；

其中：n_g为电力网络发电机台数；n_GS为天然气网络气源个数；n_EH为能源中心个数。

状态变量：

1)电力网络状态变量：节点电压向量的实部e和虚部f；

2)天然气网络状态变量：节点的压力π和通过压缩机的流量f_C；

3)能源中心状态变量：电能的输入量P_e和天然气的输入量P_g。

则本发明选取的状态变量为：

$\tilde{x}={\&lsqb;e_1,...,e_{nb},f_1,...,f_{nb},p_1,...,p_{nGB},f_{CK1},...,f_{CKnGC},{\mathrm{\&pi;}}_{e1},...,{\mathrm{\&pi;}}_{enEH},P_{g1},...,P_{gnEH}\&rsqb;}^T;$

其中：n_b为电力网络的节点个数；n_GB为天然气网络节点个数；n_GC为天然气网络中压缩机的个数。

总的变量定义为：

本发明以系统总的能源成本作为目标函数：

$f\left(x\right)=\underset{i\&Element;\{P_G,N_G\}}{\&Sigma;}(a_{2i}{P_i}^2+a_{1i}{P}_i+a_{0i});$

其中：a_2i，a_1i，a_0i为能量损耗曲线参数；P_i为广义能量，这里包括电力网中发电机发出的有功功率和天然气网络气源点提供气体的能量。

等式约束

1)电力网络等式约束:

ΔP_i＝P_Gi-P_Di-P_ei-P_i＝0；

ΔQ_i＝Q_Ri-Q_Di-Q_i＝0；

$\mathrm{\&Delta;}angle={\mathrm{tan\&theta;}}_{bal}-\frac{f_{bal}}{e_{bal}}=0;$

其中：ΔP_i，ΔQ_i为各节点有功、无功功率不平衡量；Δangle为平衡节点角度约束；tanθ_bal为平衡节点的角度正切值；P_Gi，Q_Ri分别为发电机i的有功、无功出力；P_Di，Q_Di分别为节点i的有功、无功负荷；P_ei为节点i对能源中心的电能输入。

2)天然气网络等式约束:

ΔW_i＝N_Gi+w_i-F_i-P_gi＝0；

${\mathrm{\&Delta;P}}_{Ck}={\stackrel{\&OverBar;}{R}}_k-R_k=0;$

其中：ΔW_i为天然气网络中各节点流量值的不平衡量；N_Gi为气源点对节点i的气体注入量；w_i，F_i分别为节点i的气体注入量和负荷量；P_gi为节点i对能源中心的天然气输入；ΔP_Ck为天然气网络中压缩机k状态量的不平衡量；R_k分别为压缩机k状态量的给定值和计算值。

3)能源中心等式约束:

ΔL_ei＝L_ei-η_TP_ei-ν_iη_GTeP_gi；

ΔL_hi＝L_hi-ν_iη_GThP_gi-(1-ν_i)η_FP_gi；

其中：ΔL_ei，ΔL_hi分别为能源中心i的电能负荷与热负荷的不平衡量；L_ei，L_hi分别为能源中心i的电能负荷与热负荷；P_ei，P_gi分别为能源中心i的电能与天然气输入；ν_i为能源中心i的调度因子。

不等式约束

1)电力网络不等式约束：

P_Gimin≤P_Gi≤P_Gimax；

Q_Rimin≤Q_Ri≤Q_Rimax；

$V_{imin}^2\&le;{e_i}^2+{f_i}^2\&le;V_{imax}^2;$

其中：P_Gimin，P_Gimax为发电机所发出有功功率的下限和上限；Q_Rimin，Q_Rimax为发电机所发无功功率的下限和上限；为节点电压幅值平方的下限和上限。

2)天然气网络不等式约束：

0≤N_Gi≤N_Gimax；

π_imin≤π_i≤π_imax；

其中：N_Gimax为天然气网络中各气源点气体供应量的上限；π_imin与π_imax分别为各节点压力值的下限和上限。

3)能源中心不等式约束：

0≤ν_i≤1；

通过原-对偶内点法求解，首先引入广义松弛变量l＝[l₁,…,l_r]^T,u＝[u₁,…,u_r]^T将广义不等式转化为等式：

g(x)+u-g_max＝0；

g(x)-l-g_min＝0；

应满足：u>0,l>0。

然后，为了保证目标函数在可行域内取得最小值，对目标函数引入对数障碍参数μ，当u,l靠近约束边界时，目标函数趋于无穷大。此时优化调度模型为只有等式约束的模型：

$\min .f\left(x\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}log\left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}log\left(u_j\right);$

s.t.h(x)＝0；

g(x)+u-g_max＝0；

g(x)-l-g_min＝0；

其中：扰动因子(障碍参数)μ＞0。

这时可以直接通过拉格朗日乘子法求解。引入拉格朗日乘子y＝[y₁,...,y_m]^T，z＝[z₁,…,z_r]^T，w＝[w₁,...,w_m]^T得到拉格朗日函数：

$\begin{array}{c}Lag=f\left(x\right)-y^{T}h\left(x\right)-z^T\&lsqb;g\left(x\right)-l-g_{\min }\&rsqb;\\ -w^T\&lsqb;g\left(x\right)+u-g_{\max }\&rsqb;-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(u_j\right);\end{array}$

根据拉格朗日函数极值存在的必要条件(KKT条件)可以得到非线性方程组：

$L_x=\frac{\&part;Lag}{\&part;x}\&equiv;{\&dtri;}_{x}f\left(x\right)-{\&dtri;}_{x}h\left(x\right)y-{\&dtri;}_{x}g\left(x\right)(z+w)=0;$

$L_y=\frac{\&part;Lag}{Qy}\&equiv;h\left(x\right)=0;$

$L_z=\frac{\&part;Lag}{\&part;z}\&equiv;g\left(x\right)-l-g_{min}=0;$

$L_w=\frac{\&part;Lag}{\&part;w}\&equiv;g\left(x\right)+u-g_{max}=0;$

$L_l=\frac{\&part;Lag}{\&part;l}\&equiv;z-{\mathrm{\&mu;L}}^{-1}e\&DoubleRightArrow;L_l^{\mathrm{\&mu;}}=LZe-\mathrm{\&mu;}e=0;$

其中：e为各元素为1的_r维列向量；L＝diag(l)；U＝diag(u)；Z＝diag(z)；W＝diag(w)均为对角矩阵。

计算可以得到：

$\mathrm{\&mu;}=\frac{l^{T}z-u^{T}w}{2r};$

定义互补间隙(对偶间隙)为：Gap＝l^Tz-u^Tw。实践表明，目标函数中的扰动因子μ按下式计算收敛性较好：

$\mathrm{\&mu;}=\mathrm{\&sigma;}\frac{Gap}{2r};$

其中：σ＝[0,1]为中心参数。σ＝1时，有利于提高算法的可行性，但对减少对偶间隙无作用；σ＝0减少了算法的对偶间隙，有利于解的最优性。在保证可行性的基础上尽量提高解的最优性，σ一般取0.1。

可以用牛顿-拉弗逊法求解上述非线性方程组，其一阶修正方程组的矩阵形式为：

$\left[\begin{array}{cccccc}I& L^{-1}Z& 0& 0& 0& 0\\ 0& I& 0& 0& -{\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& I& U^{-1}W& 0& 0\\ 0& 0& 0& I& {\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& 0& 0& H^{\&prime;}& {\&dtri;}_x^{T}h\left(x\right)\\ 0& 0& 0& 0& {\&dtri;}_x^{T}h\left(x\right)& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}z\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}w\\ \mathrm{\&Delta;}u\\ \mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\end{array}\right]=\left[\begin{array}{c}-L^{-1}{L}_l^{\mathrm{\&mu;}}\\ L_z\\ -U^{-1}{L}_u^{\mathrm{\&mu;}}\\ -L_w\\ L^{\&prime;}\\ -L_y\end{array}\right];$

其中： $H^{\&prime;}=-\&lsqb;{\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w)\&rsqb;-{\&dtri;}_{x}g\left(x\right)(L^{-1}Z-U^{-1}W){\&dtri;}_x^{T}g\left(x\right);$ $L_x^{\&prime;}=L_x+{\&dtri;}_{x}g\left(x\right)\&lsqb;L^{-1}(L_l^{\mathrm{\&mu;}}+{\mathrm{ZL}}_z)+U^{-1}(L_u^{\mathrm{\&mu;}}-{\mathrm{WL}}_w)\&rsqb;;$ L_x，L_y，L_z，L_w，为KKT方程上次迭代的残差；和分别为h(x)和g(x)的海森矩阵。

求解修正方程得到修正量，计算原变量和对偶变量步长：

$\left.\begin{array}{c}{\mathrm{\&alpha;}}_p=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-l_i}{{\mathrm{\&Delta;l}}_i},{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i},{\mathrm{\&Delta;u}}_i<0\right\},1\right\}\\ {\mathrm{\&alpha;}}_d=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-z_i}{{\mathrm{\&Delta;z}}_i},{\mathrm{\&Delta;z}}_i<0;\frac{-w_i}{{\mathrm{\&Delta;w}}_i},{\mathrm{\&Delta;w}}_i>0\right\},1\right\}\end{array}\right\}\left(i=1,2,\mathrm{...},r\right);$

按下式修正变量：

$\left[\begin{array}{c}{x}^{(k+1)}\\ l^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ l^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_p\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}u\end{array}\right];$

$\left[\begin{array}{c}{y}^{(k+1)}\\ z^{(k+1)}\\ w^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{y}^{\left(k\right)}\\ z^{\left(k\right)}\\ w^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_d\left[\begin{array}{c}\mathrm{\&Delta;}y\\ \mathrm{\&Delta;}z\\ \mathrm{\&Delta;}w\end{array}\right];$

依次迭代，直到对偶间隙Gap小于收敛精度。

实施例

为验证本发明效果，将Matpower中IEEE14节点电力网络与某14节点天然气网络通过图3中的能源中心耦合构造一个算例。其中电力网络与天然气网络通过5个能源中心联系，将天然气网络中原来的5个负荷点连接到能源中心，作为能源中心的天然气输入；选取电力网络中负荷比较重的5个负荷点与能源中心连接，作为能源中心的电能输入，能源中心负荷侧考虑电能负荷和热能输出，并同时考虑能源中心调度因子的影响。各能源中心与电力网络及天然气网络具体连接情况和能源中心输出情况见表1所示。

用Matlab对构造的算例进行编程并仿真，对仿真结果进行分析。

当将天然气网络中各节点压力值限制在500～1500psia内，取天然气损耗曲线系数a₂＝0,a₁＝2,a₀＝0，能源中心输出取表1中数据时，混联系统综合优化的最小成本为$5338.5239。电力网单独优化的最小成本为$8081.526，天然气网络单独优化的最小成本为a₁·N_G·GHV/baseMVA＝$89.98，故当电力网与天然气单独优化的总成本为$8081.526+$89.98＝$8171.506。可见混联系统的综合优化成本明显低于单独优化成本，说明了电力网络与天然气网混联系统的综合优化具有显著的经济效益。

为了更充分说明混联系统综合优化成本低于单独优化成本，选取与上面相同的压力范围和天然气损耗曲线系数，对不同的能源中心电能输出值与热能输出值仿真，优化结果如图4所示。从对比图可知，对于不同的能源中心输出，混联系统综合优化成本均低于独立优化的成本，进一步说明了混联系统综合优化的优势。

表1混联系统节点连接与能源中心输出

Claims (3)

1.一种基于能源中心的电气混联系统建模和优化调度方法，其特征在于：包括以下步骤：

1)获取电力网络的参数信息，包括：输电线路的首端节点和末端节点编号，支路π型等效电路的电阻、电抗，对地并联电导、电纳，变压器变比和阻抗，各节点负荷以及发电机输出有功、无功约束，各节点电压约束；

2)获取天然气网络的参数信息，包括：输气管道的首端节点和末端节点编号，管道的长度、内径和传输效率等物理特性，压缩机的状态量，各节点气负荷以及气源点天然气输出约束，各节点压力约束；

3)多输入多输出能源中心模型可以用矩阵方程描述：

其中：P和L分别为输入输出向量；矩阵C被称为耦合矩阵；

列出能源中心输入输出平衡方程：

$\left[\begin{array}{c}{L}_e\\ L_h\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\&eta;}}_T& {\mathrm{\&nu;\&eta;}}_{GTe}\\ 0& {\mathrm{\&nu;\&eta;}}_{GTh}+(1-\mathrm{\&nu;}){\mathrm{\&eta;}}_F\end{array}\right]\left[\begin{array}{c}{P}_e\\ P_g\end{array}\right];$

其中：η_T表示变压器效率；η_GTe表示燃气轮机气转电的效率；η_GTh表示燃气轮机气转热的效率；η_F表示燃气锅炉的效率；

4)建立电气混联系统优化调度模型：

minf(x)；

s.t.h(x)＝0；

g_min≤g(x)≤g_max；

其中：h表示等式约束；g表示不等式约束；g_max，g_min分别表示不等式约束的上、下限；

5)引入松弛变量和障碍参数将模型转化成只有等式约束的模型：

$min.f\left(x\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}log\left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}log\left(u_j\right);$

s.t.h(x)＝0；

g(x)+u-g_max＝0；

g(x)-l-g_min＝0；

其中：扰动因子(障碍参数)μ＞0；

6)引入拉格朗日乘子得到拉格朗日函数：

$\begin{array}{c}Lag=f\left(x\right)-y^{T}h\left(x\right)-z^T\&lsqb;g\left(x\right)-l-g_{\min }\&rsqb;\\ -w^T\&lsqb;g\left(x\right)+u-g_{\max }\&rsqb;-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(l_j\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\log \left(u_j\right)\end{array};$

7)根据拉格朗日极值存在条件(KKT条件)得到一组非线性方程组，用牛顿-拉弗逊法求解，其一阶修正方程为：

$\left[\begin{array}{cccccc}I& L^{-1}Z& 0& 0& 0& 0\\ 0& I& 0& 0& -{\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& I& U^{-1}W& 0& 0\\ 0& 0& 0& I& {\&dtri;}_x^{T}g\left(x\right)& 0\\ 0& 0& 0& 0& H^{\&prime;}& {\&dtri;}_x^{T}h\left(x\right)\\ 0& 0& 0& 0& {\&dtri;}_x^{T}h\left(x\right)& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}z\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}w\\ \mathrm{\&Delta;}u\\ \mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\end{array}\right]=\left[\begin{array}{c}-L^{-1}{L}_l^{\mathrm{\&mu;}}\\ L_z\\ -U^{-1}{L}_u^{\mathrm{\&mu;}}\\ -L_w\\ L^{\&prime;}\\ -L_y\end{array}\right];$

其中： $H^{\&prime;}=-\&lsqb;{\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w)\&rsqb;-{\&dtri;}_{x}g\left(x\right)(L^{-1}Z-U^{-1}W){\&dtri;}_x^{T}g\left(x\right);$ $L_x^{\&prime;}=L_x+{\&dtri;}_{x}g\left(x\right)\&lsqb;L^{-1}(L_l^{\mathrm{\&mu;}}+{\mathrm{ZL}}_z)+U^{-1}(L_u^{\mathrm{\&mu;}}-{\mathrm{WL}}_w)\&rsqb;;$ L_x，L_y，L_z，L_w，为KKT方程上次迭代的残差；和分别为h(x)和g(x)的海森矩阵；

8)求解修正方程得到修正量，计算原变量和对偶变量步长：

$\left.\begin{array}{c}{\mathrm{\&alpha;}}_p=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-l_i}{{\mathrm{\&Delta;l}}_i},{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i},{\mathrm{\&Delta;u}}_i<0\right\},1\right\}\\ {\mathrm{\&alpha;}}_d=0.9995\min \left\{\underset{i}{\min }\left\{\frac{-z_i}{{\mathrm{\&Delta;z}}_i},{\mathrm{\&Delta;z}}_i<0;\frac{-w_i}{{\mathrm{\&Delta;w}}_i},{\mathrm{\&Delta;w}}_i>0\right\},1\right\}\end{array}\right\}\left(i=1,2,\mathrm{...},r\right);$

9)按下式修正变量：

$\left[\begin{array}{c}{x}^{(k+1)}\\ l^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ l^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_p\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}u\end{array}\right];$

$\left[\begin{array}{c}{y}^{(k+1)}\\ z^{(k+1)}\\ w^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{y}^{\left(k\right)}\\ z^{\left(k\right)}\\ w^{\left(k\right)}\end{array}\right]+{\mathrm{\&alpha;}}_d\left[\begin{array}{c}\mathrm{\&Delta;}y\\ \mathrm{\&Delta;}z\\ \mathrm{\&Delta;}w\end{array}\right];$

10)按上述步骤迭代，直到对偶间隙小于收敛精度，得到电气混联系统优化调度的最优解。

2.根据权利要求1所述的基于能源中心的电气混联系统建模和优化调度方法，其特征在于：所述步骤1中，对于电力网络中节点i：

$P_i=e_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{e}_j-B_{ij}{f}_j)-f_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{f}_j+B_{ij}{e}_j);$

$Q_i=f_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{e}_j-B_{ij}{f}_j)+e_i\underset{j=1}{\overset{nb}{\&Sigma;}}(G_{ij}{f}_j+B_{ij}{e}_j);$

其中：P_i，Q_i分别为电力网络中节点i有功功率和无功功率；e_i，f_i分别为节点i电压向量的实部和虚部；G_ij，B_ij分别为节点导纳矩阵第i行第j列元素的实部和虚部。

3.根据权利要求1所述的基于能源中心的电气混联系统建模和优化调度方法，其特征在于：所述步骤2中，引入包括管道流量方程、压缩机流量消耗方程及流量平衡方程；

理想条件下，管道k从节点i到节点j的流量值可以用下面的方程表示：

$f_k=f_{kij}=S_{ij}\&times;6.4774\frac{T_0}{{\mathrm{\&pi;}}_0}\sqrt{S_{ij}\frac{\left({\mathrm{\&pi;}}_i^2-{\mathrm{\&pi;}}_j^2\right)D_k^5}{F_k{\mathrm{GL}}_k{T}_{ka}{Z}_a}};$

其中：f_kij为管道流量值； $S_{\mathrm{ij}}=\left\{\begin{array}{cc}+1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j\&GreaterEqual;0\\ -1& {\mathrm{\&pi;}}_i-{\mathrm{\&pi;}}_j<0\end{array}\right.;$ F_k为管道摩擦系数；D_k为管道内径；G为气体重力系数；L_k为管道长度；π_i为节点i压力值；π_j为节点j压力值；π₀为标准压力值；T₀为标准温度值；T_ka为平均气体温度；Z_a为平均气体压缩系数；

对于高压网络的全紊流形式，流量方程可以进一步简化为：

$f_k=f_{kij}=M_k{S}_{ij}\sqrt{S_{ij}({\mathrm{\&pi;}}_i^2-{\mathrm{\&pi;}}_j^2)};$

其中： $M_k=\mathrm{\&epsiv;}\frac{18.73T_0{D}_k^{8/3}}{{\mathrm{\&pi;}}_0\sqrt{{\mathrm{GL}}_k{T}_{ka}{Z}_a}},$ ε为管道效率；

理想气体条件下，压缩机的能量消耗方程可以表示为：

$H_k=H_{kij}=B_k{f}_{Ck}\&lsqb;{\left(\frac{{\mathrm{\&pi;}}_j}{{\mathrm{\&pi;}}_i}\right)}^{Z_{ki}\left(\frac{\mathrm{\&alpha;}-1}{\mathrm{\&alpha;}}\right)}-1\&rsqb;;$

其中：f_Ck为通过压缩机的气体流量；π_i为气体注入压缩机压力；π_j为气体输出压缩机压力；Z_ki为压缩机入口的气体压缩系数；T_ki为压缩机输出温度；α为热量系数；η_k为压缩机效率；

转化为消耗的流量值：

${\mathrm{\&tau;}}_k={\mathrm{\&alpha;}}_{Tk}+{\mathrm{\&beta;}}_{Tk}{H}_{kij}+{\mathrm{\&gamma;}}_{Tk}{H}_{kij}^2.$