Microsoft Word - 1301_final


Adv Syst Sci Appl 2022; 04; 79-91 
Published online at http://ijassa.ipu.ru. 

 

Bifurcation in the Model of Cargo 
Transportation Organization 

Nerses K. Khachatryan1* 
1) State Academic University for the Humanities, Moscow, Russia 

E-mail: nerses-khachatryan@yandex.ru 

Abstract: This article is devoted to the study of the model of cargo transportation organization 
between two nodal stations. The main characteristic of an arbitrary station is the degree of 
inconsistency between receiving and sending cargo, which is the difference between the volume of 
incoming and outgoing cargo per unit of time. The initial node station accepts goods depending on 
the demand for them within its technical potential, which is determined by the maximum allowable 
increase in the degree of inconsistency between the reception and dispatch of goods per unit of 
time. The movement of goods from one station to another is carried out within the framework of 
their technical potentials. The distribution of goods from the final node station is carried out in a 
certain mode. Such a model is described by a system of differential equations with a number of 
parameters that define the characteristics of the demand for cargo transportation, the degree of use 
of the technical potential of the stations and the mode of cargo distribution from the final node 
station. When changing the parameters of the model, a bifurcation effect occurs and multiple 
stationary solutions appear, among which the most acceptable from the point of view of economic 
feasibility are identified. 

Keywords: cargo transportation organization model, differential equations, stationary solutions, 
model parameters, stability 

1. INTRODUCTION 

The mathematical models used for the analysis of transport networks are diverse in terms of 
the tasks to be solved, the mathematical apparatus, the data used and the degree of detail of the 
description of traffic. Based on the functional role of models, i.e. on the tasks for which they 
are used, three main classes can be conditionally distinguished: predictive, simulation and 
optimization models [24]. 

Predictive models allow to determine what will be the traffic flows in the network with 
known geometry and properties of the transport network. The forecast of the load of the 
transport network includes the calculation of the average characteristics of traffic, such as the 
volume of inter-district movements, the intensity of the flow, the distribution of vehicles along 
the paths, etc. With the help of such models, it is possible to predict the consequences of 
changes in the transport network or in the placement of objects. 

Simulation modeling aims to reproduce all the details of the movement, while the averaged 
value of the flows and the distribution along the paths are considered known and serve as the 
initial data for these models. Thus, flow forecasting and simulation modeling are 
complementary directions [9]. 

The third class of models is aimed at optimizing the functioning of transport networks. 
With their help the problems of optimization of transport routes, development of the optimal 
network configuration, etc. are solved [10, 23, 26, 27]. 

According to the method of describing traffic flows, all models of transport networks can 
be divided into classes: models-analogues, models-following the leader, models-probabilistic. 

 
* Corresponding author: nerses-khachatryan@yandex.ru 



80  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

In analog models, the movement of a vehicle is likened to some kind of physical flow 
(hydro- and gas-dynamic models) [6, 7, 9]. 

In models of following the leader, it is important to assume that there is a connection 
between the movement of the driven and the main vehicle [5]. 

In probabilistic models, the traffic flow is considered as a result of the interaction of 
vehicles on the elements of the transport network. Due to the rigid nature of network 
restrictions and the mass nature of traffic in the traffic flow, there are clear patterns of 
formation of queues, intervals, loads on the road lanes, etc. These patterns are significantly 
stochastic [25]. 

One of the most popular modes of transport for cargo transportation in Russia is rail. 
Publications devoted to railway logistics can be divided into the following main groups 
according to the type of tasks studied: 

1) the tasks of designing the infrastructure of the railway network; 
2) tasks of managing the fleet of locomotives and wagons; 
3) tasks of railway planning. 
In the first group we can highlight the works [8, 11, 19, 22]. The task of managing the fleet 

of locomotives and wagons is described in [4]. The third group, in particular, is represented by 
the tasks of forming the schedule of freight trains and the organization of freight flows [20, 
21]. 

One of the approaches to the organization of cargo traffic is described in [1, 2, 12-16]. They 
present macroscopic dynamic models in which the process of organizing railway freight 
transportation is the formation of cargo traffic based on the interaction of neighboring stations. 
By their functional role, they are predictive, because they allow predicting the dynamics of 
loading of the stations and flows arising in the railway network, at the set procedure of the 
organization of a cargo flow. In the way of describing traffic flows, they are close to the models 
of following the leader, if the flow on a particular section of the railway network is identified 
with the congestion of the respective stations. At the same time, these models have a significant 
difference from the models of following the leader, which means that each station interacts not 
with one, but with the two nearest neighboring stations. 

In [17] and [18], another approach to the organization of cargo traffic is presented, in which 
the process of its formation is determined by the demand for cargo transportation and the 
degree of use of the technical potential of stations. 

This work is devoted to the study of the modification of the model described in [17], and 
unlike it, it is initially aimed at finding practically feasible modes of cargo transportation. Such 
models are characterized by the global stability of a stationary solution, which is sometimes 
economically impractical. At the same time, when the model parameters change, a bifurcation 
effect occurs and multiple stationary solutions appear, among which there are economically 
feasible ones. Thus, in this model, the presence of bifurcation makes it possible to obtain 
solutions that describe the optimal functioning of the system. 

2. PROBLEM STATEMENT 

Consider the movement of goods on a section of the railway network between two nodal 
stations connected by a set of intermediate stations. Denoting the number of intermediate 
stations by 𝑚 , we get the following set of station numbers {0,1, . . , 𝑚, 𝑚 + 1}, where 0 – is 
the number of the initial node station, and, а 𝑚 + 1 – is the number of the final node station. 
Denote by 𝑛  the number of roads at the station with the number 𝑖. We assume that all paths 
are used with the same degree of efficiency. 

Consider discrete moments of time 𝑡 , 𝑡 , 𝑡 , . . . ;   𝑡 = 𝑡 + Δ𝑡, 𝑘 = 1,2, . .. 
Let 𝑉 (𝑡 ) is the volume of cargo received on the 𝑗th road of the 𝑖th station for a period of 

time [𝑡 , 𝑡 ], and 𝑉 (𝑡 ) is the volume of cargo sent from the 𝑖th road of the 𝑖th station for 
a period of time [𝑡 , 𝑡 ]. Denote 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  81 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

𝑥 (𝑡 ) =

𝑉 (𝑡 ) − 𝑉 (𝑡 )

𝑉 (𝑡 )
,  𝑖𝑓        𝑉 (𝑡 ) > 𝑉 (𝑡 )

0,  𝑖𝑓        𝑉 (𝑡 ) ≤ 𝑉 (𝑡 ).

 

It's obvious that 

0 ≤ 𝑥 (𝑡 ) ≤ 1,    𝑖 = 0,1, . . . , 𝑚 + 1;     𝑗 = 1,2, . . . , 𝑛  

and characterize the degree of inconsistency between receiving and sending on the 𝑗th road 
of the 𝑖th station at time 𝑡 . Denote  

𝑧 (𝑡 ) =
1

𝑛
𝑥 (𝑡 ). 

It is also obvious that 0 ≤ 𝑧 (𝑡 ) ≤ 1  and characterizes the degree of inconsistency 
between receiving and sending at the 𝑖th station at time 𝑡 .  

The technical potential of the station is determined by the maximum allowable increase in 
the degree of inconsistency between the reception and dispatch of goods per unit of time and 
is given by a non-negative decreasing function ( )z  defined on the segment [0,1]  and 
satisfying the condition 𝜑(1) = 0. 

The initial node station accepts cargo depending on the demand for transportation within 
its technical potential and sends it to the next station within its technical potential. Each of the 
intermediate stations receives cargo within its technical potential and sends it within the 
technical potential of the next station. The final node station accepts cargo within its technical 
potential and distributes it in a certain mode. Taking into account the above, we will write 
down a system of finite-difference equations describing the change in the degree of 
inconsistency between the reception and dispatch of goods at stations.  

𝑧 (𝑡 ) = 𝑧 (𝑡 ) + min 𝑑 , 𝜑 𝑧 (𝑡 ) − 𝜆𝜑 𝑧 (𝑡 ) Δ𝑡,    𝑘 = 1, 2, 3 … (2.1) 

𝑧 (𝑡 ) = 𝑧 (𝑡 ) + 𝜆𝜑 𝑧 (𝑡 ) − 𝜆𝜑 𝑧 (𝑡 ) Δ𝑡,   𝑖 = 1, 𝑚,    𝑘 = 1, 2, 3 … (2.2) 

𝑧 (𝑡 ) = 𝑧 (𝑡 ) + 𝜆𝜑 𝑧 (𝑡 ) − 𝑑 Δ𝑡,    𝑘 = 1, 2, 3 …   (2.3) 

0 ≤ 𝑧 (𝑡 ) ≤ 1,    𝑖 = 0, 1, … , 𝑚 + 1,     𝑘 = 0, 1, 2, …     (2.4) 

Here 𝑑 > 0,    0 < 𝜆 ≤ 1,    𝑑 > 0 are the model parameters: 
𝑑  is a characteristics of the demand for transportation; 
𝜆 is a characteristics of the degree of use of the technical potential of the stations; 
𝑑  is a characteristics of the cargo distribution mode from the final node station. 
Let's move on to the continuous analog of the system of discrete-difference equations 

(2.1)–(2.4), presented below 

�̇� (𝑡) = min(𝑑 , 𝜑(𝑧 (𝑡))) − 𝜆𝜑(𝑧 (𝑡)),      𝑡 ∈ [𝑡 , +∞),   (2.5) 

�̇� (𝑡)   = 𝜆[𝜑(𝑧 (𝑡)) − 𝜑(𝑧 (𝑡))],       𝑖 = 1,2, . . . , 𝑚,      𝑡 ∈ [𝑡 , +∞),  (2.6) 

�̇� (𝑡) = 𝜆𝜑(𝑧 (𝑡)) − 𝑑  ,        𝑡 ∈ [𝑡 , +∞),    (2.7) 

0 ≤ 𝑧 (𝑡) ≤ 1 ,      𝑖 = 0, 1, . . ., 𝑚 + 1,      𝑡 ∈ [𝑡 , +∞).    (2.8) 

Next, consider the function that sets the technical potential of the stations, of the following 
type 

𝜑(𝑧) = 𝑎(1 − 𝑧),     𝑎 > 0.     (2.9) 



82  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

The parameter 𝑎 > 0 , which participates in the definition of the function 𝜑(𝑧) , is a 
characteristic of the ability of stations to increase cargo traffic. Since 𝜑(𝑧) ≤ a  for all 0 ≤
𝑧 ≤ 1 then the parameter 𝑑 , which is a characteristic of the demand for transportation and is 
involved in equation (1.1), can be represented as follows: 

𝑑 = 𝜇𝑎,    0 < 𝜇 ≤ 1.     (2.10) 

The parameter   𝑑 , which is a characteristic of the cargo distribution mode from the 
final node station, is represented as follows: 

𝑑 = 𝛾𝑎,    𝛾 > 0.      (2.11) 

Let's rewrite the system (2.5)-(2.8), in which the function 𝜑(𝑧) is defined according to 
(2.9), and the parameters 𝑑  and 𝑑  are defined according to (2.10) and (2.11), respectively. 

�̇� (𝑡) = min(𝜇𝑎, 𝑎(1 − 𝑧 (𝑡))) − 𝜆𝑎(1 − 𝑧 (𝑡)),      𝑡 ∈ [𝑡 , +∞), (2.12) 

�̇� (𝑡) = 𝜆𝑎 (𝑧 (𝑡) − 𝑧 (𝑡)),      𝑖 = 1,2, . . . , 𝑚,      𝑡 ∈ [𝑡 , +∞),  (2.13) 

�̇� (𝑡) = 𝜆𝑎(1 − 𝑧 (𝑡)) − 𝛾𝑎 ,      𝑡 ∈ [𝑡 , +∞),   (2.14) 

0 ≤ 𝑧 (𝑡) ≤ 1 ,      𝑖 = 0,1, . . . , 𝑚 + 1,      𝑡 ∈ [𝑡 , +∞).   (2.15) 

Here 𝜇,    𝑎,    𝜆,    𝛾 are the model parameters: 
𝜇       (0 < 𝜇 ≤ 1) is a characteristics of the range of demand for transportation, which can 

be satisfied with the existing technical potential of the stations; 
𝑎        (𝑎 > 0) is a characteristic of the ability of stations to increase cargo traffic; 
𝜆        (0 < 𝜆 ≤ 1) is a characteristics of the degree of use of the technical potential of the 

stations; 
𝛾       (𝛾 > 0) is a characteristics of the cargo distribution mode from the final node station. 
Here are the main objectives of the study:  
- to determine the ranges of change of parameters 𝜇,    𝑎,    𝜆,    𝛾, in which the system of 

cargo transportation can function smoothly, i.e. system (2.12)-(2.15) has a solution, 
describe the qualitative behavior of solutions depending on the parameters.  

- for a given value of the demand characteristics for cargo transportation (parameter 𝜇) 
set the most acceptable achievable levels of the degree of inconsistency between the 
reception and dispatch of goods at all stations, by controlling the values of the following 
characteristics: the ability of stations to increase cargo traffic (parameter 𝑎), the degree 
of use of the technical potential of stations (parameter 𝜆) and the mode of cargo 
distribution from the final node station (parameter 𝛾). 

3. INVESTIGATION OF SYSTEM SOLUTIONS (2.12)-(2.15) 

The study of the set of solutions of the system (2.12)–(2.15) begins with the study of all 
solutions of the system of differential equations (2.12)–(2.14). 

First of all, we will highlight the stationary solutions of the system (2.12)–(2.14). With the 
help of direct verification, you can verify the validity of the following statement 

Proposition 3.1: 

System (2.12)–(2.14) for any parameter values 

0 < 𝜇 ≤ 1,    𝑎 > 0,     0 < 𝜆 ≤ 1,     𝛾 > 0 

such that 𝛾 ≤ 𝜇 has stationary solutions: 

𝑧 (. ) ≡ 1 − 𝛾,       𝑧 (. ) ≡ 1 − ,     𝑖 = 1, . . . , 𝑚 + 1, for 𝛾 < 𝜇;  (3.1) 

𝑧 (. ) ≤ 1 − 𝜇,   𝑧 (. ) ≡ 1 − ,      𝑖 = 1, . . . , 𝑚 + 1, for 𝛾 = 𝜇.  (3.2) 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  83 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

For 𝛾 > 𝜇 the system (2.12)–(2.14) has no stationary solutions. 
Let's proceed to the study of the remaining solutions of the system (2.12)–(2.14). 
Theorem 3.1:  

An arbitrary solution of the system (2.12)–(2.14) at 𝛾 < 𝜇 eventually goes to a stationary 
solution (3.1), and at 𝛾 = 𝜇 to one of the stationary solutions (3.2). For 𝛾 > 𝜇 coordinates 
𝑧 (. ),      𝑖 = 1, . . ., 𝑚 + 1  of solutions of the system (2.12)–(2.14) eventually enter the 
stationary mode specified in (3.1), and the function 𝑧 (. ) becomes linearly decreasing. 

Proof. 
Let is find the general solution of the system (2.12)–(2.14). Let is start with the last equation, 

which contains one variable (𝑧 ). Let 's rewrite it in the following form  

�̇� (𝑡) + 𝜆𝑎𝑧 (𝑡) = (𝜆 − 𝛾)𝑎.    (3.3) 

It is not difficult to verify that linear equation (3.3) has the following general solution  

𝑧 (𝑡) = 1 − + 𝑐 𝑒 .    (3.4) 

Substituting the expression for 𝑧  from (3.4) into the penultimate equation of the system 
(2.12)–(2.14), we find its general solution. It has the following form 

𝑧 (𝑡) = 1 − + 𝑒 (𝜆𝑎𝑐 𝑡 + 𝑐 ).   (3.5) 

Similarly, we will find general solutions to all other equations of the system (2.12)-(2.14) 
except the initial one:  

𝑧 (𝑡) = 1 − + 𝑒 ( 𝑡 + 𝑐 𝑡 + 𝑐 ); 

… 

𝑧 (𝑡) = 1 − + 𝑒 (
( )!

𝑡 + 𝑐 𝑡 + ⋯ + 𝑐 𝑡 + 𝑐 );   (3.6) 

𝑧 (𝑡) = 1 − + 𝑒 (
!

𝑡 + 𝑐 𝑡 + ⋯ + 𝑐 𝑡 + 𝑐 ). 

It follows from (3.5) and (3.6) that  

lim
→

𝑧 (𝑡) = 1 − ,      𝑖 = 1, 2, … , 𝑚. 

Let's move on to solving the first equation of the system (2.12)-(2.14). Let 's rewrite it in 
the following form 

�̇� (𝑡) =
𝜇𝑎 − 𝜆𝑎 1 − 𝑧 (𝑡) ,   if  𝑧 (𝑡) < 1 − 𝜇,   𝑡 ∈ [𝑡 , +∞),

𝑎 1 − 𝑧 (𝑡) − 𝜆𝑎 1 − 𝑧 (𝑡) ,   if  𝑧 (𝑡) ≥ 1 − 𝜇 ,    𝑡 ∈ [𝑡 , + ∞).
 (3.7) 

Consider the following two equations 

�̇� (𝑡) = 𝜇𝑎 − 𝜆𝑎 1 − 𝑧 (𝑡)    , 𝑡 ∈  𝑡, +∞ ,     (3.8) 

�̇� (𝑡) = 𝑎 1 − 𝑧 (𝑡) − 𝜆𝑎 1 − 𝑧 (𝑡) ,      𝑡 ∈  𝑡, +∞ ,     (3.9) 

where 𝑡 ≥ 𝑡 . 
Using the expression for 𝑧 (𝑡) from (3.6), we obtain the solution of equations (3.8) and 

(3.9). They can be represented as follows 

𝑧 (𝑡) = 𝑎(𝜇 − 𝛾)𝑡 + 𝐹 (𝑡) + 𝑐 , where 𝐹 (𝑡) ∈ С  𝑡, +∞ ,   lim
→

𝐹 (𝑡) = 0, (3.10) 

𝑧 (𝑡) = 1 − 𝛾 + 𝐹 (𝑡), where 𝐹 (𝑡) ∈ С  𝑡, +∞ ,  lim
→

𝐹 (𝑡) = 0           (3.11) 

These solutions allow us to investigate the asymptotic behavior of the solution of equation 
(3.7). It is easy to see that when 𝛾 < 𝜇 the asymptotics of the solution of equation (3.7) is 



84  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

determined by the relation (3.11), i.e. lim
→

𝑧 (𝑡) = 1 − 𝛾 . For 𝛾 > 𝜇  the asymptotics of 

equation (3.7) is determined by the relation (3.10), i.e. 𝑧 (. )  decreases linearly and 
lim
→

𝑧 (𝑡) = − ∞ . For 𝛾 = 𝜇  the asymptotics of equation (3.7), depending on the initial 

conditions, can be determined by both equation (3.10) and equation (3.11), i.e. either 
lim
→

𝑧 (𝑡) = с , where с < 1 − 𝜇 or lim
→

𝑧 (𝑡) = 1 − 𝜇. □ 

Let us proceed to the study of solutions of system (2.12)-(2.14) satisfying constraints (2.15).  
Lemma 3.1:  

For all values of parameters 𝑎 > 0 , 𝛾 > 0  and 𝜆 , 𝛾 ≤ 𝜆 ≤ 1  components 
𝑧 (. ),   𝑧 (. ), … 𝑧 (. )  of an arbitrary solution of the system (2.12)-(2.14) satisfying the 
constraints (2.15) at the initial moment of time will satisfy them at subsequent moments of time. 

Proof. 
Let's start by considering the last component of the solution of the system (2.12)-(2.14), i.e. 

𝑧 (. ). It has the form (3.4), where с  is determined from the condition 

1 − + 𝑐 𝑒 = 𝑧 , where 0 ≤ 𝑧 ≤ 1, 

i.e.  

𝑐 = ( − 1 + 𝑧 )𝑒 , where 0 ≤ 𝑧 ≤ 1.  (3.12) 

From (3.12) follows 

( − 1)𝑒 ≤ 𝑐 ≤ 𝑒  . 

Using this estimate for 𝑐  and expression (3.4), we get an estimate for 𝑧 (. ). 
It will take the following form 

(1 − )(1 − 𝑒 𝑒 ) ≤ 𝑧 (𝑡) ≤ 1 − (1 − 𝑒 𝑒 ).  (3.13) 

It follows from (3.13) that for all 𝛾 ≤ 𝜆 there is an inequality 

0 ≤ 𝑧 (𝑡) ≤ 1,     𝑡 ∈ [𝑡 , +∞) .    (3.14) 

We show that for the other components of the solution of the system (2.12)-(2.14), 
inequalities similar to inequality (3.14) are fulfilled. Let's start with the component 𝑧 (. ). To 
do this, consider equation (2.13) for 𝑖 = 𝑚: 

�̇� (𝑡) = 𝜆𝑎 𝑧 (𝑡) − 𝑧 (𝑡) ,      𝑡 ∈ [𝑡 , +∞).      

We show that the function 𝑧 (. ) cannot take a value greater than 1. Indeed, otherwise, due 
to the continuity of the function 𝑧 (. ) there must be a point  𝑡∗ > 𝑡 , such that 𝑧 ( 𝑡∗) = 1. 
Then it follows from (3.14) that �̇� ( 𝑡∗) ≤ 0. Similarly, the function 𝑧 (. ) cannot take a value 
less than 0. Thus, it is proved that for all 𝛾 ≤ 𝜆 function 𝑧 (. ) also satisfies an inequality 
similar to inequality (3.14). Similarly, the satisfiability of all other inequalities is proved. □ 

We formulate a similar lemma for the zero component of the solution of system (2.12)-
(2.14). 

Lemma 3.2:  

For all values of parameters 𝑎 > 0, 𝛾 > 0 and 𝜇, 𝛾 ≤ 𝜇 ≤ 1 there exists 𝜆(𝜇, 𝑧 (𝑡 ),

𝑧 (𝑡 ), . . ., 𝑧 (𝑡 )), 𝜇 ≤ 𝜆 ≤ 1 such that for any the value of the parameter 𝜆 from the 
segment [𝜇, 𝜆] the zero component 𝑧 (. )  of an arbitrary solution of the system (2.12)-(2.14) 
satisfying the constraint (2.15) at the initial moment of time, satisfies it at subsequent moments 
of time.  

 

 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  85 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

Proof. 
As for the other components, we show that the function 𝑧 (. ) cannot take a value greater 

than 1. Indeed, otherwise, due to the continuity of the function 𝑧 (. )  there must be a 
point 𝑡∗∗ > 𝑡 , such that 𝑧 ( 𝑡∗∗) = 1. Then it follows from (2.12) that 

�̇� ( 𝑡∗∗) =  𝜆𝑎(𝑧 ( 𝑡∗∗) − 1), 

that is, according to lemma 1, �̇� ( 𝑡∗∗) ≤ 0. 
Let's move on to evaluating the function 𝑧 (. ) from below. To do this, we investigate the 

behavior of its derivative at 𝑧 (. ) → 0 +. According to (2.12), it is described by the equation 

�̇� (𝑡) = 𝑎(𝜇 − 𝜆(1 − 𝑧 (𝑡))). 

Let 's investigate the inequality 

𝜇 − 𝜆(1 − 𝑧 (𝑡)) ≥ 0. 

Let 's rewrite it in the form 

𝑧 (𝑡) ≥ 1 − .    (3.15) 

According to lemma 1, for arbitrary 𝛾 > 0,   𝜆 satisfying the condition   𝛾 ≤ 𝜆 ≤ 1 there is 
an inequality 

0 ≤ 𝑧 (𝑡) ≤ 1,     𝑡 ∈ [𝑡 , +∞) .   (3.16) 

 It follows from (3.16) that an arbitrary μ satisfying the condition 𝛾 ≤ 𝜇 ≤ 1 exists 𝜆, 𝜇 ≤

𝜆 ≤ 1 such that for any value of the parameter 𝜆 from the segment [𝜇, 𝜆] inequality (3.15) 
will hold for all 𝑡 ∈ [𝑡 , +∞), i.e. �̇� (𝑡) ≥ 0 for 𝑧 (𝑡) → 0 +, which shows that the function 
𝑧 (. )  is limited from below by the value 0. Obviously 𝜆  depends on both 𝜇 , and initial 
conditions, so we denote it 𝜆(𝜇, 𝑧 (𝑡 ), 𝑧 (𝑡 ), . . ., 𝑧 (𝑡 )). □ 

We formulate the main result of this paragraph. 
Theorem 3.2:  

For any initial values 0 ≤ 𝑧 (𝑡 ) ≤ 1,      𝑖 = 0, 1, . . ., 𝑚 + 1, parameters 𝑎 > 0, 𝛾 > 0 
and 𝜇, 𝛾 ≤ 𝜇 ≤ 1    there exists 𝜆(𝜇, 𝑧 (𝑡 ), 𝑧 (𝑡 ), . . ., 𝑧 (𝑡 )) , 𝜇 ≤ 𝜆 ≤ 1  such that 
for any value of the parameter 𝜆 from the segment [𝜇, 𝜆] the solution of the system (2.12)-
(2.15) exists and converges to the stationary solution (3.1) (for 𝛾 < 𝜇 ) or to one of the 
stationary solutions (3.2), the same for all γ and a (for 𝛾 = 𝜇).  

Proof. The proof follows directly from Theorem 1, Lemma 1 and Lemma 2. □ 
Corollary 3.1:  

The system of differential equations (2.12)–(2.15) has a globally stable stationary solution 
(3.1) and a family of stable solutions of the form (3.2).  

Proof. The proof follows directly from Theorem 2. □ 

4. EXAMPLES OF SOLUTIONS OF THE SYSTEM (2.12)–(2.15) 

Let's consider some examples of solutions of the system (2.12)–(2.15) that clearly 
demonstrate the results of the previous paragraph. In all the examples below, the number of 
stations is 10 (the initial node station, eight intermediate stations, the final node station), the 
values of parameters, 𝜇, 𝑎 and the initial conditions are fixed and take the following values:  

𝜇 = 0.6,    𝑎 = 1.5     (4.1) 

𝑧 (𝑡 ) = 0.4,    𝑧 (𝑡 ) = 0.2,    𝑧 (𝑡 ) = 0.1,    𝑧 (𝑡 ) = 0.4,    𝑧 (𝑡 ) = 0.7, 

𝑧 (𝑡 ) = 0.1,    𝑧 (𝑡 ) = 0.2,    𝑧 (𝑡 ) = 0.3,    𝑧 (𝑡 ) = 0.5,    𝑧 (𝑡 ) = 0.3. (4.2) 



86  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

According to theorem 3.2, for any initial values satisfying condition (2.15), parameters 𝑎 >

0, 𝛾 > 0 and 𝜇, 𝛾 ≤ 𝜇 ≤ 1   there exists 𝜆 (𝜇 ≤ 𝜆 ≤ 1) depending on the initial conditions 
and the parameter 𝜇, such that for all 𝜆 from the segment [𝜇, 𝜆] the solution of the system 
(2.12)–(2.15) exists. Calculations have shown that for the above initial conditions (4.2) and 
parameter 𝜇 (4.1), 𝜆 takes the value equal to 0.898 (𝜆 = 0.898). The nature of the solution of 
the system (2.12)–(2.15) depends on the value of the parameter 𝛾. Depending on whether it is 
less than the value of the parameter 𝜇 or equal to it, there are two types of solutions of the 
system (2.12)–(2.15).  

Let's start with the first case (𝛾 < 𝜇). In fig. 4.1, fig. 4.2 and fig. 4.3. graphs of solutions of 
the system (2.12)–(2.15) with a fixed value of parameter 𝛾 (𝛾 = 0.5) and three different values 
of parameter 𝜆 are given: two at the ends of the segment [𝜇, 𝜆] and one at an internal point.  

 

Fig. 4.1. Graph of the solution of the system (2.12)–(2.15) 
(𝛾 = 0.5,    𝜇 = 𝜆 = 0.6). 

 

Fig. 4.2. Graph of the solution of the system (2.12)–(2.15)  
(𝛾 = 0.5, 𝜇 = 0.6,    𝜆 = 0.85) 

 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  87 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

 

Fig. 4.3. Graph of the solution of the system (2.12)–(2.15)  
(𝛾 = 0.5,    𝜇 = 0.6,    𝜆 = 𝜆 = 0.898) 

Note that due to the global stability of the stationary solution (3.1), the asymptotic behavior 
of the solutions shown in figures 4.1–4.3 does not change when the initial values change. These 
figures clearly show that the steady – state mode for the zero component of the solution of the 
system (2.12)–(2.15) does not depend on the parameter 𝜆, and for the remaining components, 
an increase in the parameter 𝜆 leads to their asymptotic increase. Recall that the components 
of the solution of the system (2.12)–(2.15) determine the dynamics of the degree of 
inconsistency between the reception and dispatch of goods at stations. Therefore, in this case, 
the question of choosing a parameter 𝜆 is uniquely determined, its value should be taken equal 
to the value of the parameter 𝜇. 

Let's move on to the second case (𝛾 = 𝜇). Fig. 4.4, fig. 4.5 and fig. 4.6 show graphs of 
solutions of the system (2.12)–(2.15) at 𝛾 = 0.6 and the same values of parameter 𝜆, as for the 
previous case.  

 

Fig. 4.4. Graph of the solution of the system (2.12)–(2.15) 
𝛾 = 𝜆 = 𝜇 = 0.6 



88  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

 

Fig. 4.5. Graph of the solution of the system (2.12)–(2.15) 
𝛾 = 𝜇 = 0.6,    𝜆 = 0.85 

 

Fig. 4.6. Graph of the solution of the system (2.12)–(2.15) 
𝛾 = 𝜇 = 0.6,    𝜆 = 𝜆 = 0.898 

Comparing figures 4.4, 4.5 and 4.6, it can be seen that an increase in the value of parameter 
𝜆 leads to an asymptotic decrease in the zero component of the solution of the system (2.12)–
(2.15), in contrast to the first case. The remaining components of the solution, as in the first 
case, increase asymptotically with an increase in the value of the parameter 𝜆. This means that 
in this case, by controlling the parameter 𝜆 you can set the desired degree of inconsistency 
between the reception and dispatch of goods at all stations (including at the zero station). If we 
take 𝜆  equal to 𝜇 , then over time the degree of inconsistency between the reception and 
dispatch of goods at all stations except zero will become zero, and at zero station 1 − 𝜇. If the 
value of 𝜇 is close to one, then this choice of the parameter 𝜆 will be optimal. Otherwise, 
everything will depend on the specific value of the parameter 𝜇, as well as on how important 
it is in a given situation to reduce the degree of inconsistency between receiving and sending 
goods at the initial node station by increasing this characteristic at other stations.  

Comparing figures 4.1 and 4.4, figures 4.2 and 4.5, as well as figures 4.3 and 4.6, it can be 
concluded that in terms of minimizing the degree of inconsistency between the reception and 
dispatch of goods at stations, it is advisable to take the parameter 𝛾 equal to the parameter 𝜇. 
This means that the mode of cargo distribution from the final node station must be coordinated 
with the characteristics of the demand for transportation. 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  89 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

5. INVESTIGATION OF THE ATTRACTION AREA OF STATIONARY 
SOLUTIONS OF THE SYSTEM (2.12)–(2.15) OF THE FORM (3.2) 

According to consequence 1, the system of differential equations (2.12)–(2.15) has a 
globally stable stationary solution (3.1) and a family of stable solutions of the form (3.2). We 
select the following stationary solution from this family  

𝑧 (. ) ≡ 1 − 𝜇,    𝑧 (. ) ≡ 1 − ,     𝑖 = 1, . . . , 𝑚 + 1.  (5.1) 

As numerical experiments have shown, this stationary solution is in contrast to other 
stationary solutions of the form (3.2) presented below 

𝑧 (. ) < 1 − 𝜇,   𝑧 (. ) ≡ 1 − ,     𝑖 = 1, . . . , 𝑚 + 1   (5.2) 

and being locally stable, has a certain area of attraction.  
Denote  

с(𝜇, 𝜆) =
( )( / )

. 

The analysis of a large number of numerical experiments made it possible to describe the 
region of attraction of the stationary solution (5.1). The results of this analysis are given in the 
following proposition. 

Proposition 5.1: 

For any parameter values   0 < 𝜇 ≤ 1,    𝑎 > 0,        𝜇 ≤ 𝜆 ≤ 𝜆 ,    𝛾 = 𝜇    the solution of 
the system (2.12)–(2.15) converges to the stationary solution (5.1) if the initial values satisfy 
the condition  

𝑧 (𝑡 ) ≥ �̃�(𝜇, 𝜆),    𝑖 = 0, 1, . . . , 𝑚 + 1. 

Otherwise, i.e. if the condition𝑧 (𝑡 ) ≤ �̃�(𝜇, 𝜆),    𝑖 = 0, 1, . . . , 𝑚 + 1 and ∃ 𝚤̅ ∈ {0, 1, . . . .,
𝑚 + 1} such that 

𝑧 ̅(𝑡 ) < �̃�(𝜇, 𝜆) 

is met, then the solution of the system (2.12)–(2.15) converges to one of the stationary 
solutions (5.2). □ 

6. CONCLUSION 

This article presents a model of organization of cargo transportation between two nodal 
stations, described by a system of differential equations with a number of parameters that 
define the characteristics of demand for cargo transportation, the degree of use of the technical 
potential of the stations and the mode of cargo distribution from the final nodal station. The 
ranges of parameter changes under which the cargo transportation system can function 
smoothly are determined. For a given value of the demand characteristics for cargo 
transportation, by controlling the degree of use of the technical potential of the stations and 
the mode of cargo distribution from the final node station, the most acceptable achievable 
levels of the degree of inconsistency between the reception and dispatch of goods at all stations 
are established. 

REFERENCES 

[1] Beklaryan L.A. & Khachatryan N.K. (2006). Traveling wave type solutions in 
dynamic transport models, Functional Differential Equations, 13(12), 125–155. 

[2] Beklaryan L.A. & Khachatryan N.K. (2013). On One Class of Dynamic 
Transportation Models, Computational Mathematics and Mathematical Physics, 
53(10), 1466–1482, https://doi.org/10.1134/S0965542513100035. 



90  N.K. KHACHATRYAN 

Copyright ©0000 ASSA Adv. in Systems Science and Appl. (0000) 

[3] Beklaryan L.A., Khachatryan N.K. & Akopov A.S. (2019). Model for organization 
cargo transportation at resource restrictions. International Journal of Applied 
Mathematics, 32(4), 627–640, http://dx.doi.org/10.12732//ijam.v32i4.7. 

[4] Belousov F.A., Nevolin I.V. & Khachatryan N.K. (2020). Modeling and 
optimization of plans for railway freight transport performed by a transport operator, 
Business Informatics, 14(2), 21–35, https://doi.org/10.17323/2587-814X. 
2020.2.21.35. 

[5] Brackstone, M. & McDonald M. (1999). Car following: A historical review, 
Transportation Research Part F: Traffic Psychology and Behavior, 2, 181–196. 

[6] Daganzo, C.F. (1994). The cell transmission model: A dynamic representation of 
highway traffic consistent with the hydrodynamic theory, Transportation Research 
Part B: Methodological, 28, 269–287. 

[7] Daganzo, C.F. (1995). The cell transmission model. Part II: Network traffic, 
Transportation Research Part B: Methodological, 29(2), 79–93. 

[8] Ferreira L., Murray M.H. (1997). Modelling Rail Track Deterioration and 
Maintenance: Current Practices and Future Needs, Transport Reviews, 17(3), 207–
221. 

[9] Gasnikov A.V, Klenov S. L., Nurminskii E. A., Kholodov Ya.A. & Shamrai N.B 
(2013). Introduction to mathematical model operation of traffic flows, Moscow: 
MCCME (in Russian). 

[10] Galaburda V.G. (1985). Optimal Planning of Cargo Traffic. [Optimal’noe 
planirovanie gruzopotokov.] Moscow: Transport (in Russian).  

[11] Higgins A., Ferreira L. & Kozan E. (1995). Modeling Single-Line Train 
Operations, Transportation Research Record, 1489, 9–16. 

[12] Khachatryan, N.K. (2013). Dinamicheskaya model' organizatsii 
gruzoperevozok pri ogranichennosti emkostey peregonnykh putey [Dynamic 
model of organization of cargo transportation with a limited caracity of the 
distillation ways] Biznes-informatika, 26(4), 62–68 [in Russian]. 

[13] Khachatryan N.K., Akopov A.S. (2017). Model for Organizing Cargo 
Transportation with an Initial Station of Departure and a Final Station of Cargo 
Distribution, Business Informatics, 1, 25–35, https://doi.org/10.17323/1998-
0663.2017.1.25.35. 

[14] Khachatryan N.K., Akopov A.S. & Belousov F.A. (2018). About Quasi-
Solutions of Traveling Wave Type in Models for Organizing Cargo Transportation, 
Business Informatics, 43(1), 61–70, https://doi.org/10.17323/1998-
0663.2018.1.61.70  

[15] Khachatryan N.K., Beklaryan G.L., Borisova S.V. & Belousov F.A. (2019). 
Research into the dynamics of railway track capacities in a model for organizing 
cargo transportation between two node stations, Business Informatics, 13(1), 59–
70, https://doi.org/10.17323/1998-0663.2019.1.59.70.  

[16] Khachatryan N. & Beklaryan L. (2021). Study of flow dynamics in the model 
of cargo transportation organization along a circular chain of stations, Economics 
and Mathematical Methods, 57(1), 83–91 [in Russian], 
http://dx.doi.org/10.31857/S042473880013024-5. 



                 BIFURCATION IN THE MODEL OF CARGO TRANSPORTATION ORGANIZATION  91 

Copyright ©0000 ASSA.  Adv. in Systems Science and Appl. (0000) 

 

[17] Khachatryan Nerses K. (2020). Study of flow dynamics in the model of cargo 
transportation organization between node stations// International Journal of 
Applied Mathematics. 33(5), 937–949, http://dx.doi.org/10.12732/ijam.v33i5.14. 

[18] Khachatryan, Nerses K. (2021). Modeling the process of cargo transportation 
between node stations, International Journal of Applied Mathematics, 34(6), 1223–
1235, http://dx.doi.org/10.12732/ijam.v34i6.12. 

[19] Kraay D., Barker P.T. & Chen B.T. (1991). Optimal Pacing of Trains in Freight 
Railroads: Model Formulation and Solution, Operations Research, 39(1), 82–99. 

[20] Lazarev A., Musatova E., Grafov E. & Kvaratskheliya A. (2012). Teoriya 
raspisaniy. Zadachi zheleznodorojnogo planirovaniya [Schedule theory. Problems 
of railway planning]. Moscow, Russia: ICS RAS, [in Russian]. 

[21] Lazarev A.A. (2009). Estimates of the absolute error and a scheme for an 
approximate solution to scheduling problems, Computational Mathematics and 
Mathematical Physics, 49(2), 373–386, https://doi.org/10.1134/ 
S0965542509020158. 

[22] LeBlanc L.J. (1976). Global Solutions for a Nonconvex Nonconcave Rail 
Network Model, Management Science, 23(2), 131–139. 

[23] Leventhal T., Nemhauser G.L. & Trotter L. (1973). A Column Generation 
Algorithm for Optimal Traffic Assignment, Transportation Science, 7, 168–176.  

[24] Shvetsov V. I. (2003). Mathematical Modeling of Traffic Flows, Automation 
and Remote Control, 64(11), 1651–1689. 

[25] Solomon, H. & Wang P. (1972). Nonhomogeneous Poisson fields of Random 
Lines with Applications to Traffic Flow. Proceedings of the Sixth Berkeley 
Symposium on Mathematical Statistics and Probability, 3, 383–400. 

[26] Steenbrink P.A. (1981). Optimizaciya transportnih system. [Optimization of 
Transport Networks]. Moscow: Transport [in Russian]. [Translated from 
Steenbrink P.A. (1974). Optimization of Transport Networks. New York: Wiley.]  

[27] Vasil’eva E.M., Igudin R.V. & Livshits V.N. (1987). Optimizaciya 
planirovaniya i upravleniya transportnimi sistemami. [Optimization of planning 
and control of transport systems]. Moscow: Transport [in Russian].