A NEW MATHEMATICAL MODEL OF DYNAMIC PROCESSES IN DIRECT CURRENT TRACTION POWER SUPPLY SYSTEM

A new autonomous 4D nonlinear model with two nonlinearities describing the dynamics of change of voltage and current in the contact railway electric network is offered. This model is a connection of two 2D oscillatory circuits for current and voltage in the contact electric network. In the found system for the defined values of parameters an existence of limit cycles is proved. By introduction of new variables this system can be reduced to 5D system only with one quadratic nonlinearity. The constructed model may be used for the control by voltage stability in a direct current power supply system.


Introduction
The modern stage of functioning of railways is conditioned by the necessity of providing competitiveness with other types of transport. The decision of this problem supposes introduction of high-speed passenger transport as well as heavy freight trains. For this purpose on the railways measures on the increase of speed of movement are carried out, new electric locomotives of large power are created, different ways of strengthening traction power supply are applied (see, for example, [1]). The traction power system of the electrified section of the railway (TPS) is a set of territorially dispersed and operating electric power stations. This system can include traction substations, sectionalizing stations, parallel connection points, contact network devices and power transmission lines between them, united by a common purpose and intended for processing and transmission of qualitative electric energy to electric rolling stock (ERS). The peculiarities of electric power transmission through the traction network is the change in the position of ERS and the change in their operating modes, the restrictions imposed by trains on each other depending on their relative location, as well as the restrictions associated with the technology of the transportation process as a whole. One of the main indicators of the quality of transmitted electrical energy to ERS is the level of voltage on the current collectors of electric locomotives, and the nature of the factors affecting on this level, which are nonlinear and non-stationary.
Ensuring the stable and reliable operation of any technical system is an important task that requires its solution. Voltage resilience is the ability of the power system to maintain stable and acceptable voltage levels on all bus systems (BS) in both normal and post-emergency and repair modes. The criterion for the stability of the power system in terms of voltage is that, in the current mode, the value of reactive power on the same BS should increase at each BS with increasing voltage. Dynamic voltage stability is associated with the evaluation and support of the voltage within 1 -2 seconds immediately after a large disturbance. Static voltage stability belongs to the form of stability, determined mainly by the static characteristics of the load and network parameters.
The existent system of traction electric supply of direct current not always is able to provide the transmission electric power of necessary capacity for speed trains. In this connection there can be the following limitations: lowering of voltage in the contact network is below than normal level 2.7 kV for ordinary motion (below than level 2.9 kV for high speed motion) and heatings of wires of contact network, that will promote in the loss of their mechanical durability. Lowering of voltage in the receiver of current diminishes the speed of movement of trains. Thus, the saving of level of consumable power results in the increase of current of electric locomotive and loss of electric power in the contact network.
Since the dynamical state of any technical system is described by a system of differential equations, the study of the problem of the stability of its motion reduces to the study of the stability of solutions of differential equations. To calculate the static stability of the power system, it is necessary to compile a system of differential equations of transient processes, linearize these equations, and obtain the characteristic equation of a system of linearized equations. Together, these equations constitute a mathematical model of the energy system [2], as a result of which solutions it is possible to obtain algebraic and other stability criteria for the system under consideration.
In this connection, there is an urgent need to consider the problem of determining the stability of TPS as an initially nonlinear problem. At present, an approach based on the analysis of signals produced by the system is widely used to study the properties of complex systems, including experimental studies. This is very important in cases when it is practically impossible to describe the pro-cess under study mathematically, but we have at our disposal certain observable quantities, which allow to build a wanted model [3]. Now there is a great number of methods of diminishment of losses of electric power in the contact network [2], [3]. In the present paper the diminishment of losses of electric power and increase of power efficiency of the systems of electric supply will be attained due to introduction of new dynamic models describing the behavior of current and voltage in these networks. Due to these models new methods of calculation of parameters of contact electric networks can be offered. Such methods allow to apply new organizational measures resulting in diminishment of losses of electric energy in contact networks. Let be a finite sequence of numerical values of some scalar dynamical variable x(t) measured with the constant time step ∆t in the moments t i = t 0 +i∆t; x i = x(t i ); i = 0, 1, ..., n. Sequence (1.1) is called a time series [4] - [10]. In future the time series (1.1) will characterize a voltage (or current) in the contact network measured through the equal intervals of time.
A common practice in chaotic time series analysis has been to reconstruct the phase space by utilizing the delay-coordinate embedding technique, and then to compute the dynamical invariant magnitudes such as unstable periodic orbits, a fractal dimension of the underlying chaotic set, and its Lyapunov spectrum. As a large body of literature exists on applying of the technique of the time series to study chaotic attractors [10] - [15], a relatively unexplored issue is its applicability to dynamical systems of differential equations depending on parameters. Our focus will be concentrated on the analysis of influence of parameters of found dynamic system on the behavior of its solutions. These parameters are determined by the structure of the time series (1.1) and choice of approximating functions in right sides of the got system of differential equations. To create a model by measuring the variables characterizing any dynamic process, it is necessary to solve the following four main problems.
It is known that any dynamic process depends on many variables. Most of these variables are functions of some small number of independent variables. Identifying these independent variables leads to the first problem.
Problem 1. Determine the dimension of phase space in which the explored process takes place.
Usually, a continuous dynamic process is described using a system of differential equations. It is important to specify the class of those functions that form the right-hand sides of the differential equations. In applied modeling problems, the right-hand sides of differential equations are given by polynomials. In this connection, the second problem arises.
Problem 2. Establish degrees of monomials and their composition that form the polynomial right-hand sides of the differential equations. Problem 3. After the structure of the differential equations describing the dynamic process is established, it is necessary to determine the numerical value of coefficients of these equations. Problem 4. Using the specifics of the problem, to establish analytical formulas of the coefficients of the found system of differential equations as functions of characteristics of individual elements of the direct current traction power supply system (it can be capacitors, resistances, inductances, and so on).
We must say that of all these problems, the fourth is the most difficult. Indeed, for such a complex system as the direct current traction power supply system we can not specify all electrical elements included in its composition. Therefore, in this paper we confine ourselves to solving only the first three problems.
Note that in great numbers papers on the problem of minimization of losses of electric power in contact networks, analytical simulation techniques are used. In these papers mathematical models are linear or, at the best, linearized. Their use gives positive results. However, losses of electric power in the contact network may be diminished only on 20 -30 percents [1].
Therefore, construction of mathematical models allowing suggest practical methods to decrease the losses of electric power in the contact network is the primary purpose of the present paper. These models are constructed on the basis of the real information taken from the electric system of locomotive. As experimental information the results of measuring of currents and voltages taken on the segment of Nyzhne-Dneprovsk Knot -Pyatykhatky of Prydneprovskaya Railway (Ukraine) are used.
It should be said that we do not know publications in which the direct current traction power supply system was modeled only by measuring currents and voltages (see [16,17]). The fact is that a large number of electrical elements forming such system practically exclude the possibility of constructing a detailed model through the combination of the equations of its individual elements. (Such model would represent a set of several hundred different equations.) Therefore, the main achievements of this work are: -creation of a new model of direct current traction power supply system, which describes this system not as a set of interconnected elements, but as a single dynamic element whose behavior is determined only by changes in voltages and currents flowing through it; -development of a universal methodology for modeling of direct current traction power supply systems suitable for any segments of contact networks located anywhere in the world.

Embedding Method for Chaotic Time Series Analysis
The material of this section is well known (see [4]). The results of section are placed in the present article only for the convenience of readers.
Let sequence (1.1) be the time series. In principle, the measured time series comes from an underlying dynamical system that evolves the state variable in time according to a set of deterministic rules, which are generally represented by a set of differential equations, with or without the influence of noise. Mathematically, any such set of differential equations can be easily converted to a set of first-order, autonomous equations. The dynamical variables from all the first-order equations constitute the phase space, and the number of such variables is the dimension of the phase space, which we denote by M. The phase space dimension can in general be quite large (in some cases it may be infinite) [10,11,16,18].
However, it often occurs that the asymptotic evolution of the system lives on a dynamical invariant set of a finite dimension. The assumption here is that the details of the system equations in the phase space and of the asymptotic invariant set that determines what can be observed through experimental probes, are unknown. The task is to estimate, based solely on one or few time series, practically useful statistical quantities characterizing the invariant set, such as its dimension, its dynamical skeleton, and its degree of sensitivity on initial conditions. The delay-coordinate embedding technique established by Takens [9], in particular, his famous embedding theorem guarantees that a topological equivalence of the phase space of the intrinsic unknown dynamical system can be reconstructed from the time series, based on which characteristics of the dynamical invariant set can be estimated.

Letẏ
(t) = F(y(t)), y ∈ M ⊂ R p (2.1) be the autonomous p-dimensional system of ordinary differential equations in the phase space M.
We will consider that system (2.1) satisfies in the phase space M (an open region in R p ) to the conditions of the known Cauchy Theorem about existence and uniqueness of solutions. Then for any initial condition y(0) = y 0 ∈ M it is possible uniquely to define the solution y(t) systems (2.1) on the formula y(t) = W t (y 0 ), where W t is an evolution operator. (A domain G ⊂ M of the phase space M under action of the evolution operator passes, generally speaking, in another domain It is known [12] that it is possible to get the attractor satisfactory image of a small dimension, if instead of the phase vector y(t) to use m-dimensional vectors derived from the time series (1.1) on the following formula: . . .
where l is a positive integer. Consider the m-dimensional autonomous dynamical systeṁ for which the following conditions The magnitude x i depends on x 0 and τ , but it does not depend on t 0 . We will especially emphasize that the map Q : R m → R m , determining the right side of systems (2.3), it is not known. In addition, it is clear that the role of the number l in (2.2) plays the number τ . The magnitude τ is called a delay parameter of the time series (1.1). Introduce the evolution operator P t : R m → R m of system (2.3). For any vector x ∈ R m the action of this operator in a coordinate form looks like: Let t = τ . Consider the sequence of real numbers Introduce the new vector z k under the formula: z k = (h k , h k+1 , ..., h k+m−1 ) T . Then there must be an operator ∆ : R m → R m depending only on Q and τ such that z = ∆(x), where x = (x i , x i+1 , ..., x i+m−1 ) T is one of vectors (2.2).
Theorem 2.1. [9] Let d be a dimension of the attractor Σ generated by system (2.3). Then for almost all τ > 0 and m ≥ 2d+1 the mapping ∆ will be continuous and one-to-one. Theorem 2.1 means that if in the space R m to select the set H k such that ∀x k ∈ H k we have ∆(x k ) ∈ H k , then on this set the map ∆ is invertible and ∀k By N denote the set of natural numbers. Thus, Theorem 2.2 (which is sometimes called the Poincare recurrence theorem) asserts that in the phase space of the dissipative system any trajectory beginning from the almost liked point A of this space in some finite time (even very large) will pass as much as close to A.
Theorems 2.1 and 2.2 allowed to create the necessary research instrument which is used presently in the theory of the dynamic systems. Indeed, as the time series (1.1) has only a finite number of terms and, consequently, it is bounded, there are no justified arguments in order to assert that at the further measurements we will derive very large values of terms of this series. Further, the time series (1.1) describes the behavior of some phase variable of the explored dynamic process. If we assume that the number of such phase variables is finite, then it is possible to consider that there exists the evolution operator, which controls by the behavior of this dynamic system in some finite-dimensional space. In addition, most systems describing the dynamics of one or another processes in our world are dissipative. Thus, the use of Theorems 2.1 and 2.2 for description of dynamics of the dissipative finite-dimensional systems becomes more than justified.
Eckmann [6] have introduced tools which visualize the recurrence of states x i in the phase space. Usually, the phase state does not have a dimension (it is more than two or three) which allows it to be pictured. Higher dimensional phase spaces can only be visualized by projection into the two or three dimensional subspaces [10], [12].
., x i+m−1 ) T , which is built from the elements of the time series (1.1) describing the change of some scalar variable (or some coordinate of the vector variable, if a phase trajectory in m-dimensional space is considered); i → 0, 1l, 2l, ..., il. If il + m − 1 > n, then number i must be replaced by the number k → kl = il − [n/l]l, where [n/l] is an integer part of the number n/l. We will consider that l = τ .
Introduce in the first quadrant of the cartesian system of coordinates the graphic square matrix T ∈ R (n+1)×(n+1) , which is built on the following algorithm: if point x(i) is close enough to the point x(j) (the concept of "closeness" will be defined below) then such points are called recurrence, and in the matrix T a black point with coordinates (i, j) are put. If point x(i) is not near to the point x(j), then in the matrix T no marks is done. The matrix T is called a recurrence plot of time series (1.1) [6], [10]. Let be a real function accepting only two values: 0 and 1.
is the Euclidian norm of the vector v ∈ R m ; i is a radius of ball with a center in the point x i .) In the future, it is possible to be restricted to the situation, when ∀i, j i = j = . In this case positive number is called a recurrence threshold and we have symmetry of the recurrence plot with respect to the diagonal of the first quadrant. Indeed, if point x i is near to the point x j , then the reverse statement must be right: the point x j is near to the point x i . In the present paper we want to apply the instruments of the recurrence analysis for research of periodic trajectories in the dynamic systems described by nD autonomous systems of differential equations. In order that such research was correct it is necessary to provide the boundedness of solutions of the explored systems.

Mathematical Statement of Problem and Its Discussion
We will assume that we can measure the voltage and current, and also if it is possible other dynamic characteristics of contact electric network. We also suppose that among these characteristics can be derivatives with respect to t from the voltage and current. (If the derivatives can not be measured, it is assumed that there exist smooth enough approximations of these derivatives.) A choice of equations of model of describing the dynamics of one or another processes is a difficult task. The experiments show that the most logical approach describing dynamics of electrical engineering models is based on the use of the known physical laws. In particular there can be the conservation of energy laws.
By U (t) and I(t) denote respectively the voltage and current in contact electric network. The following laws are most known: the electric energy is accumulated in a capacitor according to the law E C = k C U 2 ; the electric energy is accumulated in an inductor according to the law E L = k L I 2 ; the electric energy is transformed into heat energy on a resistor according to the law E R = k R U I. (Here k C , k L , and k R are constants.) In addition, a speed of change of energyĖ C = 2k C UU , E L = 2k L Iİ orĖ R = k R (U I + Uİ), and magnitudes U,U , I,İ also influences on the dynamics of electric network.
Thus, the vast class of electric networks can be described by quadratic differential equations depending on the linear U,U , I,İ, and quadratic U 2 , UU , I 2 , Iİ, U I,U I, Uİ terms.
Let (c 1 , ..., c n ) T and A = (a ij ), B 1 , ..., B n ∈ R n×n be a real vector and real matrices, and let the matrices B 1 , ..., B n be symmetrical.
Principal problem. Construct the quadratic system of differential equations Further, we use the procedure for determining unknown quadratic right sides of the system of differential equations (3.2), which was suggested in [12] - [15]. This procedure is based on the least squares method and the fact that we know sufficient precision the components of x(t) and its derivativeẋ(t).
elements of which are known. Then by the least square method [12,20,21], we have Further, the following is said in work [13]. In view of the fact that number N may be chosen arbitrary large, a high precision reconstruction may be achieved. Thus, we can expect that the solution of reconstructed system will be near the purified solution x(t).
However, it should be said that one important circumstance, which can arise up at a reconstruction, remained outside the attention of the authors of article [14]. The point is that in [14] it is assumed that this interval (t 1 , t N ) is finite. If the problem of long-term prediction is considered, it is necessary to assume that t N → ∞. In this case a reconstruction must be fulfilled so that system (3.2) had the bounded solutions [22]- [24].
Finally, we mark that the presence of all quadratic elements in the right side of system (3.2) do not always result in that the built model will adequately describe the explored process. Therefore, the choice of base of quadratic part of system (3.2) must be argued by the real information about the studied process, by which it is impossible to neglect.

Simulation of Dynamics of Current and Voltage in Contact Electric Network
We consider the following scheme of placing of traction substations on the segment Pyatykhatky -Nyzhne-Dneprovsk Knot [1]: Measurings of current and voltage will be realized by a measuring laboratory which moves on the segment together with a locomotive with a constant speed. Time of passing of all segment Pyatykhatky -Nyzhne-Dneprovsk Knot is 12470 sec.
On the segment Pyatykhatky -Nyzhne-Dneprovsk Knot the temporal dependences for current and voltage in the contact network were built. In all 12470 measurings with an interval of 1 second were done. The graphs of these dependences are represented on A standard procedure for modeling of the direct current traction power supply system consists of the following steps: 1. Construct a scheme replacing the direct current traction power supply system of the considered segment (no load; see Fig. 4.3).
2. Introduce in the model obtained at the first step, an element describing a motion of locomotive (it is a load; see Fig. 4.4).
3. Calculate using the Ohm and Kirchhoff laws, the voltage and current variations along in all the segment Pyatykhatky -Nyzhne-Dneprovsk Knot.
4. If the measured values of current and voltage significantly differ from calculated in the step 3, then the topology of substitution scheme (see Fig. 4.4) and the characteristics of elements forming this scheme must be changed.  It is clear that such procedure leads to a very approximate model of the direct current traction power supply system (see, for example, [17]). Therefore, in this paper we propose another modeling method based on recurrence analysis.
A preliminary analysis of the obtained data shows that they are nonstationary. In this connection, the solution of problems of modeling and forecasting of nonstationary processes is of particular relevance. We point out that nonstationarity can manifest itself in the appearance of a deterministic or stochastic trend that varies in time with variance and covariance. There are two main purposes of analyzing time series: the determination of the nature time series and prediction (the prediction of future values of time series by present and past values). Both these goals require that the series model be identified and, more or less, formally described [25] - [28].
In order that to successfully fulfill the modeling of processes of represented on Fig. 4 We take advantage of the least squares method. In accord to the done calculations an embedding dimension space n must be not less than 4. In future we will consider that n = 4. It should be noted here that dimension 5 can also be considered.
Thus, the dynamic system assumes existence of limit cycle of dimension N = 1 ( in this case, we have n = 4 > 2 · N + 1 = 3). Exactly bifurcations of limit cycles result in an appearance of chaotic dynamics.
On all following Fig. 4.6-4.12 the voltage U (t) = x(t) and current I(t) = z(t) is measured in kilovolts (kV) and kiloamperes (kA).
(4.2) 3. The base of quadratic part of system (3.2) consists of two elements {xz, xy +zu}: . 4. The base of quadratic part of system (3.2) consists of two elements {xz, xy +zu}: 1269(x(t)u(t) + z(t)y(t)). (4.5) 6. The base of quadratic part of system (3.2) consists of four elements {xy, xu, zy, zu}: 7. The base of quadratic part of system (3.2) consists of six elements {xy, xz, xu, x 2 , y 2 , z 2 }: It should be noted that all the considered models are chaotic: an arbitrarily small change in the parameters of the model leads to a radically different behavior of this model. (Nevertheless, one can always find such parameter values for which a limit cycle appears in the system. his fact can be used to construct stabilizing control laws. Another application of the obtained limit cycle can be the search for limit tori.) Now there is a question: what from models (4.1)-(4.7) most adequately describes the behavior of process of represented on Fig. 4.2? Comparison of trajectories of the real chaotic system and its model does not enable to speak about their adequacy. Therefore, we decided to define the adequacy beginning from comparison of U − I characteristics of model and real process. In this connection we build a parallelogram ABCD in the coordinate system U − I (see Fig. 4.13): This parallelogram must possess the following features: a) taking into account that we deal with the direct current power supply system bases AB and CD of parallelogram ABCD must be parallels to axis U ; b) the vertices of parallelogram must be disposed in points: A(U min , I max ), B(U min +δ , I max ), C(U min , I min ), D(U max −δ, I min ), where I min (I max ) is a minimum (maximum) current with the exception of some random fluctuations; U min (U max ) is a minimum (maximum) voltage with the exception of some random fluctuations; δ is a magnitude of change of voltage at the fixed current.
Definition 4.1. The parallelogram ABCD is called a framework of U −I characteristic for the direct current power supply system. Let S ∈ R 2 be the set of all internal points of the parallelogram ABCD. Now we are ready to answer on the question about adequacy of model and real process.
By M 1 , ..., M k ∈ R 2 denote U − I characteristics of models, which it is possible to build beginning from the real process of represented on Fig. 4.2).
Introduce the Hausdorff distance d H (S, M i ) between the sets S and M i , i ∈ {1, ..., k} [29]. We notice that exact calculation of the Hausdorff distance d H (S, M) for the complex sets S and M is practically impossible. Therefore in the present work we will be limited to rough enough estimations of the magnitude d H (S, M).
Thus, the analysis of models (4.1)-(4.7) shows that only the behavior of model (4.1) is adequate to the real process of represented on Fig. 4

Modeling of Other Contact Networks
All the above equations modeled the experimental data of voltage and current changes (they are shown in Fig. 4.2) in the contact network. In order to dwell on specific equations modeling the dynamics, it is necessary to use other data representing the dynamics of processes in other contact networks. These data, describing changes in voltage and current on other railway lines, are presented in Fig. 5.17. (Note that the above data, generally speaking, is non-stationary. Therefore, for good modeling it is necessary to have such data as much as possible.) Now we use the methodology of Sections 3 and 4. Then, we get the following equations and the corresponding behavior of their chaotic solutions (see Fig. 5.14- Fig. 5.16): 1. The base of quadratic part of system (3.2) consists of two elements {xz, xu+ zy}: +0.0005x(t)z(t) + 1 (x(t)u(t) + z(t)y(t)), z(t) = u(t), u(t) = 0.01016 − 0.0029x(t) + 0.0006y(t) + 0.00228z(t) + 1.14108u(t) −0.001481x(t)z(t) − 2 (x(t)u(t) + z(t)y(t)). 2. The base of quadratic part of system (3.2) consists of two elements {z 2 , zu}: 3. The base of quadratic part of system (3.2) consists of two elements {x 2 , xy}: Using the methodology of Sections 3 and 4, as well as Fig. 4.13, we arrive at the following conclusion: the most adequate description of processes in the contact network is achieved when pairs of quadratic nonlinearities {x 2 , xy} or {z 2 , zu} or {xz, zy + ux} are used as nonlinear terms in system (3.2).

Research of System (4.1)
For verification of conditions Theorem 2.2 we will carry beginning of coordinates of system (4.1) in the point (3.3657, 0, −0.0388, 0). In the total we get such system The condition of dissipativity for system (6.1) ∂ẋ ∂x + ∂ẏ ∂y + ∂ż ∂z + ∂u ∂u = 0 + 0.0086 + 0 − 0.0095 = −0.0009 < 0 (and for system (4.1)) is fulfilled. On Fig. 6.18 possible solutions of system (4.1) are shown. Fig. 4.2 and 5.17 show that processes in contact networks are chaotic. It is known that any chaotic processes begin from bifurcations of limit cycles. Therefore, we show the existence of such cycles using the example of system (4.1).

Consider the trigonometric polynomial
where a number V > 0 must be chosen so that the real trigonometric periodic function H 1 (cos φ, sin φ, V ) has three real distinct roots on a period. (The number V must be large enough.) By V * denote the minimal value V > 0 at which the periodic function H 1 (cos φ, sin φ, V ) has three real roots on the period.
Further, we have that if λ * 1 < tan φ < λ * 2 , then H 1 (cos φ, sin φ, V * ) > 0, and if Since (0, 0) T there is a repellent point, then from here it follows that the maximal value of function V must be bounded.
Let S be a bounded set containing the point (0, 0). (It can be a circle of small enough radius.) Now we will organize the iterative process (6.5) from any point of the set S. In this case, in acording to LaSall's Theorem [30] all limit points of the process (6.5) form a positively invariant set C.
Since system (6.3) has also the saddle equilibrium, then her separatrix is a restriction for any trajectory of the system of beginning in the small neighbouring of the point (0, 0) T . It means that sequence (6.4) (it is (x 0 , y 0 ) T , (x 1 , y 1 ) T , ...) must converge to the set C. Thus, this is a stable limit cycle. Consider the following system: In system (6.6) it is possible to select two subsystems:   ẋ (t) = y(t), y(t) = a 10 + a 11 x(t) + a 12 y(t) + a 13 z(t) + a 14 u(t) +b 12 x(t)y(t) + b 11 x 2 (t).
Now we can apply Theorem 6.1 to system (6.12) (and (6.9)). To do this, we show the sequence of bifurcations leading to the appearance of chaotic dynamics from the limit cycle (see Fig. 6.19).
Let A, ω be an amplitude and frequency of external perturbation A sin(ωt). (If A = 0, then the perturbation there doesn't exist.) For system (6.7) the perturbed system is   ẋ (t) = y(t), y(t) = 0.0193 − 0.0072x(t) + 0.0218y(t) − 0.0039x(t)y(t) +0.000422x 2 (t) + A sin(ωt). (6.13) All conditions of Theorem 6.1 for system (6.12) are satisfied. (It can be confirmed directly.) Appearance in system (6.13) of a stable limit cycle at the corresponding external perturbations is shown on Fig. 6.19. Further, on system (6.8) it is also possible to look as on a linear system with external perturbation G(t) ≡ a 21 x(t) + a 22 y(t) + b 21 x(t)y(t) + b 22 x 2 (t). In order that the solution of this linear system was bounded it is necessary that the equilibrium (−a 20 /a 23 , 0) unperturbed system was stable and the function G(t) was bounded [30]. For system ż(t) = u(t), u(t) = 0.0294 − 0.0019z(t) − 0.0095u(t) + G(t) (6.14) both these conditions are fulfilled. The equilibrium (z * = −15.4737, u * = 0) (at G(t) ≡ 0) is a stable focus. In system (6.14) the boundedness of function G(t) is guaranteed by the boundedness of solutions of system (6.13) [23,24]. Thus, on system (4.1) it is also possible to look as on a system of two connected 2D circuits of describing oscillations of current and voltage in the contact electric network.

Remarks on Design of Voltage Regulator
At certain values of parameters system (6.6) describes the dynamics of changes in voltage and current in a contact network. Note that voltage U (t) = x(t) and current I(t) = z(t) are measured using a mobile laboratory that moves at constant velocity v along a contact network [31]. Note that in order to study processes in the contact network, it is more convenient to employ a dynamic model in which voltage U (t) = x(t) and current I(t) = z(t) are represented as functions U (s) = x(s), I(s) = z(s) of distance s from some starting point. Such representation is very convenient when the stabilization of voltage in contact network is implemented from some fixed points along the route of the train, for example, at traction substations or gain points. Furthermore, we shall assume that voltage control is executed using the regulator U input(s) = f (U (s),U (s), I(s),İ(s)), where f (...) is the real function of its arguments. Thus, the transition from model (6.6), where x(t) and z(t) are represented as a function of time t, has been achieved through the replacement of independent variable t with independent variable s, according to formula s = v s t. In this case, y(t) → v s y(s), u(t) → v s u(s), and system (6.6) transforms to the following system: (For simplicity, we kept the former designations for dependent variables x and z in the newly derived system. The variables x(t) and z(t) are replaced with variables x(s) and z(s). Now we will introduce a control law. Model (7.1) is designed to study the stability of the voltage in the contact network. This model was obtained by observing the chaotic behavior of voltage and current. In future, for simplicity we will consider that v s = 1.

(7.2)
Hence it is already possible to establish parameters under which the voltage in system (7.2) (or (6.6)) can be stabilized [30,32]. (The corresponding characteristics of real behavior of the voltage and current are given in Fig. 5.17.) Now we will do an attempt yet to simplify system (7.2) and to do this system of more suitable for further researches.
We will consider that system (6.6) has an equilibrium P = (x * , y * , z * , u * ). Then we can transfer the origin of coordinates in the point P . In this case in system (7.2) we have a 10 = a 20 = 0. Therefore, we can consider that the condition a 10 = a 20 = 0 is satisfied.
Introduce the matrices Then by the algorithm of indicated in [33], the following result can be got: 2) can be reduced to the following canonical form: The canonical form (7.3) is generalization of the known Bezout's states column model [33]. (The system (7.3) is interesting to those that unlike the system (7.2), this system has only one nonlinearity.) The last equation can be generalized in the following way. We will assume that the model of direct current power supply system is association of several oscillatory circuits for description of the voltage, current, electromagnetic induction, and so on. Then model (7.3) may be generalized in the following form: (7.4) Here, .., d n , p i , q i ∈ R; i = 1, ..., n + 1. This equation can be useful at the further study of the model of direct current power supply system. Now we show a simple method of construction of stabilizing linear feedback for system (6.6).
Introduce in system (7.6) a linear feedback by the formula where k 1 , ..., k 4 are indeterminate coefficients. Then a linear part of the closed by the feedback system has a characteristic polynomial h(λ) = λ 4 + ( Choose the coefficients (k i − d i ) of polynomial h(λ); i = 1, ..., 4, so that this polynomial became the Hurwitz polynomial [34]. In this case the origin of the closed by feedback system will be stable (see Fig. 7 Let S be a linear transformation reducing system (7.5) to form (7.4): Then in order that the origin of system (7.5) (a 10 = a 20 = 0) would be stable it is enough to give the control law by the formula F (t) = (k 1 , k 2 , k 3 , k 4 ) · S −1 (x(t), y(t), z(t), u(t)) T .
An area of stability of power system is the set of its modes, in which static stability is provided for a certain composition of the generators and a fixed circuit of the electric network. A surface bounding a set of stable regimes is called a boundary of region of static stability [32]. The stability regions are constructed in the coordinates of the parameters that affect the stability of the regime. The calculated and experimentally determined areas of stability are used to set the dispatch restrictions on the regime of the power system (in the form of dispatch instructions) and to configure automation facilities to prevent possible violations of static stability. Obviously, reliable and stable operation of the power system in modes directly adjacent to the boundary of the stability region is impossible. In these modes, any, even weak disturbances in the power system or spontaneous minor weighting of the regime will lead to a violation of stability. Changes in the regime of the power system (active and reactive overflows, voltage and frequency) are primarily associated with load fluctuations in load nodes: the inclusion and disconnection of individual electrical installations, start up and shutdown of enterprises, changes in their operating mode according to technology conditions, etc. Part of these changes are regular character, due to daily, weekly, seasonal regimes. Such changes are described by the corresponding load schedules and predictable enough.
Thus, the estimate of stability regions is another problem for future research. The basis for such studies is system (7.3).
The problem of voltage regulation is the topic of future work. For this regulation, external controls can be introduced into system (7.2). Then the system (7.2) can be taken in one of the following forms: either Ü (t)

Conclusion and Analysis of Results
One of the main problems arising in modeling any dynamic process is the problem of determination of dimension of phase space in which this process occurs. In article [35], which is devoted to the study of chaotic processes in the self-exciting homopolar disc dynamo, for modeling of the dynamics of this system three and five dimensional systems of differential equations were used.
It should be said that researches fulfilled in [35] are based on the known models, for which Problems 1 -4 were already solved (see Section 1). A purpose of these publications it is the search of hidden attractors and establishment of their properties.
Note that the problems considered in [35] can be raised and for system (6.6). However, it is possible only when the results of verifications and tests fully will confirm adequacy of system (6.6) and the direct current traction power supply system. This adequacy can be set by the recurrence analysis methods [10].
At the recurrence analysis of recurrence plots an important role play lengths of diagonal lines (we will emphasize that on recurrence plots the length of line characterizes a response time of trajectory in some region of phase space) [5], [10], [13], [18]. We performed such analysis of diagonal lines, but its results were rather rough.
It should be said that in the general case it is impossible to achieve a good correspondence between the model and the complex process that this model describes. Therefore, in this work, the adequacy of the model and the process was evaluated by the deviations of the current and voltage obtained in the simulation and experiment (a current-voltage characteristic U − I).
Comparison of the experimental information on Fig. 4.2 and solutions of system equations (4.1) shows that there is satisfactory description of dynamics of the direct current power supply system on the interval 4000 seconds -10000 seconds. (However, it is necessary to notice that among all systems of equations (4.1) -(4.7) only the phase portrait of system (4.1) most adequate to the real phase portrait on Fig. 4.2. Thus, there is good quality coincidence of the experimental U − I characteristic with U − I characteristic of system (4.1).) In order to attain a greater accuracy it is necessary to specify the coefficients of system (4.1) (or (6.6)). It is possible to do by the use of artificial neural networks. In the future authors hope to get back to this problem.