- Email: [email protected]
Mathematical modeling and computational principles for the analysis and simulation of long-distance energy systems
WEMATICS AND coh4lwTERS N SIMULATION ELSEVIER
Mathematics and Computers in Simulation 39 (1995) 565-572
Mathematical modeling and computational principles for the analysis and simulation of long-distance energy systems Kurt Schlacher
*, Andreas
Kugi
Johannes Kepler University Linz, Institute for Automatic Control and Electrical Drives, Altenbergerstrafle 69, A-4040 Linz, Auhof; Austria
Abstract Methods for the modeling and simulation of long-distance energy systems are considered. Due to Kirchhoff’s laws a special type of non linear equation is crucial for the steady state and the transient analysis. A reliable algorithm to solve these equations based on Newton’s method is presented. The uniqueness of the solution and the convergence of the method are proved by the stability theory of Liapunov. To apply this method, the calculation of some derivatives is necessary. An approach for C+ + to do this in an automatic way without rewriting a program is presented.
1. Summary This paper presents methods for the analysis and simulation of long-distance energy systems. Essentially two problems make the analysis and the modeling of such systems very difficult. The first one is the nonlinearity of the equations, which describe the steady state and the dynamical model. The second one follows from the size of the mathematical model, which can only be set up and solved by computer aided methods. In order to get an efficient solution to this type of systems of nonlinear equations a reliable algorithm based on a modified Newton method is introduced. The unambiguity of the solution and the convergence of the procedure is proved by the help of Liapunov functions. Since realistic models of long-distance energy systems consist of many equations, the computer is not only used to solve these nonlinear equations, but also to derive them automatically. The system describing a network can also change during simulation according to switching terminals. Therefore a simple computer aided algorithm based on graph theory is used to set up the network equations.
* Corresponding
author.
037%4754/95/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDI 0378-4754(94)00119-7
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
566
The partial derivatives which are needed for solving the equations are calculated in a symbolic way. This method is based on finite taylor series and special features of the programming language C + + .
2. Introduction Long-distance energy systems are growing fast because of ecological improvements. The size, the complexity and the nonlinearity of the hydraulic equations command the modeling and simulation of these systems. Since a computer must be used to integrate the differential equations, it is natural to use it for setting up the equations, too. The propagation of temperature and pressure are the most important phenomena in these systems. One builds the network for the first effect but a big part of the control equipment is needed to handle the second one. The wave speed of a water hammer is about 1000 [ms-l]. The delay due to the finite propagation velocity is not negligible any more in energy systems for whole cities. Simulation based on mathematical modeling and reliable algorithms offer a good way in understanding the phenomena of those systems.
3. Hydraulic networks A hydraulic network is formed by connecting together terminals like pumps, valves, pipes, etc.. The connection points are called nodes and the terminals connecting two nodes are called branches. This type of network obeys the two laws of Kirchhoff. The current law says that the total mass flow to a node is zero. The density of water may be assumed to be constant in the interesting range of temperature and pressure. Therefore one can replace the mass flow by the volume flow. The voltage law says that the difference of the pressure of two nodes connected by a branch is equal to the difference pressure of this branch. Now the most important branches are presented.
3.1. Pipes The one-dimensional transient flow in pipes can be described by a pair of quasi linear hyperbolic partial differential equations of first order 0
1 0
f
-~xh(X, t) +-~-;q(x, t) + 2gDA---------~lq(x,t)lq(x, t)=O 0
a2
and
0
-~h(x, t) + -~ O---xO(X,t)=0. h(x, t)[m] denotes the pressure head, q(x, t)[mas -1] the volume flow, A the area of the pipe, D the pipe diameter, g the gravitational acceleration, f the friction factor and a the equivalent wave speed [7]. The boundary conditions of a pipe of length l at a time t are of the form a
h(l,t)=--~q(l,t)+R1,
a
h(O,t)=--~q'(O,t)+R 2
and
q'(O,t)=-q(O,t)
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
567
with R 1 and R z depending only on values for "r < t. These conditions describe a 1-port which consists of a linear port and a controlled pressure source in series. The minus sign for q' is used to have the same flow reference direction on both end points. The boundary conditions of a pipe can be embedded in the theory of 1-port networks [4].
3.2. Pumps In long-distance energy systems only rotary pumps are used. As a good approximation the transient behavior of a pump is described by the steady state characteristic curves of the total head increase dh = dh(q, to) and the shaft torque re(q, to) with to as the angular velocity. The torque equation takes the form d
0-~oo = md(to) - m(q, w) with 0 as the moment of inertia and md(to) the torque of the driving motor [7]. If necessary the dynamics of the driving motor is taken into consideration, too.
3.3. Valves The difference pressure of the valve is given by
Iqlq dh =
-r
and
d -~a = f ( q , a, u)
with 0 ~
with r as the resistance factor, a as the grade of opening and u as the control input. The topology of the hydraulic network changes at each time, when a switching component like a valve is opening or closing [7].
3.4. Steady state and transient analysis A large system of nonlinear equations has to be solved in the case of steady state analysis. The transient behavior is described by systems of partial and ordinary differential equations with nonlinear equations as restrictions. The latter ones are a consequence of Kirchhoff's circuit laws. Investigations about the computing time have shown, that the solution of this system of nonlinear equations takes much more time than the integration of the differential equations. These circuit equations are described by the theory of 1-port networks most efficiently. They are of the same type for the steady state and the transient analysis, therefore only the first one is treated. It is very important to use efficient algorithms to solve this problem in order to obtain a short computing time. A fast algorithm will be presented after the uniqueness of the solution of the hydraulic network equations is shown.
4. Facts of networks
Only the steady state equations are considered according to the above-mentioned observations. The circuit is formed by connecting together special terminals. A graph G = {N, B} is
568
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
related to the network. N is the finite set of nodes of cardinality n and B is the finite set of branches. A branch has exactly two end points which must be nodes. T h e n u m b e r of branches is b. T h e flow qi and the difference pressure d h i are assigned to the branch i. T h e pressure hj is related to a n o d e j. Now the laws of Kirchhoff can be expressed as hid i+hkd~=dh with
i{i
dt -
,
and
~-'~d~qJ=O J
if n o d e i and branch 1 are connected, the direction of the flow to i is positive, if n o d e i and branch l are connected, the direction of the flow off i is positive, otherwise.
T h e flow q of a network is a point q = ( q l , . . . , qb) ~ Rb. q is said to be admissible if q satisfies Kirchhoff's current law. A difference pressure dh is a point dh = ( d h l , . . . , dh b) ~ (Rb) * . dh is said to be admissible if dh satisfies Kirchhoff's voltage law. Let p be the natural bilinear m a p p . ( ~ b ) * X ~b ~
p ( d h , q ) = ~_,dhjq j J then Tellegen's t h e o r e m can be expressed in the following form [5]. Theorem 1. For all admissible q and dh is p(dh, q)=0. T h e linear m a p D" Rb .__,Rn X i=
Ed~q j J
is well defined. For a basis
{a i} of
Ker D there exists a linear m a p A : R c ~ Ker D
qi= Ea x j J which is bijective for some c. Tellegen's t h e o r e m leads now to the following results [2]. Theorem 2. A flow q is admissible iff Dq=O
or
q=Ax.
(1)
A difference pressure dh is admissible iff dh ~ Im D*
or A *dh = O.
(2)
It is a well-known fact that one can always find a basis of Ker D that x j is the flow qi of a branch i.
K. Schlacher,A. Kugi/ Mathematicsand Computersin Simulation39 "(1995)565-572
569
5. A potential function Only 1-ports are considered as branches for this type of network. Then there is a branch equation of the form: (a) dh i =fi(qi), (b) dhi = const., or (c) qZ = const. for each branch i. For the sake of simplicity only equations of type a) are considered, the extensive case is presented in [4]. The short form of the branch equations is
dh = f ( q ) .
(3) A flow q is said to be a solution of the network equations if q and dh(q) are admissible (equations (1), (2)) and satisfy the equation above or q=Ax
A*f(q)=O
and
(4)
hold. Based on Tellegen's theorem we define the function q
qi
.
For an admissible flow the derivative of v with respect to x ~ is given by
ai OX
v(Ax)= ~O ax ) ( x ) qU(AX)~x/(A
~--- --
~ fj ( ~ , a j x l ) a [ = - p ( f ( A x ) , •
ai).
l
Each solution of the network satisfies the equation a
7x(Ax)=O.
Now it is possible to give the main result of this paper. Theorem 3. If the branch equations fi (3) of the network satisfy the conditions f i e c l ( - ~ , oo)
and d then the solution q of the network equations (4) is unique. To prove this theorem first the functions qi
gi(qi) = -- fo f i ( q ) d x ' are considered. It can be easily shown that gi is strictly convex [3] and that the condition limlqil _,®gi(qi) = oo is satisfied, v(q) is strictly convex, since v(q) is the sum of strictly convex functions and the condition limllql I ~ v ( q ) =
570
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
is satisfied, too. Because of Ker A = {0} the restriction of v(q) to v(Ax) is strictly convex and lim4rxjI __,oov(Ax)= oo
holds. It is a well-known fact of the theory of convex functions that v(Ax) has a global minimum and x is the unique solution of ~ v ( A x ) = 0. This completes the proof. It must be emphasized that the conditions of Theorem 3 are no restriction for hydraulic networks.
6. A modified Newton method
Due to the theory of gradient systems and the present considerations a simple algorithm to solve the network equations is given by d - ~ x i = o ( f ( A x ) , ai). In general a faster algorithm can be constructed using Newton's method [3]. This algorithm is applicable since v(Ax) ~ C 2 and v has a unique minimum. Let Jv be the Jacobian and Hv the Hessian of v. Then the sequence
x(i + 1 ) = x ( i ) - a ( x ( i ) ) H ~ l ( x ( i ) ) J f ( x ( i ) )
(5)
with the tuning function a(x(i)) is well defined. The Hessian of v(Ax) is positive definite since v(Ax) is strictly convex. In this special case the Hessian H~ is given by
H~(Ax) = A * ( - J f ) A with JI as Jacobian of f of equation (3). Expanding v into a Taylor series the equation
v(x(i + 1 ) ) = v ( x ( i ) ) + Jv(x(i))(x(i + 1 ) - x ( i ) ) + O( [[x(i + 1 ) - x ( i ) [[ z) follows for sufficiently small [[ x(i + 1) -x(i)[[. With the help of
v(x(i + 1)) - v ( x ( i ) ) = - a ( x ( i ) ) J . ( x ( i ) ) H ~ ' ( x ( i ) ) J f (x(i)) + O(aZ(x(i))) one can ensure that the inequality
v(x(i+ 1))-v(x(i)) 0. If the process of minimizing is considered as a (time) discrete process the theory of Liapunov gives an easy proof of the next theorem [6]. Theorem 4. The Newton series (5) converges for each initial value to the unique solution of the network equations (4) for a sufficiently small choice of the tuning function a. If v ~ C 3 the order of convergence is at least two provided that the initial value is sufficiently close to the solution.
K. Schlacher,A. Kugi/ Mathematics and Computersin Simulation 39 (1995)565-572
571
This theorem ensures not only the convergence it offers a way to choose a by checking the difference v(x(i + 1)) - v(x(i)), too.
7. An algorithm for calculating derivatives The set of the network equations is too large to be set up by hand. It changes also during a simulation according to the action of switching terminals. It can be shown that the graph searching algorithm "depth first" offers an efficient way setting up the equations in an automatic way [4]. Although this algorithm can be used to calculate all required derivatives for the Newton procedure, too, a more general method is presented. This approach needs a modern programming language like C + + which allows overloading of functions of binary operations and dynamic memory management [1]. First a new type of variable x = (x v, (xi)) is defined as a pair of a float variable x~ and an array of float variables (x i) of dimension one. Next the binary operations x+y=(xv+Yv,(xi+Yi)),
x'y=(xvY,,(YvXi+XvYi))
and the rule
f(x)= (f(xv), ( d for the value of the real function f are introduced. Let x b e a finite Taylor series of the form 0
x
x(z),
The value of the variable x at z is defined as
x= x(z),--x(z)l. az~ ] It is easy to verify that these operations implement the rules to handle this special kind of series. If these operations are used in a program all required derivatives are calculated in an automatic way. The advantage of this approach is that only a new type of variable and overloading of some operations is required. It is not necessary to write a new program.
8. Results The above introduced methods are the basis of the programming system F L U I D I X for the simulation of long-distance energy systems. This package consists of a graphic editor, a data base of mathematical models, an interactive simulator and a graphic post processor. The user only draws the plan of the circuit assisted by mouse and menus. Parameters of the models are added by the user or by the data base. The system derives the equations for the steady state or
572
K. Schlacher, A. Kugi /Mathematics and Computers in Simulation 39 (1995) 565-572
transient analysis automatically. During the simulation the user can trigger events with the mouse or keyboard. The results of a simulation are d o c u m e n t e d by the graphic post processor. The programming system F L U I D I X has been used for the analysis of long-distance energy systems of cities with more than 250.000 inhabitants. Therefore one can say that the methods presented in this p a p e r are applicable for real world problems.
References [1] K.E. Gorlen, S.M. Orlow and P.S. Plexico, Data Abstraction and Object-Oriented Programming in C+ + (Wiley, 1991). [2] M. Hirsch and S. Smale, Differential Equations, Dynamical Systems and Linear Algebra (Academic Press, 1974). [3] D.G. Luenberger, Linear and Nonlinear Programming (Addison Wesley, 1989). [4] K. Schlacher, Systemanalytische Methoden f'tir hydraulische Netze, Habilitationschrift an der TU-Graz, 1990. [5] P. Penfield, R. Spence Jr and S. Diunker, Tellegen's Theorem and Electrical Networks (MIT Press, 1970). [6] M. Vidyasagar, Nonlinear Systems Analysis (Prentice Hall, 1993). [7] E.B. Wylie, V.L. Streeter and L. Suo, Fluid Transients in Systems (Prentice Hall, 1993).
Mathematics and Computers in Simulation 39 (1995) 565-572
Mathematical modeling and computational principles for the analysis and simulation of long-distance energy systems Kurt Schlacher
*, Andreas
Kugi
Johannes Kepler University Linz, Institute for Automatic Control and Electrical Drives, Altenbergerstrafle 69, A-4040 Linz, Auhof; Austria
Abstract Methods for the modeling and simulation of long-distance energy systems are considered. Due to Kirchhoff’s laws a special type of non linear equation is crucial for the steady state and the transient analysis. A reliable algorithm to solve these equations based on Newton’s method is presented. The uniqueness of the solution and the convergence of the method are proved by the stability theory of Liapunov. To apply this method, the calculation of some derivatives is necessary. An approach for C+ + to do this in an automatic way without rewriting a program is presented.
1. Summary This paper presents methods for the analysis and simulation of long-distance energy systems. Essentially two problems make the analysis and the modeling of such systems very difficult. The first one is the nonlinearity of the equations, which describe the steady state and the dynamical model. The second one follows from the size of the mathematical model, which can only be set up and solved by computer aided methods. In order to get an efficient solution to this type of systems of nonlinear equations a reliable algorithm based on a modified Newton method is introduced. The unambiguity of the solution and the convergence of the procedure is proved by the help of Liapunov functions. Since realistic models of long-distance energy systems consist of many equations, the computer is not only used to solve these nonlinear equations, but also to derive them automatically. The system describing a network can also change during simulation according to switching terminals. Therefore a simple computer aided algorithm based on graph theory is used to set up the network equations.
* Corresponding
author.
037%4754/95/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDI 0378-4754(94)00119-7
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
566
The partial derivatives which are needed for solving the equations are calculated in a symbolic way. This method is based on finite taylor series and special features of the programming language C + + .
2. Introduction Long-distance energy systems are growing fast because of ecological improvements. The size, the complexity and the nonlinearity of the hydraulic equations command the modeling and simulation of these systems. Since a computer must be used to integrate the differential equations, it is natural to use it for setting up the equations, too. The propagation of temperature and pressure are the most important phenomena in these systems. One builds the network for the first effect but a big part of the control equipment is needed to handle the second one. The wave speed of a water hammer is about 1000 [ms-l]. The delay due to the finite propagation velocity is not negligible any more in energy systems for whole cities. Simulation based on mathematical modeling and reliable algorithms offer a good way in understanding the phenomena of those systems.
3. Hydraulic networks A hydraulic network is formed by connecting together terminals like pumps, valves, pipes, etc.. The connection points are called nodes and the terminals connecting two nodes are called branches. This type of network obeys the two laws of Kirchhoff. The current law says that the total mass flow to a node is zero. The density of water may be assumed to be constant in the interesting range of temperature and pressure. Therefore one can replace the mass flow by the volume flow. The voltage law says that the difference of the pressure of two nodes connected by a branch is equal to the difference pressure of this branch. Now the most important branches are presented.
3.1. Pipes The one-dimensional transient flow in pipes can be described by a pair of quasi linear hyperbolic partial differential equations of first order 0
1 0
f
-~xh(X, t) +-~-;q(x, t) + 2gDA---------~lq(x,t)lq(x, t)=O 0
a2
and
0
-~h(x, t) + -~ O---xO(X,t)=0. h(x, t)[m] denotes the pressure head, q(x, t)[mas -1] the volume flow, A the area of the pipe, D the pipe diameter, g the gravitational acceleration, f the friction factor and a the equivalent wave speed [7]. The boundary conditions of a pipe of length l at a time t are of the form a
h(l,t)=--~q(l,t)+R1,
a
h(O,t)=--~q'(O,t)+R 2
and
q'(O,t)=-q(O,t)
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
567
with R 1 and R z depending only on values for "r < t. These conditions describe a 1-port which consists of a linear port and a controlled pressure source in series. The minus sign for q' is used to have the same flow reference direction on both end points. The boundary conditions of a pipe can be embedded in the theory of 1-port networks [4].
3.2. Pumps In long-distance energy systems only rotary pumps are used. As a good approximation the transient behavior of a pump is described by the steady state characteristic curves of the total head increase dh = dh(q, to) and the shaft torque re(q, to) with to as the angular velocity. The torque equation takes the form d
0-~oo = md(to) - m(q, w) with 0 as the moment of inertia and md(to) the torque of the driving motor [7]. If necessary the dynamics of the driving motor is taken into consideration, too.
3.3. Valves The difference pressure of the valve is given by
Iqlq dh =
-r
and
d -~a = f ( q , a, u)
with 0 ~
with r as the resistance factor, a as the grade of opening and u as the control input. The topology of the hydraulic network changes at each time, when a switching component like a valve is opening or closing [7].
3.4. Steady state and transient analysis A large system of nonlinear equations has to be solved in the case of steady state analysis. The transient behavior is described by systems of partial and ordinary differential equations with nonlinear equations as restrictions. The latter ones are a consequence of Kirchhoff's circuit laws. Investigations about the computing time have shown, that the solution of this system of nonlinear equations takes much more time than the integration of the differential equations. These circuit equations are described by the theory of 1-port networks most efficiently. They are of the same type for the steady state and the transient analysis, therefore only the first one is treated. It is very important to use efficient algorithms to solve this problem in order to obtain a short computing time. A fast algorithm will be presented after the uniqueness of the solution of the hydraulic network equations is shown.
4. Facts of networks
Only the steady state equations are considered according to the above-mentioned observations. The circuit is formed by connecting together special terminals. A graph G = {N, B} is
568
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
related to the network. N is the finite set of nodes of cardinality n and B is the finite set of branches. A branch has exactly two end points which must be nodes. T h e n u m b e r of branches is b. T h e flow qi and the difference pressure d h i are assigned to the branch i. T h e pressure hj is related to a n o d e j. Now the laws of Kirchhoff can be expressed as hid i+hkd~=dh with
i{i
dt -
,
and
~-'~d~qJ=O J
if n o d e i and branch 1 are connected, the direction of the flow to i is positive, if n o d e i and branch l are connected, the direction of the flow off i is positive, otherwise.
T h e flow q of a network is a point q = ( q l , . . . , qb) ~ Rb. q is said to be admissible if q satisfies Kirchhoff's current law. A difference pressure dh is a point dh = ( d h l , . . . , dh b) ~ (Rb) * . dh is said to be admissible if dh satisfies Kirchhoff's voltage law. Let p be the natural bilinear m a p p . ( ~ b ) * X ~b ~
p ( d h , q ) = ~_,dhjq j J then Tellegen's t h e o r e m can be expressed in the following form [5]. Theorem 1. For all admissible q and dh is p(dh, q)=0. T h e linear m a p D" Rb .__,Rn X i=
Ed~q j J
is well defined. For a basis
{a i} of
Ker D there exists a linear m a p A : R c ~ Ker D
qi= Ea x j J which is bijective for some c. Tellegen's t h e o r e m leads now to the following results [2]. Theorem 2. A flow q is admissible iff Dq=O
or
q=Ax.
(1)
A difference pressure dh is admissible iff dh ~ Im D*
or A *dh = O.
(2)
It is a well-known fact that one can always find a basis of Ker D that x j is the flow qi of a branch i.
K. Schlacher,A. Kugi/ Mathematicsand Computersin Simulation39 "(1995)565-572
569
5. A potential function Only 1-ports are considered as branches for this type of network. Then there is a branch equation of the form: (a) dh i =fi(qi), (b) dhi = const., or (c) qZ = const. for each branch i. For the sake of simplicity only equations of type a) are considered, the extensive case is presented in [4]. The short form of the branch equations is
dh = f ( q ) .
(3) A flow q is said to be a solution of the network equations if q and dh(q) are admissible (equations (1), (2)) and satisfy the equation above or q=Ax
A*f(q)=O
and
(4)
hold. Based on Tellegen's theorem we define the function q
qi
.
For an admissible flow the derivative of v with respect to x ~ is given by
ai OX
v(Ax)= ~O ax ) ( x ) qU(AX)~x/(A
~--- --
~ fj ( ~ , a j x l ) a [ = - p ( f ( A x ) , •
ai).
l
Each solution of the network satisfies the equation a
7x(Ax)=O.
Now it is possible to give the main result of this paper. Theorem 3. If the branch equations fi (3) of the network satisfy the conditions f i e c l ( - ~ , oo)
and d then the solution q of the network equations (4) is unique. To prove this theorem first the functions qi
gi(qi) = -- fo f i ( q ) d x ' are considered. It can be easily shown that gi is strictly convex [3] and that the condition limlqil _,®gi(qi) = oo is satisfied, v(q) is strictly convex, since v(q) is the sum of strictly convex functions and the condition limllql I ~ v ( q ) =
570
K. Schlacher, A. Kugi / Mathematics and Computers in Simulation 39 (1995) 565-572
is satisfied, too. Because of Ker A = {0} the restriction of v(q) to v(Ax) is strictly convex and lim4rxjI __,oov(Ax)= oo
holds. It is a well-known fact of the theory of convex functions that v(Ax) has a global minimum and x is the unique solution of ~ v ( A x ) = 0. This completes the proof. It must be emphasized that the conditions of Theorem 3 are no restriction for hydraulic networks.
6. A modified Newton method
Due to the theory of gradient systems and the present considerations a simple algorithm to solve the network equations is given by d - ~ x i = o ( f ( A x ) , ai). In general a faster algorithm can be constructed using Newton's method [3]. This algorithm is applicable since v(Ax) ~ C 2 and v has a unique minimum. Let Jv be the Jacobian and Hv the Hessian of v. Then the sequence
x(i + 1 ) = x ( i ) - a ( x ( i ) ) H ~ l ( x ( i ) ) J f ( x ( i ) )
(5)
with the tuning function a(x(i)) is well defined. The Hessian of v(Ax) is positive definite since v(Ax) is strictly convex. In this special case the Hessian H~ is given by
H~(Ax) = A * ( - J f ) A with JI as Jacobian of f of equation (3). Expanding v into a Taylor series the equation
v(x(i + 1 ) ) = v ( x ( i ) ) + Jv(x(i))(x(i + 1 ) - x ( i ) ) + O( [[x(i + 1 ) - x ( i ) [[ z) follows for sufficiently small [[ x(i + 1) -x(i)[[. With the help of
v(x(i + 1)) - v ( x ( i ) ) = - a ( x ( i ) ) J . ( x ( i ) ) H ~ ' ( x ( i ) ) J f (x(i)) + O(aZ(x(i))) one can ensure that the inequality
v(x(i+ 1))-v(x(i))
K. Schlacher,A. Kugi/ Mathematics and Computersin Simulation 39 (1995)565-572
571
This theorem ensures not only the convergence it offers a way to choose a by checking the difference v(x(i + 1)) - v(x(i)), too.
7. An algorithm for calculating derivatives The set of the network equations is too large to be set up by hand. It changes also during a simulation according to the action of switching terminals. It can be shown that the graph searching algorithm "depth first" offers an efficient way setting up the equations in an automatic way [4]. Although this algorithm can be used to calculate all required derivatives for the Newton procedure, too, a more general method is presented. This approach needs a modern programming language like C + + which allows overloading of functions of binary operations and dynamic memory management [1]. First a new type of variable x = (x v, (xi)) is defined as a pair of a float variable x~ and an array of float variables (x i) of dimension one. Next the binary operations x+y=(xv+Yv,(xi+Yi)),
x'y=(xvY,,(YvXi+XvYi))
and the rule
f(x)= (f(xv), ( d for the value of the real function f are introduced. Let x b e a finite Taylor series of the form 0
x
x(z),
The value of the variable x at z is defined as
x= x(z),--x(z)l. az~ ] It is easy to verify that these operations implement the rules to handle this special kind of series. If these operations are used in a program all required derivatives are calculated in an automatic way. The advantage of this approach is that only a new type of variable and overloading of some operations is required. It is not necessary to write a new program.
8. Results The above introduced methods are the basis of the programming system F L U I D I X for the simulation of long-distance energy systems. This package consists of a graphic editor, a data base of mathematical models, an interactive simulator and a graphic post processor. The user only draws the plan of the circuit assisted by mouse and menus. Parameters of the models are added by the user or by the data base. The system derives the equations for the steady state or
572
K. Schlacher, A. Kugi /Mathematics and Computers in Simulation 39 (1995) 565-572
transient analysis automatically. During the simulation the user can trigger events with the mouse or keyboard. The results of a simulation are d o c u m e n t e d by the graphic post processor. The programming system F L U I D I X has been used for the analysis of long-distance energy systems of cities with more than 250.000 inhabitants. Therefore one can say that the methods presented in this p a p e r are applicable for real world problems.
References [1] K.E. Gorlen, S.M. Orlow and P.S. Plexico, Data Abstraction and Object-Oriented Programming in C+ + (Wiley, 1991). [2] M. Hirsch and S. Smale, Differential Equations, Dynamical Systems and Linear Algebra (Academic Press, 1974). [3] D.G. Luenberger, Linear and Nonlinear Programming (Addison Wesley, 1989). [4] K. Schlacher, Systemanalytische Methoden f'tir hydraulische Netze, Habilitationschrift an der TU-Graz, 1990. [5] P. Penfield, R. Spence Jr and S. Diunker, Tellegen's Theorem and Electrical Networks (MIT Press, 1970). [6] M. Vidyasagar, Nonlinear Systems Analysis (Prentice Hall, 1993). [7] E.B. Wylie, V.L. Streeter and L. Suo, Fluid Transients in Systems (Prentice Hall, 1993).