Controlled Positive Dynamic Systems with an Entropy Operator: Fundamentals of the Theory and Applications

: Controlled dynamic systems with an entropy operator (DSEO) are considered. Mathematical models of such systems were used to study the dynamic properties in demo-economic systems, the spatiotemporal evolution of trafﬁc ﬂows, recurrent procedures for restoring images from projections, etc. Three problems of the study of DSEO are considered: the existence and uniqueness of singular points and the inﬂuence of control on them; stability in “large” of the singular points; and optimization of program control with linear feedback. The theorems of existence, uniqueness, and localization of singular points are proved using the properties of equations with monotone operators and the method of linear majorants of the entropy operator. The theorem on asymptotic stability of the DSEO in “large” is proven using differential inequalities. Methods for the synthesis of quasi-optimal program control and linear feedback control with integral quadratic quality functional, and ensuring the existence of a nonzero equilibrium, were developed. A recursive method for solving the integral equations of the DSEO using the multidimensional functional power series and the multidimensional Laplace transform was developed. The problem of managing regional foreign direct investment is considered, the distribution of ﬂows is modeled by the corresponding DSEO. It is shown that linear feedback control is a more effective tool than program control.


Introduction
Dynamic systems with an entropy operator (DSEOs) form a class of nonlinear systems in which the nonlinearity is described by a perturbed mathematical programming problem with an entropy objective function. DSEOs turned out to be a useful framework for studying a fairly wide range of applied problems: modeling of demo-economic systems [1], traffic flows [2,3], the labor market [4], and urban agglomerations [5]; the development of diagnostics algorithms for atmospheric waves [6] and oil-bearing strata [7]; analysis of stationary modes in dynamic image reconstruction procedures from projections in computerized tomography [8]; identification of dynamic systems [9], and others.
The mathematical class of DSEOs is based on a physical model of a dynamic system with self-reproduction of matter, energy, or information exchanged stochastically between the subsystems. Moreover, self-reproduction processes have an evolutionary nature ("slow"), whereas exchange processes occur much more intensively ("fast"). The latter processes relax in rather short periods (in the time scale of "slow" processes) to local equilibrium states. Since exchange processes are assumed to be stochastic, their locally stationary states can be described by the corresponding entropy operator.
Research into DSEOs has some history and problems examined in related areas. Apparently, the monograph [10], devoted to studying processes in urban and regional systems, was the first publication on models with an entropy operator.
A feature of DSEOs is a specific nonlinearity description in the form of the argmax function. In [11], forced periodic modes in a dynamic system with the argmin function without constraints were investigated. The paper [12] considered the convergence of iterative algorithms for finding the argmin of a continuous function.
The entropy operator in DSEOs is also described by the argmax function, but the constrained maximum of the entropy function. Constraints can be specified as a system of equalities or inequalities. In the latter case, the mathematical model of the entropy operator reduces to a perturbed mathematical programming problem, which significantly complicates the study of the operator itself and the corresponding dynamic system.
The general properties of an entropy operator with equality constraints (existence, continuity, boundedness, differentiability, the Lipschitz condition) were studied in [13].
This paper considers the dynamic properties of controlled positive DSEOs with the constrained argmax entropy function. Control in such systems affects the entropy operator. Two classes of controls are studied: program and linear feedback controls. We formulate a non-negativity theorem for the solutions of equations describing the DSEO under nonnegative control and initial conditions. The problem of equilibrium states is investigated: the existence of a unique singular point is established using the method of monotonic operators [14] and bilinear majorants [15], and the sizes of a multidimensional rectangular box (a vector interval) containing this point are estimated. The sizes of the asymptotic stability domain of the singular point are estimated by the state vector norm using the method of differential inequalities [16]. Optimal control design problems in the two classes mentioned are considered using an integral quadratic criterion and constraints associated with the existence, uniqueness, and stability of the equilibrium of positive DSEOs.
The developed methods are applied to study the qualitative properties and optimize a stochastic regional exchange system for foreign direct investment.

Mathematical Model of Positive Controlled DSEOs
Consider a DSEO with some control u(t) ∈ R r + applied to the entropy operator: where the controlled entropy operator is given by Note that in the absence control, the DSEO has zero singular point. It is stable if w < 0, and unstable if w > 0.
The other notations are as follows: ⊗ means the coordinate-wise product of vectors; w ∈ R n + is a vector with constant components; S > 0 is a non-negative matrix of dimensions (n × m); T is a non-negative matrix of dimensions (r × m) that has full rank and normalized columns: where e (r) and e (m) denote unit column vectors of dimensions r and m, respectively. The entropy has the form H(y) = − y, ln y e a , a ∈ A = [0, 1] for all x ∈ R n + , e = 2.718.
The entropy operator is a vector y[u(x)] with the components where the vector z consists of the exponential Lagrange multipliers for problem (2). (1), and let w 0, S ≥ 0, and y[u(x)] 0 for all u(x) 0.

Theorem 1. In Equation
Then for x(0) 0, all solutions x(t) of (1) are non-negative: x * (t) ∈ R n + , for all t ≥ 0, and the system (1) belongs to the class of positive DSEOs.
Proof. We introduce three subsets in the state space of the DSEO (1) as follows: where and s (i) is the ith row of the matrix S; 1. Consider an initial point x(0) ∈ R n + whose components with numbers i 1 , . . . , i k belong to the set X + , and those with numbers i k+1 , . . . , i n to the set X − . From (1) it follows that the right-hand sides of the equations with numbers i 1 , . . . , i k are positive. Hence, the components x i 1 (t), . . . , x i n (t) evolving from the point x(0) will increase monotonically, remaining positive. Now, we analyze the components x i k+1 , . . . , x i n . Since they belong to the set X − , the right-hand sides of the corresponding equations in (1) are negative. As a result, the components x i k+1 (t), . . . , x i n (t) evolving from the point x(0) will gradually vanish, remaining positive.
2. Consider the opposite case: in the initial point x(0) ∈ R n + , the components with numbers i 1 , . . . , i k belong to the set X − , and those with numbers i k+1 , . . . , i n to the set X + . This case is a "mirror image" of the previous one: the components x i 1 (t), . . . , x i n (t) will decrease, remaining non-negative, whereas the components x i k+1 , . . . , x i n (t) will increase, remaining positive.
3. The set X * is the set of nonzero singular points, and their existence is determined by equality (10). The proof of Theorem 1 is complete.

Singular Points and Their Localization
The system (1) has the trivial singular point x 0 = 0 and nonzero singular points satisfying the equations where the entropy operator y and the function Φ have the components (5). The relationship between the exponential Lagrange multipliers z and the state vector x, i.e., an implicit function z(x), is given by (5). We transform Equation (11) tõ are the components of the entropy operator y[z], and the functions Φ(x, z) are given by (5). Next, we multiply component-wise the left-and right-hand sides of Equation (12) by the vectors x and z, respectively. As a result, where (G, x,ỹ) ∈ R n + and Ψ,Φ, z ∈ R r + . We introduce the following notations: Then Equation (15) can be written as The operator A(v) maps the vector v into itself. Therefore, under several useful properties (contraction, monotonicity), the equation will have a unique solution, and we may estimate the domain containing it.

Theorem 2.
Let the control function have the following properties: for all x ∈ R n + . Then the operator A(v) is monotonically increasing for all v ∈ R n+r + .
Proof. Consider the Jacobian of (17). Due to (16), According to (16), the elements of the Jacobian J A are: By conditions (18) of this theorem, all component matrices of (19) have non-negative elements, and the conclusion follows. The proof of Theorem 2 is complete.
Thus, the operator A(v) is monotonically increasing on the set v 0. For monotonically increasing operators in Equation (17), there exists a pair of vectors v (l f ) < v (rh) such that: A(v (l f ) ) > v (l f ) and A(v (rh) ) < v (rh) for a concave operator, and A(v (l f ) ) < v (l f ) and A(v (rh) ) > v (rh) for a convex operator. Then the solution v * of Equation (17) belongs to the corresponding vector interval V = [v (l f ) , v (rh) ]; see Theorem 3.1 in [14].
The problem is estimating the boundaries of this vector interval. The null vector 0 can be adopted as the left boundary. For estimating the right boundary (the vector v (1) ), we define the majorant of the operator A(v) (15).
Recall that the majorant is a vector function C(v), such that where Lemma 1. Assume that condition (3) holds, and the control function u(x) is monotonically decreasing, i.e., its Jacobian satisfies the inequality Then the majorant C(v) (24) exists, and its components are given by where the vector x (2) and z / u(x) consist of x 2 i , i = 1, n, and z k / u k (x), k = 1, r, respectively.
Replacing the operator A(v) in (17) with its majorant (27) and (28), we obtain equations for the vector v (rh) = (x (rh) , z (rh) ): Hence, the right boundary of the domain locating the unique singular point of the system (1) satisfies the equation Here, the vector 1 x consists of 1 x i , i = 1, n, and the control function u(x) is monotonically decreasing on R r + . As an example, consider a feedback law with the exponential control function: We construct an approximate solution to localize the right boundary of the vector interval. The linear approximation of the exponential control function has the form Then Equation (31) reduces to where the vector p(x) consists of Consider a parameterized family of the following systems of equations: The case ε = 0 corresponds to the so-called generating system: Its solution is For ε = 1, we obtain the original system (36). Expanding its solution into a formal power series [17] in the parameter ε yields where x (k) is the kth approximation to the solution. Next, we substitute the series (39) into Equation (36) and equate the terms with the same powers of ε. As a result, the approximations will satisfy the following chain of relations: • the zeroth approximation, • the first approximation, • · · · · · · Therefore, the approximate right boundary of the vector interval locating the singular point is Well, the vector interval I = [0,x (rh) ] contains the singular point of the positive DSEO (1) and (2).

Stability of Nonzero Singular Point
Consider the deviation from the singular point x * , z * : The differential equation for this deviation has the form It taken into account here that Theorem 3. Let the following condition be satisfied: where L > 0 is the Lipschitz constant for the operator y[x].
Then the singular point x * is stable in large, i.e., stable under the initial deviations where σ = S .

Proof. Introducing the matrix
we consider the linear differential equation The corresponding matricant is defined as and its norm satisfies the upper bound where w min = min i w i . Using the matricant (51), we pass to the integral equation Denote the Euclidean norm of the vectors ξ(τ) in the form: Then in term norms and take into account (46), we have the following inequality: where Consider the equality: Differentiating equality (57) yields According to the theorem on differential inequalities [16], The solution of the first-order Equation (59) will asymptotically vanish if the initial conditions for v(t) (the norm of the deviation from the singular point x * ) belong to the domain, where the right side of the Equation (59) is negative.
Hence, the singular point x * will be asymptotically stable in the domain O (48).

Optimal Control for a Class of Positive DSEOs
It was previously seen that a control is necessary for the occurrence of nonzero singular point. However, in the process of reaching the point, the selected indicators of its quality can be performed. We will consider the integral indicators.
Consider the system (1) in which the entropy operator is described by the constrained argmax problem n ∑ s=1 y ks = u k (x), k = 1, n.
This problem has the analytical solution y * ks = a ks u k (x), y * sk = a sk u s (x), (k, s) = 1, n.
The entropy operator is given by the vector where A = [a sk | (s, k) = 1, n]. Then Equation (1) can be written as Here u ∈ R r + for all x ∈ R n + . Consider the case of program control: where c 0 > 0 is a vector of dimension r. Then Equation (64) reduces to Using the matrizer (51), we express (66) in the integral form: Equations (64) and (66) belong to the class of equations with the polynomial right-hand side. To find the solution, we employ the functional Volterra series [18][19][20]: where: -The vectors have dimension n.
Clearly, all elements of the weight function matrices in (69) are explicitly expressed through the matrix W(t) = diag[exp(−w i t) | i = 1, n] and the control characteristics (vector)c 0 .
Substituting the weight function matrices (73) into (69) gives where: All matrices in these equalities are square and have dimensions (n × n), and all vectors are n-dimensional; see (70).
Note that the vector x(t) linearly depends on the program controlc 0 , and the parameters of this dependence are positive.
Consider a quadratic optimization problem: Here the terminal time T is known.
Next, we substitute the expression for x(t) into (77), integrate, and perform some trivial transformations. As a result, we obtain: where the quadratic and linear form matrices are given by and respectively. In the equalities above, Clearly, the matrix K[x(0)] (78) has non-negative elements. The gradient of the function (77) is given by For calculating the minimum of (77), we apply the gradient projection recursive formula: The parameter γ is chosen to guarantee the convergence of this algorithm [22]. Now consider the case of linear feedback control: The equation of the controlled DSEO (68) reduces to Using the procedure described above, we get the following images of the weight function matrices of the series (69): The inverse Laplace transform finally yields Thus, the solution of (62) has the form For choosing the control matrix C, we adopt the criterion (77). In the case under consideration, it depends on the matrix C: where x(t) is given by (87). Therefore, where Within the constant term, the gradient of the criterion (90) with respect to the matrix C is calculated as The minimum of (90) can be found numerically by the gradient projection recursive formula: The parameter β is chosen to guarantee the convergence of this algorithm [22].

Control of Stochastic Flows of Regional Foreign Direct Investment (FDI)
Consider a regional system consisting of n regions exchanging investments [23]. Each region i has an investment appeal x i (t), i = 1, n, measured in the total cost of all projects implemented in it. Each region invests in its own projects and projects implemented in other regions (foreign direct investment, FDI).
The region's investment appeal changes due to two processes. On the one hand, it decreases over time since the projects proposed for investment are "aging." On the other, it grows as the result of FDI. Assume that the regions exchange their FDI stochastically and in portions.
The relative rate of change in the region's investment appeal, v i =ẋ i /x i , is proportional to the difference between the aging Φ i (t) and renewal Ψ i (t) functions of investment projects: v The aging function is linear one: where w i ; s i are positive constants. The renewal function is characterize by distribution of the FDI flows. As the FDI exchange is stochastic and partial, the locally stationary distribution of the FDI flows Y = [y ij | (i, j) = 1, n] is described by the controlled entropy operator where A = [a ij | (i, j) = 1, n] denotes the prior probability matrix, e = 2.71, and u = {u 1 , . . . , u n } is the control vector. The entropy-optimal distribution of the FDI flows has the form subject to the constraint n ∑ j=1 a ij = 1 for all i = 1, n.
The integral equation corresponding to (100) has the form where matricant

Equilibria in the System Stochastic FDI-Exchange
Consider the case of program control u = c 0 , where c 0 is a real positive vector. In FDI terms, this vector specifies a fixed distribution of FDI among the regions.
The equation of the system stochastic FDI-exchange with program control takes the form A unique nonzero singular point of (104) is calculated as Hence, the optimal program control satisfies the inequality Now consider the case of linear feedback control: In this case, the investment appeal equation reduces to Its unique nonzero singular point is given by where I is unique matrix. Let, in terms of the norm, we have According to (109) and (111) the control matrix take the form: The linear inequality (112) determine the set of admissible control matrices. Notice that an elements of control matrix C may be negative.

Optimization of Stochastic FDI Exchange
We apply the general optimization procedure of controlled DSEOs (see above) to the stochastic FDI exchange process. We optimize the control function on a time horizon T = [0, T] by the integral quadratic criterion J(u | y(0)) = T 0 y (t, u(t) | y(0)) y(t, u(t) | y(0)) dt. (113) 1. Program control u(t) = c 0 . In this case, the criterion takes the form and the integral equation describing the stochastic FDI exchange process (see (104)) reduces to and the matricant is given by We find the solution by expanding into a functional power series: In terms of the multidimensional Laplace transform, the images of the weight function matrices are written as Consider the first approximation to the solution, characterized by the matrix K (1) (s). Due to (119), this matrix is diagonal with the elements where a i denotes the ith row of the matrix A . The weight functions of the first approximation have the form In view of (119), the weight matrix of the second approximation is given by Since all weight function matrices are square and diagonal, we arrive at expressing the first approximation of J 1 depending on the variable c 0 (program control). An optimization problem for the program control c 0 can be stated as follows: 2. Linear feedback control u = C y. In this case, the criterion takes the form and the integral equation describing the stochastic FDI exchange process (see (104)) reduces to where the matrizer M t τ (W) is given by (117). The solution of this equation can be constructed in the form (118). Similar to the previous case, we apply the multidimensional Laplace transform and equate the terms with the same powers of y(0). As a result, we obtain the following chain of recursive equations for the images of the weight function matrices of the series (118): · · · · · · · · · Within the quadratic term, the solution of (127) has the image The approximate solution of (127) is calculated as Obviously, the control matrix C appears in the second term only (and even linearly). Therefore, within the constant terms, the criterion J(C | y(0)) (126) is a quadratic form of the control matrix: where:W The optimal linear feedback control problem can be written as This is a quadratic programming problem with linear constraints. It can be solved by standard algorithms using the matrix gradient of the criterion (132) [24].

Simulation of FDI-Exchange
Consider the FDI exchange process between three countries: China (1), France (2), USA (3). Assume that the exchange is adequately described by the stochastic model from Section 6, and the linear part of the aging function is absent. The Table 1 shows the data on the investment attractiveness of these countries in terms of shares of the total number of investment projects.
The notations are as follows: c 1 , c 2 , and c 3 are the components of the program control vector; A = [a ij | (i, j) = 1, 3] is the prior probability matrix; the vector w > 0 characterizes the aging function of investment projects. We choose the following numerical values of the parameters: The row sums of the prior probability matrix are equal to 1.
Without program control, the regional FDI flows quickly fades away, approaching the zero singular point. At the zero singular point, the investment attractiveness of the region is zero.
Therefore, when optimizing program control, it is necessary to take into account the conditions for the existence of a nonzero singular point: The optimal program control problem is written as where T = 4 and The values of the initial conditions y i (0), (i = 1, 2, 3) are taken from Table 1, 2016. The solution of the problem (139) takes the form: As we can see, the solution to this problem lies in the feasible region. Optimal programmed control provides the following nonzero points of the regional equilibrium attractions: y * = {0.07; 0.50; 0.14}.
(142) Table 2 shows the calculated data on the investment attractiveness of these countries in the interval 2016-2019 (initial conditions from 2016) for optimal program (141). The last row of the table contains the values of the relative mean square error between the data in Tables 1 (z) and 2 (y): So, we can see that the use of program control leads to rather significant errors when compared with real data.
2. Linear feedback control u = C y. In this case, the system equations take the form Here c ik , (i, k) = 1, 3 are elements of the control matrix. The numerical values of the parameters are given by (136).
The non-negativity condition for a singular nonzero point has the form: The solution to the Equation (144) in terms of (130) is: Diagonal matrices are (3 × 3). The optimal linear feedback control problem is written as J(C, y(0)) = 3 t=1 y (t,C) y(t,C) dt ⇒ min, To solve this problem, a random search algorithm was used with inversion at each step of the matrix and checking the boundary of the feasible region. The calculation results are presented in Table 3.  Tables 2 and 3, it can be seen that the mathematical model of a positive DSEO with feedback control more adequately describes the dynamics of investment attractiveness in the three-regional system.

Conclusions
A positive dynamical system with an entropy operator was considered. A theorem was proved on the stay of the trajectories of the system in the non-negative orthant of the state space near the initial perturbations from this subspace. Conditions for the existence, uniqueness ,and localization of a nonzero singular point were obtained under program control and feedback control. Applying the linear majorant of the entropy operator, the boundaries of the localization region were obtained.
The problem of stability of singular points was investigated. It was shown that, in the absence of control, the zero singular point was unstable, and not zero did not exist. Conditions for the existence and stability of a nonzero singular point were obtained under program control and feedback control.
A method for the synthesis of quasi-optimal control using the integral quadratic functional of quality and conditions for the existence of a singular non-zero point was developed.
The methods used in the article could be applied to study the dynamic properties and optimize the distribution system of regional flows of foreign direct investment. The advantages of using feedback control were shown.
It should be noted that the proposed method for studying program control and feedback control is based on sequential approximation of control using functional series. Therefore, the question of estimates of the accuracy of the approximation to the optimal solution remains open.