Patent ID: 11972182
Assignee: TSINGHUA UNIVERSITY
Field: Computer technology (Electrical engineering)
Classification: CPC G | IPC G

Claim 0:
1. A method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines to obtain a real-time, reliable and complete operating state of the natural gas network, the method comprising:
step 1 of establishing a time-domain window and a frequency-domain window for the dynamic state estimation of the natural gas network, the step 1 comprising:
sub-step 1-1 of defining a time-domain window width as It, where It is a positive integer, and a value of It is determined by a dispatcher of the natural gas network; defining a u-th sampling time point in the time-domain window as τu=τ−uΔt, u=0, 1, . . . , It−1, where τ represents a current time point of the natural gas network, and Δt represents a sampling interval of the natural gas network; defining a current time-domain window width as It,e, where It,e is a positive integer, and a value of It,e is determined by the dispatcher of the natural gas network; and defining a historical time-domain window width as It,h where It,h is a positive integer, and a value of It,h is determined by the dispatcher of the natural gas network, wherein It, It,e and It,h satisfy the following relational expression:

It=It,e+It,h; and

sub-step 1-2 of defining a frequency-domain window width as If, where a value of If is determined by the dispatcher of the natural gas network; and defining a d-th frequency component in the frequency-domain window as ωd, d=0, 1, . . . , If−1, where ωd is calculated by the following formula:, ω
      d
    
    =
    
      d
      
        
          
            I
            t
          
          ·
          Δ
        
        ⁢
        
          
        
        ⁢
        t
      
    
  
  .

step 2 of constructing a measurement vector for the dynamic state estimation of the natural gas network, the step 2 comprising:
sub-step 2-1 of acquiring, from a data acquisition and monitoring control system of the natural gas network, all operation data of the natural gas network at a sampling time point τu in the time-domain window where the current time point τ of the natural gas network belongs, wherein the all operation data of the natural gas network comprises: a measurement value zG+,uip of a natural gas flow at a head end of each pipeline in the natural gas network, and a measurement value zG−,uip of a natural gas flow at a tail end of each pipeline in the natural gas network, where ip represents a serial number of a pipeline in the natural gas network; a measurement value zG+,uic of a natural gas flow at a head end of each compressor, and a measurement value zG−,uic of a natural gas flow at a tail end of each compressor, where ic represents a serial number of a compressor; a pressure measurement value zpr,uin of each node of the natural gas network, where in represents a serial number of a node of the natural gas network; a measurement value zgs,uis of a natural gas flow of each natural gas source, where is represents a serial number of a natural gas source; and a measurement value zgl,uil of a natural gas flow of each natural gas load, where il represents a serial number of a natural gas load; and
sub-step 2-2 of constructing a measurement vector zu for the dynamic state estimation of the natural gas network at the sampling time point τu:, z
      u
    
    =
    
      [
      
        
          
            
              Z
              
                
                  G
                  +
                
                ,
                u
              
            
          
        
        
          
            
              Z
              
                
                  G
                  -
                
                ,
                u
              
            
          
        
        
          
            
              Z
              
                pr
                ,
                u
              
            
          
        
        
          
            
              Z
              
                gs
                ,
                u
              
            
          
        
        
          
            
              Z
              
                gl
                ,
                u
              
            
          
        
      
      ]
    
  
  ,

where zG+,u represents a column vector consisting of all the measurement values zG+,uip of natural gas flows at head ends of respective pipelines in the natural gas network and all the measurement values zG+,uic of natural gas flows at head ends of respective compressors at the sampling time point τu; zG−,u represents a column vector consisting of all the measurement values zG−,uip of natural gas flows at tail ends of respective pipelines in the natural gas network and all the measurement values zG−,uic of natural gas flows at tail ends of respective compressors at the sampling time point τu; zpr,u represents a column vector consisting of all the pressure measurement values zpr,uin of respective nodes of the natural gas network at the sampling time point τu; zgs,u represents a column vector consisting of all the measurement values zgs,uis of natural gas flows of respective natural gas sources in the natural gas network at the sampling time point τu; and zgl,u represents a column vector consisting of all the measurement values zgl,uil of natural gas flows of respective natural gas loads in the natural gas network at the sampling time point τu;
step 3 of constructing a state vector xu for the dynamic state estimation of the natural gas network at the sampling time point τu:, x
      u
    
    =
    
      [
      
        
          
            
              x
              
                
                  G
                  +
                
                ,
                u
              
            
          
        
        
          
            
              x
              
                
                  G
                  -
                
                ,
                u
              
            
          
        
        
          
            
              x
              
                pr
                ,
                u
              
            
          
        
        
          
            
              x
              
                gs
                ,
                u
              
            
          
        
        
          
            
              x
              
                gl
                ,
                u
              
            
          
        
      
      ]
    
  
  ,

where xG+,u represents a column vector consisting of all natural gas flows Gip,u+ at the head ends of the respective pipelines in the natural gas network and all natural gas flows Gic,u+ at the head ends of the respective compressors at the sampling time point τu; xG−,u represents a column vector consisting of all natural gas flows Gip,u− at the tail ends of the respective pipelines in the natural gas network and all natural gas flows Gic,u− at the tail ends of the respective compressors at the sampling time point τu; xpr,u represents a column vector consisting of all pressures hin of the respective nodes of the natural gas network at the sampling time point τu; xgs,u represents a column vector consisting of all natural gas flows Gis,ugs of the respective natural gas sources in the natural gas network at the sampling time point τu; and xgl,u represents a column vector consisting of all natural gas flows Gilugl of the respective natural gas loads in the natural gas network at the sampling time point τu;
step 4 of establishing, based on the measurement vector constructed in the step 2 and the state vector constructed in the step 3, an objective function of the dynamic state estimation of the natural gas network as follows:

min J=Σu=0It,e−1{[zu−xu]W−1[zu−xu]T}+Σu=It,eIt−1{[zu−xu]W−1δu-It,e[zu−xu]T},

where J represents an expression of the objective function; W represents a covariance matrix of a measurement error and is determined by the dispatcher of the natural gas network; a superscript T represents a matrix transpose; and δ represents a decay factor of a historical time window and is determined by the dispatcher of the natural gas network;
step 5 of establishing constraint conditions for the dynamic state estimation of the natural gas network, the step 5 comprising:
sub-step 5-1 of establishing constraints related to a flow and a pressure of the compressor in the natural gas network, the sub-step 5-1 comprising:
establishing a flow constraint of the head end and the tail end of a compressor:

Gic,u+=Gic,u−,∀ic∈Ωc,∀u=0,1, . . . ,It−1,

where Ωc represents a set of serial numbers of all the compressors in the natural gas network; and
establishing a pressure constraint at the head end and the tail end of the compressor, wherein, for the compressor with a constant tail end pressure, the pressure constraint of the head end and the tail end of the compressor is as follows:

hic,u−=hic,con−,∀ic∈Ωc,1,∀u=0,1, . . . ,It−1,

where hic,u− represents a tail end pressure of a compressor ic at the sampling time point τu; hic,con− represents a set value of a tail end pressure of the compressor ic and is a constant determined by the dispatcher of the natural gas network; and Ωc,1 represents a set of serial numbers of all the compressors with the constant tail end pressure in the natural gas network;
wherein, for a compressor with a constant compression ratio, the pressure constraint of the head end and the tail end of the compressor is as follows:

hic,u−=ric,con·hic,u+,∀ic∈Ωc,2,∀u=0,1, . . . It−1,

where hic,u+ represents a head end pressure of the compressor ic at the sampling time point τu; ric,con represents a set value of a compression ratio of the compressor ic and is a constant determined by the dispatcher of the natural gas network; and Ωc,2 represents a set of serial numbers of all the compressors with the constant compression ratio in the natural gas network; and
wherein, for a compressor with a constant pressure difference, the pressure constraint of the head end and the tail end of the compressor is as follows:

hic,u−−hic,u+=Δhic,con·,∀ic∈Ωc,3,∀u=0,1, . . . ,It−1,

where Δhic,con represents a set value of a pressure difference between the tail end and the head end of the compressor ic and is a constant determined by the dispatcher of the natural gas network; and ψc,3 represents a set of serial numbers of all the compressors with the constant pressure difference in the natural gas network;
sub-step 5-2 of establishing a flow constraint and a pressure constraint of the natural gas in the pipeline in the natural gas network, the sub-step 5-2 comprising:
establishing a two-port constraint of the pipeline in the natural gas network of each frequency component ωd in the frequency-domain window:, [
      
        
          
            
              h
              
                
                  i
                  p
                
                ,
                d
              
              -
            
          
        
        
          
            
              G
              
                
                  i
                  p
                
                ,
                d
              
              -
            
          
        
      
      ]
    
    =
    
      
        [
        
          
            
              
                A
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
            
            
              
                B
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
            
          
          
            
              
                C
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
            
            
              
                D
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
            
          
        
        ]
      
      ⁡
      
        [
        
          
            
              
                h
                
                  
                    i
                    p
                  
                  ,
                  d
                
                +
              
            
          
          
            
              
                G
                
                  
                    i
                    p
                  
                  ,
                  d
                
                +
              
            
          
        
        ]
      
    
  
  ,
  
    

  
  ⁢
  
    
      
        ∀
      
      ⁢
      
        i
        p
      
    
    ∈
    
      Ω
      p
    
  
  ,
  
    

  
  ⁢
  
    
      
        ∀
      
      ⁢
      d
    
    =
    0
  
  ,
  1
  ,
  …
  ⁢
  
    
  
  ,
  
    
      I
      f
    
    -
    1
  
  ,

where hip,d− represents a value of a d-th component of a tail end pressure of a pipeline ip in the natural gas network in a frequency-domain window If, and hip,d− is a complex variable to be solved; hip,d+ represents a value of a d-th frequency component of a head end pressure of the pipeline ip in the natural gas network, and hip,d+ is a complex variable to be solved; Gip,d− represents a value of a d-th component of a natural gas flow at the tail end of the pipeline ip in the natural gas network in the frequency-domain window If, and Gip,d− is a complex variable to be solved; Gip,d+ represents a value of a d-th component of a natural gas flow at the head end of the pipeline ip in the natural gas network in the frequency-domain window, and Gip,d+ is a complex variable to be solved; and Aip,d, Bip,d, Cip,d and Dip,d represent two-port parameters of a d-th component of the pipeline ip in the natural gas network in the frequency-domain window, and values of Aip,d, Bip,d, Cip,d and DiP,d are respectively expressed as:, A
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        [
        
          
            cosh
            (
            
              
                l
                
                  i
                  p
                
              
              ·
              
                a
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
            
            )
          
          -
          
            
              b
              
                
                  i
                  p
                
                ,
                d
              
            
            ·
            
              sinh
              (
              
                
                  l
                  
                    i
                    p
                  
                
                ·
                
                  a
                  
                    
                      i
                      p
                    
                    ,
                    d
                  
                
              
              )
            
          
        
        ]
      
      ·
      
        e
        
          -
          
            
              
                k
                
                  i
                  p
                
              
              ·
              
                l
                
                  i
                  p
                
              
            
            2
          
        
      
    
  

  
    
      B
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        -
        
          
            Z
            
              
                i
                p
              
              ,
              d
            
          
          
            
              
                
                  (
                  
                    k
                    
                      i
                      p
                    
                  
                  )
                
                2
              
              +
              
                
                  4
                  ·
                  
                    Z
                    
                      
                        i
                        p
                      
                      ,
                      d
                    
                  
                
                ⁢
                
                  Y
                  
                    
                      i
                      p
                    
                    ,
                    d
                  
                
              
            
          
        
      
      ⁢
      
        
          sinh
          (
          
            
              l
              
                i
                p
              
            
            ·
            
              a
              
                
                  i
                  p
                
                ,
                d
              
            
          
          )
        
        ·
        
          e
          
            -
            
              
                
                  k
                  
                    i
                    p
                  
                
                ·
                
                  l
                  
                    i
                    p
                  
                
              
              2
            
          
        
      
    
  

  
    
      C
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        -
        
          
            Y
            
              
                i
                p
              
              ,
              d
            
          
          
            
              
                
                  (
                  
                    k
                    
                      i
                      p
                    
                  
                  )
                
                2
              
              +
              
                
                  4
                  ·
                  
                    Z
                    
                      
                        i
                        p
                      
                      ,
                      d
                    
                  
                
                ⁢
                
                  Y
                  
                    
                      i
                      p
                    
                    ,
                    d
                  
                
              
            
          
        
      
      ⁢
      
        
          sinh
          (
          
            
              l
              
                i
                p
              
            
            ·
            
              a
              
                
                  i
                  p
                
                ,
                d
              
            
          
          )
        
        ·
        
          e
          
            -
            
              
                
                  k
                  
                    i
                    p
                  
                
                ·
                
                  l
                  
                    i
                    p
                  
                
              
              2
            
          
        
      
    
  

  
    
      
        D
        
          
            i
            p
          
          ,
          d
        
      
      =
      
        
          [
          
            
              cosh
              (
              
                
                  l
                  
                    i
                    p
                  
                
                ·
                
                  a
                  
                    
                      i
                      p
                    
                    ,
                    d
                  
                
              
              )
            
            -
            
              
                b
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
              ·
              
                sinh
                (
                
                  
                    l
                    
                      i
                      p
                    
                  
                  ·
                  
                    a
                    
                      
                        i
                        p
                      
                      ,
                      d
                    
                  
                
                )
              
            
          
          ]
        
        ·
        
          e
          
            -
            
              
                
                  k
                  
                    i
                    p
                  
                
                ·
                
                  l
                  
                    i
                    p
                  
                
              
              2
            
          
        
      
    
    ,
  

where lip represents a length of the pipeline ip in the natural gas network, kip, aip,d, bip,d, Zip,d and Yip,d represent values of the d-th frequency component of pipeline parameters of the natural gas network, and values of kip, aip,d, bip,d, Zip,d and Yip,d are respectively expressed as:, k
      
        i
        p
      
    
    =
    
      -
      
        
          
            2
            ⁢
            g
            ⁢
            
              D
              
                i
                p
              
            
            ⁢
            sin
            ⁢
            
              
            
            ⁢
            
              θ
              
                i
                p
              
            
          
          -
          
            
              
                λ
                
                  i
                  p
                
              
              (
              
                v
                
                  base
                  ,
                  
                    i
                    p
                  
                
              
              )
            
            2
          
        
        
          2
          ⁢
          
            RTD
            
              i
              p
            
          
        
      
    
  

  
    
      a
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        1
        2
      
      ⁢
      
        
          
            
              (
              
                k
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
              )
            
            2
          
          +
          
            4
            ⁢
            
              Z
              
                
                  i
                  p
                
                ,
                d
              
            
            ⁢
            
              Y
              
                
                  i
                  p
                
                ,
                d
              
            
          
        
      
    
  

  
    
      b
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        k
        
          
            i
            p
          
          ,
          d
        
      
      
        
          
            
              (
              
                k
                
                  
                    i
                    p
                  
                  ,
                  d
                
              
              )
            
            2
          
          +
          
            4
            ⁢
            
              Z
              
                
                  i
                  p
                
                ,
                d
              
            
            ⁢
            
              Y
              
                
                  i
                  p
                
                ,
                d
              
            
          
        
      
    
  

  
    
      Z
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      
        R
        
          i
          p
        
      
      +
      
        j
        ⁢
        
          ω
          d
        
        ⁢
        
          L
          
            i
            p
          
        
      
    
  

  
    
      Y
      
        
          i
          p
        
        ,
        d
      
    
    =
    
      j
      ⁢
      
        
      
      ⁢
      
        ω
        d
      
      ⁢
      
        C
        
          i
          p
        
      
    
  

where g represents an acceleration of gravity; Dip represents an inner diameter of the pipeline ip in the natural gas network; θip represents an angle of inclination of the pipeline ip in the natural gas network; λip represents a friction coefficient of the pipeline ip in the natural gas network; vbase,ip represents a basic value of a flow velocity in the pipeline ip in the natural gas network; R represents a gas constant of natural gas; T represents a temperature of the natural gas; j represents a complex number unit; and Rip, Lip and Cip are parameters of the natural gas network, and values of Rip, Lip and Cip are respectively expressed as:

Rip=λipvbase,ip/(AipDip)

Lip=1/Aip 

Cip=Aip/(RT),

where Aip represents a cross-sectional area of the pipeline ip in the natural gas network;
establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the head end of the pipeline of the natural gas network:

Gip,u+=Σd=0If−1[Re(Gip,d+)·cos(θd−ωd·uΔt)−Im(Gip,d+)·sin(θd−ωd·uΔt)],

where Re( ) represents valuing a real part of a complex number; Im( ) represents valuing an imaginary part of the complex number; and Od represents a parameter calculated with ωd as follows:

θd=Ifωd−ωd;

establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the tail end of the pipeline of the natural gas network:

Gip,u−=Σd=0If−1[Re(Gip,d−)·cos(θd−ωd·uΔt)−Im(Gip,d−)·sin(θd−ωd·uΔt)]; and

establishing a time domain-frequency-domain mapping constraint of the node of the natural gas network:, h
    
      
        i
        n
      
      ,
      u
    
  
  =
  
    
      ∑
      
        d
        =
        0
      
      
        
          I
          f
        
        -
        1
      
    
    ⁢
    
      [
      
        
          
            Re
            (
            
              h
              
                
                  i
                  n
                
                ,
                d
              
            
            )
          
          ·
          
            cos
            ⁡
            
              (
              
                
                  θ
                  d
                
                -
                
                  
                    
                      ω
                      d
                    
                    ·
                    u
                  
                  ⁢
                  
                    
                  
                  ⁢
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
              
              )
            
          
        
        -
        
          
            Im
            (
            
              h
              
                
                  i
                  n
                
                ,
                d
              
            
            )
          
          ·
          
            sin
            (
            
              
                θ
                d
              
              -
              
                
                  
                    ω
                    d
                  
                  ·
                  u
                
                ⁢
                
                  
                
                ⁢
                Δ
                ⁢
                
                  
                
                ⁢
                t
              
            
            ]
          
        
      
    
  

where hin,u represents a value of a d-th component of a pressure of a node in in the frequency-domain window; and hin,u represents a complex variable to be solved; and
sub-step 5-3 of establishing a topological constraint of the natural gas network, the sub-step 5-3 comprising:
establishing a flow balance constraint of a node of the natural gas network:, ∑
        
          
            i
            p
          
          ∈
          
            Ω
            p
            
              +
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            p
          
          ,
          u
        
        +
      
    
    -
    
      
        ∑
        
          
            i
            p
          
          ∈
          
            Ω
            p
            
              -
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            p
          
          ,
          u
        
        -
      
    
    +
    
      
        ∑
        
          
            i
            c
          
          ∈
          
            Ω
            c
            
              +
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            c
          
          ,
          u
        
        +
      
    
    -
    
      
        ∑
        
          
            i
            c
          
          ∈
          
            Ω
            p
            
              -
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            p
          
          ,
          u
        
        -
      
    
    +
    
      
        ∑
        
          
            i
            l
          
          ∈
          
            Ω
            l
            
              -
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            l
          
          ,
          u
        
        gl
      
    
    -
    
      
        ∑
        
          
            i
            s
          
          ∈
          
            Ω
            s
            
              +
              
                ,
                
                  i
                  n
                
              
            
          
        
      
      ⁢
      
        G
        
          
            i
            s
          
          ,
          u
        
        gs
      
    
  
  =
  0

where Ωp+,in represents a set of serial numbers of pipelines connected to the node in at head ends; Ωp−,in represents a set of serial numbers of pipelines connected to the node in at tail ends; Ωc+,in represents a set of serial numbers of compressors connected to the node in at head ends; Ωc−,in represents a set of serial numbers of compressors connected to the node in at tail ends; Ωs+,in represents a set of serial numbers of natural gas sources connected to the node in; and Ωl−,in represents a set of serial numbers of natural gas loads connected to the node in;
establishing constraints of a pipeline-compressor-node time-domain pressure relationship in the natural gas network:

hip,u+=hin,u,∀ip∈Ωp+,in 

hip,u−=hin,u,∀ip∈Ωp−,in 

hic,u+=hin,u,∀ic∈Ωc−,in 

hic,u−=hin,u,∀ic∈Ωc−,in; and

establishing constraints of a pipeline-node frequency-domain pressure relationship in the natural gas network:

hip,d+=hin,d,∀ip∈Ωp+,in 

hip,d−=hin,d,∀ip∈Ωp−,in; and

step 6 of forming a dynamic state estimation model of the natural gas network by using the objective function of the dynamic state estimation of the natural gas network established in the step 4 and the constraint conditions for the dynamic state estimation of the natural gas network established in the step 5; solving, by using a Lagrange method or an interior point method, the dynamic state estimation model of the natural gas network, to obtain the state vector xu for the dynamic state estimation of the natural gas network at the sampling time point τu; and performing the dynamic state estimation of the natural gas network by considering the dynamic characteristics of the natural gas pipelines to provide sufficient data support for operation and control of an integrated energy system.