Design and computer research of a nonlinear stochastic models describing the dynamics of interacting populations Anastasia V. Demidova1 , Olga V. Druzhinina2 , Olga N. Masina3 and Alexander V. Shcherbakov3 1 Peoples’ Friendship University of Russia (RUDN University), 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation 2 Federal Research Center «Computer Science and Control» of Russian Academy of Sciences, 44, building 2, Vavilov St., Moscow, 119333, Russian Federation 3 Bunin Yelets State University, 28, Kommunarov St., Yelets, 399770, Russian Federation Abstract The construction of nonlinear three-dimensional models of interconnected communities number dy- namics is considered, taking into account competition in populations of victims. A qualitative research of the systems is carried out, equilibrium states are found, the species number dynamics graphs are constructed. For these models, an estimate of the model parameters is given and local phase portraits are constructed. The transition to the corresponding stochastic models is made. In stochastic cases, the method of constructing self-consistent stochastic models is used. A comparative analysis of deterministic and stochastic models is carried out. Effects typical for three-dimensional models with regard to compe- tition in prey populations are revealed. A software package for the numerical solution of differential equations systems by modified Runge–Kutta methods is used as a software tool for researching the model. The software package allows performing numerical experiments based on the implementation of algorithms for generating trajectories of multidimensional Wiener processes and multipoint distributions and algorithms for solving stochastic differential equations. The formulation of the optimal control problem is proposed. Computer research of the models makes it possible to obtain the results of numerical experiments on the search for trajectories and the estimation of parameters. The results obtained can find application in problems of ecological systems computer modeling, as well as in problems of synthesis, optimal control and analysis of the multidimensional stochastic models stability describing the dynamics of interacting populations. Keywords computer modeling, nonlinear model of population dynamics, optimal control, stochastization of one-step processes, symbolic computation libraries, software package Woodstock’21: Symposium on the irreproducible science, June 07–11, 2021, Woodstock, NY Envelope-Open demidova-av@rudn.ru (A. V. Demidova); ovdruzh@mail.ru (O. V. Druzhinina); olga121@inbox.ru (O. N. Masina); shcherbakov_al.vl@mail.ru (A. V. Shcherbakov) Orcid 0000-0003-1000-9650 (A. V. Demidova); 0000-0002-9242-9730 (O. V. Druzhinina); 0000-0002-0934-7217 (O. N. Masina); 0000-0002-4752-8422 (A. V. Shcherbakov) © 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) 19 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 1. Introduction An important tool in solving problems of predicting the state of natural systems and managing them is mathematical modeling. To solve these problems, both traditional and new methods and approaches are used. Ecological systems with various interconnections between subsystems and a change in the structure of these interconnections in the process of functioning leads to the mathematical models construction, the analytical research of which is very difficult. Mathematical models of the interacting communities dynamics taking into account competi- tion and with food chains are considered in [1, 2, 3, 4, 5, 6, 7]. In [2], the stability conditions of the «predator–two preys» system are investigated. In [3], a mathematical model of a system with two competing prey and one predator is analyzed, the influence of predation on the species coexistence is described. A model of two competitors’ prey dynamics with the addition of a predator species to change the competition results is studied in [4]. In [5], the deterministic stability of the three-dimensional model «predator–two preys» limit cycles is investigated. In [6, 7], three-dimensional models of population dynamics with competition and with trophic chains are considered. In the process of stochastic modeling for various dynamical systems a method for constructing self-consistent one-step models [8] is proposed and a software package [9] is developed. Some systems of population dynamics based on the construction of stochastic self-consistent models are considered in [10, 11, 12, 13, 14]. In [12, 13, 14, 15, 16], a number of control problems for the models of population dynamics are considered. When modeling population systems, various software tools are used that provide ample opportunities for conducting computational experiments. In [17, 18], the models research is carried out using the Python language and symbolic computation libraries. In this paper we considered several types of three-dimensional models, taking into account the competition among prey species and predation are studied. The stability of these models is investigated, the equilibrium states are calculated. The transition from deterministic models to stochastic ones is performed. A computer research is carried out to study the stability. The estimation of the model parameters is carried out, the phase portraits of the system are constructed, as well as the graphs of the population size dynamics in the deterministic and stochastic cases. The research is carried out using a software package for constructing stochastic dynamic models and searching for appropriate trajectories, as well as the Jupyter application package. The obtained effects are analyzed. The formulations of optimal control problems for the models with trophic chains are proposed. 2. The deterministic models description We consider a model described by differential equations of the form 𝑥1̇ = 𝑥1 (𝑎1 − 𝜀𝑥1 − 𝛿𝑥2 − 𝑏1 𝑦), 𝑥2̇ = 𝑥2 (𝑎2 − 𝛿𝑥1 − 𝜀𝑥2 − 𝑏2 𝑦), (1) 𝑦 ̇ = 𝑦(−𝑐 + 𝑑1 𝑥1 + 𝑑2 𝑥2 − 𝛾 𝑦), 20 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 where 𝑥1 is the population density of the first competitor, 𝑥2 is the population density of the second competitor, 𝑦 is the population density of the predator, 𝑎𝑖 is the reproduction rate of the competitor’s population in the absence of a predator, 𝑏𝑖 is the specific rate of consumption of the prey population by the predator population, di is the conversion factor of the prey biomass consumed by the predator in own biomass, 𝑖 = 1, 2, 𝜀 is the coefficient of intraspecific competition, 𝛿 is the coefficient of interspecific competition, 𝑐 is the natural mortality of the predator, 𝛾 is the intraspecific competition of the predator, 𝑦(0) ≥ 0, 𝑥𝑖 (0) ≥ 0, 𝑖 = 1, 2. According to the ecological sense, the constraints on the coefficients are: 𝜀 > 0, 𝛿 ≥ 0, 𝛾 ≥ 0, 𝑎𝑖 > 0, 𝑏𝑖 > 0, 𝑑𝑖 > 0, 𝑐 > 0, 𝑖 = 1, 2. Model (1) is a modification of the model described in [1] with the following notation: 𝜀11 = 𝜀22 = 𝜀, 𝜀12 = 𝜀21 = 𝛿. We introduce the following notation: 𝜌 = 𝑏1 𝑑2 + 𝑏2 𝑑1 , 𝜙 = 𝑏1 𝑑1 + 𝑏2 𝑑2 , 𝐷 = 𝛾 (𝜀 2 − 𝛿 2 ) − 𝛿𝜌 + 𝜀𝜙, (𝑎2 𝑏1 − 𝑎1 𝑏2 )𝑑2 + (𝑏2 𝛿 − 𝑏1 𝜀)𝑐 + (𝑎2 𝛿 − 𝑎1 𝜀)𝛾 𝑥1∗ = , 𝛾 (𝛿 2 − 𝜀 2 ) + 𝛿𝜌 − 𝜀𝜙 (𝑎1 𝑏2 − 𝑎2 𝑏1 )𝑑1 + (𝑏1 𝛿 − 𝑏2 𝜀)𝑐 + (𝑎1 𝛿 − 𝑎2 𝜀)𝛾 𝑥2∗ = , 𝛾 (𝛿 2 − 𝜀 2 ) + 𝛿𝜌 − 𝜀𝜙 𝑐(𝜀 2 − 𝛿 2 ) + 𝛿(𝑎1 𝑑2 + 𝑎2 𝑑1 ) − 𝜀(𝑎1 𝑑1 + 𝑎2 𝑑2 ) 𝑦∗ = . 𝛾 (𝛿 2 − 𝜀 2 ) + 𝛿𝜌 − 𝜀𝜙 Equilibrium states of the (1) in general form are found. The indicated equilibrium states are as follows: −𝑐 𝑎 𝑎 𝐸0 (0, 0, 0), 𝐸1 (0, 0, ) , 𝐸2 (0, 2 , 0) , 𝐸3 ( 1 , 0, 0) , 𝛾 𝜀 𝜀 𝑎2 𝛾 + 𝑏2 𝑐 𝑎2 𝑑2 − 𝑐𝜀 𝑎1 𝛾 + 𝑏1 𝑐 𝑎 𝑑 − 𝑐𝜀 𝐸4 (0, , ), 𝐸5 ( , 0, 1 1 ), 𝛾 𝜀 + 𝑏2 𝑑2 𝛾 𝜀 + 𝑏2 𝑑2 𝛾 𝜀 + 𝑏1 𝑑1 𝛾 𝜀 + 𝑏1 𝑑1 𝑎 𝛿 − 𝑎1 𝜀 𝑎1 𝛿 − 𝑎2 𝜀 𝐸6 ( 2 2 , , 0) , 𝐸7 (𝑥1∗ , 𝑥2∗ , 𝑦 ∗ ) . 𝛿 − 𝜀2 𝛿 2 − 𝜀2 The state of equilibrium 𝐸7 is an internal state of equilibrium for which the condition of positivity is satisfied. Permanent coexistence of populations in the model (1) is established under the following conditions: 1) 𝜀 > 0, 𝛿 ≥ 0,𝛾 ≥ 0, 𝑎𝑖 > 0, 𝑏𝑖 > 0, 𝑑𝑖 > 0, 𝑐 > 0, 𝑖 = 1, 2; 2) there is a unique internal equilibrium state 𝐸7 such that 𝐷 ≠ 0 and 𝑥1∗ > 0, 𝑥2∗ > 0, 𝑦 ∗ > 0; 3) 𝐷 > 0; 4) at least one of the following inequalities holds 𝑎1 𝜀 > 𝑎2 𝛿 or 𝑎2 𝜀 > 𝑎1 𝛿. Next, we consider the model: 𝑥1̇ = 𝑥1 (𝑎 − 𝜀11 𝑥1 − 𝜀12 𝑥2 − 𝑏𝑦), 𝑥2̇ = 𝑥2 (𝑎 − 𝜀21 𝑥1 − 𝜀22 𝑥2 − 𝑏𝑦), (2) 𝑦 ̇ = 𝑦(𝑐 + 𝑑𝑥1 + 𝑑𝑥2 − 𝛾 𝑦), 21 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 where 𝜀𝑖𝑗 for 𝑖 = 𝑗 = 1 and 𝑖 = 𝑗 = 2 are the coefficients of intraspecific competition, 𝜀𝑖𝑗 for 𝑖 not equal to 𝑗 are the coefficients of interspecies competition. According to the ecological sense, the constraints on the coefficients are: 𝜀𝑖𝑗 ≥ 0, 𝑖 ≠ 𝑗, 𝜀𝑖𝑗 > 0, 𝑖 = 𝑗, 𝑖, 𝑗 = 1, 2, 𝛾 ≥ 0, 𝑎 > 0, 𝑏 > 0, 𝑑 > 0, 𝑐 > 0. The meaning of other parameters included in the system (2) is similar to the model (1). Model (2) is a modification of the model described in [1] with the following values: 𝑏1 = 𝑏2 = 𝑏, 𝑎1 = 𝑎2 = 𝑎, 𝑑1 = 𝑑2 = 𝑑. Next, we introduce the following notation: 𝜌 = (𝜀11 − 𝜀12 − 𝜀21 + 𝜀22 ), 𝜙 = (𝜀11 𝜀22 − 𝜀12 𝜀21 ), 𝐷 = 𝑏𝑑𝜌 + 𝛾 𝜙, (𝑎𝛾 + 𝑏𝑐)(𝜀22 − 𝜀12 ) (𝑎𝛾 + 𝑏𝑐)(𝜀11 − 𝜀21 ) 𝑎𝑑𝜌 − 𝑐𝜙 𝑥̂1∗ = , 𝑥̂2∗ = , 𝑦̂ ∗ = . 𝐷 𝐷 𝐷 Equilibrium states of the model (2) in general form are found. The indicated equilibrium states are as follows: −𝑐 𝑎 𝑎 𝐸0 (0, 0, 0), 𝐸1 (0, 0, ), 𝐸2 (0, , 0) , 𝐸3 ( , 0, 0) , 𝛾 𝜀22 𝜀11 𝑎𝛾 + 𝑏𝑐 𝑎𝑑 − 𝑐𝜀 𝑎(𝜀22 − 𝜀12) 𝑎(𝜀11 − 𝜀21) 𝐸4 (0, , ), 𝐸5 ( , , 0) , 𝛾 𝜀22 + 𝑏𝑑 𝛾 𝜀22 + 𝑏𝑑 𝜙 𝜙 𝑎𝛾 + 𝑏𝑐 𝑎𝑑 − 𝑐𝜀11 𝐸6 ( , 0, ,), 𝐸7 (𝑥̂1∗ , 𝑥̂2∗ , 𝑦̂ ∗ ) . 𝛾 𝜀11 + 𝑏𝑑 𝛾 𝜀11 + 𝑏𝑑 The state of equilibrium 𝐸7 is an internal state of equilibrium for which the condition of positivity is satisfied. Permanent coexistence of populations in the model (1) is established under the following conditions: 1) 𝜀12 ≥ 0, 𝜀21 ≥ 0, 𝜀11 > 0, 𝜀22 > 0, 𝛾 ≥ 0, 𝑎 > 0, 𝑏 > 0, 𝑑 > 0, 𝑐 > 0; 2) there is a unique internal equilibrium state 𝐸7 such that 𝐷 ≠ 0 and 𝑥̂1∗ > 0, 𝑥̂2∗ > 0, 𝑦̂ ∗ > 0; 3) 𝐷 > 0; 4) at least one of the following inequalities holds 𝜀22 > 𝜀12 or 𝜀11 > 𝜀21 . In [1], a theoretical research of the generalized model «predator–two preys» stability is carried out. However, the numerical research of this model in order to identify the conditions of oscillating modes causes a number of difficulties. We carry out a computer research of the models (1) and (2), which are special cases of the model [1]. Subsequently, the transition to the stochastic case is made based on the method of self-consistent stochastic models. 3. Transition to stochastic models We carry out the transition to stochastization of the models (1) and (2) using the method of constructing self-consistent stochastic models. This method is based on a combinatorial methodology. We write down the scheme of interaction between elements for the «predator–two preys» system in general form 22 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 𝑎𝑖 𝑋𝑖 −→ 2𝑋𝑖 , 𝑖 = 1, 2, 𝜀𝑖𝑗 𝑋𝑖 + 𝑋𝑗 −−→ 𝑋𝑗 , 𝑖, 𝑗 = 1, 2, 𝑑𝑖 𝑏𝑖 𝑋𝑖 + 𝑌 −→ 2𝑌 , 𝑋𝑖 + 𝑌 −→ 𝑌 , 𝑖 = 1, 2, (3) 𝛾 𝑌 +𝑌− → 𝑌, 𝑐 𝑌→ − 0. In the interaction scheme (3), the first line corresponds to the natural reproduction of prey in the absence of other factors, the second line symbolizes intraspecific (at 𝑖 = 𝑗) interspecific (at 𝑖 ≠ 𝑗) competition, and the third describes the predator-prey relationship. The fourth line is responsible for intraspecific competition among predators, and the fifth describes their natural mortality. The method of constructing self-consistent stochastic models assumes, in the course of math- ematical transformations, a transition from the interaction scheme to obtaining the coefficients of the Fokker–Planck equation. This transition is carried out using the upgraded software package described in [9]. This software package is implemented in the Python programming language using the NumPy and SyPy libraries. The software package consists of the following modules: IS_to_SDE.py and stochastic.py. The IS_to_SDE.py module is designed to obtain the coefficients of the Fokker–Planck equation from the interaction scheme. The stochastic.py module is a module for obtaining solutions for the stochastic model. The algorithm of the software package is shown in Diagram 1. The IS_to_SDE.py module takes as input the matrices M and N of the system states before and after interaction, vectors K_plus and K_minus interaction coefficients) and a vector X (the system state vector). As a result, we obtain a symbolic representation of the Fokker–Planck equation coefficients. Using the SymPy library allows you to get the code of these coefficients in TeX, which makes them easier to read. The IS_to_SDE.py module consists of several functions. Hereafter there is a description of the main ones. The S_plus function for obtaining the forward interaction coefficients has the following description: def S_plus(X, K_plus, M): """ interaction coefficient [s^{-}_{1}(x),...,s^{-}_{s}(x)]""" res = [] for i in range(len(K_plus)): Prod_s = [Prod_(x, int(n)) for (x, n) in zip(X, M[i, :])] res.append(K_plus[i]* sp.prod(Prod_s)) return res The derivation function drift_vector for obtaining the drift vector A in the Fokker–Planck equation is described as follows: 23 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 Figure 1: Algorithm of the software package. def drift_vector(X, K_plus, K_minus, N, M): """ drift vector""" res = sp.zeros(rows=len(X), cols=1) R = M.T - N.T for i in range(len(K_plus)): res += R[:, i] * S(X, K_plus, K_minus, N, M)[i] return res In addition, we use the following description of the diffusion_matrix function to obtain the diffusion matrix B in the Fokker–Planck equation: def diffusion_matrix(X, K_plus, K_minus, N, M): """ diffusion matrix""" res = sp.zeros(rows=len(X), cols=len(X)) R = M.T - N.T R = sp.Matrix(R) for i in range(len(K_plus)): 24 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 res += R[:, i] * R[:, i].T * S(X, K_plus, K_minus, N, M)[i] return res In fig. 2. the result of the functions for obtaining the Fokker–Planck equation coefficients for the model (1) is presented. Figure 2: The Fokker-Planck equation coefficients for the model (1). Figure 3 shows the derivation and the result of the functions for obtaining the Fokker–Planck equation coefficients for the model (2). Figure 3: The Fokker-Planck equation coefficients for the model (2). Further, the obtained coefficients are transferred to the stochastic.py module for the numerical solution of the generated stochastic differential equations and drawing the graphs of these solutions. 25 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 The stochastic.py module can be used to study and numerically solve systems of ordinary differential equations and their corresponding stochastic differential equations based on the Runge–Kutta method and its modifications. A detailed description of this module is performed in [8, 9]. 4. Results of computer experiments For the models (1) and (2), we carry out a series of computational experiments using the above-described software package designed to study and numerically solve systems of ordinary differential equations and the corresponding stochastic differential equations. Computational experiments are carried out for both the deterministic case and the stochastic case. For the model (1), calculations are carried out at 𝜀 = 𝛿. We consider the following sets of pa- rameters: (𝑥1 , 𝑥2 , 𝑦) = (2, 1.6, 1.2), (𝑎1 , 𝑎2 , 𝑐, 𝜀, 𝑏1 , 𝑏2 , 𝑑1 , 𝑑2 , 𝛾 ) = (0.2, 0.4, 0.8, 0, 0.2, 0.4, 0.3, 0.6, 1). With this set of parameters, the approximate equilibrium states are obtained. Further, for the model (1), we consider the following sets of parameters: (𝑥1 , 𝑥2 , 𝑦) = (0.5, 0.4, 0.3), (𝑎1 , 𝑎2 , 𝑐, 𝜀, 𝑏1 , 𝑏2 , 𝑑1 , 𝑑2 , 𝛾 ) = (0.2, 0.4, 0.8, 0.2, 0.2, 0.4, 1.3, 2.4, 0.1). With this set of parameters, we obtained the approximate equilibrium states. Figures 4 and 5 show the dynamics of population density for given sets of initial conditions. The dashed line indicates the fluctuations of species number of the deterministic model (1), the solid line indicates the stochastic one. Taking into account the results presented in Fig. 4, we note that the trajectories have an oscillating character with conservation of amplitudes. For the corresponding stochastic model, damping of oscillations takes place with an approximation to the stationary mode. Taking into account the results presented in Fig. 5, we note the proximity of the trajectories of the deterministic and stochastic models. In the case under consideration, stochastization does not affect the of the system behavior, which is characterized by damping of oscillations. For the model (2), the following sets of parameters are considered: (𝑥1 , 𝑥2 , 𝑦) = (2, 1, 6), (𝑎, 𝑐, 𝜀11 , 𝜀12 , 𝜀21 , 𝜀22 , 𝑏, 𝑑, 𝛾 ) = (8, 2.5/3, 8, 4, 4.125, 1, 1, 1, 0). With this set of parameters, the approximate equilibrium states are obtained. Fig. 6 shows the population density dynamics for two sets of initial conditions indicated above. For stochastic and deterministic cases, the trajectories are located close to each other. As in the deterministic case, the mean values graphs of various realizations of the stochastic differential equation solutions reach a stationary mode. Computational experiments show that the character of the systems (1) and (2) stability is significantly influenced by the coefficients of intraspecific and interspecific competition. Oscillations are formed if 𝜀 = 𝛿 = 0 and if 𝑎1 = 𝑘𝑎2 ,𝑑1 = 𝑙𝑑2 at 𝑙 = 𝑘. At 𝜀 = 𝛿 ≠ 0, the oscillations have a damping character. At 𝑎1 = 𝑘𝑎2 , 𝑑1 = 𝑙𝑑2 at 𝑙 ≠ 𝑘 one of the prey populations dies out. Next, we proceed to the consideration of controlled models and formulate the optimal control problem. 26 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 Figure 4: The populations dynamics 𝑥1 , 𝑥2 , 𝑦 of the model (1) under initial conditions: (𝑥1 , 𝑥2 , 𝑦) = (2, 1.6, 1.2), (𝑎1 , 𝑎2 , 𝑐, 𝜀, 𝑏1 , 𝑏2 , 𝑑1 , 𝑑2 , 𝛾 ) = (0.2, 0.4, 0.8, 0, 0.2, 0.4, 0.3, 0.6, 1). Figure 5: The populations dynamics 𝑥1 , 𝑥2 , 𝑦 of the model (1) under initial conditions: (𝑥1 , 𝑥2 , 𝑦) = (0.5, 0.4, 0.3), (𝑎1 , 𝑎2 , 𝑐, 𝜀, 𝑏1 , 𝑏2 , 𝑑1 , 𝑑2 , 𝛾 ) = (0.2, 0.4, 0.8, 0.2, 0.2, 0.4, 1.3, 2.4, 0.1). 5. The problem of optimal control Let us formulate an optimal control problem for a three-dimensional model of the interconnected communities number dynamics, taking into account competition in populations of preys. For a three-dimensional model (1), the controlled model is given by a system of differential equations 27 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 Figure 6: The populations dynamics 𝑥1 , 𝑥2 , 𝑦 of the model (2) under initial conditions: (𝑥1 , 𝑥2 , 𝑦) = (2, 1, 6), (𝑎, 𝑐, 𝜀11 , 𝜀12 , 𝜀21 , 𝜀22 , 𝑏, 𝑑, 𝛾 ) = (8, 2.5/3, 8, 4, 4.125, 1, 1, 1, 0). 𝑥1̇ = 𝑥1 (𝑎1 − 𝜀𝑥1 − 𝜀𝑥2 − 𝑏1 𝑦) − 𝑢1 𝑥1 , 𝑥2̇ = 𝑥2 (𝑎2 − 𝜀𝑥1 − 𝜀𝑥2 − 𝑏2 𝑦) − 𝑢2 𝑥2 , (4) 𝑦 ̇ = 𝑦(−𝑐 + 𝑑1 𝑥1 + 𝑑2 𝑥2 − 𝛾 𝑦) − 𝑢3 𝑦, where 𝑢𝑖 = 𝑢𝑖 (𝑡) are control functions. The input parameters is explained in section 2. Constraints for model (4) are specified in the form 𝑥1 (0) = 𝑥10 , 𝑥2 (0) = 𝑥20 , 𝑥3 (0) = 𝑥30 , 𝑥1 (𝑇 ) = 𝑥11 , 𝑥2 (𝑇 ) = 𝑥21 , 𝑥3 (𝑇 ) = 𝑥31 , 𝑡 ∈ [0, 𝑇 ], (5) 0 ≤ 𝑢1 ≤ 𝑢11 , 0 ≤ 𝑢2 ≤ 𝑢21 , 0 ≤ 𝑢3 ≤ 𝑢31 , 𝑡 ∈ [0, 𝑇 ]. (6) With regard to problem (4)–(6), we consider the functional to be minimized in the form 𝑇 3 𝐽 (𝑢) = ∫ ∑ 𝑘𝑖 𝑢𝑖 (𝑡)𝑑𝑡. (7) 0 𝑖=1 Control quality criterion (4) corresponds to minimizing losses from regulation of the popula- tion, and in this case, the positive coefficients are denoted by 𝑘𝑖 in (7). For the model (4), the optimal control problem can be formulated as follows: find the minimum of functional (7) under conditions (5), (6) taking into account 𝑥𝑖 ≥ 0. We also generalize model (2) for the controlled case and formulate the corresponding optimal control problem. At the same time, we consider the criterion of control quality, which also consists in minimizing losses from regulation of the population. 28 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 It is possible to construct the control laws 𝑢1 , 𝑢2 , 𝑢3 by different methods. For example, it is possible to use PID controllers [19] or controllers using sliding mode [20]. For a population model with competition and migration flows in [13], the authors use a polynomial control of the form 𝑢𝑖 (𝑡) = ‖𝑅𝑇‖ , 𝑅 = (𝑟𝑖1 , 𝑟𝑖2 , ..., 𝑟𝑖𝑛 )𝑇 , 𝑇 = (𝑡 0 , 𝑡 1 , ..., 𝑡 𝑚 ), 𝑖 = 1, 2, 3. In this case, the model parameters are the coefficients 𝑟𝑖1 , 𝑟𝑖2 , ..., 𝑟𝑖𝑛 of polynomial functions. Methods of global parametric optimization [21, 22] are usually used to calculate the parameters. We propose to use control methods based on machine learning and regulators using artificial intelligence. In particular, it is possible to use machine learning in combination with controllers based on fuzzy logic or artificial neural networks [23, 24]. The generalized algorithm for constructing the optimal trajectory of the model (4) based on machine learning is shown in Fig. 7. Thus, the optimal control problem is to find 𝑢𝑖 = 𝑢𝑖 (𝑡) those that satisfy, firstly, the phase constraints of problem (5), and secondly, the optimality criterion (7). To solve the problem, it is necessary: a) to construct the loss function, b) to build a parametric control model, c) to use the global parametric optimization algorithm for search the coefficients of the control model with the minimum loss function. We propose the construction of stochastic controlled models taking into account the «predator–two preys» interaction. To study such models, it is advisable to consider the control laws 𝑢1 , 𝑢2 , 𝑢3 using the algorithm for constructing the optimal trajectory of the model (4). For population models with model migration, a number of computer experiments are carried out in [12, 14]. In [14], the «predator–prey» model is studied taking into account migration flows. In [12], multidimensional models of competing populations with migration, without trophic chains, are studied. In [13, 14], the indicated approach to the construction of stochastic controlled models is effective. A similar approach can be used for the models (1) and (2). 6. Conclusion Computer research of two competing individuals interaction of prey with a predator population nonlinear models makes it possible to study the proposed models stability. The obtained results of the models research at different variables sets and initial conditions make it possible to estimate the influence of the predator species on the result of the prey competition. It is established that the presence of trophic chains has a positive character on the result of competition, and this, in turn, contributes to the coexistence of species. By the aid of developed software package, graphs of population dynamics are constructed. For the systems (1) and (2), the transition to the corresponding stochastic differential equations is carried out. The introduction of stochastics makes it possible to take into account the probabilistic nature of the reproduction processes and species death, as well as random fluctuations occurring in the environment in time and leading to random fluctuations of the model parameters. We formulated an optimal control problem, 29 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 Figure 7: Generalized algorithm for solving the optimal control problem. proposed a criterion for the quality of control and developed the corresponding generalized algorithm. To solve the problem of optimal control, it is proposed to use numerical optimization methods and intelligent algorithms for symbolic computations. The obtained results can find application in the research of the proposed models, taking into account the requirements of control and optimization. 30 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 References [1] V. Hutson, G. T. Vickers, A criterion for permanent coexistence of species, with an application to a two-prey one-predator system, Mathematical Biosciences 63 (1983) 253–269. doi:10.1016/0025- 5564(82)90042- 6 . [2] K. Fujii, Complexity-stability relationship of two-prey-one-predator species system model: Local and global stability, Journal of Theoretical Biology 69 (1977) 613–623. doi:10.1016/ 0022- 5193(77)90370- 8 . [3] S. B. Hsu, Predator-mediated coexistence and extinction, Mathematical Biosciences 54 (1981) 231–248. doi:10.1016/0025- 5564(81)90088- 2 . [4] M. van Baalen, V. Křivan, P. C. J. van Rijn, M. W. Sabelis, Alternative food, switching predators, and the persistence of predator‐prey systems, The American Naturalist 157 (2001) 512–524. doi:10.1086/319933 . [5] L. V. Bashkirceva, I. A. andKarpenko, L. B. Ryashko, Stochastic sensitivity of limit cycles for «predator–two preys» model, Izvestiya VUZ. Applied Nonlinear Dynamics 18 (2010) 42–64. doi:10.18500/0869- 6632- 2010- 18- 6- 42- 64 . [6] A. V. Epifanov, V. G. Tsibulin, On the dynamics of symmetric systems, Computer research and modeling 9 (2017) 799–813. [7] E. A. Aponina, Y. M. Aponin, A. D. Bazykin, Analysis of complex dynamic behavior in the predator–two prey model, Problems of ecological monitoring and modeling of ecosystems 5 (1982) 163–180. [8] M. N. Gevorkyan, T. R. Velieva, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, Stochas- tic runge–kutta software package for stochastic differential equations, in: W. Zamojski, J. Mazurkiewicz, J. Sugier, T. Walkowiak, J. Kacprzyk (Eds.), Dependability Engineer- ing and Complex Systems, Springer International Publishing, Cham, 2016, pp. 169–179. doi:10.1007/978- 3- 319- 39639- 2_15 . [9] M. Gevorkyan, A. Demidova, T. Velieva, A. Korolkova, D. Kulyabov, L. Sevastyanov, Imple- menting a method for stochastization of one-step processes in a computer algebra system, Programming and Computer Software 44 (2018) 86–93. doi:10.1134/S0361768818020044 . [10] A. V. Demidova, Equations of population dynamics in the form of stochastic differential equations, RUDN Journal of Mathematics, Information Sciences and Physics (2013) 67–76. [11] O. V. Druzhinina, O. N. Masina, Methods of stability research and controllability of fuzzy and stochastic dynamical systems, Dorodnitsyn Computing Center of RAS, Moscow, 2009. [12] A. V. Demidova, O. V. Druzhinina, O. N. Masina, E. D. Tarova, Computer research of nonlinear stochastic models with migration flows, in: Proceedings of the Selected Papers of the 9th International Conference “Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems” (ITTMM-2019). CEUR Workshop Proceedings, volume 2407, 2019, pp. 26–37. [13] A. V. Demidova, O. V. Druzhinina, O. N. Masina, A. A. Petrov, Computer research of the controlled models with migration flows, in: Proceedings of the Selected Papers of the 10th International Conference “Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems” (ITTMM-2020). CEUR Workshop Proceedings, volume 2639, 2020, pp. 117–129. [14] A. V. Demidova, O. V. Druzhinina, M. Jacimovic, O. N. Masina, N. Mijajlovic, N. Olenev, A. A. 31 Anastasia V. Demidova et al. CEUR Workshop Proceedings 19–32 Petrov, The generalized algorithms of global parametric optimization and stochastization for dynamical models of interconnected populations, in: N. Olenev, Y. Evtushenko, M. Khachay, V. Malkova (Eds.), Optimization and Applications, volume 12422, Springer International Publishing, Cham, 2020, pp. 40–54. doi:10.1007/978- 3- 030- 62867- 3_4 . [15] O. Kuzenkov, G. Kuzenkova, Optimal control of auto-reproduction systems, Bulletin of the Russian Academy of Sciences, Theory and control systems (2012) 26–37. [16] I. Danilova, The problem of optimal behavior of the population in the area taking into account intraspecific competition and migration, in: Proceedings of the XIII All-Russian Meeting on Control Problems of the VSPU-2019, ICS of RAS, Moscow, 2019, pp. 1476–1479. [17] C. Fuhrer, J. Solem, O. Verdier, Scientific Computing with Python 3, Packt Publishing, 2016. [18] R. Lamy, Instant SymPy Starter, Packt Publishing, 2013. [19] E. Nikulin, Fundamentals of the theory of automatic control. frequency methods of analysis and synthesis of systems, in: Textbook for universities, BHV–Petersburg, Saint-Petersburg, 2004. [20] S. Ulyanov, N. Nefedov, A fuzzy controller with a sliding mode based on soft calculations, Software products and systems, Software & Systems 4 (2014) 174–177. doi:10.15827/ 0236- 235X.108 . [21] R. Storn, K. Price, Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces, Journal of Global Optimization 11 (1997) 341–359. doi:10.1023/A:1008202821328 . [22] J. Lampinen, A constraint handling approach for the differential evolution algorithm, in: Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No.02TH8600), volume 2, 2002, pp. 1468–1473. doi:10.1109/CEC.2002.1004459 . [23] A. Uskov, A. Kuzmin, Intelligent control technologies, Artificial neural networks and fuzzy logic, Hot line–Telecom, Moscow, 2004. [24] D. Rutkovskaya, M. Pilinsky, L. Rutkovsky, Neural networks, genetic algorithms and fuzzy systems, Hot line–Telecom, Moscow, 2006. 32