Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University)2220610.22363/2658-4670-2019-27-2-154-164Research ArticleNumerical analysis of ecology-economic model for forest fire fighting in Baikal regionSukhodolovAlexander P<p>Doctor of Sciences (Economics), Professor, First Vice Rector - Vice Rector for Science</p>3952_2015@mail.ruSorokinaPolina G<p>Senior Lecturer of Department of Mathematics and Computer Science</p>ermolaeva_polina@mail.ruFedotovAndrey P<p>Doctor of Sciences (Geology and Mineralogy), Director of Limnological Institute SB RAS</p>info@lin.irk.ruBaikal State UniversityLimnological Institute SB RAS1512201927215416422112019Copyright © 2019, Sukhodolov A.P., Sorokina P.G., Fedotov A.P.2019<p>Forest fires lead to the serious damage of ecological state and national economy of the country. This problem is especially relevant for Siberians. According to Greenpeace, Siberian forest fires in 2019 reached record levels in the entire history of observation in terms of burning area and the amount of carbon dioxide emitted into the atmosphere. It leads not only to a deterioration in the health of Siberians, but also to environmental problems of the region. Note that the large-scale fire-prevention measures entails enormous financial costs. Therefore, economical, ecological and mathematical modeling of the situations, arose in forest fires countering, becomes actual. The paper is devoted to optimal control problem of forest fires fighting. Its prototype is the well-known Parks model. To investigate the model, we apply the modern programming language Julia, which is designed to mathematical calculations and numerical studies. We made an extensive computational experiment in this model and a numerical analysis of corresponding optimal control problems. The obtained results were examined both on the adequacy of the model, and on the possibility of using the Julia language and the included solvers of mathematical problems.</p>forest fires fightingmathematical modelingoptimal controlnumerical analysisJulia programming languageборьба с лесными пожарами, математическое моделирование, оптимальное управление, численный анализ, язык программирования Julia<p>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 peoples 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 DCMACS. 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 models 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 DCMACS. 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 DCMACS. 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 DCMACS. 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.</p>[A. P. Sukhodolov, A. A. Izmestiev, Economic accessibility of forest resources as rent-forming factor and assessment basis for forest raw material potential [Ekonomicheskaya dostupnost’ lesnykh resursov kak rentoobrazuyushchiy faktor i osnova otsenki lesosyr’yevogo potentsiala], Bulletin of Baikal State University (6 (86)) (2012) 32-35, in Russian.][A. V. Panov, A. S. Prokushkin, A. V. Bryukhanov, M. A. Korets, E. I. Ponomarev, N. V. Sidenko, G. K. Zrazhevskaya, A. V. Timokhina, M. O. Andreae, A complex approach for the estimation of carbonaceous emissions from wildfires in Siberia, Russian Meteorology and Hydrology 43 (5) (2018) 295-301.][I. V. Tikhonova, N. V. Efimova, Prevalence of the chronic respiratory tract pathology in teenagers: role of some factors [Chastota khronicheskoy patologii verkhnikh dykhatel’nykh putey u podrostkov: rol’ nekotorykh faktorov], Hygiene and Sanitation (6) (2012) 51-53, in Russian.][G. Dorrer, I. Buslov, S. Yarovoy, Conception of managing system for wild fire struggle [Kontseptsiya sistemy upravleniya bor’boy s prirodnymi pozharami], Siberian Fire and Rescue Bulletin (1 (1)) (2016) 38-44, in Russian.][A. V. Kolyada, Optimization of the process of extinguishing a forest fire using simulation [Optimizatsiya protsessa tusheniya lesnogo pozhara s ispol’zovaniyem imitatsionnogo modelirovaniya], Scientific Notes of the RSSU (8) (2010) 89-94, in Russian.][G. M. Parks, Development and application of a model for suppression of forest fires, Management Science 10 (4) (1964) 760-766. doi:10.1287/mnsc.10.4.760.][M. Parlar, R. G. Vicson, Optimal forest fire control: an extension of Park’s model, Forest Science 28 (2) (1982) 345-355. doi:10.1093/forestscience/28.2.345.][M. Parlar, Optimal forest fire control with limited reinforcements, Optimal Control Applications Methods 4 (1983) 185-191. doi:10.1002/oca.4660040208.][J. Rodriguez-Veiga, M. J. Ginzo-Villamayor, B. Casas-Mendez, An integer linear programming model to select and temporally allocate resources for fighting forest fires, Forests 9 (583) (2018) 2-18. doi:10.3390/f9100583.][C. Artigues, E. Hébrard, Y. Pencolé, A. Schutt, P. J. Stuckey, A study of evacuation planning for wildfires, in: The Seventeenth International Workshop on Constraint Modelling and Reformulation (ModRef 2018), Lille, France, 2018. URL https://hal.archives-ouvertes.fr/hal-01814083][J. Rodriguez-Veiga, I. Gomez-Costa, M. J. Ginzo-Villamayor, B. CasasMendez, J. L. Saiz-Diaz, Assignment problems in wildfire suppression: models for optimization of aerial resource logistics, Forest Science 64 (5) (2018) 504-514. doi:10.1093/forsci/fxy012.][E. V. Bogdanova, G. V. Davydova, Methodical approaches for forecasting forest fire in the Irkutsk region [Metodicheskiye podkhody k prognozirovaniyu lesnykh pozharov v Irkutskoy oblasti], in: S. V. Chuprova, N. N. Danilenko (Eds.), Revitalization of the Intellectual and Resource Potential of the Regions: New Challenges for the Management of Companies Materials of the 3rd All-Russian Conference, Irkutsk, 18 May 2017, BGU, Irkutsk, 2017, pp. 45-51, in Russian. URL http://vseup.ru/static/files/IRKUTSK_-_2017.pdf][E. A. Pyanova, A. A. Faleychik, L. M. Faleychik, Forest fires in Transbaikalia: numerical simulation [Lesnyye pozhary v Zabaykal’ye: chislennoye modelirovaniye], in: Kulagin readings: technology and technology of production processes XV International Scientific and Practical Conference: a collection of articles in 3 parts. Part 2, 2015, pp. 125-129, in Russian.][O. A. Belykh, G. D. Rusetskaya, Instruments effectiveness assessment to implement sustainable management principles of forest systems in Eastern Siberia [Otsenka effektivnosti instrumentov realizatsii printsipov ustoychivogo upravleniya lesnymi sistemami v vostochnoy Sibiri], Forestry Bulletin 23 (1) (2019) 5-13, in Russian.][J. B. Ditsevich, O. A. Belykh, G. D. Rusetskaya, Counteracting crimes in the sphere of forest resources’ use: problems and perspectives [Protivodeystviye prestupnosti v sfere lesopol’zovaniya: problemy i perspektivy], Russian Journal of Criminology 11 (2) (2017) 308-317, in Russian. doi:10.17150/2500-4255.2017.11(2).308-317.][I. Y. Kulagina, Problems of reforestation in Russia [Problemy vosstanovleniya lesov v Rossii], in: The economic development of society in modern crisis conditions. Collection of articles of the International Scientific and Practical Conference. In 3 parts. Samara, May 13, 2017, 2017, pp. 128-131 [in Russian].][M. Sherrington, Learn the language of Julia, Packt Publishing, 2016.][L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrilidze, E. F. Mishenko, Mathematical theory of optimal processes [Matematicheskaya teoriya optimal’nykh protsessov], Fizmatlit, Moscow, 1961, in Russian.][V. I. Gurman, Turnpike solutions in the procedures seeking optimal controls, Automation and Remote Control 64 (3) (2003) 399-408. doi:10.1023/A:1023209524049.]