KINETICS OF BIOLOGICAL PROCESSES

How can the dynamics of biological systems be described? At each moment in time, a biological system has a set of certain characteristics. For example, observing a population of a species, you can record its size, area occupied by the territory, the amount of available food, temperature environment etc. Leakage chemical reaction can be characterized by the concentrations of the participating substances, pressure, temperature, and the level of acidity of the medium. The set of values ​​of all characteristics that the researcher has chosen to describe the system is the state of the system at each moment of time. When creating a model, variables and parameters are selected in the specified set. Variables are those quantities whose changes are primarily of interest to the researcher, parameters are the conditions of the "external environment". It is for the selected variables that equations are compiled that reflect the patterns of change in the system over time. For example, when creating a model for the growth of a culture of microorganisms, its number usually acts as a variable, and the reproduction rate as a parameter. It is possible that the temperature at which growth occurs will turn out to be significant, then this indicator is also included in the model as a parameter. And if, for example, the level of aeration is always sufficient and does not have any effect on growth processes, then it is not included in the model at all. As a rule, the parameters remain unchanged during the experiment, however, it is worth noting that this is not always the case.

It is possible to describe the dynamics of a biological system (that is, the change in its state over time) using both discrete and continuous models. Discrete models assume that time is discrete quantity. This corresponds to recording the values ​​of variables at certain fixed intervals (for example, once an hour or once a year). In continuous models, a biological variable is a continuous function of time, denoted, for example, x(t).

Often great importance have initial conditions models - the state of the characteristic under study at the initial moment of time, i.e. at t = 0.

When studying the continuous change of some characteristic x(t) we may know information about the rate of its change. This information can generally be written as a differential equation:

Such a formal notation means that the rate of change of some characteristic under study is a function of time and the magnitude of this characteristic.

If the right side of a differential equation of the form does not explicitly depend on time, i.e. fair:

then this equation is called autonomous(the system described by such an equation is called autonomous). The state of autonomous systems at each moment of time is characterized by one single value - the value of the variable x at the moment t.

Let us ask ourselves a question: let a differential equation be given for x(t), is it possible to find all functions x(t) satisfying this equation? Or: if the initial value of a certain variable is known (for example, the initial size of a population, the concentration of a substance, the electrical conductivity of the medium, etc.) and there is information about the nature of the change in this variable, is it possible to predict what its value will be at all subsequent points in time? The answer to the question posed is as follows: if the initial conditions are given for and the conditions of the Cauchy theorem are satisfied for the equation (the function given in a certain region and its partial derivative are continuous in this region), then there is a unique solution of the equation that satisfies the given initial conditions. (Recall that any continuous function that satisfies a differential equation is called a solution to that equation.) This means that we can uniquely predict the behavior of a biological system if the characteristics of its initial state are known and the model equation satisfies the conditions of Cauchy's theorem.

Stationary state. Sustainability

We will consider the autonomous differential equation

In a stationary state, the values ​​of variables in the system do not change with time, that is, the rate of change in the values ​​of variables is 0: . If the left side of equation (1.2) is equal to zero, then the right one is also equal to zero: . The roots of this algebraic equation are stationary states differential equation (1.2).

Example1.1: Find the stationary states of the equation.

Solution: Let's move the term, which does not contain the derivative, to the right side of the equality: . By definition, in a stationary state, the following equality holds: . So the equality must hold . We solve the equation:

,

So, the equation has 3 stationary states: , .

Biological systems constantly experience various external influences and numerous fluctuations. At the same time, they (biological systems) have homeostasis, i.e. resistant. In mathematical language, this means that the variables, with small deviations, return to their stationary values. Will this behavior of a biological system be reflected in its mathematical model? Are the stationary states of the model stable?

The steady state is sustainable, if, for a sufficiently small deviation from the equilibrium position, the system will never go far from the singular point. The stable state corresponds to the stable mode of the system operation.

An equilibrium state of an equation is Lyapunov stable if for any one can always find such that if , then for all .

There is an analytical method for studying stability steady state– Lyapunov method. To substantiate it, we recall Taylor formula.

Speaking loosely, the Taylor formula shows the behavior of a function in the vicinity of a certain point. Let a function have derivatives at a point of all orders up to n- th inclusive. Then the Taylor formula is valid for:

Discarding the remainder term, which represents itself an infinitesimal of a higher order than , we obtain the approximate Taylor formula:

The right side of the approximate formula is called Taylor polynomial functions, it is denoted as .

Example 1.2: Expand the function in a Taylor series in a neighborhood of a point up to 4th order.

Solution: We write the Taylor series up to the 4th order in general form:

Let's find derivatives given function at point :

,

Substitute the obtained values ​​into the original formula:

An analytical method for studying the stability of a stationary state ( Lyapunov method) is as follows. Let be the stationary state of the equation . Let's set a small deviation of the variable x from its stationary value: , where . Substitute the expression for the point x into the original equation: . The left side of the equation will take the form: , since in the stationary state the rate of change of the value of the variable is equal to zero: . We expand the right side into a Taylor series in the vicinity of the stationary state, taking into account that , we leave only the linear term on the right side of the equation:

Got linearized equation or first approximation equation. There is some value constant, let's denote it a: . The general solution of the linearized equation has the form: . This expression describes the law according to which the deviation from the stationary state given by us will change in time. The deviation will decay over time, i.e. at , if the exponent in the exponent is negative, i.e. . By definition, the steady state will be sustainable. If , then with increasing time, the deviation will only increase, the stationary state is unstable. In the case when the equation of the first approximation cannot give an answer to the question of the stability of the stationary state. It is necessary to consider higher-order terms in the Taylor series expansion.

In addition to the analytical method of studying the stability of a stationary state, there is also a graphical one.

Example 1.3. Let . Find the stationary states of the equation and determine their type of stability using the graph of the function .

Solution: Let's find special points:

,

,

We build a graph of the function (Fig. 1.1).

Rice. 1.1. Function graph (example 1.3).

Let us determine from the graph whether each of the found stationary states is stable. Let's set a small deviation of the representative point from the singular point to the left: . At a point with a coordinate, the function takes a positive value: or . The last inequality means that over time the coordinate should increase, that is, the representative point should return to the point . Now let's set a small deviation of the representative point from the singular point to the right: . In this region, the function retains a positive value, therefore, over time, the coordinate x also increases, that is, the representative point will move away from the point. Thus, a small deviation takes the system out of the stationary state, therefore, by definition, the singular point is unstable. Similar reasoning leads to the fact that any deviation from the singular point decays with time, the stationary state is stable. The deviation of the representative point in any direction from the stationary state leads to its removal from the point , this is an unstable stationary state.

Solving a system of linear differential equations

Let us turn to the study of systems of equations, first linear. In general, the system of linear differential equations can be represented as:

The analysis of the system of equations begins with finding the stationary states. For systems of the form (1.3), the singular point is unique, its coordinates are (0,0). The exception is the degenerate case, when the equations can be represented as:

(1.3*)

In this case, all pairs satisfying the relation are stationary points of the system (1.3*). In particular, the point (0,0) is also stationary for the system (1.3*). On the phase plane, in this case, we have a straight line with a slope coefficient passing through the origin, each point of which is a singular point of the system (1.3 *) (see Table 1.1, item 6).

The main question that should be answered by the result of studying a system of equations is whether the stationary state of the system is stable, and what character this solution has (monotonic or nonmonotonic).

Common decision system of two linear equations has the form:

characteristic numbers can be expressed in terms of the coefficients of linear equations as follows:

Characteristic numbers can be 1) real of different signs, 2) real of the same sign, 3) complex conjugate, and also, in degenerate cases, 4) purely imaginary, 5) real coinciding, 6) real, one of which (or both) zero. These cases determine the type of behavior of the solution to a system of ordinary differential equations. The corresponding phase portraits are presented in Table 1.1.

Table 1.1. Types of stationary states of a system of two linear differential equations and the corresponding phase portraits. The arrows show the direction of movement of the representative point

Construction of phase and kinetic portraits of a system of two linear differential equations

phase plane called a plane with coordinate axes on which the values ​​of variables are plotted x and y, each point of the plane corresponds to a certain state of the system. The set of points on the phase plane, the position of which corresponds to the states of the system in the process of changing variables in time, according to the given equations of the system under study, is called phase trajectory. The set of phase trajectories for different initial values ​​of the variables gives a portrait of the system. Building phase portrait allows you to draw conclusions about the nature of changes in variables x and y without knowing the analytical solutions of the original system of equations.

Consider a system of linear differential equations:

The construction of the phase portrait begins with the construction main isoclines(an isocline is a line throughout which the slope of the phase curve (trajectory), determined by the equation, remains constant). For a system of two linear differential equations, these are always straight lines passing through the origin. The equation isoclines of horizontal tangents: . Equation of the isocline of vertical tangents: . For further construction of the phase portrait, it is useful to construct an isocline of tangents passing at an angle . To find the corresponding isocline equation, it is necessary to solve the equation . You can also find the isoclines of tangents of other angles, using the approximate values ​​of the tangents of the angles. The answer to the question at what angle the phase trajectories should intersect the coordinate axes can also help in constructing the phase portrait. To do this, in the isocline equation we substitute the corresponding equalities (to determine the angle of intersection with the OY axis) and (to determine the angle of intersection with the OX axis).

Example 1.4. Determine the type of singular point of the system of linear equations:

Construct a phase and kinetic portrait of the system.

Solution: The singular point coordinates are (0,0). The coefficients of the linear equations are: , , , . Let us define the type of stationary state (see the section on characteristic numbers):

Thus, the characteristic roots are imaginary: therefore, the singular point of the considered linear system has the center type (Fig. 1.2a).

Equation of the isocline of horizontal tangents: , Equation of the isocline of vertical tangents: . At an angle of 45°, the trajectories of the system intersect a straight line .

After constructing the phase portrait, it is necessary to determine the direction of movement along the found trajectories. This can be done in the following way. Take an arbitrary point on any trajectory. For example, on the isocline of horizontal tangents (1,1). Let us substitute the coordinates of this point into the system of equations. We obtain expressions for the rates of change of variables x,y at this point:

The obtained values ​​show that the rate of change of the variable x- negative, that is, its value should decrease, and the variable y does not change. We mark the received direction with an arrow. Thus, in the example under consideration, the motion along the phase trajectories is directed counterclockwise. Substituting the coordinates of different points into the system, you can get a "map" of the directions of velocities, the so-called vector field.

Fig 1.2. Phase (a) and kinetic (b) portrait of the system, example 1.4

Note that on the isocline of the horizontal tangents the variable y reaches its maximum or minimum value on a given trajectory. On the contrary, on the isocline of vertical tangents, the variable x.

To build a kinetic portrait of the system means to plot the dependence of the values ​​of the variables x,y from time. A phase portrait can be used to build a kinetic one and vice versa. One phase trajectory corresponds to one pair of kinetic curves. Let us choose an arbitrary point on the phase portrait on an arbitrary phase trajectory. This is the starting point corresponding to time . Depending on the direction of movement in the system under consideration, the values ​​of the variables x,y either decrease or increase. Let the coordinates of the starting point be (1,1). According to the constructed phase portrait, starting from this point, we must move counterclockwise, the coordinates x and y while they will decrease. Over time, the coordinate x passes through 0, value y while remaining positive. Further coordinates x and y continue to decrease, the coordinate y goes through 0 (value x while it is negative). Value x reaches its minimum value on the isocline of vertical tangents, then begins to increase. Value y reaches its minimum value on the isocline of horizontal tangents (value x at this point in time is negative). Next, and the value x, and the value y increase, returning to the initial values ​​(Fig. 1.2b).

Investigation of the Stability of Stationary States of Second-Order Nonlinear Systems

Let a biological system be described by a system of two autonomous differential equations of the second order general view:

The stationary values ​​of the system variables are determined from the algebraic equations:

In the neighborhood of each stationary state, one can consider first approximation system(linearized system), the study of which can make it possible to answer the question of the stability of the singular point and the nature of the phase trajectories in its small neighborhood.

outside

We have , , the singular point is rough. The characteristic roots of the system of the first approximation are equal to , both are real and negative, therefore, in the vicinity of the zero singular point, the behavior of the phase trajectories of the system will correspond to the type of a stable knot.

Introduction

Since the concept of a nonlinear dynamical system is rich enough to cover an extremely wide range of processes in which the future behavior of the system is determined by the past, the methods of analysis developed in this field are useful in a huge variety of contexts.

Nonlinear dynamics enters the literature in at least three ways. First, there are cases where experimental data on the change over time of one or more quantities are collected and analyzed using techniques based on non-linear dynamic theory, with minimal assumptions about the underlying equations that govern the process that produces the data. That is, it is a case in which one seeks to find correlations in the data that can guide the development of a mathematical model, rather than first guessing the model and then comparing it to the data.

Secondly, there are cases where nonlinear dynamical theory can be used to state that some simplified model should exhibit important features of this system, from which it follows that the describing model can be built and studied in a wide range of parameters. This often results in models that behave qualitatively differently under different parameters and demonstrate that one region exhibits behavior that is very similar to that observed in the real system. In many cases, the behavior of the model is quite sensitive to changes in parameters, so if the parameters of the model can be measured in a real system, the model exhibits realistic behavior at these values, and one can be sure that the model captures the essential features of the system.

Thirdly, there are cases when the model equations are built on the basis of detailed descriptions known physics. Numerical experiments can then provide information about variables not available to physical experiments.

Based on the second path, this work is an extension of my previous work “Nonlinear dynamic model of interdependent industries”, as well as another work (Dmitriev, 2015)

All the necessary definitions and other theoretical information needed in the work will appear in the first chapter, as needed. Two definitions will be given here, which are necessary for the disclosure of the research topic itself.

First, let's define system dynamics. According to one definition, system dynamics is a simulation modeling approach that, through its methods and tools, helps to evaluate the structure complex systems and their dynamics (Shterman). It is worth adding that system dynamics is also a modeling technique that is used to recreate correct (in terms of accuracy) computer models for complex systems for their future use in order to create a more effective company / organization, as well as improve methods of interaction with this system. Predominantly, the need for system dynamics arises when confronted with long-term, strategic models, and it is also worth noting that it is rather abstract.

Speaking of non-linear differential dynamics, we will consider a non-linear system, which, by definition, is a system in which the change in the result is not proportional to the change in the input parameters, and in which the function describes the dependence of change in time and the position of a point in space (Boeing, 2016).

Based on the above definitions, it becomes clear that this work will consider various non-linear differential systems that describe the interaction of companies, as well as simulation models built on their basis. Based on this, the purpose of the work will be determined.

Thus, the purpose of this work is to conduct a qualitative analysis of dynamic systems that describe the interaction of companies in the first approximation and build a simulation model based on them.

To achieve this goal, the following tasks were identified:

Determining the stability of the system.

Construction of phase portraits.

Finding integral trajectories of systems.

Construction of simulation models.

Each of these tasks will be devoted to one of the sections of each of the chapters of the work.

Based on practice, the construction of fundamental mathematical structures that effectively model the dynamics in various physical systems and processes, indicates that the corresponding mathematical model to some extent reflects the proximity to the original under study, when its characteristics can be derived from the properties and structure of the type of motion that forms the dynamics of the system. To date, economic science is at a stage of its development, in which new, and in many cases, non-standard methods and methods of physical and mathematical modeling are especially effectively used in it. economic processes. This is where the conclusion follows about the need to create, study and build models that can somehow describe the economic situation.

As for the reason for choosing qualitative rather than quantitative analysis, it is worth noting that in the vast majority of cases, the results and conclusions from a qualitative analysis of dynamic systems are more significant than the results of their quantitative analysis. In such a situation, it is appropriate to point to the statements of V.P. Milovanov, in which he states that they traditionally believe that the results expected when applying mathematical methods to the analysis of real objects should be reduced to a numerical result. In this sense, qualitative methods have a somewhat different task. It focuses on achieving a result that describes the quality of the system, on the search for characteristic features of all phenomena as a whole, on forecasting. Of course, it is important to understand how demand will change when prices for a certain type of goods change, but do not forget that it is much more important to understand whether there will be a shortage or surplus of these goods under such conditions (Dmitriev, 2016).

The object of this study is nonlinear differential and system dynamics.

In this case, the subject of research is the description of the process of interaction between companies through non-linear differential and system dynamics.

Speaking about the practical application of the study, it is worth immediately dividing it into two parts. Namely, theoretical, that is, a qualitative analysis of systems, and practical, in which the construction of simulation models will be considered.

The theoretical part of this study provides basic concepts and phenomena. It considers simple differential systems, as in the works of many other authors (Teschl, 2012; Nolte, 2015), but at the same time allows describing the interaction between companies. Based on this, in the future it will be possible to conduct more in-depth studies, or else begin your acquaintance with what constitutes a qualitative analysis of systems.

The practical part of the work can be used to create a decision support system. Decision Support System - Automated Information system, aimed at supporting business or decision making in an organization, allowing you to choose between many different alternatives (Keen, 1980). Even if the models are not very accurate at the moment, but by changing them for a specific company, you can achieve more accurate results. Thus, when changing in them various parameters and conditions that can arise in the market, you can get a forecast for the future and make a more profitable decision in advance.

1. Interaction of companies in the conditions of mutualism

The paper will present two-dimensional systems that are quite simple in comparison with systems of a higher order, but at the same time allow us to demonstrate the relationships between organizations that we need.

It is worth starting work with the choice of the type of interaction, which will be described in the future, since for each of the types the systems describing them, albeit slightly, are different. Figure 1.1 shows Eujima Odum's classification for population interaction modified for economic interaction (Odum, 1968), based on which we will further consider the interaction of companies.

Figure 1.1. Types of interaction between enterprises

Based on Figure 1.1, we single out 4 types of interaction and present for each of them a system of equations describing them based on the Malthus model (Malthus, 1798). According to it, the growth rate is proportional to the current abundance of the species, in other words, it can be described by the following differential equation:

where a is a parameter that depends on the natural growth of the population. It is also worth adding that in the systems considered below, all parameters, as well as variables, take non-negative values.

The production of raw materials is the production of products, which is similar to the predator-prey model. The predator-prey model, also known as the Lotka-Volterra model, is a pair of non-linear first-order differential equations describing the dynamics of a biological system with two species, one of which is predator and the other is prey (Llibre, 2007). The change in the abundance of these species is described by the following system of equations:

(1.2)

where - characterizes the growth of the production of the first enterprise without the influence of the second one (in the case of the predator-prey model, the growth of the prey population without predators),

It characterizes the growth of production of the second enterprise without the influence of the first (growth of the population of predators without prey),

It characterizes the growth of the production of the first enterprise, taking into account the influence of the second enterprise on it (an increase in the number of prey when interacting with predators),

It characterizes the growth of production of the second enterprise, taking into account the influence of the first enterprise on it (an increase in the number of predators during their interaction with victims).

For one, the predator, as can be seen from the system, as well as Odum's classification, their interaction imposes a favorable effect. On the other unfavorable. If considered in economic realities, then, as can be seen in the figure, the simplest analogue is the manufacturer and its supplier of resources, which correspond to the predator and prey, respectively. Thus, in the absence of raw materials, output decreases exponentially.

Competition is rivalry between two or more (in our case, we are considering two-dimensional systems, so we take exactly two-species competition) species, economic groups for territories, limited resources, or other values ​​(Elton, 1968). Changes in the number of species, or the number of products in our case, are described by the system below:

(1.3)

In this case, species or companies that produce one product adversely affect each other. That is, in the absence of a competitor, product growth will increase exponentially.

Now let's move on to a symbiotic interaction, in which both enterprises have a positive influence on each other. Let's start with mutualism. Mutualism is a type of relationship between different species in which each of them benefits from the actions of the other, and it is worth noting that the presence of a partner is a necessary condition for existence (Thompson, 2005). This type of relationship is described by the system:

(1.4)

Since interaction between companies is necessary for their existence, in the absence of the product of one company, the output of the goods of another decreases exponentially. This is possible when companies simply have no other alternatives for procurement.

Consider another type of symbiotic interaction, protocooperation. Proto-cooperation is similar to mutualism, with the only exception that there is no need for a partner to exist, since, for example, there are other alternatives. Since they are similar, their systems look almost identical to each other:

(1.5)

Thus, the absence of one company's product does not hinder the growth of another company's product.

Of course, in addition to those listed in paragraphs 3 and 4, other types of symbiotic relationships can be noted: commensalism and amensalism (Hanski, 1999). But they will not be mentioned further, since in commensalism one of the partners is indifferent to his interaction with the other, but we still consider cases where there is influence. And amensalism is not considered, because from an economic point of view, such relations, when their interaction harms one, and the other is indifferent, simply cannot exist.

Based on the influence of companies on each other, namely the fact that symbiotic relationships lead to a stable coexistence of companies, in this paper we will consider only cases of mutualism and proto-cooperation, since in both cases the interaction is beneficial to everyone.

This chapter is devoted to the interaction of companies in the conditions of mutualism. It will consider two systems that are a further development of systems based on the Malthus model, namely systems with imposed restrictions on the increase in production.

The dynamics of a pair connected by mutualistic relationships, as mentioned above, can be described in the first approximation by the system:

(1.6)

It can be seen that with a large initial amount of production, the system grows indefinitely, and with a small amount, production falls. This is where the incorrectness of the bilinear description of the effect arising from mutualism lies. To try to correct the picture, let us introduce a factor that resembles the saturation of a predator, that is, a factor that will reduce the growth rate of production, if it is in excess. In this case, we arrive at the following system:

(1.7)

where is the growth in the production of the product of the first company in its interaction with the second, taking into account saturation,

Growth in the production of the product of the second company in its interaction with the first, taking into account saturation,

Saturation coefficients.

Thus, we got two systems: the Malthusian model of growth with and without saturation.

1.1 Stability of systems in the first approximation

The stability of systems in the first approximation is considered in many foreign (Hairer, 1993; Bhatia, 2002; Khalil, 2001; Strogatz, 2001 and others) and Russian-language works (Akhromeyeva, 1992; Bellman, 1954; Demidovich, 1967; Krasovsky, 1959 and others), and its definition is a basic step for analyzing the processes occurring in the system. To do this, perform the following necessary steps:

Let's find equilibrium points.

Let us find the Jacobian matrix of the system.

Find the eigenvalues ​​of the Jacobian matrix.

We classify the equilibrium points according to the Lyapunov theorem.

Having considered the steps, it is worth dwelling on their explanation in more detail, so I will give definitions and describe the methods that we will use in each of these steps.

The first step, the search for equilibrium points. To find them, we equate each function to zero. That is, we solve the system:

where a and b mean all parameters of the equation.

The next step is to find the Jacobian matrix. In our case, this will be a 2-by-2 matrix with first derivatives at some point, as shown below:

After completing the first two steps, we proceed to finding the roots of the following characteristic equation:

Where the point corresponds to the equilibrium points found in the first step.

Having found and , we pass to the fourth step and use the following Lyapunov theorems (Parks, 1992):

Theorem 1: If all roots of the characteristic equation have a negative real part, then the equilibrium point corresponding to the original and linearized systems is asymptotically stable.

Theorem 2: If at least one of the roots of the characteristic equation has a positive real part, then the equilibrium point corresponding to the original and linearized systems is asymptotically unstable.

Also, looking at and it is possible to determine the type of stability more precisely, based on the division shown in Figures 1.2 (Lamar University).

Figure 1.2. Types of stability of equilibrium points

Having considered the necessary theoretical information, we turn to the analysis of systems.

Consider a system without saturation:

It is very simple and not suitable for practical use, since it has no restrictions. But as a first example of system analysis is suitable for consideration.

First, let's find the equilibrium points by equating the right-hand sides of the equations to zero. Thus, we find two equilibrium points, let's call them A and B: .

Let's combine the step with the search for the Jacobian matrix, the roots of the characteristic equation, and the determination of the type of stability. Since they are elementary, we immediately get the answer:

1. At the point , , there is a stable knot.

At point: . . saddle.

As I already wrote, this system is too trivial, so no explanation was required.

Now let's analyze the system from saturation:

(1.9)

The appearance of a restriction on the mutual saturation of products by enterprises brings us closer to the real picture of what is happening, and also slightly complicates the system.

As before, we equate the right parts of the system to zero and solve the resulting system. The point remained unchanged, but the other point in this case contains more parameters than before: .

In this case, the Jacobi matrix takes the following form:

Subtract from it the identity matrix multiplied by , and equate the determinant of the resulting matrix at points A and B to zero.

At the point of a similar early picture:

stable node.

But at the point everything is somewhat more complicated, and although the mathematics is still quite simple, the complexity causes the inconvenience of working with long literal expressions. Since the values ​​turn out to be rather long and inconveniently written down, they are not given, it is enough to say that in this case, as with the previous system, the type of stability obtained is a saddle.

2 Phase portraits of systems

The vast majority of nonlinear dynamic models are complex differential equations that either cannot be solved, or this is some kind of complexity. An example is the system from the previous section. Despite the apparent simplicity, finding the type of stability at the second equilibrium point was not an easy task (albeit not from a mathematical point of view), and with an increase in parameters, restrictions and equations to increase the number of interacting enterprises, the complexity will only increase. Of course, if the parameters are numerical expressions, then everything will become incredibly simple, but then the analysis will somehow lose all meaning, because in the end, we will be able to find equilibrium points and find out their stability types only for a specific case, and not a general one.

In such cases, it is worth remembering the phase plane and phase portraits. In applied mathematics, in particular the context of nonlinear systems analysis, the phase plane is a visual representation of certain characteristics of certain types of differential equations (Nolte, 2015). Coordinate plane with axes of values ​​of any pair of variables characterizing the state of the system - a two-dimensional case of a common n-dimensional phase space.

Thanks to the phase plane, it is possible to graphically determine the existence of limit cycles in solutions of a differential equation.

The solutions of a differential equation are a family of functions. Graphically, this can be plotted in the phase plane as a two-dimensional vector field. Vectors are drawn on the plane, representing derivatives at characteristic points with respect to some parameter, in our case, with respect to time, that is (). With enough of these arrows in one area, the behavior of the system can be visualized and limit cycles can be easily identified (Boeing, 2016).

The vector field is a phase portrait, a particular path along the flow line (that is, a path always tangent to the vectors) is a phase path. Flows in a vector field indicate the change in the system over time, described by a differential equation (Jordan, 2007).

It is worth noting that a phase portrait can be built even without solving a differential equation, and at the same time, good visualization can provide a lot useful information. In addition, at the present time there are many programs that can help with the construction of phase diagrams.

Thus, phase planes are useful for visualizing the behavior of physical systems. In particular, oscillatory systems, such as the predator-prey model already mentioned above. In these models, phase trajectories can "twist" towards zero, "go out of a spiral" to infinity, or reach a neutral stable situation called centers. This is useful in determining whether dynamics are stable or not (Jordan, 2007).

The phase portraits presented in this section will be built using WolframAlpha tools, or provided from other sources. Malthusian growth model without saturation.

Let us build a phase portrait of the first system with three sets of parameters to compare their behavior. Set A ((1,1), (1,1)), which will be referred to as a single set, set B ((10,0.1), (2,2)), when selected, the system experiences a sharp decline in production , and the set C ((1,10), (1,10)) for which, on the contrary, a sharp and unlimited growth occurs. It should be noted that the values ​​along the axes in all cases will be in the same intervals from -10 to 10, for the convenience of comparing phase diagrams with each other. Of course, this does not apply to a qualitative portrait of the system, whose axes are dimensionless.

Figure 1.3 Phase portrait with parameters A

mutualism differential limit equation

Figure 1.3 above shows the phase portraits of the system for the three specified sets of parameters, as well as the phase portrait describing the qualitative behavior of the system. Do not forget that the most important from a practical point of view is the first quarter, since the amount of production, which can only be non-negative, is our axes.

In each of the figures, stability at the equilibrium point (0,0) is clearly visible. And in the first figure, the “saddle point” is also noticeable at the point (1,1), in other words, if we substitute the values ​​of the set of parameters into the system, then at the equilibrium point B. When the boundaries of the model construction change, the saddle point is also found on other phase portraits.

Malthusian model of growth from saturation.

Let us construct phase diagrams for the second system, in which there is saturation, with three new sets of parameter values. Set A, ((0.1,15,100), (0.1,15,100)), set B ((1,1,0.5), (1, 1,0.5)) and set C ((20,1,100), (20,1,100 )).

Figure 1.4. Phase portrait with parameters A

As you can see, for any set of parameters, the point (0,0) is equilibrium, and also stable. Also in some figures, you can see a saddle point.

In this case, different scales were considered in order to more clearly demonstrate that even when a saturation factor is added to the system, the qualitative picture does not change, that is, saturation alone is not enough. It should be taken into account that in practice, companies need stability, that is, if we consider non-linear differential equations, then we are most interested in stable equilibrium points, and in these systems, only zero points are such points, which means that such mathematical models are clearly not suitable for enterprises. . After all, this means that only with zero production, companies are in stability, which is clearly different from the real picture of the world.

In mathematics, an integral curve is a parametric curve that represents a particular solution to an ordinary differential equation or system of equations (Lang, 1972). If the differential equation is represented as a vector field, then the corresponding integral curves are tangent to the field at every point.

Integral curves are also known by other names, depending on the nature and interpretation of the differential equation or vector field. In physics, integral curves for an electric field or magnetic field are known as field lines, and the integral curves for the fluid velocity field are known as streamlines. In dynamic systems, integral curves for a differential equation are called trajectories.

Figure 1.5. Integral curves

Solutions of any of the systems can also be considered as equations of integral curves. Obviously, each phase trajectory is a projection of some integral curve into space x,y,t to the phase plane.

There are several ways to construct integral curves.

One of them is the isocline method. An isocline is a curve passing through points at which the slope of the function under consideration will always be the same, regardless of the initial conditions (Hanski, 1999).

It is often used as a graphical method for solving ordinary differential equations. For example, in an equation of the form y "= f (x, y), isoclines are lines in the (x, y) plane obtained by equating f (x, y) to a constant. This gives a series of lines (for different constants) along which the curves solutions have the same gradient.By calculating this gradient for each isocline, the slope field can be visualized, making it relatively easy to draw approximate solution curves.The figure below shows an example of using the isocline method.

Figure 1.6. Isocline method

This method does not require computer calculations, and was very popular in the past. Now there are software solutions that will build integral curves on computers extremely accurately and quickly. However, even so, the isocline method has shown itself well as a tool for studying the behavior of solutions, since it allows you to show the areas of typical behavior of integral curves.

Malthusian growth model without saturation.

Let's start with the fact that despite the existence of different construction methods, it is not so easy to show the integral curves of a system of equations. The isocline method mentioned earlier is not suitable because it works for first order differential equations. And software tools that have the ability to plot such curves are not in the public domain. For example, Wolfram Mathematica, which is capable of this, is paid. Therefore, we will try to use the capabilities of Wolfram Alpha as much as possible, the work with which is described in various articles and works (Orca, 2009). Even despite the fact that the picture will be clearly not entirely reliable, but at least it will allow you to show the dependence in the planes (x, t), (y, t). First, let's solve each of the equations for t. That is, we derive the dependence of each of the variables with respect to time. For this system we get:

(1.10)

(1.11)

The equations are symmetric, so we consider only one of them, namely x(t). Let the constant be equal to 1. In this case, we will use the plotting function.

Figure 1.7. Three-dimensional model for equation (1.10)

Malthusian model of growth from saturation.

Let's do the same for the other model. Ultimately, we obtain two equations that demonstrate the dependence of variables on time.

(1.12)

(1.13)

Let's build a three-dimensional model and level lines again.

Figure 1.8. Three-dimensional model for equation (1.12)

Since the values ​​of the variables are non-negative, then in the fraction with the exponent we get a negative number. Thus, the integral curve decreases with time.

Previously, a definition of system dynamics was given to understand the essence of the work, but now let's dwell on this in more detail.

System dynamics - methodology and method mathematical modeling for the formation, understanding and discussion of complex problems, originally developed in the 1950s by Jay Forrester and described in his work (Forrester, 1961).

System dynamics is one aspect of systems theory as a method for understanding the dynamic behavior of complex systems. The basis of the method is the recognition that the structure of any system consists of numerous relationships between its components, which are often as important in determining its behavior as the individual components themselves. Examples are chaos theory and social dynamics, described in the works of various authors (Grebogi, 1987; Sontag, 1998; Kuznetsov, 2001; Tabor, 2001). It is also argued that since whole-properties can often not be found in element properties, in some cases the behavior of the whole cannot be explained in terms of the behavior of the parts.

Simulation can really show the whole practical significance dynamic system. Although it is possible in spreadsheets, there are many software packages that have been optimized specifically for this purpose.

Modeling itself is the process of creating and analyzing a prototype of a physical model in order to predict its performance in the real world. Simulation modeling is used to help designers and engineers understand under what conditions and in what cases a process can fail and what loads it can withstand (Khemdy, 2007). Simulation can also help predict the behavior of fluid flows and other physical phenomena. The model analyzes the approximate working conditions due to the applied simulation software (Strogalev, 2008).

Limitations on the possibilities of simulation modeling have a common cause. The construction and numerical calculation of an exact model guarantees success only in those areas where there is an exact quantitative theory, i.e., when the equations describing certain phenomena are known, and the task is only to solve these equations with the required accuracy. In those areas where there is no quantitative theory, the construction of an exact model is of limited value (Bazykin, 2003).

However, the modeling possibilities are not unlimited. First of all, this is due to the fact that it is difficult to assess the scope of the simulation model, in particular, the period of time for which the forecast can be built with the required accuracy (Law, 2006). In addition, by its nature, the simulation model is tied to a specific object, and when trying to apply it to another, even similar object, it requires a radical adjustment or, at least, a significant modification.

There is a general reason for the existence of limitations on simulation. The construction and numerical calculation of an “exact” model is successful only if a quantitative theory exists, that is, only if all equations are known, and the problem is reduced only to solving these equations with a certain accuracy (Bazykin, 2003).

But even despite this, simulation modeling is an excellent means of visualizing dynamic processes, allowing, with a more or less correct model, to make decisions based on its results.

In this work, system models will be built using the system dynamics tools offered by the AnyLogic program.

Malthusian growth model without saturation/

Before building a model, it is necessary to consider the elements of system dynamics that we will use and relate them to our system. The following definitions have been taken from the help information of the AnyLogic program.

The drive is the main element of system dynamics diagrams. They are used to represent objects of the real world, in which certain resources accumulate: money, substances, numbers of groups of people, some material objects, etc. Accumulators reflect the static state of the simulated system, and their values ​​change over time in accordance with the flows existing in the system. It follows that the dynamics of the system are determined by the flows. Flows entering and leaving the accumulator increase or decrease the values ​​of the accumulator.

The flow, as well as the aforementioned drive, is the main element of system-dynamic diagrams.

While the bins define the static part of the system, the flows determine the rate of change of the bins, that is, how stocks change over time, and thus determine the dynamics of the system.

The agent may contain variables. Variables are typically used to model the changing characteristics of an agent or to store the results of the model. Typically, dynamic variables consist of accumulator functions.

The agent may have parameters. Parameters are often used to represent some of the characteristics of the modeled object. They are useful when object instances have the same behavior as described in the class, but differ in some parameter values. There is a clear difference between variables and parameters. The variable represents the state of the model and can change during simulation. The parameter is usually used to describe objects statically. During one "run" of the model, the parameter is usually a constant and is changed only when the behavior of the model needs to be reconfigured.

A link is an element of system dynamics that is used to determine the relationship between elements of a flow diagram and accumulators. It does not automatically create links, but forces the user to explicitly draw them in the graphical editor (however, it is worth noting that AnyLogic also supports a mechanism for quickly setting missing links). As an example, if any element of A is mentioned in the equation or the initial value of element B, then you first need to connect these elements with a link going from A to B, and only then enter the expression in the properties of B.

There are some other elements of system dynamics, but they will not be involved in the course of work, so we will omit them.

To begin with, let us consider what the model of system (1.4) will consist of.

First, we immediately mark two drives, which will contain the values ​​of the quantity of production of each of the enterprises.

Secondly, since we have two terms in each equation, we get two flows to each of the drives, one incoming, the other outgoing.

Thirdly, we pass to variables and parameters. There are only two variables. X and Y, responsible for the growth of production. We also have four options.

Fourth, with regard to connections, each of the flows must be associated with the variables and parameters included in the flow equation, and both variables must be associated with accumulators to change the value over time.

We will leave a detailed description of building a model, as an example of working in the AnyLogic modeling environment, for the next system, since it is somewhat more complicated and uses more parameters, and we will immediately proceed to consider the finished version of the system.

Figure 1.9 below shows the constructed model:

Figure 1.9. System dynamics model for system (1.4)

All elements of system dynamics correspond to those described above, i.e. two drives, four streams (two incoming, two outgoing), four parameters, two dynamic variables, and necessary links.

The figure shows that the more products, the stronger its growth, which leads to a sharp increase in the number of goods, which corresponds to our system. But as mentioned earlier, the absence of restrictions on this growth does not allow the application of this model in practice.

Malthusian growth model from saturation/

Considering this system, let us dwell on the construction of the model in more detail.

The first step is to add two drives, let's call them X_stock and Y_stock. We assign an initial value of 1 to each of them. Note that in the absence of flows in the classical given equation there is nothing in storage.

Figure 1.10. Building a system model (1.9)

The next step is adding threads. Let's build an incoming and outgoing stream for each drive using a graphical editor. We must not forget that one of the edges of the flow must be in the drive, otherwise, they will not be connected.

You can see that the equation for the drive was set automatically, of course, the user can write it himself by choosing the “arbitrary” equation mode, but the easiest way is to leave this action to the program.

Our third step is to add six parameters and two dynamic variables. Let's give each element a name in accordance with its literal expression in the system, and also set the initial values ​​of the parameters as follows: e1=e2=1, a12=a21=3, n1=n2=0.2.

All elements of the equations are present, it remains only to write the equations for the flows, but for this you first need to add connections between the elements. For example, the outgoing stream responsible for the term must be associated with e1 and x. And each dynamic variable must be associated with its corresponding stock (X_stock x, Y_stock y). Creating links is similar to adding threads.

After creating the necessary connections, you can proceed to writing equations for the flows, which is shown in the right figure. Of course, you can go in the reverse order, but if there are connections, when writing equations, hints appear for substituting the necessary parameters / variables, which makes the task easier in complex models.

After completing all the steps, you can run the simulation model and look at its result.

Having considered the systems of non-linear differential equations for the interaction of companies in the conditions of mutualism, we can draw several conclusions.

There are two states of the system: a sharp unlimited growth, or the tendency of the quantity of production to zero. Which of the two states the system will assume depends on the parameters.

None of the proposed models, including the model taking into account saturation, is not suitable for practical use, due to the lack of a non-zero stable position, as well as the reasons described in paragraph 1.

In the case of an attempt to further study this type of symbiotic interaction in order to create a model applicable by companies in practice, it is necessary to further complicate the system and introduce new parameters. For example, Bazykin in his book gives an example of the dynamics of two mutualistic populations with the introduction of an additional factor of intraspecific competition. Due to which the system takes the form:

(1.15)

And in this case, a non-zero stable position of the system appears, separated from the zero by a “saddle”, which brings it closer to the real picture of what is happening.

2. Interaction of companies in the conditions of proto-cooperation

All the basic theoretical information was presented in the previous chapter, so in the analysis of the models considered in this chapter, for the most part, the theory will be omitted, with the exception of a few points that we did not encounter in the previous chapter, and there may also be a reduction in calculations. The model of interaction between organizations considered in this chapter under conditions of proto-cooperation, which consists of systems of two equations based on the Malthusian model, looks like system (1.5). The systems analyzed in the previous chapter showed that for their maximum approximation to the existing models, it is necessary to complicate the systems. Based on these findings, we will immediately add a growth constraint to the model. Unlike the previous type of interaction, when growth that does not depend on another company is negative, in this case all signs are positive, which means that we have constant growth. Avoiding the shortcomings described earlier, we will try to limit it to the logistic equation, also known as the Verhulst equation (Gershenfeld, 1999), which has the following form:

, (2.1)

where P is the population size, r is the parameter showing the growth rate, K is the parameter responsible for the maximum possible population size. That is, over time, the population size (in our case, production) will tend to a certain parameter K.

This equation will help curb the rampant output growth we have seen so far. Thus, the system takes the following form:

(2.2)

Do not forget that the volume of goods stored in the warehouse for each company is different, so the parameters that limit growth are different. Let's call this system "", and in the future we will use this name when we consider it.

The second system we will consider is further development models with the Verhulst constraint. As in the previous chapter, we introduce a saturation constraint, then the system will take the form:

(2.3)

Now each of the terms has its own limitation, so without further analysis it can be seen that there will be no unlimited growth, as in the models of the previous chapter. And since each of the terms demonstrates positive growth, then the quantity of production will not fall to zero. Let's call this model the “two-constrained proto-operation model”.

These two models are discussed in various sources on biological populations. Now we will try to expand the systems somewhat. To do this, consider the following figure.

The figure shows an example of the processes of two companies: the steel and coal industries. In both enterprises there is an increase in production that is independent of the other, and also there is an increase in production, which is obtained due to their interaction. We have already taken this into account in earlier models. Now it is worth paying attention to the fact that companies not only produce products, they also sell them, for example, to the market or to a company interacting with it. Those. based on logical conclusions, there is a need for negative growth of companies due to the sale of products (in the figure, parameters β1 and β2 are responsible for this), as well as due to the transfer of part of the production to another enterprise. Previously, we took this into account only with a positive sign for another company, but did not consider the fact that the number of products decreases for the first enterprise when transferring products. In this case, we get the system:

(2.4)

And if it can be said about the term that if it was indicated in the previous models that , characterize the natural increase, and the parameter can be negative, then there is practically no difference, then about the term this cannot be said. In addition, in the future, when considering such a system with a restriction imposed on it, it is more correct to use the terms of positive and negative growth, since in this case different restrictions may be imposed on them, which is impossible for natural growth. Let's call it the "extended proto-cooperation model".

Finally, the fourth model under consideration is the extended proto-cooperation model with the previously mentioned logistic growth constraint. And the system for this model is as follows:

, (2.5)

where is the increase in the production of the first enterprise, independent of the second, taking into account the logistical constraint, - the increase in the production of the first company, depending on the second, taking into account the logistical constraint, - the increase in the production of the second enterprise, independent of the first, taking into account the logistical constraint, - increase in the production of the second company, depending on the first, taking into account the logistical constraint, - consumption of the goods of the first company, not related to another, - consumption of goods of the second company, not related to another, - consumption of goods of the first industry by the second industry, - consumption of goods of the second industry first industry.

In the future, this model will be referred to as the "extended proto-operation model with a logistic constraint".

1 Stability of systems in the first approximation

Proto-operation model with Verhulst constraint

Methods for analyzing the stability of the system were indicated in a similar section of the previous chapter. First of all, we find equilibrium points. One of them, as always, is zero. The other is a point with coordinates .

For the zero point λ1 = , λ2 = , since both parameters are non-negative, we obtain an unstable node.

Since it is not very convenient to work with the second point, due to the lack of the ability to shorten the expression, we will leave the definition of the type of stability to the phase diagrams, since they clearly show whether the equilibrium point is stable or not.

The analysis of this system is more complicated than the previous one due to the fact that the saturation factor is added, thus new parameters appear, and when finding equilibrium points, it will be necessary to solve not a linear, but a bilinear equation due to the variable in the denominator. Therefore, as in the previous case, we leave the definition of the stability type to phase diagrams.

Despite the appearance of new parameters, the Jacobian at the zero point, as well as the roots of the characteristic equation, looks similar to the previous model. Thus, at the zero point, an unstable node.

Let's move on to advanced models. The first of them does not contain any restrictions and takes the form of system (2.4)

Let's make a change of variables, , and . New system:

(2.6)

In this case, we get two equilibrium points, point A(0,0), B(). Point B lies in the first quarter because the variables have a non-negative value.

For the equilibrium point A we get:

. - unstable knot

. - saddle,

. - saddle,

. - stable knot

At point B, the roots of the characteristic equation are complex numbers: λ1 = , λ2 = . We cannot determine the type of stability relying on Lyapunov's theorems, so we will carry out numerical simulations that will not show all possible states, but will allow us to find out at least some of them.

Figure 2.2. Numerical simulation of the search for the type of stability

Considering this model, one will have to face computational difficulties, since it has a large number of various parameters, as well as two limitations.

Without going into details of calculations, we arrive at the following equilibrium points. Point A(0,0) and point B with the following coordinates:

(), where a =

For point A, determining the type of stability is a trivial task. The roots of the characteristic equation are λ1 = , λ2 = . Thus we get four options:

1. λ1 > 0, λ2 > 0 - unstable node.

2.λ1< 0, λ2 >0 - saddle.

3. λ1 ​​> 0, λ2< 0 - седло.

4.λ1< 0, λ2 < 0 - устойчивый узел.

Speaking about point B, it is worth agreeing that substituting abbreviations into the expression for it will complicate work with the Jacobian and finding the roots of the characteristic equation. For example, after trying to find them using WolframAlpha computing tools, the output of the roots took about five lines, which does not allow working with them in literal terms. Of course, if there are already existing parameters, it seems possible to quickly find an equilibrium point, but this is a special case, since we will find the equilibrium state, if any, only for these parameters, which is not suitable for the decision support system for which the model is planned to be created. .

Due to the complexity of working with the roots of the characteristic equation, we construct the mutual arrangement of zero-isoclines by analogy with the system analyzed in Bazykin's work (Bazykin, 2003). This will allow us to consider the possible states of the system, and in the future, when constructing phase portraits, to find equilibrium points and types of their stability.

After some calculations, the zero-isoclinic equations take the following form:

(2.7)

Thus, isoclines have the form of parabolas.

Figure 2.3. Possible null-isoclinic location

In total, there are four possible cases of their mutual arrangement according to the number of common points between the parabolas. Each of them has its own sets of parameters, and hence the phase portraits of the system.

2 Phase portraits of systems

Let us construct a phase portrait of the system, provided that and the remaining parameters are equal to 1. In this case, one set of variables is sufficient, since the quality will not change.

As can be seen from the figures below, the zero point is an unstable node, and the second point, if we substitute the numerical values ​​of the parameters, we get (-1.5, -1.5) - a saddle.

Figure 2.4. Phase portrait for the system (2.2)

Thus, since no changes should occur, then for this system there are only unstable states, which is most likely due to the possibility of unlimited growth.

A proto-operation model with two restrictions.

In this system, there is an additional limiting factor, so the phase diagrams must differ from the previous case, as can be seen in the figure. The zero point is also an unstable node, but a stable position appears in this system, namely a stable node. With these parameters, its coordinates (5.5,5.5), it is shown in the figure.

Figure 2.5. Phase portrait for the system (2.3)

Thus, the restriction on each term made it possible to obtain a stable position of the system.

Extended proto-operation model.

Let's build phase portraits for the extended model, but immediately using its modified form:

Let us consider four sets of parameters, such as to consider all cases with a zero equilibrium point, and also to demonstrate the phase diagrams of the numerical simulation used for a non-zero equilibrium point: the set A(1,0.5,0.5) corresponds to the state , set B(1,0.5,-0.5) corresponds to set C(-1.0.5,0.5) and set D(-1.0.5,-0.5) , that is, a stable node at the zero point. The first two sets will demonstrate the phase portraits for the parameters that we considered in the numerical simulation.

Figure 2.6. Phase portrait for system (2.4) with parameters А-D.

In the figures, it is necessary to pay attention to the points (-1,2) and (1,-2), respectively, a “saddle” appears in them. For a more detailed representation, the figure shows a different scale of the figure with a saddle point (1,-2). In the figure, at points (1,2) and (-1,-2), a stable center is visible. As for the zero point, starting from figure to figure on the phase diagrams, we can clearly distinguish an unstable node, a saddle, a saddle and a stable node.

Extended proto-cooperation model with logistic constraint.

As in the previous model, we will demonstrate phase portraits for four cases of a zero point, and also try to note nonzero solutions in these diagrams. To do this, take the following sets of parameters with the parameters specified in the following order (): A (2,1,2,1), B (2,1,1,2), C (1,2,2,1) and D (1,2,1,2). The remaining parameters for all sets will be as follows: , .

In the figures presented below, one can observe the four equilibrium states of the zero point described in the previous section for this dynamical system. And also in the figures, the stable position of a point with one non-zero coordinate.

Figure 2.7. Phase portrait for system (2.5) with parameters A-B

3 Integral trajectories of systems

Proto-operation model with Verhulst constraint

As in the previous chapter, we solve each of the differential equations separately and explicitly express the dependence of the variables on the time parameter.

(2.8)

(2.9)

It can be seen from the obtained equations that the value of each of the variables increases, which is demonstrated in the three-dimensional model below.

Figure 2.8. Three-dimensional model for equation (2.8)

This type of plot initially resembles the unsaturated 3D Malthusian model discussed in Chapter 1 in that it has a similar rapid growth, but later on you can see a decrease in the growth rate as the output limit is reached. Thus the final appearance of the integral curves is similar to the plot of the logistic equation that was used to limit one of the terms.

A proto-operation model with two restrictions.

We solve each of the equations using Wolfram Alpha tools. Thus, the dependence of the function x(t) is reduced to the following form:

(2.10)

For the second function, the situation is similar, so we omit its solution. Numerical values ​​appeared due to the replacement of the parameters by certain appropriate values, which does not affect the qualitative behavior of the integral curves. The charts below show the use of limits on growth as exponential growth becomes logarithmic over time.

Figure 2.9. Three-dimensional model for equation (2.10)

Extended proto-operation model

Almost similar to the models with mutualism. The only difference is faster growth relative to those models, which can be seen from the equations below (if you look at the degree of the exponent) and graphs. The integral curve must take the form of an exponent.

(2.11)

(2.12)

Extended proto-cooperation model with logistic constraint

Dependence x(t) looks like this:

Without a graph, it is difficult to evaluate the behavior of the function, so using the tools already known to us, we will build it.

Figure 2.10 3D Model for Equation

The value of the function decreases for non-small values ​​of another variable, which is due to the absence of restrictions on the negative bilinear term, and is an obvious result

4 System dynamics of interacting companies

Proto-operation model with Verhulst constraint.

Let us construct system (2.2). Using the tools already known to us, we build a simulation model. This time, unlike the mutualistic models, the model will have a logistical constraint.

Figure 2.11. System dynamics model for system (2.2)

Let's run the model. In this model, it is worth noting the fact that the growth from the relationship is not limited by anything, and the growth of output without the influence of the other has a specific limitation. If you look at the expression of the logistic function itself, you can see that in the case when the variable (number of goods) exceeds the maximum possible storage volume, the term becomes negative. In the case when there is only a logistic function, this is impossible, but with an additional always positive growth factor, this is possible. And now it is important to understand that the logistics function will cope with the situation of not too fast growth in the number of products, for example, linear. Let's take a look at the pictures below.

Figure 2.12. An example of the operation of the system dynamics model for system (2.2)

The left figure shows the 5th step of the program corresponding to the proposed model. But at the moment it is worth paying attention to the right figure.

First, for one of the incoming streams for Y_stock, the link to x, expressed in terms of , has been removed. This is done in order to show the difference in the performance of the model with a linear always positive flow, and bilinear growth, which is presented for X_stock. With linear unlimited flows, after exceeding the parameter K, the system at some point comes to equilibrium (in this model, the equilibrium state is 200 thousand units of goods). But much earlier, bilinear growth leads to a sharp increase in the quantity of goods, passing into infinity. If we leave both the right and left constantly positive flows bilinear, then already at about 20-30 steps, the value of the accumulator comes to the difference of two infinities.

Based on the above, it is safe to say that in the case of further use of such models, it is necessary to limit any positive growth.

A proto-operation model with two restrictions.

Having found out the shortcomings of the previous model and introducing a restriction on the second term by the saturation factor, we will build and run a new model.

Figure 2.13. Model of system dynamics and an example of its operation for system (2.3)

This model, in the end, brings the long-awaited results. It turned out to limit the growth of accumulator values. As can be seen from the right figure, for both enterprises, the equilibrium is reached with a slight excess of storage volume.

Extended proto-operation model.

When considering the system dynamics of this model, the capabilities of the AnyLogic software environment for colorful visualization of models will be demonstrated. All previous models were built using only elements of system dynamics. Therefore, the models themselves looked inconspicuous, they did not allow tracking the dynamics of changes in the amount of production over time and changing the parameters while the program was running. When working with this and the next models, we will try to use a wider range of program capabilities to change the three above disadvantages.

Firstly, in addition to the “system dynamics” section, the program also contains the “pictures”, “3D objects” sections, which make it possible to diversify the model, which is useful for its further presentation, since it makes the model look “more pleasant”.

Secondly, to track the dynamics of changes in the values ​​of the model, there is a "statistics" section that allows you to add charts and various data collection tools by linking them to the model.

Thirdly, for changing parameters and other objects during the execution of the model, there is a section "controls". The objects in this section allow you to change parameters while the model is running (for example, “slider”), select different states of the object (for example, “switch”) and perform other actions that change the initially specified data during work.

The model is suitable for teaching acquaintance with the dynamics of changes in the production of enterprises, but the lack of restrictions on growth does not allow using it in practice.

Extended proto-cooperation model with logistic constraint.

Using the already prepared previous model, we add to it the parameters from the logistic equation for limiting growth.

We omit the construction of the model, since the previous five models presented in the work have already demonstrated all the necessary tools and principles for working with them. It is only worth noting that its behavior is similar to the proto-cooperation model with the Verhulst constraint. Those. the lack of saturation hinders its practical application.

After analyzing the models in terms of proto-cooperation, we define several main points:

The models considered in this chapter are better suited in practice than mutualistic ones, since they have non-zero stable equilibrium positions even with two terms. Let me remind you that in the models of mutualism we were able to achieve this only by adding a third term.

Suitable models must have restrictions on each of the terms, because otherwise, a sharp increase in bilinear factors "destroys" the entire simulation model.

Proceeding from paragraph 2, when adding a proto-operation with the Verhulst restriction of the saturation factor to the extended model, as well as adding a lower critical amount of production, the model should become as close as possible to the real state of affairs. But do not forget that such manipulations of the system will complicate its analysis.

Conclusion

As a result of the study, an analysis was made of six systems that describe the dynamics of production by enterprises that mutually influence each other. As a result, the equilibrium points and types of their stability were determined in one of the following ways: analytically, or thanks to the constructed phase portraits in cases where an analytical solution is not possible for some reason. For each of the systems, phase diagrams were built, as well as three-dimensional models were built, on which, when projecting, it is possible to obtain integral curves in the planes (x, t), (y, t). After that, using the AnyLogic modeling environment, all models were built and their behavior options were considered under certain parameters.

After analyzing the systems and building their simulation models, it becomes obvious that these models can only be considered as training, or for describing macroscopic systems, but not as a decision support system for individual companies, because of their low accuracy and in in some places not quite a reliable representation of the ongoing processes. But also do not forget that no matter how true the dynamic system describing the model is, each company / organization / industry has its own processes and limitations, so it is not possible to create and describe a general model. In each specific case, it will be modified: to become more complicated or, on the contrary, to be simplified for further work.

Making a conclusion from the conclusions for each chapter, it is worth focusing on the revealed fact that the introduction of restrictions on each of the terms of the equation, although it complicates the system, but also allows you to detect stable positions of the system, as well as bring it closer to what is happening in reality. And it is worth noting that proto-cooperation models are more suitable for study, since they have non-zero stable positions, in contrast to the two mutualistic models we have considered.

Thus, the purpose of this study was achieved, and the tasks were completed. In the future, as a continuation of this work, an extended model of the interaction of the type of proto-operation with three restrictions introduced on it will be considered: logistic, saturation factor, lower critical number, which should allow creating a more accurate model for a decision support system, as well as a model with three companies. As an extension of the work, we can consider two other types of interaction besides symbiosis, which were mentioned in the work.

Literature

1. Bhatia Nam Parshad; Szegh Giorgio P. (2002). Stability theory of dynamical systems. Springer.

2. Blanchard P.; Devaney, R. L.; Hall, G. R. (2006). Differential Equations. London: Thompson. pp. 96-111.

Boeing, G. (2016). Visual Analysis of Nonlinear Dynamical Systems: Chaos, Fractals, Self-Similarity and the Limits of Prediction. systems. 4(4):37.

4. Campbell, David K. (2004). Nonlinear physics: Fresh breather. Nature. 432 (7016): 455-456.

Elton C.S. (1968) reprint. animal ecology. Great Britain: William Clowes and Sons Ltd.

7. Forrester Jay W. (1961). Industrial Dynamics. MIT Press.

8. Gandolfo, Giancarlo (1996). Economic Dynamics (Third ed.). Berlin: Springer. pp. 407-428.

9. Gershenfeld Neil A. (1999). The Nature of Mathematical Modeling. Cambridge, UK: Cambridge University Press.

10 Goodman M. (1989). Study Notes in System Dynamics. Pegasus.

Grebogi C, Ott E, and Yorke J. (1987). Chaos, Strange Attractors, and Fractal Basin Boundaries in Nonlinear Dynamics. Science 238 (4827), pp 632-638.

12 Hairer Ernst; Nørsett Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York

Hanski I. (1999) Metapopulation Ecology. Oxford University Press, Oxford, pp. 43-46.

Hughes-Hallett Deborah; McCallum, William G.; Gleason, Andrew M. (2013). Calculus: Single and Multivariable (6 ed.). John Wiley.

15. Llibre J., Valls C. (2007). Global analytic first integrals for the real planar Lotka-Volterra system, J. Math. Phys.

16. Jordan D.W.; Smith P. (2007). Non-Linear Ordinary Differential Equations: Introduction for Scientists and Engineers (4th ed.). Oxford University Press.

Khalil Hassan K. (2001). nonlinear systems. Prentice Hall.

Lamar University, Online Math Notes - Phase Plane, P. Dawkins.

Lamar University, Online Math Notes - Systems of Differential Equations, P. Dawkins.

Lang Serge (1972). Differential manifolds. Reading, Mass.-London-Don Mills, Ont.: Addison-Wesley Publishing Co., Inc.

Law Averill M. (2006). Simulation Modeling and Analysis with Expertfit Software. McGraw-Hill Science.

Lazard D. (2009). Thirty years of Polynomial System Solving, and now? Journal of Symbolic Computation. 44(3):222-231.

24 Lewis Mark D. (2000). The Promise of Dynamic Systems Approaches for an Integrated Account of Human Development. child development. 71(1): 36-43.

25. Malthus T.R. (1798). An Essay on the Principle of Population, in Oxford World's Classics reprint. p 61, end of Chapter VII

26. Morecroft John (2007). Strategic Modeling and Business Dynamics: A Feedback Systems Approach. John Wiley & Sons.

27. Nolte D.D. (2015), Introduction to Modern Dynamics: Chaos, Networks, Space and Time, Oxford University Press.

Automation and telemechanics, L-1, 2007

RAS B 02.70.-c, 47.ll.-j

© 2007 Yu.S. POPKOV, Dr. tech. Sci. (Institute for System Analysis RAS, Moscow)

QUALITATIVE ANALYSIS OF DYNAMIC SYSTEMS WITH Vd-ENTROPY OPERATOR

A method is proposed for studying the existence, uniqueness, and localization of singular points of the considered class of DSEE. Conditions for stability "in the small" and "in the large" are obtained. Examples of application of the obtained conditions are given.

1. Introduction

Many problems of mathematical modeling of dynamic processes can be solved on the basis of the concept of dynamic systems with an entropy operator (DEOS). DSEE is a dynamic system in which the nonlinearity is described by the parametric problem of entropy maximization. Feio-moyologically, DSEO is a model of a macrosystem with "slow" self-reproduction and "fast" resource allocation. Some properties of DSEO were studied in. This work continues the cycle of studies of the qualitative properties of DSEOs.

We consider a dynamical system with a Vd-entropy operator :

^ = £(x, y(x)), x e En:

y(x) = a^max(Hv(y) | Ty = u(x), y e E^) > 0.

In these expressions:

C(x, y), u(x) are continuously differentiable vector functions;

Entropy

(1.2) Hv (y) = uz 1n as > 0, s = T~m;

T - (r x w)-matrix with elements ^ 0 has a total rank equal to r;

The vector function u(x) is assumed to be continuously differentiable, the set

(1.3) Q = (q: 0<оТ ^ ц ^ а+} С Е+,

where a- and a + are vectors from E+, where a- is a vector with small components.

Using the well-known representation of the entropy operator in terms of Lagrange multipliers. we transform system (1.1) to the following form:

- = £(x, y(z)), x e Kn, y(z) e K?, r e Er+

Uz (r) \u003d az\\ ^, 3 \u003d 1, m-

O(x, z) = Ty(z) = q(x),

where rk = exp(-Ak) > 0 are the exponential Lagrange multipliers.

Along with the DSEE of the general form (1.1), we will consider, following the classification given in .

DSEE with separable flow:

(1-5) ^ = I (x) + Vy (z),

where B (n x m)-matrix;

DSEO with multiplicative flow:

(1.6) ^ = x ® (a - x ® Xu(r)), ab

where W is a (n x m)-matrix with non-negative elements, a is a vector with positive components, ® is the sign of coordinate-wise multiplication.

The aim of this paper is to study the existence, uniqueness, and localization of singular points of the DSEE and their stability.

2. Singular points

2.1. Existence

Consider system (1.4). The singular points of this dynamical system are determined by the following equations:

(2.1) C^(x, y(z))=0, r = TP;

(2.2) uz(r) = a^ r^, 3 = T^:

(2.3) bk(r) = ^as r^ = dk(x), k = 1,r.

Consider first the auxiliary system of equations:

(2.4) C(q, z) = r, q e R,

where the set R is defined by equality (1.3) and C(q, r) is a vector function with components

(2.5) Sk(d, r) = - Ok(r), a-< дк < а+, к =1,г.

Equation (2.4) has a unique solution r* for each fixed vector q, which follows from the properties of the Vg-entropy operator (see ).

From the definition of the components of the vector function С(g, z), the obvious estimate takes place:

(2.6) C(a+, r)< С(д, г) < С(а-,г), г в Е+. Рассмотрим два уравнения:

Let us denote the solution of the first equation by r+ and the second - by r-. Let's define

(2.7) C (a+,z) = z, C(a

(2.8) zmaX = max z+, zmin = mm zk

and r-dimensional vectors

(2.9) z(zmax, zmax), z(zmin , zmin).

Lemma 2.1. For all q G Q (1 . 3) solutions z*(q) of equation (2.4) belong to the vector 1 to the segment

zmin< z*(q) < zmax,

where the vectors zmin and zmax are defined by expressions (2.7)-(2.9).

The proof of the theorem is given in the Appendix. Qq

qk(x) (1.3) for x G Rn, then we have

Corollary 2.1. Let the conditions of Lemma 2.1 be satisfied and the functions qk(x) satisfy conditions (1.3) for all ex x G Rn. Then for all x G Rm the solutions z* of Eq. (2.3) belong to the vector segment

zmin< z* < zmax

Let us now return to equations (2.2). which determine the components of the vector function y(z). The elements of its Jacobian have the form

(2.10) jb aj zk JJ & > 0

for all z G R+ except for 0 and g. Therefore, the vector function y(z) is strictly monotonically increasing. According to Lemma 2.1, it is bounded from below and above, i.e., for all z G Rr (hence for all x G Rn) its values ​​belong to the set

(2.11) Y = (y: y-< y < y+},

where the components of the vectors yk, y+ are determined by the expressions:

(2.12) yk = aj y+ = aj znlax, j = h™.

(2.13) bj = Y, tsj, 3 =1,

Consider the first equation in (2.1) and rewrite it as:

(2.14) L(x, y) = 0 for all y e Y ⊂ E^.

This equation determines the dependence of the variable x on the variable y belonging to Y

we (1.4) reduces to the existence of an implicit function x(y) defined by equation (2.14).

Lemma 2.2. Let the following conditions be satisfied:

a) the vector function L(x, y) is continuous in the set of variables;

b) lim L(x, y) = ±<ж для любого фиксированного у е Y;

c) det J (x, y) = 0 for all ex x e En for any fixed y e Y.

Then there is a unique implicit function x*(y) defined on Y. In this lemma, J(x, y) is the Jacobian with elements

(2.15) Ji,i (x,y) = --i, i,l = l,n.

The proof is given in the Appendix. From the above lemmas it follows

Theorem 2.1. Let the conditions of Lemmas 2.1 and 2.2 be satisfied. Then there exists a unique singular point of the DSEE (1.4) and, accordingly, (1.1).

2.2. Localization

The study of the localization of a singular point is understood as the possibility of establishing the interval in which it is located. This task is not very simple, but for some class of DSEE such an interval can be established.

Let us turn to the first group of equations in (2.1) and represent them in the form

(2.16) L(x,y)=0, y- y y y+,

where y- and y+ are defined by equalities (2.12), (2.13).

Theorem 2.2. Let the vector function L(x,y) be continuously differentiable and monotonically increasing in both variables, i.e.

--> 0, --> 0; i,l = 1, n; j = 1,m. dxi dyj

Then the solution of system (2.16) with respect to the variable x belongs to the interval (2.17) xmin x x x xmax,

a) the vectors xmin, xmax have the form

Min \u003d i x 1 xmax \u003d r x t;

\xmin: . .., xminlxmax, . . ., xmax) :

xmin - ^Qin ^ ■ , xmax - ^QaX ^ ;

6) x- and x+ - components of the solution of the following equations

(2.19) L(x,y-)=0, L(x,y+) = 0

with oo m of course.

The proof of the theorem is given in the Appendix.

3. Sustainability of DSEA "in the small"

3.1. DSEE with a separable flow Let us turn to the DSEE equations with a separable flow, presenting them in the form:

- \u003d / (x) + Bu (r (x)), x e Kp ab

Y- (r (X)) \u003d azP (X) Y33, 3 \u003d 1, "~ 8 \u003d 1

0(x, r(x)) = Ty(r(x)) = q(x), r e Hr.

Here the values ​​of the components of the vector function q(x) belong to the set Q (1.3), the (n × w)-matrix B has a total rank equal to n (n< ш). Вектор-функция / (х) непрерывно дифференцируемая.

Let the system under consideration have a singular point x. To study the stability of this singular point "in the small" we construct a linearized system

where A is an (n x n)-matrix, the elements of which are calculated at the point x, and the vector t = x - x. According to the first equation in (3.1), the matrix of the linearized system has

A \u003d 7 (x) + BUg (g) Xx (x), x \u003d g (x),

| 3 \u003d 1, w, k \u003d 1,

I k \u003d 1, g, I \u003d 1, p

From (3.1) the elements of the matrix Yr: dy are determined.

"bkz P" 8=1

3, r8 x8, 5 1, r.

To determine the elements of the matrix Zx, we turn to the last group of equations in (3.1). B shows that these equations define an implicit vector function r(x), which is continuously differentiable if the vector function g(x) is continuously differentiable. The Jacobian Zx of the vector function z(x) is defined by the equation

<Эг (z)Zx(Х) = Qx(Х),

vg (X) \u003d T Ug (X),

ddk, -t-, - "- k \u003d 1, r, I \u003d 1, n dx \

From this equation we have (3.9) Zx(x) = s-1(z)Qx(x).

Substituting this result into equality (3.3). we get:

A \u003d 1 (x) + P (x), P (x) \u003d VUg (g) [Tg (g)] -1 Qx (x).

Thus, the equation of the linearized system takes the form

(c.i) | = (j+p)e

Here, the elements of the matrices J, P are calculated at a singular point. Sufficient stability conditions "in the small" DSEE (3.1) are determined by the following

Theorem 3.1. The DSEE (3.1) has a singular point x that is stable "in the small" if the following conditions are satisfied:

a) the matrices J, P (3.10) of the linearized system (3.11) have real and different eigenvalues, and the matrix J has the maximum eigenvalue

Pmax = max Pg > 0,

Wmax = maxUi< 0;

Umax + Ptah<

It follows from this theorem and equality (3.10) that for singular points for which Qx(x) = 0 and (or) for X, = 0 and tkj ^ 1 for all ex k,j, the sufficient conditions of the theorem are not satisfied.

3.2. DSEE with multiplicative flow Consider equations (1.6). presenting them in the form:

X ® (a - x ® Wy(z(x))), x e Rn;

yj(z(x)) = aj ПZs(x)]isi" j = 1,m;

(ZL2) yj (z(x)) = a^<~"ts

Q(x, z(x)) = Ty(z(x)) = q(x), z e R++.

systems. Will have:

(3.13)

In this expression, diag C] is a diagonal matrix with positive elements a1,..., an, Yr, Zx are matrices defined by equalities (3.4)-(3.7).

We represent the matrix A in the form

(3.14) A = diag + P (x),

(3.15) P(x) = -2xWYz(z)Zx(x).

Denote: maxi ai = nmax and wmax is the maximum eigenvalue of the matrix P(x) (3.15). Then Theorem 3.1 is also valid for the DSEE (1.6). (3.12).

4. Sustainability of DSEA "in the big"

Let us turn to the DESO equations (1.4), in which the values ​​of the components of the vector function q(x) belong to the set Q (1.3). In the system under consideration, there is a singular point Z, to which the vectors z(x) = z ^ z-> 0 and

y(x) = y(z) = y > y- > 0.

Let us introduce the deviation vectors £, C, П from the singular point: (4.1) £ = x - x, (= y - y, n = z - z.

ZHEZHERUN A.A., POKROVSKY A.V. - 2009

1

The aim of the study is to develop a supercomputer-oriented logical method (Boolean constraint method) and a service-oriented technology for creating and using a computer system for a qualitative study of the dynamics of the behavior of trajectories of autonomous binary dynamic systems over a finite time interval. The relevance of the topic is confirmed by the continuously increasing range of applications of binary models in scientific and applied research, as well as the need for a qualitative analysis of such models with a large state vector dimension. A mathematical model of an autonomous binary system on a finite time interval and a Boolean equation equivalent to this system are presented. The specification of a dynamic property is proposed to be written in the language of predicate logic using limited existential and universal quantifiers. Boolean equations for the search for equilibrium states and cycles of a binary system and conditions for their isolation are obtained. The main properties of the reachability type are specified (reachability, safety, simultaneous reachability, reachability under phase constraints, attraction, connectivity, total reachability). For each property, its model is built in the form of a Boolean constraint (a Boolean equation or a quantified Boolean formula) that satisfies the logical specification of the property and the equations of the system dynamics. Thus, checking the feasibility of various properties of the behavior of trajectories of autonomous binary dynamic systems over a finite time interval is reduced to the problem of the feasibility of Boolean constraints using modern SAT and TQBF solvers. A demonstration example of using this technology to test the feasibility of some of the previously given properties is given. In conclusion, the main advantages of the Boolean constraint method are listed, the features of its software implementation within the framework of a service-oriented approach, and the directions for further development of the method for other classes of binary dynamic systems are indicated.

binary dynamic system

dynamic property

qualitative analysis

boolean constraints

boolean satisfiability problem

1. Biere A., Ganesh V., Grohe M., Nordstrom J., Williams R. Theory and Practice of SAT Solving. Dagstuhl Reports. 2015.vol. 5. no. 4. R. 98–122.

2. Marin P., Pulina L., Giunchiglia E., Narizzano M., Tacchella A. Twelve Years of QBF Evaluations: QSAT Is PSPACE-Hard and It Shows. fundam. inform. 2016.vol. 149. R. 133–58.

3. Bohman D., Posthof H. Binary dynamical systems. M.: Energoatomizdat, 1986. 400 p.

4. Maslov S.Yu. The theory of deductive systems and its application. Moscow: Radio and communication, 1986. 133 p.

5. Jhala R., Majumdar R. Software model checking. ACM Computing Surveys. 2009.vol. 41 no. 4 R. 21:1–21:54.

6. Vasiliev S.N. Reduction method and qualitative analysis of dynamical systems. I–II // Izvestiya RAN. Theory and control systems. 2006. No. 1. S. 21–29. No. 2, pp. 5–17.

7. DIMACS format [Electronic resource]. Access mode: http://www.cs.utexas.edu/users/moore/acl2/manuals/current/manual/index-seo.php/SATLINK____DIMACS (accessed 24.07.2018).

8. QDIMACS standard [Electronic resource]. Access mode: http://qbflib.org/qdimacs.html (accessed 07/24/2018).

9. Delgado-Eckert E., Reger J., Schmidt K. Discrete Time Systems with Event-Based Dynamics: Recent Developments in Analysis and Synthesis Methods. Mario Alberto Jordan (Ed.). Discrete Time Systems. intech. 2011. R. 447–476.

10. Vasiliev S.N. Reachability and Connectivity in an Automaton Network with a General Switching Rule // Differential Equations. 2002. V. 38. No. 11. S. 1533–1539.

11. Bychkov I.V., Oparin G.A., Bogdanova V.G., Gorsky S.A., Pashinin A.A. Multi-agent technology for automating the parallel solution of Boolean equations in a distributed computing environment // Computational technologies. 2016. V. 21. No. 3. S. 5–17.

12. Lonsing F., Biere A. DepQBF. A Dependency-Aware QBF solver. Journal on Satisfiability. Boolean Modeling and Computing. 2010. vol. 9. R. 71–76.

13. Oparin G.A., Bogdanova V.G., Pashinin A.A., Gorsky S.A. Distributed Solvers of Applied Problems Based on Microservices and Agent Networks. Proc. Of the 41st Intern. Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO-2018). R. 1643–1648.

14. Bogdanova V.G., Gorsky S.A. Scalable Parallel Solver of Boolean Satisfiability Problems. Proc. Of the 41st Intern. Convention on Information and Communication Technology. Electronics and Microelectronics (MIPRO-2018). R. 244–249.

15. Bychkov I.V., Oparin G.A., Bogdanova V.G., Pashinin A.A. The Applied Problems Solving Technology Based on Distributed Computational Subject Domain Model: a Decentralized Approach // Parallel Computing Technologies XII International Conference, PaVT’2018, Rostov-on-Don, April 2–6, 2018. Short articles and poster descriptions. Chelyabinsk: Publishing Center of SUSU, 2018. P. 34–48.

The range of applications of binary dynamic models is extremely wide, and every year the number of objects and tasks where their use is required only increases. A classic example is a binary synchronous automaton, which is a model of many discrete devices in control systems, computer technology, telemechanics. Modern applications of binary dynamic models include the problems of bioinformatics, economics, sociology, and a number of other areas that seem far from the use of two-valued variables. In this regard, the relevance of developing new and improving existing methods for qualitative analysis of the behavior of the trajectories of binary dynamical systems (DDS) increases significantly.

As is known, the goal of a qualitative analysis of a dynamic system (not just a binary one) is to obtain a positive or negative answer to the question: Does the required dynamic property hold in a given system? Let's rephrase this question as follows: Does the behavior of the trajectories of a dynamical system satisfy a certain set of constraints that characterize the property? Further, we will use this interpretation of the goal of a qualitative analysis of the dynamic properties of the system.

For DDS, whose operation is considered on a finite time interval, such constraints are Boolean and are written in the language of Boolean equations or Boolean formulas with quantifiers. The first type of constraints leads to the need to solve the SAT problem (boolean satisfiability problem); the second type of constraints is associated with the solution of the TQBF problem (checking the truth of quantified Boolean formulas). The first problem is a typical representative of the NP complexity class, and the second problem is the PSPACE complexity class. As is known, PSPACE-completeness of a discrete problem gives stronger evidence of its intractability than NP-completeness. Because of this, the reduction of the problem of qualitative analysis of DDS to the SAT problem is more preferable than the reduction to the TQBF problem. In the general case, the study of not every property of the DDS can be represented in the language of Boolean equations.

The theoretical possibility of using Boolean constraints (namely, Boolean equations) in the qualitative analysis of DDS was first demonstrated in . However, it should be noted that the application of this approach in practice at that time was constrained by the lack of efficient algorithms and programs for solving Boolean equations (especially with a large number of unknown variables), which would significantly reduce the search space. In the last decade, as a result of intensive research in this area, a sufficient number of various efficient Boolean equation solvers (SAT solvers) have appeared that use modern achievements (new heuristics, fast data structures, parallel computing, etc.) in solving the Boolean satisfiability problem. Similar processes (but with some delay) are also observed in the field of creating more and more efficient algorithms and programs for solving the TQBF problem. Thus, to date, there are all the necessary prerequisites for the systematic development of the method of Boolean constraints in the qualitative analysis of DDS, its software implementation and application in solving scientific and applied problems.

In addition to the Boolean constraint method, other methods of qualitative analysis are also applicable to DDS, which include deductive analysis, model checking, and the reduction method. Each of these methods (including the boolean constraint method) has its limitations, advantages, and disadvantages. A common disadvantage is that all methods are enumerative in nature and the problem of reduction of enumeration is fundamental for these methods.

The importance of deductive analysis, which involves the application of axioms and inference rules to prove the correct functioning of a system, is recognized by a wide range of specialists, but this is a laborious and therefore rarely used method. In the model checking method, the required property specification language uses the language of temporal logics, which is unusual for specialists in automata dynamics. The reduction method is associated with the construction of a simplified (in a certain sense) model of the original system, the study of its properties and the conditions for transferring these properties to the original complex system. The conditions for the transferability of properties are only sufficient in this case. The simplicity of the idea of ​​the reduction method in the qualitative analysis of DDS faces the problem of choosing a simplified system that satisfies all the conditions of the method.

The practical use of the Boolean constraint method involves the algorithmization and automation of the following processes:

1) development of a logical language for the specification of dynamic properties focused on a specialist in system dynamics;

2) building a model of a dynamic property in the form of a Boolean constraint of one type or another that satisfies the logical specification of the property and the equations of dynamics of a binary system;

3) presentation of the resulting model in the international format DIMACS or QDIMACS;

4) selection (development) of an efficient parallel (distributed) solver of the problem of satisfiability of Boolean constraints (SAT or TQBF solver);

5) development of tools for creating software services;

6) development of services for qualitative research of various dynamic properties of DDS.

aim of the present study is the solution of only the first two problems in relation to the algorithmization of qualitative studies of autonomous (without control inputs) synchronous DDS. Such systems in English-language publications are called synchronous Boolean networks (Boolean network). Other aspects of applying the Boolean constraint method (including for DDS with control inputs) are the subject of the following publications.

Mathematical model of autonomous DDS

Let X = Bn (B = (0, 1) be the set of binary vectors of dimension n (the DDS state space). Let t∈T = (1,…,k) denote the discrete time (cycle number).

For each state x0∈X, called the initial state, we define the trajectory x(t, x0) as a finite sequence of states x0, x1,…, xk from the set X. Further, we will consider DDS in which each pair of adjacent states xt, x(t - 1) (t∈T) trajectories are related by the relation

xt = F(xt - 1). (one)

Here F:X>X is a logic algebra vector function, called the transition function. Thus, for any x0∈X, the system of Boolean equations (1) represents a model of the dynamics of the behavior of DDS trajectories in the state space X on a finite time interval T = (1, 2,…,k). Here and below, the value k in the definition of the set T is assumed to be a predetermined constant. This limitation is quite natural. The point is that in a qualitative analysis of the behavior of DDS trajectories, the question of what can be said about the feasibility of some dynamic property for a fixed, not too large k is of practical interest. The choice of the value of k in each specific case is based on a priori information about the duration of the processes in the simulated discrete system.

It is known that the system of Boolean equations (1) with the initial state x0∈X for T = (1, 2,…,k) is equivalent to one Boolean equation of the form

For k = 1 (only one-step transitions are considered), Eq. (2) takes the form

(3)

Solutions to this equation define a directed graph consisting of 2n vertices marked by one of the 2n states of the set X. The vertices x0 and x1 of the graph are connected by an arc directed from state x0 to state x1. Such a graph in the theory of binary automata is called a transition diagram. The representation of the DDS behavior in the form of a transition diagram is very clear both when constructing trajectories and studying their properties, but is practically realizable only for small dimensions n of the state vector x∈X.

Language means for specifying dynamic properties

It is most convenient to specify a dynamic property specification in the language of formal logic. Following the paper , we denote by X0∈X, X1∈X, X*∈X the sets of initial, admissible and target states.

The main syntactic elements of the dynamic property logical formula are: 1) subject variables (components of the vectors x0, x1,…, xk, time t); 2) limited quantifiers of existence and universality; 3) logical connectives v, &; final formulas. The final formula represents the assertion that some states of the set of trajectories x(t, x0) (x0∈X0) belong to the evaluation sets X* and X1.

It should be noted that the use of limited existential and universal quantifiers provides a way of writing a dynamic property that is familiar to a specialist in dynamics. In the process of constructing a Boolean model, the properties for system (1) are replaced by restricted quantifiers by ordinary ones according to the following definitions:

where A(y) is a predicate that limits the value of the variable y.

Due to the finiteness of the range of the variable t, the limited quantifiers of existence and universality with respect to this variable are replaced by equivalent formulas that do not contain quantifiers

In what follows, we will assume that the elements of the sets X0, X1, X* are determined, respectively, by the zeros of the following Boolean equations

or characteristic functions of these sets , .

Taking into account the restriction on the initial states G0(x) = 0, along with equations (2, 3), we will use the following Boolean equations to shorten the notation:

(4)

Preliminary qualitative analysis of autonomous DDS

At the stage of preliminary analysis, the branching of the state (the set of its immediate predecessors), the presence of equilibrium states and closed trajectories (cycles) can be revealed (if necessary).

The state x1 in (3) will be called the successor of the state x0, and x0 the predecessor of the state x1. In an autonomous DDS, each state has only one successor, and the number of predecessors of a given state can vary from zero to 2n - 1. All immediate predecessors x0 of a state s∈X are zeros of the Boolean equation

If equation (6) has no solutions, then there are no predecessors of the state s.

The equilibrium states (if they exist) are solutions to the Boolean equation

The trajectory x0, x1,…, xk is called a cycle of length k if the states x0, x1,…, xk-1 are pairwise different from each other and xk = x0. A cyclic sequence of length k (if it exists) is a solution to the Boolean equation

where = 0 ( ) - pairwise difference conditions for the set of states C of a cycle of length k. If none of the cycle states has predecessors that do not belong to the set C, then such a cycle is called isolated. Let the elements s of the set C be determined by the solution of the Boolean equation Gc(s) = 0. Then it is easy to show that the cycle isolation condition is equivalent to the absence of zeros in the following Boolean equation:

Solutions to equation (7) (if they exist) determine the cycle states that have predecessors that do not belong to the set C.

Since the equilibrium state is a cycle of length k = 1, its isolation condition is similar to the isolation condition with k ≥ 2, with the difference that Gc(s) has the form of a complete disjunction that determines this equilibrium state.

In what follows, non-isolated equilibrium states and cycles will be called attractors.

Specification of dynamic properties of a reachability type

The main property of DDS, the need to verify which most often arises in practice, is the reachability property traditionally studied in graph theory (in our case, such a graph is a transition diagram) and its various variations. The reachability is defined as the classical problem of analyzing the behavior of DDS trajectories.

The definition of this property is connected with the assignment of the previously introduced sets X0, X*, X1 (corresponding to these sets of Boolean equations). It is assumed that the sets X0, X*, X1 satisfy the constraint

Since the set T is finite, the reachability property and its variations will be further understood as the property of practical reachability (reachability in a finite number of cycles). The following reachability type properties are considered:

1. The main reachability property of a set X* from a set X0 is formulated as follows: any trajectory launched from the set of initial states X0 reaches the target set X*. Using the restricted existential and universal quantifiers, the formula for this property is:

2. The security property ensures that for any trajectory launched from X0 the set X* is unreachable:

3. Property of simultaneous reachability. In some cases, a more “strict requirement” may be set, which consists in that each trajectory reaches the target set in exactly k cycles (k∈T):

4. Reachability property under phase constraints:

This property guarantees that all trajectories emitted from the set X0, until they hit the target set X*, are in the set X1.

5. Property of attraction. Let X* be an attractor. Then the logical formula of the attraction property coincides with the formula of the main reachability property:

those. for each trajectory released from the set X0, there is a time t∈T, starting from which the trajectory does not go beyond the set X*. The set X0 in this case belongs to a part of the area of ​​attraction of the set X*(X0∈Xa, where Xа is the full area of ​​attraction (pool) of the attractor).

Note that all variables in the above formulas of properties are actually connected, since the trajectory x0, x1,…, xk is completely determined by the initial state. Since the quantifiers with respect to the variable t are replaced by operations of multiplace disjunction or conjunction of the corresponding predicates, in each of the formulas there remains a single bounded universal quantifier (), which allows us to write the conditions for the feasibility of these properties in the language of Boolean equations (in the form of a SAT problem).

We present two properties, the verification of which leads to the necessity of solving the TQBF problem.

6. Connectivity property of the target set:

those. there is an initial state x0∈X0 such that each target state x*⊆X* is reachable at some time t∈T, which means that there is a trajectory corresponding to this state, such that all target states x*∈X* belong to this trajectory.

7. Property of total reachability of a set X* from X0:

those. each target state is reachable from X0.

Checking the feasibility of dynamic properties

For properties (1-5), checking their feasibility is reduced to finding the zeros of the Boolean equation, the technology of formation of which is of a standardized nature and is considered in detail only for the main achievability property. Properties (6, 7) lead to the problem of checking the truth of a quantified Boolean formula.

1. The main property of reachability. Its logical formula is

Taking into account (4), we write formula (8) as

where is the characteristic function of the set of states of the trajectory released from the initial state x0∈X0. Let's get rid of the existential quantifier in (9). Then we will have

where is the characteristic function of the set X*. We replace the restricted universal quantifiers with ordinary quantifiers. As a result, we get

Formula (10) is true if and only if the subquantifier expression is identically true, i.e.

The identical truth of the implication means that the Boolean function is a logical consequence of the function , i.e. any trajectory with initial state x0∈X0 reaches the target set X*.

Satisfaction of identity (11) is equivalent to the absence of zeros in the Boolean equation

In deriving (12), we got rid of the implication and replaced ϕ*(x0, x1,..., xk) with . If equation (12) has at least one solution, then the reachability property does not hold. Such a solution represents (in a certain sense) a counterexample for the property being checked and can help the researcher to identify the cause of the error.

Further, for brevity, for each property (2-4) we write out only an equation of type (12), suggesting the reader to independently reproduce the necessary arguments close to those given for the main reachability property.

2. Safety property

3. Property of simultaneous reachability

4. Reachability property under phase constraints

5. Property of attraction. The feasibility of this property is checked in two stages. At the first stage, it is found out whether the set X* is an attractor. If the answer is yes, then the main reachability property is checked at the second stage. If X* is reachable from X0, then all conditions of the attraction property are satisfied.

6. Connectivity property

7. Property of total reachability`

For properties (6, 7), the scalar form of the equality of two Boolean vectors xt = x* has the form

Let us demonstrate the above technology for qualitative analysis of autonomous DDS using the Boolean constraint method when checking the feasibility of some of the above properties for model 3.2 from work:

Denote by x0∈X = B3 the initial state of model (13). Let T = (1, 2). Let us write out the functions of one-step and two-step transitions of model (13) required for the specification of properties:

(14)

where the sign "." denotes the operation of conjunction.

To check the satisfiability of each property, the initial (X0) and target (X*) sets are specified, which are determined by the zeros of the equations G0(x) = 0, G*(x) = 0 or by the characteristic functions of these sets (see Section 2). As a SAT solver, the REBUS instrumental complex (IC) solver is used, and the TQBF solver is DepQBF . The coding of variables in Boolean models of the properties considered below for these solvers is given in Table. 1, Boolean models of these properties in DIMACS and QDIMACS formats are located in Table. 2.

Table 1

Variable encoding

Variable number in Boolean model

Property 1

Property 2

Property 3

Property 4

Property 5

table 2

Boolean property models

Property 1	Property 2	Property 3	Property 4 (A)	Property 4 (B)	Property 5
					e 1 2 3 4 5 6 7 8 9 0 4 -5 -6 7 -8 -9 -10 11 12 0 4 5 6 -7 8 9 10 -11 -12 0

1. The main reachability property (k = 2). Let X0 = (x∈X: x1 = 0), X*=(x∈X: x1 = 1). The initial and target sets are defined respectively by the equations G0(x) = x1 = 0 and . Boolean equation (12) in this case takes the form

where the function ϕ(x0, x1, x2) is defined in (14). The IR REBUS solver gives the answer "unsat" (the equation has no zeros), thus the reachability property X* from X0 is satisfied, which is clearly seen from the next transition diagram shown in the figure.

2. Cycles of length k = 2. A cyclic sequence of length 2 (if it exists) is a solution to the Boolean equation

The function looks like

The expression R(x0, x1) was not included in the equation when the cycle was found, since there are no cycles of length k = 1 (equilibrium states) in model (13). Using the IR REBUS solver, two answers were obtained (in the DIMACS output format): 1 2 3 4 5 -6 0 and 1 2 -3 4 5 6 0, corresponding to cyclic sequences (figure): ((1 1 1), (1 1 0)) and ((1 1 0), (1 1 1)). The sets of states of both cycles coincide, which means that model (13) has one cycle of length k = 2.

System Transition Diagram (13)

3. The property of cycle isolation. If the elements s of the set of states C of a cycle of length k = 2 are determined by the solution of the Boolean equation Gc(s) = 0, then the cycle isolation condition is equivalent to the absence of zeros in the following Boolean equation:

Since C = ((1 1 1), (1 1 0)), we have

For this equation, the IR REBUS solver finds two solutions: -1 2 3 4 5 -6 0 and -1 2 -3 4 5 -6 0 (in binary representation, according to the coding of variables in Table 1, these are pairs of states (0 1 1), (1 1 0) and ((0 1 0), (1 1 0)) Thus, the cycle state (1 1 0) has two predecessors, (0 1 1) and (0 1 0), which do not belong to the state set cycle This means that the isolation property of the cycle is not satisfied, i.e. this cycle is an attractor.

4. Property of attraction. Let X* = C be an attractor. The logical formula of the attraction property is the same as the formula of the main reachability property

and the corresponding Boolean equation for our case has the form

Let us write out the functions G0(x0), ϕ(x0, x1, x2) and . The function ϕ(x0, x1, x2) is given in (14). For X* = C, the expression is . Consider two options for setting the set of initial states X0, for the cases of fulfillment (A) and non-fulfilment (B) of the attraction property for k = 2 cycles.

A. Let . Then

In this case, for the Boolean equation (15), the answer is "unsat". The attraction property for a given set X0 is satisfied.

B. Let . Then

In this case, IR REBUS for equation (15) finds a solution: 1 -2 3 4 -5 -6 -7 8 9 0, which corresponds to the trajectory ((1 0 1),(1 0 0),(0 1 1)) . This trajectory with the initial state x0 = (1 0 1) does not reach the set X* = C in two cycles, which means that the attraction property cannot be satisfied for the given X0.

5. Connectivity property. The logical formula of the connectivity property has the form of the following statement:

For k = 2 ϕ*(x0, x1, x2) = G0(x0)∨ϕ(x0, x1, x2), where the function ϕ(x0, x1, x2) is given in (14). Let us choose the state (1 0 1) as the initial state. Then . Let the target set X* = ((0 1 1), (1 0 0)). In this case, the function G*(x*) has the form

Let's write G*(x*) in CNF format:

Using De-Morgan's law, we find the negation of the function ϕ*(x0, x1, x2). Substituting all the obtained functions into (16) and taking into account the encoding of Boolean variables (Table 1), we obtain a Boolean model in the QDIMACS format (Table 2). The DepQBF solver gives the answer "sat", which means the truth of the statement (16). The connectedness property for the given X0, X*, T = (1, 2) is satisfied.

Conclusion

The main advantages of the Boolean constraint method in the qualitative study of DDS include:

1. The logical language used by a specialist in automata dynamics to specify a dynamic property through the use of limited existence and universality quantifiers.

2. Based on the property formula and dynamic equations, the construction of the corresponding Boolean equation or a quantified Boolean formula is automatically performed.

3. It is quite simple to automate the process of converting the resulting Boolean expressions into conjunctive normal form with further generation of a file in the DIMAX and QDIMAX formats, which are input for SAT solvers and QBF solvers.

4. The problem of reducing the enumeration is to some extent solved by the developers of these solvers and is shielded from specialists in the qualitative analysis of DDS.

5. The possibility of solving the problem of qualitative analysis of DDS for large dimensions of the state vector n on a sufficiently long time interval T is provided. In terms of the number of states, the Boolean constraint method is quantitatively commensurate with the model checking method. Due to the fact that in recent years there has been a significant increase in the performance of specialized algorithms for solving SAT and TQBF problems, the total number of variables in the Boolean property model for modern solvers can be measured in thousands.

The software for the qualitative analysis of DDS based on the Boolean constraint method is implemented within the framework of a service-oriented approach using specialized Boolean equation solvers. The paper presents an example of the implementation of the Boolean constraint method based on a service-oriented approach for searching for cycles and equilibrium states in gene regulatory networks.

It should be noted that the Boolean constraint method is a fairly general method for qualitative analysis of DDS over a finite time interval. It is applicable not only to autonomous systems, but also to systems with control inputs, to systems with a memory depth greater than one, to general DDS, when the transition function is unsolvable with respect to the state xt and has the form F(xt, xt-1) = 0. For DDS with inputs, the controllability property and its various variations are of particular importance. In addition to the problems of DDS analysis, the Boolean constraint method is applicable to the problems of feedback synthesis (static or dynamic, by state or by input), which ensure the fulfillment of the required dynamic property in the synthesized system.

The study was supported by the Russian Foundation for Basic Research, project No. 18-07-00596/18.

Bibliographic link

Oparin G.A., Bogdanova V.G., Pashinin A.A. BOOLEAN CONSTRAINTS IN THE QUALITATIVE ANALYSIS OF BINARY DYNAMIC SYSTEMS // International Journal of Applied and Fundamental Research. - 2018. - No. 9. - P. 19-29;
URL: https://applied-research.ru/ru/article/view?id=12381 (date of access: 03/18/2020). We bring to your attention the journals published by the publishing house "Academy of Natural History"