Numerical analysis of ecology-economic model for forest fire fighting in Baikal region

Full Text

1. Introduction The Irkutsk Region is one of the largest constituent entities of the Russian Federation. Its area is about 774.846 square km., it is slightly less than the area of the Republic of Turkey (780.580 sq. km.), and also exceeds the territory of France, Germany and many other European countries. Most of the territory of the Irkutsk region, namely 71.5 million hectares, or 92% of its territory, is occupied by the forest. About 12% of timber reserves of ripe forests of the country are concentrated here, and the share of especially valuable coniferous species is significant even on the scale of the planet. According to the press service of the Ministry of Emergency Situations of Russia, on July 28, the 143 forest fires with a total area of 597,298 hectares are active in the Irkutsk Region. Totally, in Siberia and the Far East, forest fires are burning in the area comparable to Belgium. Of course, such an essential loss of natural resources leads to negative environmental, economic and social consequences [1-3]. For instance, in the work [2] it was established that during intensive burning of the taiga, the concentration of carbon monoxide increases by almost 30 times in comparison with the background content in the air, methane by 2 times, carbon dioxide in 8%. Such exceedences lead to the health deterioration of the inhabitants of the Irkutsk region [3]. In addition, due to the annual, large-scale forest fires blazing near Lake Baikal, chemical components, such as ammonium, expedite the reproduction of various microorganisms that destroy the aquatic ecosystem of the Baikal region. The most important problem in forest fires fighting, besides the protecting of people’s lives, is a quick and effective fire suppression, planned to minimize the total damage. Controlling of the process of suppression, transportation of forces to the place of fire is made by employees of forest protection organizations. In the most cases they make decisions based on their personal experience. But even with experience, defining an optimal fire fighting plan is, often, a quite difficult task. For many years, scientists have been studying models that allow them to find optimal solutions of fire fighting forces control under an active forest fire. These attempts are being made to take into account the characteristics of the spread of fire, the capabilities of the available fire-fighting forces and equipment, topographic features and other factors. Here we mention works [4-16]. Note that the papers [4, 5] continues the research originated in [7, 8]. The analysis of some foreign works on the subject [9-11] shows that more investigations use modern programming languages, such as Python, R, Java, and so on. These software products are applied to analyze and visualize the data. Actually, it improves the quality of the research. However, despite a sufficient number of software tools, in explorations of Russian scientists such products are not so widely used as abroad. In this paper, we apply rather new universal programming language Julia [17]. Its development was begun by scientists in 2009, and its first version was published in 2012. The research uses the latest version of Julia, presented in 2018. Julia is a modern high-level programming language with dynamic typing for mathematical calculations, which is used to develop research software, approbation and test of new problem-solving methods. The essential advantages of this language are simple syntax and speed of program execution. It is often chosen by astronomers, robotics and financiers. In the computational experiment, the IPOPT (Interior Point OPTimizer) package is used. This package is meant to numerical solving of optimization problems of a large dimension. 2. Problem statement The paper focuses on the optimization model for the dispatching and withdrawal of fire-fighting forces to the place of a forest fire [6, 8]. After 156 DCM&ACS. 2019, 27 (2) 154-164 certain transformations, it can be written as the following optimal control problem ( ): = ∫ (22() - ()() + ()) + 2() → min, ̇ = 1 - 2, () = , - ⩽ 1() - 2() ⩽ , 1() ⩾ 0, 2() ⩾ 0, 0 ⩽ () ⩽ , () (, ()) = () - ( - = 0, ) where () = ( - ) - and () = + (). The state variable () denotes the size of the fire fighting force at the moment . The pair of control variables 1() and 2() represent dispatching and withdrawal rates of the reinforcements at the time , respectively. The trajectory (⋅) is supposed to be a piecewise smooth function, while the control functions 1(⋅), 2(⋅) are piecewise continuous. Let us give an economical interpretation of model’s parameters: § and are the time moment of initial attack and the final time moment, when the fire is brought under control, respectively; , in general, is supposed to be non-fixed; § () is a function of the fire spread rate in the absence of fire fighting forces; § is the cost to transportation (i.e., dispatching or withdrawal) of fire fighting forces (currency unit per force unit); § is the parameter characterizing the loss of the forest during uncontrolled burning per time unit (currency unit per time); § is the cost per unit area of forest damaged by fire (currency unit per area unit); § is the cost of the fire-fighting (currency unit per unit of forces - time); § is the maximal rate of fire-fighters withdrawal (force unit per time); § is the maximal rate of fire-fighters dispatching (force unit per time); § is the ratio of the effectiveness of fire fighting in this area (unit of forces per time); § is the initial attack force (force unit); § is the maximal limit of fire fighting forces (force unit). We point out some features of problem ( ). This linear optimal control problem contains a) the terminal state constraint in the form of equality, and b) the pointwise state constraint. These features significantly complicate the analytical investigation of problem ( ), even under the linearity of the dynamical system [8]. In the mentioned work, the problem was analyzed using the Pontryagin Maximum Principle [18]. As a result, the author of [8] gave an explicit formula for the optimal control, A. P. Sukhodolov, P. G. Sorokina, A. P. Fedotov, Numerical analysis… 157 which depends on values of two unspecified parameters: the final time and the scalar Lagrange multiplier which corresponds to the terminal state constraint. For a further specification of the optimal control, it was proposed to use a grid search of proper values of the unknown parameters. Then, the optimality conditions should be checked for each choice of and . Our paper is devoted to the numerical investigation of problem ( ). We use the so-called “direct approach”, i.e., an approach to the study of optimal control problems, when the total discretization of the dynamic optimization problem is applied. Further, the obtained problem is solved by methods and tools of mathematical programming. This approach is often criticized by experts in the field of optimal control, but it often turns out to be effective in solving practical optimization problems. At the discretization stage of problem ( ), we apply the explicit Euler scheme, firstly reducing problem ( ) to the Mayer form. We introduce an unessential state variable , which derivative coincides with the integrand of the cost functional (the initial condition for is trivial). The problem may be rewritten as follows (problem (1)): = () + 2() → min, ̇ = 1 - 2, ̇ = 22 - () + (), () = , () = 0, - ⩽ 1() - 2() ⩽ , 1() ⩾ 0, 2() ⩾ 0, 0 ⩽ () ⩽ , () (, ()) = () - ( - = 0, ) () = ( - ) - , () = + (). Note that problems ( ) and (1) are equivalent to each other. Meanwhile, problem (1) contains the terminal and pointwise state conditions as well. 3. Numerical analysis Let us consider a discrete analogue of problem (1) using the direct Euler scheme. Here we suppose that the final time moment is given. Introduce the -point grid of the time interval [, ]: = 0 < 1 < … < -1 < = . As usual, the time lag is calculated as ℎ = - 0 . Discrete problem ( ) takes the following form: = () + 2() → min; (+1) = () + ℎ[1() - 2()], (+1) = () + ℎ[22() - ()() + ()], 158 DCM&ACS. 2019, 27 (2) 154-164 (0) = , (0) = 0, - ⩽ 1() - 2() ⩽ , 1() ⩾ 0, 2() ⩾ 0, = 0, 1, … , - 1; 0 ⩽ () ⩽ , = 0, 1, … , ; () where (, ()) = () - ( - 0 = 0, ) () = ( - 0) - , () = + (). Here, we use the previous notations for state and control variables and parameters of the problem. Furthermore, note that we have the sequences {1, 2}, {, } thought as control and state, respectively: 1 = {1()}, 2 = {2()}, = 0, 1, … , -1, = {()}, = {()}, = 0, 1, … , . Notice that () is a linear programming problem. It contains 4 variables and 8 + 1 conditions. We have tested more than 30 variants of the parameters of problem (). Some of them were unsuccessful. Such outcomes we associate with the absence of admissible plans of the problem. We think that the choice of values of the final time moment was improper in certain trials. Let us show the results of certain numerical experiments. Here, we present some interesting examples. Each of them is accompanied by a table indicating the values of the parameters. Illustrations are arranged as follows: on the upper graph, the trajectory component () is depicted, and on the lower graphs we show controls 1() and 2(), from left to right, respectively. 1. The first group of examples: the rate of fire spread is constant Example 1. The parameters of problem () are presented in Table 1. Parameters for Example 1 Table 1 Parameter 0 0 Value 0 22 1000 10 1 5 30 30 1 50 2 1000 Let us give an interpretation of the found solution (see Figure 1). The corresponding value of the cost functional ≈ 75. The initial attack force at the moment is characterized by 0 = 2. The high rate of fire spread ( ≡ 50) compels us to use the maximum speed dispatching of new fire fighting forces. Then, the value of forces takes a turnpike state ≈ 32, and we keep it up to the time of withdrawal. The initial and final time intervals are characterized by maximum values of dispatching (1) and withdrawal (2) speeds. A. P. Sukhodolov, P. G. Sorokina, A. P. Fedotov, Numerical analysis… 159 Figure 1. Graphs of optimal trajectory (a) and optimal controls (b, c) in Example 1 Note that the state constraints on () are inactive in this example (the contrary case is shown by the Example 3). Example 2. We decrease some parameters: the rate of fire spread , the maximal rate of fire-fighters withdrawal , the final time . The updated data is shown in Table 2. Parameters for Example 2 Table 2 Parameter 0 0 Value 0 10 1000 5 1 5 1 30 1 5 2 100 The solution is presented on Figure 2. It refers to the high speed of withdrawal fire forces at the final period. The most effective attack needs only forces in the place (at the initial time). In this case the optimal value of the cost functional is ≈ 51. Figure 2. Graphs of optimal trajectory (a) and optimal controls (b, c) in Example 2 Note that control 1 is rather close to zero, and we look at its graph as computational error. Although the value of fire-fighting forces in the place is decreased to = 100, the pointwise state constraints remain inactive. 160 DCM&ACS. 2019, 27 (2) 154-164 Parameters for Example 3 Table 3 Parameter 0 0 Value 0 100 1000 10 10 1 30 30 1 5 2 1000 Example 3. We essentially increase the time interval and change values of and , and some other parameters (see Table 3). Figure 3 shows that the pointwise state constraint becomes active. The trajectory graph means increasing of the initial value of forces 0 = 2 to the maximal level = 1000. Herewith, both controls 1 and 2 take maximum values on the initial and final time intervals, respectively. Figure 3. Graphs of optimal trajectory (a) and optimal controls (b, c) in Example 3 2. The second group of examples: the variable rate of fire spread Here, we suppose that the rate of fire spread increases by the following law: 1, ∈ [1, 6], () = { 15, ∈ (6, 12]. Also, we fix a number of time-grid points = 300. Note that problem ( ) was investigated in [5], where some features (advantages and shortcomings) of the model were indicated. Particularly, the author said that the optimal solution is characterized by three stages of control. The first of them corresponds to the maximum speed of fire fighting forces dispatching. On the second stage all involved forces fight with the fire. And then, fire fighters are withdrew with the maximal rate. The previous examples correspond to this consequence. However, even linear problems of dynamic optimization admit “singular” intervals of control, where the last one can take intermediate values. In particular, the considered model is able to contain such intervals. This feature is illustrated by the following examples. Example 4. The data of the problem is presented in Table 4. Consider Figure 4. A. P. Sukhodolov, P. G. Sorokina, A. P. Fedotov, Numerical analysis… 161 Parameters for Example 4 Table 4 Parameter 0 0 Value 0 12 1 0 0 0 3 3 1 0 10 Figure 4. Graphs of optimal trajectory (a) and optimal controls (b, c) in Example 4 One can see that control 1 takes intermediate values of admissible set [0; 3]. Note that the “jump” of control 2 is inessential, and its values are close to zero, in fact. The found solution admits the following conclusions. On the first stage, when a fire spreading rate = 1, the most effective fire extinguishing strategy was achieved at the time point = 6. The fire was localized. However, a significant increase in the rate of fire spread (up to the level = 15) required the use of additional forces (such a situation, for example, is due to the weather deterioration). The short time of the pause in dispatching of fire-fighters is associated with an excess of the cost of transporting fire forces in comparison with the damage of the action of fire. In this example, the pointwise state constraint was again inactive. Example 5. Minor changes of parameters and 0 entail certain changes in the controls (see Table 5 and Figure 5). Note, there are also time intervals with intermediate values of the controls. Parameters for Example 5 Table 5 Parameter 0 0 Value 0 12 1 0 0 0 3 3 1,1 1 10 In order to localize the fire, the fire prevention forces increased gradually. It is noteworthy that the found solution contains two “turnpike” intervals [19]. These intervals are characterized by the constancy of the trajectory on the plot of the function . Each of them corresponds to the fire spread rate levels = 1 and = 15, respectively. The graph of control of 2 is close to zero (we 162 DCM&ACS. 2019, 27 (2) 154-164 assume, as before, the depicted “blow-ups” are computational errors). At the same time, the descending control peaks of 1 are currently not explicable by the authors and require additional analysis. Figure 5. Graphs of optimal trajectory (a) and optimal controls (b, c) in Example 5 4. Conclusion Modern mathematical methods and tools are currently very accessible for applied research. Of course, their use requires certain skills and understanding of the field of research. The article shows how using the Julia programming language one can numerically investigate some applied mathematical models and solve the corresponding optimization problems. The results of the experiments are very interesting. Examples 1-3 illustrate rather obvious strategies of forest fire fighting, which applied on practice. Note that the results of calculations for example 2 correspond to the behavior of decision makers in the north of the Irkutsk region. Such tactics, of the non-attraction of additional countervailing forces, is substantiated by the economic inefficiency of them. Examples 4, 5 also interesting from an applied point of view, are entertaining by mathematics view as well. Apparently, their solutions contain the trunk modes of dynamic systems (in the last example one can see two turnpike intervals) [19]. The noted finding requires further study with the involvement of the corresponding mathematical apparatus. Further research will be associated with a more detailed analytical study of the obtained results, consideration of nonlinear modifications of the presented model and more complicated statements of optimization problems.

About the authors

Alexander P Sukhodolov

Baikal State University; Limnological Institute SB RAS

Author for correspondence.
Email: 3952_2015@mail.ru

Doctor of Sciences (Economics), Professor, First Vice Rector - Vice Rector for Science

11, Lenin str., Irkutsk, 664003, Russian Federation; 3, Ulan-Batorskaya St., Irkutsk, 664033, Russian Federation

Polina G Sorokina

Baikal State University; Limnological Institute SB RAS

Email: ermolaeva_polina@mail.ru

Senior Lecturer of Department of Mathematics and Computer Science

11, Lenin str., Irkutsk, 664003, Russian Federation; 3, Ulan-Batorskaya St., Irkutsk, 664033, Russian Federation

Andrey P Fedotov

Limnological Institute SB RAS

Email: info@lin.irk.ru

Doctor of Sciences (Geology and Mineralogy), Director of Limnological Institute SB RAS

3, Ulan-Batorskaya St., Irkutsk, 664033, Russian Federation

References