Hierarchical control of transient flow in natural gas pipeline systems

PII:

Int. Trans. Opl Res. Vol. 5, No. 4, pp. 285±302, 1998 # 1998 IFORS. Elsevier Science Ltd All rights reserved. Printed in Great Britain S0969-6016(97)00035-X 0969-6016/98 $19.00 + 0.00

Hierarchical Control of Transient Flow in Natural Gas Pipeline Systems J. OSIADACZ ANDRZEJ Warsaw University of Technology, Poland An algorithm for optimal control of a gas network with any con®guration based upon hierarchical control and decomposition of the network is described. Local problems are solved using a gradient technique. The subsystems are coordinated using the `good coordination' method to ®nd the overall optimum. Discrete state equations for the case in which output pressures are treated as elements of the control vector has been formulated. Results of investigations are included. # 1998 IFORS. Published by Elsevier Science Ltd. All rights reserved Key words: Flow control, large scale systems, networks, optimal control, transient analysis.

1. INTRODUCTION The growth of the complexity of gas transmission systems is accompanied by increasing opportunities for more ecient management. Central Control, whose main task is the overall management of the national system has, as the number of compressors increased, recognized the rising importance of fuel usage. Gas compressors stations from a major part of the operational plant on each Transmission System. Their function is to restore the gas pressure reduction caused by frictional pressure losses. The compressors are driven mostly by gas turbines which use natural gas as fuel, taken directly from the transmission pipelines. The compressor unit comprises three main components, a gas generator, a power turbine and a centrifugal gas compressor. The maximum shaft powers of the units range from 5.5 MW to more than 20 MW, the associated fuel consumptions are 2.5 to 5.5 MSCFD equivalent to 8600 to 19,000 pounds per day fuel cost. The value of compressor fuel used in the UK is currently in excess of £30 million per annum; this represents 80% of the total energy costs used by British Gas. According to AGA sources, the operating cost of running the compressor stations represents anywhere between 25% and 50% of the company's total operating budget. Minimizing this fuel usage is a major objective in the control of gas transmission costs. Although the fuel usage of compressors is signi®cant and should be optimized there are other objectives. Above all, central Control must operate the system so that has is supplied where needed, in the quantities needed and at the appropriate pressure. So there would be no point in minimizing compressor fuel usage if we then failed to supply the gas. In a nutshell this is the basic problem of running the transmission system; security of supply versus costs. The security needed must be carefully judged. It is usually expressed as a margin of pressure above the minimum essential at any o€take. This brings us to the topic of linepack. This is the volume of gas held in the transmission system. Whenever pressure exceed the minimum some of this gas can be used as storage. This can have massive bene®ts. Linepack also helps to accommodate the unanticipated variations of demand. Correspondence to: Andrzej J. Osiadacz, Warsaw University of Technology, Institute of Heating and Ventilation, ul. Nowowiejska 20, 00-653 Warsaw, Poland 285

286

A. J. OsiadaczÐControl of Flow 2. OBJECTIVE FUNCTION

This work is concerned with the minimization of operating costs for high pressure gas networks under transient conditions. Depending on the character of gas ¯ow in the system we distinguish steady and unsteady states. The steady state in a gas network is described by systems of algebraic non-linear equations. In steady-state problems, since loads and supplies are not functions of time, an algorithm for optimization determines, once, the structure of the network (i.e. the number of sources, compressors, values and regulators called-units which must be on). In addition, the algorithm must determine the optimal parameters of the operation, namely nodal pressures and ¯ows through branches (pipes). For these reasons, the problem of optimization is formulated as a mixed integer linear programming problem. Alternatively, assuming that the structure is known, an algorithm for steady-state optimization of large gas networks is treated as a non-linear problem with non-linear constraints. The transient optimization is more dicult mathematically than the steady state, but can achieve higher savings. It is not always recognized that the high pressure gas transmission system is never in a `steady state', on the contrary conditions are always transient. Every day, and within each day, demand changes. Change of weather requires tuning of input and compression. Sometimes even the location of demand is altered by switching load from one o€take to another. To optimize such a continuously changing situation is complex. Dynamic optimization requires the use of distributed-parameter models: a partial di€erential equation or a system of such equations. The form of these equations varies with the assumptions made as regards the conditions of operation of the gas pipeline. The necessity of taking into account transient e€ect makes the problem of dynamic optimization very complicated from a computational point of view. For such a problem the optimal parameters of the operation of the system (structure and pressures and ¯ows) are functions of time. The goal of optimization is to minimize the following expression ( )  … tf X K pdj …t† Rj Aj
Qj …t†
ÿ1
dt …1† Iˆ psj …t† t0 j ˆ 1 where: KÐthe number of operating compressors is based on the assumption (Francis, 1981) that the cost of running a compressor is proportional to the integral of horsepower over the control interval ( )  … tf pdj …t† Rj Aj
Qj …t†
ÿ1
dt …2† Ij ˆ ps …t† t0 pdj Ð discharge pressure for jth compressor [Pa] psj Ð suction pressure for jth compressor [Pa] Qj Ð ¯ow through jth compressor [m3/s] Aj, Rj Ð constants for the jth compressor, (Aj=mlpsj/(Zm(mlÿ1)), Rj=(mlÿ1)/ml) t0, tf Ð beginning and ®nal time for function evaluation ml Ð isentropic exponent ZmÐ the compressor (overall) eciency.

3. MATHEMATICAL MODEL OF GAS NETWORK UNDER TRANSIENT CONDITIONS The basic equations describing the transient ¯ow of gas in pipes are derived from an equation of motion (or momentum) and an equation of continuity. In practice the form of the mathematical models varies with the assumptions made as with regards to the conditions of operation of the gas in pipes or networks. Models for gas ¯ow are described by partial di€erential equations

International Transactions in Operational Research Vol. 5, No. 4

287

or a system of such equations. Depending on the degree of simpli®cation with respect to the set of basic equations, the equations may be linear or quite generally non-linear. They may be parabolic or hyperbolic of the 1st or 2nd order. Transient ¯ow through a pipe is described by the following linear di€usion equation (Osiadacz, 1996): @p c2
@2 p ˆ @t A
l
@x2

…3†

where: l ˆ l…p; M† ˆ

2
fF
c2
jMj D
A2
p

1 @p Mˆÿ
l @x M Ð mass ¯ow fF Ð friction factor. A mathematical model of the dynamic properties of a gas network has been elaborated using equation (3) and based upon the generalization of the idea of a node including grid points along pipes, multi-junctions and pipe-ends being treated as o€-takes with demand equal to zero. For the jth node we obtain the following equation (Osiadacz, 1987):    X k dpj 1 dpi dpj ‡
ÿ Sij
…pi ÿ pj † ÿ Mj …4† ˆ dt 3 dt dt iˆ1 where: k Ð number of nodes incident to node jth Sij=1/(liDxi) Vji=Vij=AiDxi/2 Mj Ð (mass ¯ow-load) at jth node. Using the trapezoidal rule, integration of equation (4) between t and t + Dt, (where Dt is the time step) yields a linear equation in pn+1 of the form: ˆ hjj
pn‡1 j

k X iˆ1

hji
pn‡1 ˆ rj ÿ M 0j j

…5†

where: k k X 2 X Vji ‡ Sji hjj ˆ
3 i ˆ 1 Dt
c2 i ˆ 1

Vji ÿ Sij 3
Dt
c2   k X Vji 2 n 1 n
p
p rj ˆ
‡ Dt
c2 3 j 3 i iˆ1

hji ˆ

M 0j ˆ

1
Dt

… t‡Dt t

Mj
dt

For the whole network, equation (5) can be written in the matrix form as H
pn‡1 ˆ R
pn ÿ M dim H ˆn  n; dim R ˆ n  n; dim M ˆ n  1 nÐnumber of nodes.

…6†

288

A. J. OsiadaczÐControl of Flow

Matrices H and R are sparse and symmetric. Adding in sources, compressors (compressor stations), regulators and valves, and assuming the ¯ow through each unit is a positive demand at the inlet node and a negative demand at the outlet, equation (6) can be written as H
pn‡1 ˆ ÿ K
f n ‡ R
pn ÿ M dim K ˆ n  u…u ÿ number of units† 8 < ‡1; kij ˆ ÿ1; : 0;

…7†

if the jth unit has its inlet at node i if the jth unit has its outlet at node i otherwise dim f ˆ u  1

fi Ð the ¯ow through the ith unit. Equation (7) can be written in the form: pn‡1 ˆ Hÿ1
R
pn ÿ Hÿ1
K
f n ÿ Hÿ1
M

…8†

x…k ‡ 1† ˆ A…k†
x…k† ‡ B…k†
m…k† ‡ C…k†
z…k†

…9†

or where: x=p (state vector) m=f (control vector) z=M (input vector) A=Hÿ1
R, B=ÿHÿ1
K, C=ÿHÿ1. Equation (9) is a discrete state equation for a dynamic gas network with the assumption that nodal pressures are elements of the state vector and ¯ows through units are elements of the control vector. Writing equation (7) in the form H1
pn‡1 ‡ H2
pn‡1 ‡ K
fn ˆ R1
pn1 ‡ R2
pn2 ÿ Mn 1 2

…10†

where: p1Ðthe vector of non-outlet node pressures p2Ðthe vector of outlet node pressures dim p1=(n ÿ u)  1, dim p2=u  1 ®nally, we have . ‰H1 ..KŠ ˆ

"

pn‡1 1 f n‡1

#

 n . p ˆ ‰R1 ..0Š
n1 ‡ …R2 ÿ H2 †
pn2 ÿ Mn f

. . ‰H1 ..KŠ ˆ A1 ; ‰R1 ..0Š ˆ B1 ; …R2 ÿ H2 † ˆ C1 dim A1 ˆ n  n; dim B1 ˆ n  n; dim C1 ˆ n  u

…11†

or x…k ‡ 1† ˆ A…k†
x…k† ‡ B…k†
m…k† ‡ C…k†
z…k† where x is partitioned:



 p1 ; m ˆ p2 ; z ˆ M ‡ K
f xˆ f ÿ1 ÿ1 A ˆ Aÿ1 1
B1 ; B ˆ A1
C1 ; C ˆ ÿA1 :

…12†

International Transactions in Operational Research Vol. 5, No. 4

289

equation (12) is a discrete state equation for the case in which output pressures are treated as elements of the control vector. Elements of the state vector are non-outlet node pressures and ¯ows through units. 4. OPERATIONAL CONSTRAINTS

. envelope The centrifugal compressors used on the transmission system have particular characteristics. The operating regime can be expressed by what is known as an `envelope' (Fig. 1). If the increase in pressure produced by a machine is plotted graphically against the ¯ow, four limits emerge which enclose an area in which the compressor can properly run. The limits are de®ned as: . `surge': this is the point at which ¯ow through the compressor becomes so low that reversal of ¯ow can occur which can be damaging to the compressor by causing high stress in the bearings or impeller. If all the points of minimum ¯ow and corresponding pressure ratio are plotted we get the surge line. This the line at the left of the envelope. . `choke': at the opposite end of the diagram, a compressor can reach choke. When the pressure ratio is low there comes a point at which no further increase in ¯ow through the compressor is possible. . `maximum and minimum speed': obviously a compressor can run up to some given maximum speed consistent with machine safety and equally there is a minimum speed line. These four limits de®ne the boundaries of the `envelope' of any compressor or combination of compressors. In the algorithm which is presented, constraints resulting from the boundaries of the `envelope' were replaced by maximum compressor power (compression ratio and ¯ow through compressor). . linepack (pressures at selected nodes of the network cannot drop below a certain value). The structure of the high pressure gas network was represented by a directed graph Gg=(V, E) which consists of a set of nodes V and another set E whose elements are called branches. The nodes represent the pipe connections while the branches represent the pipes themselves, a de®nite direction being assigned to each pipe.

Fig. 1. Typical compressor envelope.

290

A. J. OsiadaczÐControl of Flow

For network analysis, it was necessary to select the following nodes and branches: c supply nodes; sources and supplying storages Vz WV, c pressure nodes; nodes at which constraints on the pressures were imposed Vw W V, c units; Ep WE; Ep 0 (Vs, Vd); Vs Ð suction nodes, Vd Ð discharge nodes, c units of the system; U 0 Vz[Ep. Finally, inequality constraints were imposed on: . maximum source ¯ow qj E qmax j where: j $ Vz, . maximum compressor ratio pdj  Emax j psj . minimum compressor ratio pdj  Emin j psj where: j $ Ep, . minimum pressure at selected nodes pjepmin j where: j $ Vw, . maximum horsepower of the compressor (  ) pd j R j ÿ1  Wmax Wj ˆ Aj Qj j psj where: j $ Vw. (In general compressors are also constrained to operate within an envelope described by four non-linear constraints. These are complicated to program, especially with derivatives and have not yet been implemented.) . range of the pressure variations at control nodes max pmin zj Epzj Epzj

where: j $ Vc. The operational constraints together with state equation [equation (12)] form a complete set of constraints imposed on the high pressure gas network.

5. METHOD OF SOLUTION To solve the above problem, an hierarchical method has been used. An algorithm for the optimal control of gas network with any con®guration based upon hierarchical control and decomposition of the network has been developed. Local problems were solved using a gradient technique. Implicit in all of hierarchical system theory is the idea that it is generally easier to deal with several low order systems than with one system of high order. Since the gas transmission system has no apparent natural hierarchical structure, the decomposition should be carried out in whatever manner seems appropriate. One possibility is to employ spatial decomposition according to which the network is divided into physically small subsystems. On this decomposition is imposed one constraint, namely each subsystem has to contain at least one operating compressor. A sec-

International Transactions in Operational Research Vol. 5, No. 4

291

ond possibility is decomposition in time. In this case a subsystem would be the complete system in one time step of optimization. Time decomposition has several disadvantages: 1. for most normal demand conditions on a gas network the control interval is subdivided into time steps which do not exceed 1 hour. For an overall time interval of 24 hours this could give at least 24 subsystems. 2. the coordination will involve a large number of variables and hence much computation time. 3. local optimization of each subsystem is a dicult task. For these reasons a spatial decomposition has been chosen to solve the given problem. After decomposition each subsystem is interconnected with some other subsystems using simple permutations: ai …k† ˆ

N X jˆ1

Lij
bj …k†

…13†

aiÐinput vector into ith subsystem. (Nodal pressures at nodes incident to the other subsystems) LÐmatrix of interconnections biÐoutput vector from jth subsystem. (Nodal pressures at nodes of jth subsystem incident to the ith subsystem) NÐnumber of subsystems. Vector b is calculated from the following matrix equation: bi …k† ˆ

N X jˆ1

Gij
xj …k†

…14†

GÐmatrix of interconnections. Discretization of the cost function yields: Iˆ

kX f ÿ1

F1 ‰x…k†; m…k†; kŠ

…15†

k ˆ k0

where: xÐstate vector (non-outlet node pressures and ¯ow through units), mÐcontrol vector (outlet node pressures), kfÐdesired ®nal stage. Constraints are taken by introducing an Augmented Lagrangian function. Thus Iˆ

kX f ÿ1 k ˆ k0

F‰x…k†; m…k†; x1 …k†; kŠ

…16†

where: F ˆ F1 ‡ m…k†T
u…x1 …k†† ‡

c…k†
u…x1 …k††T
S…k†
u…x1 …k†† 2

…17†

mÐLagrange multiplier u(x1(k))Ðcontains the constraints that are active at x1 c(k)Ðpenalty factor. The Hamiltonian for the integrated system is: H…k† ˆ F‰x…k†; m…k†; x1 …k†; kŠ ‡ lT …k ‡ 1†
f‰x…k†; m…k†; z…k†; kŠ: In terms of the subsystems, the Hamiltonian may be written as (Singh et al., 1978)

…18†

292

A. J. OsiadaczÐControl of Flow N X

H…k† ˆ

iˆ1

Fi ‰xi …k†; mi …k†; xli …k†; kŠ ‡ lT …k ‡ 1†
fi ‰xi …k†; mi …k†; zi …k†; kŠ:

‡pTi …k†

N X jˆ1

!! Lij
bj …k† ÿ ai …k†

…19†

where: li, piÐmultipliers. Because pTi …k†


N X j ÿ1

Lij
bj …k† ˆ

N X N X iˆ1 jˆ1

pTi …k†
Lji
bj …k†:

…20†

equation (19) can be written as follows: H…k† ˆ

N X iˆ1

…Fi ‰xi …k†; mi …k†; xli …k†; kŠ ‡ lT …k ‡ 1†
fi ‰xi …k†; mi …k†; zi …k†; kŠ

‡ KTi …k†
bi …k† ‡ bTi …k†
ai …k††

…21†

where: Ki …t† ˆ

N X jˆ1

Lij
pj …k†:

Thus for the ith subsystem we have: Hi …k† ˆ Fi ‰xi …k†; mi …k†; xli …k†; kŠ ‡ lT …k ‡ 1†
fi ‰xi …k†; mi …k†; zi …k†; kŠ ‡ KTi …k†
bi …k† ÿ pTi …k†
ai …k†

…22†

Then the optimization problem for the ith subsystem is the following: minimize with respect to mi Ii ˆ

kX f ÿ1 iˆ1

‰Fi …xi …k†† ‡ KTi …k†
bi …k† ÿ pTi …k†
ai …k†Š

…23†

subject to the constraints. We set dHi(k)/dxi(k)=li(k) to obtain the adjoint equation li …k† ˆ

      @bi …k† T @Hi …k† @Fi …k† @fi T @ai …k† T ˆ ‡
li …k† ‡
Ki ÿ
pi @xi …k† @x…k† @xi …k† @xi …k† @xi …k†

…24†

where the terminal condition of the adjoint equation is: li …kf † ˆ 0:

…25†

6. LOCAL OPTIMIZATION The computation process is as follows. *Guess or determine mi(k), k = k0, . . ., kf and then determine xi(k + 1) from xi ˆ fi ‰xi …k†; mi …k†; zi …k†; kŠ

…26†

xi …k0 † ˆ xi;0 :

…27†

International Transactions in Operational Research Vol. 5, No. 4

293

Hierarchical simulation of gas networks based on decomposition±coordination techniques has been used to solve the set of equations. Basically, in this method the original set of equations is decomposed into l independent subsystems referred to as local problems which can be solved independently given an initial guess value for the coordinator problem. Then, using the Newton±Raphson method, the coordinating variables are updated and the local problems are solved based on these values until accuracy criteria are satis®ed. cThe problem of simulation to be solved at each stage consist of a set of algebriac equations of the following form: F…Z† ˆ 0

…20†

where Z $ Rnz and F: Rnz 4Rnz. Now consider the case when the above system of equations is decomposed into l subsystems jointed together via interconnecting variables Yi, viz. Fi …Xi ; Yi † ˆ 0; i ˆ 1; 2; 3; . . . ; l

…29†

It is assumed (Osiadacz et al., 1995) that, given Yi, the local problem governed by (29) may be solved for Xi(Yi) independently for i = 1, 2, 3, . . . , l. Then the coordinator ®nds Y=(Y1, Y2, . . ., Yl) such that Yi ˆ Ai
Xiÿ1 …Yiÿ1 † ‡ Bi
Xi‡1 …Yi‡1 †

…30†

YˆC
X

…31†

or, for all i's,

where the `structural' matrix C is 2 0 6 A2 6 6 A3 6 Cˆ6 6... 6... 6 4 0

a block-tridiagonal matrix with the following structure: 3 B1 7 0 B2 0 7 0 B3 ...7 7 ... ... ... ... ... ... ...7 7; 7 ... ... ... ... ... ... 7 5 Alÿ1 0 Blÿ1 A1 0 tx

where A1=B1=0. cThe solution of the coordinator problem as given by the set of equation (31) may be obtained in a di€erent way by de®ning D…Y†  Y ÿ C
X…Y†:

…32†

Now, assuming that D(Y) is a di€erentiable function of Y and that (@D(Y)/@Y) is nonsingular, at any iteration step the coordinator updates the value of Y based on a Newton±Raphson procedure:

where

Yk‡1 ˆ Yk ÿ …D…Yk †=D 0 …Yk ††

…33†

@D…Yk † @X…Y† D …Y † ˆ ˆIÿC
: @Y YˆYk @Y YˆYk

…34†

0

k

In order to compute the Jacobian (@D(Y)/@Y), it is required to ®nd @X…Y† ˆ block diag …@Xi …Yi †=@Yi †: @Y By di€erentiating equation (29) with respect to Yi it can be shown that:

…35†

294

A. J. OsiadaczÐControl of Flow   @XXi …Yi † @Fi p 41 @Fi ˆÿ Xi …Yi †; Yi
: @Xi @Yi Xi …Yi †; Yi @Yi

…36†

The matrices (@Fi/@Xi)ÿ1 and (@Fi/@Yi) are both available at the local level and therefore can be used to compute the (@Xi/@Yi) matrices. It is important to note that it is not necessary to transfer the entire (@Xi/@Yi) matrices from the subsystems to the coordinator, but only the required elements of the product (C@X/@Y) contributed by each system. * Solve the adjoint equation, equation (24) backwards from stage kf with the terminal condition of equation (25) to stage k0. * Calculate   @Hi …k† @Fi …k† @fi …k† T ˆ ‡
li …k ‡ 1† …37† @mi …k† @mi …k† @mi …k† which will not normally be zero. The perturbation in the cost function with perturbations in xi(k) and mi(k) may be obtained by rewriting equation (23) as (Sage and White, 1977) Ii ˆ

kˆ kf ÿ1 X k ˆ k0

‰Hi …k† ÿ lTi …kf †
fi …k†Š ˆ

kˆ kf ÿ1 X

ˆ lTi …k0 †
wi …k0 † ÿ lTi …kf †
xi …kf † ‡

k ˆ k0

‰Hi …k† ÿ lTi …k ‡ 1†
fi …k ‡ 1†Š

kˆ kf ÿ1 X k ˆ k0

‰Hi …k† ÿ lTi …k†
xi …k†Š:

…38†

Thus we obtain: DIi ˆ lTi …k0 †
Dxi …k0 † ÿ lTi …kf †
Dxi …kf † ( ) T   kX ˆ kf  @Hi …k† @Hi …k† T T ÿ li …k†
Dxi …k† ‡
Dmi …k† : ‡ @xi …k† @mi …k† kˆk

…39†

0

Finally, we have: DIi ˆ

kˆk Xf  k ˆ k0

 @Hi …k† T
Dmi …k†: @mi …k†

…40†

If we want to make the longest change in DI we have to calculate @H(k)/@m(k) and then make Dm(k) directed opposite to the gradient, i.e. the steepest descent method. In this case a conjugate gradient method to calculate direction is used because, we believe, it should given a better rate of convergence. A search direction is computed by a conjugate gradient technique according to the formulas: . step 1 at M0 calculate I where: M=[mik]i = 1, 2, . . ., q; k = 1, 2, . . ., r qÐ the number of compressor stations rÐ the number of discrete stages of optimal control process mikÐ discharge pressure of ith compressor station at kth stage. Let s0=ÿgk where: gk ˆ @Hi …k†=@mi …k†

International Transactions in Operational Research Vol. 5, No. 4

295

. step 2 calculate m1k=m0k+t
s0k by minimizing equation (16) with respect too t in the s0 direction. To ®nd a minimum of objective function in that search direction, quadratic interpolation is applied. For the jth iteration, the relation is the following (Bazarra et al., 1993): j‡1 j s j‡1 k ˆ ÿg k ‡ s k


j‡1 j T …g j‡1 k †
…g k ÿ g k †

…g jk †T
g jk

:

…41†

The procedure is repeated until either the control vector or the cost function does not change `too much' from iteration to iteration. Two sided constraints of the form gj  gj …x†  dj

…42†

where gj and dj are some scalars equivalent to the inequalities gjEgj(x) ÿ ujEdj, uj=0, j = 1, . . ., r (Bertsekas, 1982). The multipliers m jk are updated by means of the iteration: m jk‡1 ˆ m jk ‡ ck u jk

…43†

where u1k, . . ., urk together with a vector xk solve problem (16). Minimization in problem (16) can be carried out ®rst with respect to uj yielding the equivalent problem: minimize F1 ‡

r X jˆ1

pj ‰uj …x†; mj …k†; c…k†Š

…44†

subject to x $ Rn where:  pj ‰uj …x†; mj …k†; c…k†Š ˆ

min

al  gj …x†  bl

mj …k†
uj …x…k†† ‡

 c…k†
uj …x…k††
Sj …k†
uj …x…k†† : 2

…45†

A straightforward calculation shows that the minimum above is attained at the point uj(x(k)) given by: 8 if mj …k† ‡ c…k†
‰gj …x† ÿ dj Š > 0 < gj …x† ÿ dj if mj …k† ‡ c…k†
‰gj …x† ÿ gj Š > 0 gj …x† ÿ gj …46† uj …x…k†† ˆ : ÿm …k†=c…k† otherwise j and pj is given by: pj ‰uj …x…k††; mj …k†; c…k†Š ˆ 8 m …k†
‰gj …x† ÿ dj Š ‡ 12
c…k†
j gj …x† ÿ dj j2 > > < j mj …k†
‰gj …x† ÿ gj ‡ 12
c…k†
j gj …x† ÿ gj j2 > > : ÿ…mj …k††2 =2
c…k†

if mj …k† ‡ c…k†
‰gj …x† ÿ dj Š > 0 if mj …k† ‡ c…k†
‰gj …x† ÿ gj Š < 0

…47†

otherwise:

It is easy to see, that if gj $ C1, then pj[uj(x(k)), mj(k), c(k)] is continuously di€erentiable in x. If gj $ C2, then pj(uj(x(k)), mj(k), c(k)] is twice continuously di€erentiable on the {xvmj(k) + c(k)
[gj(x) ÿ dj]$ 0, mj(k) + c(k)
[gj(x) ÿ gj] $0}. The form of the function pj is shown in Fig. 2. The preceding analysis shows that a method of multipliers for problem (16) consists of sequential minimizations of the form (44), (47), which do not involve the variables u1, . . ., ur.

296

A. J. OsiadaczÐControl of Flow

Fig. 2. ?

The multiplier iteration is given by: 8 < mj …k† ‡ c…k†
‰gj …x…k†† ÿ dj Š mj …k ‡ 1† ˆ mj …k† ‡ c…k†
‰gj …x…k†† ÿ gj Š : 0

if mj …k† ‡ c…k†
‰gj …x…k†† ÿ dj Š > 0 if mj …k† ‡ c…k†
‰gj …x…k†† ÿ gj Š < 0 otherwise

…48†

where x(k) solves problem (44).

7. COORDINATION PROBLEM To ®nd the overall optimum the coordination variables (p) for optimization are introduced. The subproblems treat the coordination variables as known inputs which remain ®xed until the coordinator supplies new values. The following theorem relates the subproblems to the integrated problem. TheoremÐIf a solution exists to the integrated problem and to each of the subproblems, then there exists a p*(k) such that solutions satisfying the necessary conditions for the subproblems also satisfy the necessary conditions for the integrated problem. The purpose of the coordination variables is to ensure satisfaction of the interconnection constraints, for only when these constraints are satis®ed can it be claimed that the subproblem solutions also solve the overall problem. It can be seen that the main problem of coordinating subproblem solutions so that they solve the integrated problem in the determination of l*. Since adjusting p adjusts the subproblem objective functions or goals this is a goal coordination method. If p$ p*, then the interconnection constraints are not satis®ed, so that it is rational to examine the e€ect of p on the interconnection error, which will be de®ned as: 3 2 3 j j 7 6 j 7 j 7 6 6 7 7 6 6 j 7 6 j 7 6 6 7 7 6 6 j 7 j 7 6 PN 7 6 6 j ˆ 1 Lij
b …P, k† ÿ ai …P, k† 7 6 7 ˆ e …P; k† 7 6 i 6 j 7 7 6 6 7 j 7 6 6 j 7 7 6 6 7 j 7 6 6 j 7 7 4 6 5 j 5 4 j j j 2

…error vector†

International Transactions in Operational Research Vol. 5, No. 4

297

For updating p(k) to minimize ei(p, k) a gradient procedure to produce a new p is applied: pn‡1 …k† ˆ pnj …k† ‡ knj
j

@H…k† : @Pj …k†

…36†

Thus, the original problem has been exchanged for N ( j = 1, 2, . . ., N) subproblems.

8. TEST NETWORK The correctness of the elaborated algorithm for optimization was checked using part of the National Grid of UK shown in Fig. 3. The algorithm has been tested using a network containing 23 nodes, 13 pipes, 3 compressor stations, 2 storage supply nodes and 1 source. There are also 15 demands in the network (Fig. 3).

Fig. 3. Structure of the gas network.

298

A. J. OsiadaczÐControl of Flow

Fig. 4. Changes of load at selected nodes.

International Transactions in Operational Research Vol. 5, No. 4

299

1. Pipe data: Sending node 1 2 4 5 7 9 8 8 11 12 14 13 15 1 17 18 20 21 22

Receiving node 2 3 5 6 8 8 10 11 12 13 13 15 16 17 18 19 21 22 23

Length (m) 32,991 34,279 88,675 30,095 27,037 13,197 10,461 24,784 1,770 1,770 19,795 14,001 19,312 64,856 53,591 13,679 13,358 28,163 6,437

Fig. 5. The variations of ¯ow through compressor station 3±4.

Fig. 6. The variations of ¯ow through compressor station 5±7.

Diameter (mm) 884 884 884 584 884 439 584 884 732 732 584 732 732 884 884 884 884 884 732

300

A. J. OsiadaczÐControl of Flow

Fig. 7. The variations of ¯ow through compressor station 19±20.

Fig. 8. The variations of total power of compressor stations.

2. Unit data: Compressor stations: 3±4 Pmin=3.45 [MPa] Pmax=7.36 [MPa] Emax=1.8 5±7 Pmin=3.45 [MPa] Pmax=5.83 [MPa] Emax=1.8 19±20 Pmin=3.45 [MPa] Pmax=6.91 [MPa] Emax=1.8 Supply nodes (node 1) Pmin=3.45 [MPa] Pmax=6.72 [MPa]

Qmax=985 [Nm3/s]

Storage supply nodes: Node 9 Q = 2.6 [Nm3/s] = const. Node 14 Q = 32.5 [Nm3/s] = const. ÐPeriod of optimization: TM=25 h Ðdiscretization step: Dt = 7200 s, Dx = 16,000 m The loads at nodes 2, 6, 8, 9, 10, 13, 14, 15, 16, 17, 18, 21, 22 and 23 varied in accordance with the curves presented in Fig. 4. The inequality constraints are imposed on:

International Transactions in Operational Research Vol. 5, No. 4

301

Fig. 9. The variations of pressure at node 7. 1Ðlower constraint; 2Ðupper constraint; 3Ðoptimal pro®le of the pressure at node 7.

pressure at selected nodes: p6 E3.5 [MPa] p16 E4.1 [MPa] p23 E4.62 [MPa] and 3.45 E p1 E6.72 [MPa] 3.45 E p4 E7.36 [MPa] 3.45 E p7 E5.83 [MPa] 3.45 E p20 E6.91 [MPa] Selected results of optimization are given in Figs 5±9.

9. CONCLUSIONS The algorithm described above has been implemented in double precision FORTRAN on PC486 computer and some preliminary testing has been done. Variables and constraints are scaled to be of order unity. To deal with the hierarchical algorithm the network was split into three subsystems. Investigations have shown that the developed algorithm works properly although certain modi®cations should be done to speed up calculations and to improve ¯exibility. The validity of the solution was checked by composing solutions with those from simulation programs run with appropriate control values, it was veri®ed that the equations of gas ¯ow were satis®ed. An algorithm described in (Furey, 1993) has been developed for dynamic optimal control of complex gas networks. The method is based on Sequential Quadratic Programming and takes account of the structure of the pipe¯ow equations by means of a reduced gradient technique which eliminates most of the variables from the quadratic subproblems. The latter involve only simple bound constraints, which are handled eciently by a conjugate gradient active set algorithm. Trust region techniques permit use of the exact Hessian, preventing sparsity. More general constraints are handled at an outer level by a truncated augmented Lagrangian method. An algorithm has been tested using the same input data (southern part of the National Grid of UKÐFig. 3 and Fig. 4). The results are similar. Maximum discrepancy does not exceed 15%. It is clear, however, that large dynamic problems over 24 hours cannot be solved exactly on one processor computer in reasonable time. One approach to obtaining fast real time solutions to large problems involves the use of parallel processing computers. Another approach is to reduce the scope of the problem, so that one aims for a good solution rather than an optimal one. It is very important to notice that

302

A. J. OsiadaczÐControl of Flow

the method based on decomposition-coordination is very suitable for parallel computation and hence implementation on parallel processes. REFERENCES Bazaraa, M. S., Sherali, H. D. and Shetty, C. M. (1993) Nonlinear Programming (second edition). John Wiley & Sons, Inc. Bertsekas, D. P. (1982) Constrained Optimization and Lagrange Multiplier Methods. Academic Press. Francis, R. F. (1981) The Ecient Management of the Bulk Transmission of Gas. In The Institution of Gas Engineers. Communication 1144, Birmingham. Furey, B. P. (1993) A sequential quadratic programming based algorithm for optimization of gas networks. Automatics 29 (6). Osiadacz, A. J. (1987) Simulation and Analysis of Gas Networks. E. & F. N. Spon, London. Osiadacz, A. J. and Swierczewski, S. (1994) Optimal control of gas transportation systems. In The 3rd IEEE Conference on Control Applications, pp. 795±798. The University of Strathclyde, Glasgow. Osiadacz, A. J. and Bell, D. J. (1995) Dynamic simulation of gas networks by decomposition and coordination. Math. Engng. Ind. 5 (3), 235±254. Osiadacz, A. J. (1996) Di€erent transient modelsÐlimitations, advantages and disadvantages. In PSIG Conference, 1996. San Francisco. Sage, A. P. and White, C. C. 1977. Optimal Systems Control. Prentice-Hall, New Jersey. Singh, M. G. and Titli, A. (1978) SystemsÐDecomposition, Optimisation and Control. Pergamon Press, Oxford.