CC BY
9
MSC 00A72
DOI: 10.14529/ mmp230306
PARAMETERS IDENTIFICATION ALGORITHM
FOR THE SUSUPLUME AIR POLLUTION PROPAGATION MODEL
S.M. Elsakov1, D.A. Drozin1, A.A. Zamyshlyaeva1, A.P. Basmanov2, V.A. Surin1, A.V. Herreinstein1, S.G. Nitskaya1
1 South Ural State University, Chelyabinsk, Russian Federation
2LLC MMK-Informservice, Magnitogorsk, Russian Federation
E-mail: elsakovsm@susu.ru, drozinda@susu.ru, zamyshliaevaaa@susu.ru,
petrovich1997basmanov@mail.ru, surinva@susu.ru, gerenshteinav@susu.ru,
nitckaiasg@susu.ru
The article presents the method of identifying the parameters of a dynamic dispersion calculation model SUSUPLUME. It is supposed that the model parameters contain not only the characteristics of the atmosphere and the pollutant, but also information about the influence of other particular conditions such as terrain, building, background, etc. The model parameters are configured based on instrumental measurements of concentrations of pollutants in the atmospheric air in the surface layer (2 meters above ground level). Three identification strategies are considered: identification of parameters by all measurements, identification of parameters by measurements of a given source and identification of parameters using another approved model. A method for weighing these strategies is proposed in the issue. The article also provides objective functions for optimization criteria, an acceptable set of parameters, an algorithm for solving an optimization problem, a decision tree of a feasible set and a global optimization algorithm.
Keywords: SUSUPLUME model; global optimization; ecology; propagation of pollutants in the atmosphere.
Introduction
The economic activity of industrial enterprises is associated with the probability of damaging the environment. Currently, most of them do not have operational analysis systems to control environmental risks.
The authors have created a software system "Ecomonitor" for industrial enterprises, which allows you to visualize the current pollution propagation, predict the environmental risks and identify the sources that have the greatest influence on control points in a given area. This system is based on the developed mathematical model SUSUPLUME [1]. The model is created for practical use, which takes into account the limited computing resources of enterprises. The model also provides online monitoring of the spread of pollutants in the atmosphere in real time.
Mathematical models for calculating the spread of emissions in the atmosphere are divided into Gaussian, Lagrangian and Eulerian models [2-4]. The developed mathematical model SUSUPLUME is a trajectory Gaussian model. There is a large number of trajectory models around the world, the most known of which are CALPUFF [5,6], RAPTAD [7,8], HYSPLIT [9], RIMPUFF [10] and FLEXPART [11]. The question of setting up mathematical models is considered rather complicated in the practical part and the methodology of determining the model parameters depends on their accuracy and the scope of their use.
The most common approach to find model parameters is to determine physical quantities based on instrumental measurements. For example, meteorological parameters and the dimensions of the calculated grid can be involved. Another method of finding parameters is to find the coefficients characterizing the operation of emission sources.
In this article, the authors describe the implementation of the algorithm of determining the parameters of the SUSUPLUME model based on measurements of concentrations of pollutants at control points in the atmospheric surface layer. It is supposed that the model parameters contain not only the characteristics of the atmosphere and the pollutant, but also information about the influence of other particular conditions.
1. A Pollution Propagation Model from Point Sources
The SUSUPLUME model [1] is presented briefly in this section. Let the concentration of pollutants at point (x, y, z) at time t be written as C (x, y, z,t) = iL1 C (x, y, z, t), where Ci is the concentration created by the i-th source, N is the total number of point sources, k is the coefficient of the model, which describes the conditions of pollutants accumulation. The concentration of pollution generated by the i-th source can be described by the relation
Ci (x,y,z,t)
mg,i (to)e
(xjjtQ^-x)2 (yj(t0,t)-y)2
2a£(x)
2аУ(х)
t--
(2n)3/2 ax (X)&y (X)az (X)
_ (ч(г0>г)-г)2
e 2<xi(i)
dt0
where D is the diameter of the concentration calculation area, is
the minimum wind
speed, mg>i(to) is the pollutants intensity created by the i-th source at time t0, ox,ay,az are functions which describes the diffusion rate of pollutants depending on time. These functions can be calculated as
x
X
a
VsA+sBx' VsC+sDx'
x=
\J u(s)2+v(s)2ds
to
where (u,v) is the vector of wind direction depending on time.
The movement of center of mass of the single emission in the horizontal plane is described by the relations
t t Xi (to,t) =xi + j u (t) dt, yi (to,t) =yi+ J v (t) dt,
to to
where (xi,yi) are the coordinates of the i-th source. The height of the single emission can be written as
zi (t0, t) =max
7(a+7)
c f b c \ ( ea(t-to)-1\ ce-Y(t-t0) b (t-to) n + ( ) ( - ) +—-----^-
a (a+Y)
a
Y (a+Y) a
h \Ma 11 a= —ii. b= ——1 a, c=
1 Mg 1
Ma (Tj-Tg) g MgTa '
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2023. Т. 16, № 3. С. 74-82
t
t
ax ay
where hi is the height of the i-th source, wi is the initial vertical speed of single emission from the mouth of the i-th source, mg is the mass of single emission, ma is the mass of displaced air, ^ is the viscosity coefficient, which describes the dependence of the medium resistance to the movement of the single emission, Ta is the air temperature, Tg is the temperature of the single emission, Ma is the molar mass of the air, Mg is the molar mass of the single emission.
2. Problem of the SUSUPLUME Model Parameter Identification
There are many models of pollution propagation in the atmosphere. Most models are calibrated on a large amount of raw data. To improve the accuracy of the models, researchers introduce additional factors into consideration in order to take into account the relief, building and other characteristics. Such models work well, but at some points they will give an essential error, since the calibration of the model is carried out for the entire set of data.
Three strategies can be used to calibrate the SUSUPLUME model. The first strategy is standard calibration for the entire set of data. This implies the presence of a large number of real measurements at various locations and conditions for measuring concentrations. The second strategy involves calibrating the model based on data from only one location and involves a large number of real measurements for one specific point/area. The third strategy involves calibrating the SUSUPLUME model based on the results of another model simulation.
Using a combination of these strategies allows you to determine a set of parameters for the SUSUPLUME model, which will either be as accurate as possible at some point/area, or simulate all available measurements (requirements of normative models).
Simultaneous consideration of all three strategies leads us to the problem of multi-criteria optimization, for this reason we will use the following weighted criterion to determine the appropriate strategy
F (x) =w/ (x) +W2/2 (x) +W3/3 (x) ,
where w1, w2, w3 G [0; 1] are the weights of strategies and /1 (x), /2 (x) , /3 (x) : Rn ^ R are the objective functions for each strategy.
3. Objective Functions for Optimization Criteria
Let (xj,yj, Zj) be the point coordinates, where the measurements were made, tj is the time of the j-th measurement, x is the model parameter vector. Let cmeas(xj,yj,Zj,tj) be the measurements of pollutants and cmod (x,xj ,yj, Zj ,tj) be the model calculations. Suppose we have M measurements. Measuring concentrations devices of pollutants have a number of features:
1. The measurement error is about 20-25%.
2. Each pollutant has its own approved threshold limit value (TLV). Concentration values of less than 1% TLV, as a rule, go beyond the measuring range of devices.
3. The measuring range of the device gains concentration values within the same order, but it can vary significantly more.
To solve the calibration problem, it is proposed to use several criteria that may be relevant in various situations.
1. Using the standard deviation criterion which can be written as
1 M
/и=мЕ(сГ-г)2 (i)
M ^ 3 J
j=1
This criterion may be relevant if all measurements are in the measuring range of the device or they are received from another model.
2. Using the average absolute deviation criterion which can be calculated as
1 M
/^)=MElcrd"SmeaS|- (2)
M j=1
This criterion improves the previous one and can be used for the same conditions if there are significant pollutants.
3. Using the average relative deviation criterion which can be described by the relation
i M I „mod_„meas I
f{x)=±yn—(3)
■> y ' M _/ „meas V >
j=1 j
This criterion can be applied only when all measurements are within the measuring range of the device.
4. As for the SUSUPLUME model, the objective function of the standard deviation can be improved by finding the coefficient k analytically. Let us describe the objective function with the coefficient
1 M / N \ 2
/ =m E fcE <3 fo'z^)~crs • (4)
j= 1 V i=1 J
Note that the function depends quadratically on k. This function reaches its minimum
meas V^N C (x y z t ) j Z^ii=i Ci (xj, yj, z3 ,b3)
when m N
£ Ej = 1 Cj"eaS Ei= 1 Ci (Xj > Уз > Z3 > tj)
E^i (E ill Ci (x ,yj, zj ,tj))
Then the objective function (4) can be calculated as
/ \ 2
WEM=1 j^EN=1 Ci (Xj ,yj,Zj,tj) N X
M j=i
f (x) =TT > -----9" > Ci {Xj.yi, Zj , tj) — c™
v Ej=1 ^Ei=iCi(x>yj>zj,tj)) = у
(5)
5. Using the quadratic residual. Let us describe the sets of measures: A = {j\djieas < a,j}, where a,j is the lower measure limit of the device for the j-th measurement, B = {j\cmmeas < bj}, where bj is the lower measure limit of the device for the j-th measurement, C = {j\aj < djieas < bj} which is a set of measures in range of device limits. While using the third strategy, aj can be set as 0.01TLV and bj can be infinite. Then the objective function can be written as
/(*)=if (e f(^-crl)2+E(crd-rs)2+E((crd-^)+)2] > (6)
\j e^ V + j ec jeB J
where (x), = max(x; 0) .
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2023. Т. 16, № 3. С. 74-82
6. Using the average absolute residual. The function can be described as
/ (*) =if (E («i-<rv E lsmod-casl+E (crd-^)+J • (7)
j eA jec jeB J
7. Using the average absolute residual. This can be estimated as
4
^ fa7—Cm°d^ I „mod cmeas\ b^
y -
cmeas
meas
yeA aj jec cj j eB bj j
(8)
The quality of the model is commonly evaluated by criterion (8), however, the use of this criterion for optimization algorithms is rather complicated. Initially, the function is non-smooth in the parameter space. There are also sums in the objective function that can be equal to zero.
Based on these considerations, the search for a solution is done in three stages: fitting parameters using a differentiable function (5), solving the problem using function (6) and solving the problem using function (8).
4. The Set of Admissible Parameters
It is assumed that the following identification parameters will be used for the SUSUPLUME model: molar mass of the gas-air mixture Mg G [3; 1550] and the parameters of a single emission dynamics G [10-5; 10], sA, sC G [10-6; 10], sB,sD G
[10-8; 10-3], k G [10-5; 105]. Some values have a linear scope, and some values are logarithmic. Therefore, variables with a logarithmic scope are optimized linearly, which is achieved by using the appropriate functional transformation.
The found optimal parameters may have no physical meaning and only describe the dynamics of the emission. If, after solving the problem, the parameters were on the border of the acceptable domain, then this domain should be expanded taking into account the physical meaning.
5. Solving the Optimization Problem Algorithm
Since the problem of selecting the parameters of the SUSUPLUME model is multidimensional and multi-extreme, it requires the use of multi-extreme optimization algorithms. It takes low time to count the objective function, as the calculation can be effectively parallelized. A parallel multistart algorithm in a hyperparallelepiped will be used to solve this problem.
The multistart algorithm is applied as a global optimization algorithm and the BOBYQA algorithm is used as a local optimization algorithm [12]. For guaranteed convergence of the global optimization algorithm, it is required to guarantee that the set of starting points for the BOBYQA algorithm is everywhere dense in X. A decision tree of feasible set will be used to apply the convergence.
6. A Dicision Tree of Admissible Set
Let us describe a decision tree of admissible set of parameters and the basic operations with this tree.
A decision tree of admissible set is a binary tree which consists of vertices that represent a hyperparallelepiped in a d-dimensional space. The vertex of the tree represents the entire admissible set X. The descendants of a vertex always represent two hyperparallelepipeds of dimension h mod d+1 (h is the level of the vertex in the tree (counting from the root)) which are obtained by splitting the hyperparallelepiped in half. The entire admissible set is always represented by the union of all the leaf nodes.
The operation of saving a point (a point where the calculations are done) consists of finding a leaf in the current tree in which the calculated point lies and saving it in this leaf. Such a leaf always exists and is the only one. The execution time of this operation is O(log(S)), where S is the total number of leaves in the tree.
The operation of selecting the next point is more complicated. A hyperparallepiped with the largest diagonal is found. If this hyperparallepiped does not contain saved points, then the center of this hyperparallepiped is returned as an answer.
If the hyperparallepiped contains saved points, then it is divided into two descendants according to the rules described above and all the points saved in it come down into its descendants. The selection of the leaf with the maximum diagonal is repeated again. The search time of the leaf with the maximum diagonal is O(log(S)).
7. Global Optimization Algorithm
The parallel multistart algorithm consists of the following steps: Step 0. Set the number of parallel calculations. Let this number be P. Step 1. Select P starting points for the BOBYQA algorithm from the decision tree of an admissible set.
Step 2. Run the BOBYQA algorithm for each selected point.
Step 3. Save all the points at which the BOBYQA algorithm calculated the value of the objective function in the decision tree of an admissible set.
Step 4. If the calculation time has exceeded the limit for calculation, then stop the algorithm.
The proposed algorithm makes it possible to construct an everywhere dense sequence of initial points. This is achieved because the sizes of hyperparallelepipeds decrease uniformly and an empty hyperparallelepiped with the maximum diagonal is always selected.
References
1. Elsakov S.M., Drozin D.A., Herreinstein A.V., Krupnova T.G., Nitskaya S.G., Olenchikova T.Yu., Zamyshlyaeva A.A. Numerical Study of the Susuplume Air Pollution Model. Bulletin of the South Ural State University. Series: Mathematical Modelling, Programming and Computer Software, 2020, vol. 13, no. 4, pp. 5-18. DOI: 10.14529/mmp200401
2. Leelossy A., Molnar F.Jr., Izsak F., Havasi A., Lagzi I., Meszaros R. Dispersion Modeling of Air Pollutants in the Atmosphere: a Review. Central European Journal of Geosciences, 2014, vol. 6, no. 3, pp. 257-278. DOI: 10.2478/s13533-012-0188-6
3. Yao Rentai. Atmospheric Dispersion of Radioactive Material in Radiological Risk Assessment and Emergency Response. Progress in NUCLEAR SCIENCE and TECHNOLOGY, 2011, vol. 1, pp. 7-13. DOI: 10.15669/pnst.1.7
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2023. Т. 16, № 3. С. 74-82
4. Belov I.V., Bespalov M.S., Klochkova L.V., Pavlova N.K., Suzan D.V., Jishkin V.F. Comparative Analysis of Models of Pollutions Spreading in Atmosphere. Mathematical Modeling, 1999, vol. 11, no. 8, pp. 52-64.
5. Ghannam K., El-Fadel M. Emissions Characterization and Regulatory Compliance at an Industrial Complex: An Integrated MM5/CALPUFF Approach. Atmospheric Environment, 2013, vol. 69, pp. 156-169. DOI: 10.1016/j.atmosenv.2012.12.022
6. Rzeszutek M. Parameterization and Evaluation of the CALMET/CALPUFF Model System in Near-Field and Complex Terrain - Terrain Data, Grid Resolution and Terrain Adjustment Method. Science of the Total Environment, 2019, vol. 689, pp. 31-46. DOI: 10.1016/j.scitotenv.2019.06.379
7. Ge Wang., Ostoja-Starzewski M. Influence of Topography on the Phoenix CO2 Dome: a Computational Study. Atmospheric Science Letters, 2004, vol. 5, no. 5, pp. 103-107. DOI: 10.1002/asl.67
8. Jay Wu., Chung-Hsin Lu., Shu-Jun Chang., Yung-Muh Yang., Bor-Jing Chang., Jen-Hsin Teng. Three-Dimensional Dose Evaluation System using Real-Time Wind Field Information for Nuclear Accidents in Taiwan. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2006, vol. 565, no. 2, pp. 812-820. DOI: 10.1016/j.nima.2006.05.197
9. McGowan H., Clark A. Identification of Dust Transport Pathways from Lake Eyre, Australia using Hysplit. Atmospheric Environment, 2008, vol. 42, pp. 6915-6925. DOI: 10.1016/j.atmosenv.2008.05.053
10. Lepicard S., Heling R., Maderich V. POSEIDON/RODOS Models for Radiological Assessment of Marine Environment after Accidental Releases: Application to Coastal Areas of the Baltic, Black and North Seas. Journal of Environmental Radioactivity, 2004, vol. 72, no. 1-2, pp. 153-161. DOI: 10.1016/S0265-931X(03)00197-8
11. Stohl A., Hittenberger M., Wotawa G. Validation of the Lagrangian Particle Dispersion Model Flexpart Against Large-Scale Tracer Experiment Data. Atmospheric Environment, 1998, vol. 32, no. 24, pp. 4245-4264. DOI: 10.1016/S1352-2310(98)00184-8
12. Powell M.J.D. The BOBYQA Algorithm for Bound Constrained Optimization without Derivatives. Centre for Mathematical Sciences. Report DAMTP, 2009, no. NA06, 39 p.
Received May 11, 2023
УДК 004.942 DOI: 10.14529/mmp230306
АЛГОРИТМ ИДЕНТИФИКАЦИИ ПАРАМЕТРОВ МОДЕЛИ РАСПРОСТРАНЕНИЯ ЗАГРЯЗНЕНИЯ В АТМОСФЕРНОМ ВОЗДУХЕ SUSUPLUME
С.М. Елсаков1, Д.А. Дрозин1, А.А. Замышляева1, А.П. Басманов2, В.А. Сурин1, А.В. Геренштейн1, С.Г. Ницкая1
1 Южно-Уральский государственный университет, г. Челябинск, Российская Федерация
2ООО «ММК-Информсервис», г. Магнитогорск, Российская Федерация
В статье представлен метод идентификации параметров динамической модели расчета рассеиваний SUSUPLUME. Предполагается, что параметры модели содержат не
только характеристики атмосферы и загрязняющих веществ, но и информацию о влиянии других условий, таких как тип местности, застройка, фон и т.д. Настройка параметров модели осуществляется на основе инструментальных замеров концентраций загрязняющих веществ в атмосферном воздухе в приземном слое (2 метра над уровнем земли). Рассматриваются три стратегии идентификации: идентификация параметров по всем измерениям, идентификация параметров по измерениям заданного источника и идентификация параметров на основе другой утвержденной модели. В третьей стратегии предлагается идентифицировать параметры модели так, чтобы модель была максимально похожа на утвержденную модель. Также предложен метод взвешивания этих стратегий. В статье также представлены целевые функции для критериев оптимизации, множество допустимых наборов параметров, алгоритм решения оптимизационной задачи, дерево решений для допустимого множества и алгоритм глобальной оптимизации.
Ключевые слова: модель SUSUPLUME; глобальная оптимизация; экология; распространение загрязняющих веществ в атмосферном воздухе.
Литература
1. Elsakov, S.M. Numerical Study of the Susuplume Air Pollution Model / S.M. Elsakov, D.A. Drozin, A.V. Herreinstein, T.G. Krupnova, S.G. Nitskaya, T.Yu. Olenchikova, A.A. Zamyshlyaeva // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2020. - Т. 13, № 4. - С. 5-18.
2. Leelossy, A. Dispersion Modeling of Air Pollutants in the Atmosphere: a Review / A. Leelossy, F.Jr. Molnar, F. Izsak, A. Havasi, I. Lagzi, R. Meszaros // Central European Journal of Geosciences. - 2014. - V. 6, № 3. - P. 257-278.
3. Rentai Yao, Atmospheric Dispersion of Radioactive Material in Radiological Risk Assessment and Emergency Response / Rentai Yao // Progress in NUCLEAR SCIENCE and TECHNOLOGY. - 2011. - V. 1. - P. 7-13.
4. Белов, И.В. Cравнение моделей распространения загрязнений в атмосфере / И.В. Белов, М.С. Беспалов, Л.В. Клочкова, Н.К. Павлова, Д.В. Сузан, В.Ф. Тишкин // Математическое моделирование. - 1999. - Т. 11, № 8. - С. 52-64.
5. Ghannam, K. Emissions Characterization and Regulatory Compliance at an Industrial Complex: An integrated MM5/CALPUFF approach / K. Ghannam, M. El-Fadel // Atmospheric Environment. - 2013. - V. 69. - P. 156-169.
6. Rzeszutek, M. Parameterization and Evaluation of the CALMET/CALPUFF Model System in Near-Field and Complex Terrain - Terrain Data, Grid Resolution and Terrain Adjustment Method / M. Rzeszutek // Science of the Total Environment. - 2019. - V. 689. - P. 31-46.
7. Wang Ge. Influence of topography on the Phoenix CO2 dome: a computational study / Wang Ge, M. Ostoja-Starzewski // Atmospheric science letters. - 2004. - V. 5, № 5. - P. 103-107.
8. Jay Wu. Three-Dimensional Dose Evaluation System using Real-Time Wind Field Information for Nuclear Accidents in Taiwan / Jay Wu, Chung-Hsin Lu, Shu-Jun Chang, Yung-Muh Yang, Bor-Jing Chang, Jen-Hsin Teng // Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. - 2006. - V. 565, № 2. - P. 812-820.
9. McGowan, H. Identification of Dust Transport Pathways from Lake Eyre, Australia Using Hysplit / H. McGowan, A. Clark // Atmospheric Environment. - 2008. - V. 42. -P. 6915-6925.
10. Lepicard, S. POSEIDON/RODOS Models for Radiological Assessment of Marine Environment after Accidental Releases: Application to Coastal Areas of the Baltic, Black
Вестник ЮЮУрГУ. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2023. Т. 16, № 3. С. 74-82
and North Seas / S. Lepicard, R. Heling, V. Maderich // Journal of Environmental
Radioactivity. - 2004. - V. 72, № 1-2. - P. 153-161.
11. Stohl, A. Validation of the Lagrangian Particle Dispersion Model Flexpart Against Large-
Scale Tracer Experiment Data / A. Stohl, M. Hittenberger, G. Wotawa // Atmospheric
Environment. - 1998. - V. 32, № 24. - P. 4245-4264.
12. Powell, M.J.D. The BOBYQA Algorithm for Bound Constrained Optimization without
Derivatives / M.J.D. Powell // Centre for Mathematical Sciences. Report DAMTP. - 2009. -
№ NA06. - 39 p.
Сергей Михайлович Елсаков, кандидат физико-математических наук, кафедра «Прикладная математика и программирование», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), elsakovsm@susu.ru.
Дмитрий Александрович Дрозин, кандидат экономических наук, кафедра «Прикладная математика и программирование», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), drozinda@susu.ru.
Алена Александровна Замышляева, доктор физико-математических наук, кафедра «Прикладная математика и программирование», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), zamyshliaevaaa@susu.ru.
Алексей Петрович Басманов, инженер-программист первой категории ООО «ММК-Информсервис» (г. Магнитогорск, Российская Федерация), petrovich1997basmanov@mail.ru.
Владимир Анатольевич Сурин, старший преподаватель, кафедра «Прикладная математика и программирование», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), surinva@susu.ru.
Аркадий Васильевич Геренштейн, кандидат физико-математических наук, кафедра «Прикладная математика и программирование», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), gerenshteinav@susu.ru.
Светлана Георгиевна Ницкая, кандидат технических наук, доцент, кафедра «Экология и химическая технология», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), nitckaiasg@susu.ru.
Поступила в редакцию 11 мая 2023 г.
