CC BY
35
Оригинальная статья / Original article УДК 621.311.004-021.387
DOI: http://dx.doi.org/10.21285/1814-3520-2020-3-582-595
Алгоритм поиска адаптивной траектории утяжеления для энергосистем большой размерности
© Н.Л. Бацева, В.А. Сухоруков
Национальный исследовательский Томский политехнический университет, г. Томск, Россия
Резюме: Цель исследования - создание и апробация алгоритма поиска адаптивной траектории утяжеления, соответствующей наименьшему предельному по статической апериодической устойчивости перетоку активной мощности для текущей схемно-режимной ситуации в энергосистеме. Для разработки алгоритма использовано значение определителя матрицы Якоби, методы кластерного анализа, а также программный комплекс RastrWin, пакет Microsoft Excel, функции Matlab и методы математической статистики. Предложен, разработан (на языках программирования Visual Basic Script, Visual Basic for Application) и апробирован алгоритм, базирующийся на изменении численного значения определителя матрицы Якоби в зависимости от изменения активной мощности генерирующего оборудования и активной мощности нагрузок энергосистемы (на схеме части Иркутской энергосистемы). Определены сенсорные нагрузочные узлы для исследуемого контролируемого сечения № 1, связывающего энергорайон с остальной частью энергосистемы. Установлено, что увеличение величины мощности (на 0,197 МВт - для узлов 220 кВ, на 0,076 МВт - для узлов 110 кВ, на 0,112 МВт - для узлов 6-10 кВ) в этих узлах приводит к значительному изменению значений определителя матрицы Якоби (на 0,103 о.е. - для узлов 220 кВ, на 0,926 о.е. - для узлов 110 кВ, на 0,33 о.е. - для узлов 6-10 кВ). Установлено, что разница в 27 МВт между значениями предельного перетока, определенного по статической апериодической устойчивости, и предельного перетока, рассчитанного с использованием траектории утяжеления, объясняется отличием исходных расчетных моделей, которые использованы при утяжелении режима. Апробация алгоритма на контролируемом сечении № 1 Иркутской энергосистемы позволила сделать вывод о его пригодности для поиска адаптивной траектории утяжеления с ее последующим применением при расчете и анализе предельных по статической апериодической устойчивости перетоков активной мощности в контролируемом сечении.
Ключевые слова: траектория утяжеления, определитель матрицы Якоби, статическая апериодическая устойчивость, схема большой размерности, сенсорные узлы, контролируемое сечение
Информация о статье: Дата поступления 18 декабря 2019 г.; дата принятия к печати 29 мая 2020 г.; дата он-лайн-размещения 30 июня 2020 г.
Для цитирования: Бацева Н.Л., Сухоруков В.А. Алгоритм поиска адаптивной траектории утяжеления для энергосистем большой размерности. Вестник Иркутского государственного технического университета. 2020. Т. 24. № 3. С. 582-595. https://doi.org/10.21285/1814-3520-2020-3-582-595
Search algorithm for adaptive weighting trajectory in large-scale power systems
Natalia L. Batseva, Vasiliy A. Sukhorukov
National Research Tomsk Polytechnic University, Tomsk, Russia
Abstract: The study is aimed at the creation and testing of a search algorithm for the adaptive weighting trajectory corresponding to the smallest ultimate active-power flow for the current circuit-mode situation in terms of aperiodic stability in a power system. To develop the algorithm, a cluster analysis approach was taken using the Jacobi matrix determinant in RastrWin, Microsoft Excel and Matlab software environments alongside various mathematical statistical methods. As a result, the proposed algorithm was developed (in the programming languages of Visual Basic Script and Visual Basic for Applications) and tested based on a change in the numerical value of the Jacobi matrix determinant depending on changes in the active power of both the generating equipment and the power system loads (in the diagram for a part of the Irkutsk power system). Sensory load nodes were determined for the investigated controlled section No. 1 connecting the energy district with the rest of the energy system. An increase in the power (by 0.197, 0.076 and 0.112 MW for nodes of 220, 110 and 6-10 kV, respectively) in these nodes was shown to result in significant changes in the values of the Jacobi matrix determinant (by 0.103, 0.926 and 0.33 p.u. for nodes of 220, 110 and 6-10 kV, respectively). The difference of 27 MW between the values of the limiting power flow both determined by the static aperiodic stability and calculated
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
using the weighting trajectory is explained in terms of the difference between the initial calculation models used for mode weighting. Testing the algorithm on the controlled section No. 1 of the Irkutsk power system led to a conclusion concerning its applicability in searching for an adaptive weighting trajectory with subsequent application in calculating and analysing ultimate active-power flows in the controlled section in terms of static aperiodic stability.
Keywords: weighting trajectory, Jacobian matrix determinant, aperiodic stability, large-scale power system, sensor nodes, controlled section
Information about the article: Received December 18, 2019; accepted for publication May 29, 2020; available online June 30, 2020.
For citation: Batseva NL, Sukhorukov VA. Search algorithm for adaptive weighting trajectory in large-scale power systems. Vestnik Irkutskogo gosudarstvennogo tehnicheskogo universiteta = Proceedings of Irkutsk State Technical University. 2020;24(3):582-595. (In Russ.) https://doi.org/10.21285/1814-3520-2020-3-582-595
1. ВВЕДЕНИЕ
Траектория утяжеления представляет собой набор изменений значений параметров режима энергосистемы (ЭС) на каждом шаге утяжеления, которые приводят к увеличению перетока активной мощности в контролируемом сечении (КС)1. В соответствии с принятой для определения величин допустимых перетоков мощности траекторией утяжеления должен быть получен наименьший предельный по статической апериодической устойчивости (САУ) переток активной мощности в КС, поэтому возникает задача поиска кратчайшей и «наиболее опасной» траектории утяжеления, которая будет приводить исследуемый режим ЭС к границе устойчивости.
В практических расчетах применяются два типа траекторий утяжеления -«генерация - генерация» и «генерация -нагрузка». Первый тип характерен для частей ЭС, в которых имеется достаточная мощность генерации, при этом поиск траектории утяжеления сводится к решению задачи оптимальной загрузки (или разгрузки) генераторов для получения наименьшего предельного перетока в КС. Второй тип используется, когда в передающей части ЭС в качестве изменяемого режимного параметра используется мощность генерации, а в приемной части ЭС - мощность нагрузки. Проблема поиска траектории решается путем определения опорных нагрузочных уз-
лов, в которых следует увеличить нагрузку. Опорными являются узлы, где прогнозируется увеличение мощности нагрузки или подключение новых потребителей в соответствии со схемой и программой развития (СиПР) субъекта. Выполнение задач СиПР не всегда может быть гарантировано в полном объеме, поэтому появляется неопределенность в правильности найденной траектории утяжеления. Кроме того, режимные параметры постоянно изменяются, поэтому для получения значений допустимых перетоков, соответствующих текущему режиму ЭС, необходимо искать траекторию утяжеления с учетом актуальных значений параметров режима, т.е. адаптировать ее к изменяющимся режимным параметрам.
Вопросы поиска адаптивной траектории утяжеления рассмотрены в работах [1-10]. Автор работы [1] предлагает использовать модифицированные уравнения предельных режимов для того, чтобы найти кратчайшую нормаль к границе устойчивости. Для ЭС, которая может быть замещена многомашинным эквивалентом, граница устойчивости представляется в п-мерном пространстве, и поэтому существует множество нормалей к границе устойчивости. В работах [2, 3] описан алгоритм определения кратчайшего расстояния до границы области устойчивости режима ЭС путем минимизации целевой функции. Целевая функция состоит из сумм квадратов невязок уравнений, описывающих установив-
1
СТО 59012820.27.010.001 - 2013. Правила определения максимально допустимых и аварийно допустимых перетоков активной мощности в контролируемых сечениях диспетчерского центра ОАО «СО ЕЭС». Введ. приказом ОАО «СО ЕЭС» № 10 от 18.01.2013.
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):682-Б9Б
шийся режим ЭС. В первом случае они задаются вектором регулируемых параметров режима, а во втором - узлами утяжеления, мощности которых будут изменяться. Вектор регулируемых параметров режима представляет собой совокупность узлов расчетной модели, в которых изменяют мощности генерации и нагрузки для достижения границы области САУ. При задании траектории утяжеления типа «генерация -нагрузка» и выборе вектора регулируемых параметров режима вначале необходимо выбрать опорные узлы. От выбранного состава опорных узлов будет зависеть расчетный предельный переток активной мощности. Однако перед выполнением расчетов точный состав опорных узлов неизвестен, поэтому поиск критического направления утяжеления затруднен. Выбор опорных узлов на основании экспертных оценок неоднозначен и приведет к различным значениям предельных перетоков при разных совокупностях опорных узлов, входящих в вектор регулируемых параметров режима. В работе [2] сделан акцент на то, что для реальных схем ЭС вопрос поиска направления утяжеления является достаточно сложной задачей, и чаще всего направление утяжеления выбирается с учетом знаний персонала о режимах работы ЭС и опыта выполнения расчетов, что не позволяет применять данный подход в качестве универсального.
Для того чтобы избежать многомерности задачи утяжеления режима, авторы работы [4] предлагают использовать деформацию области САУ равномерным образом, изменяя параметры схемы сети в ^ раз. Данный подход планируется применять и при выборе управляющих воздействий (УВ) с использованием алгоритма I-ДО в централизованной системе противо-аварийной автоматики (ЦСПА). С его помощью можно загрузить все КС схемы, а затем найти наиболее загруженное (опасное) сечение и относительно определенного опасного сечения выбирать УВ. Однако
использование такого подхода при определении предельных величин по САУ перетоков активной мощности и назначении максимально допустимых (МДП) и аварийно допустимых (АДП) перетоков мощности в КС затруднительно, поскольку при выполнении расчетов необходимо добиваться нарушения устойчивости именно в исследуемом КС, а не в смежных с ним сечениях.
В работах [5-7] авторы предлагают оригинальный подход определения величины МДП в КС с использованием текущей телеметрической информации на основе методов оценивания состояния ЭС. Идея подхода состоит в задании псевдоизмерения величины МДП (ПИ МДП), которое заведомо превышает существующий МДП, и далее в определении МДП при изменении заданных параметров режима в пределах регулирования. Регулируемыми параметрами режима выступают активные и реактивные мощности, а также модули напряжения в узлах ЭЭС. Определение МДП происходит путем максимального приближения перетока активной мощности в КС к величине ПИ МДП путем минимизации целевой функции взвешенных наименьших квадратов.
В работе [8] рассматривается поиск сенсорных или чувствительных узлов. Режимные параметры (напряжение и угол при напряжении) сенсорных или чувствительных узлов изменяются в большей степени при изменениях мощности нагрузки и генерации. Сенсорные узлы предполагается использовать в качестве опорных узлов в адаптивной траектории утяжеления режима. При этом не рассматривается, какие именно узлы следует отнести к сенсорным после применения алгоритмов поиска (с использованием обратной матрицы Якоби, спектрального или сингулярного разложения матрицы Якоби) и их ранжирования.
В работе2 описывается общая тенденция снижения границы области САУ при увеличении количества узлов, которые входят в траекторию утяжеления. Данный
Гайсин Б.М. Разработка методик определения влияния неоднородности электроэнергетических систем на возникновение и развитие аварийных каскадных процессов: дис. ... канд. техн. наук: 05.14.02. Уфа, 2019. 156 с.
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
факт выводит задачу отделения наиболее сенсорных узлов от несенсорных узлов на достаточно важный уровень при поиске адаптивной траектории утяжеления.
В работах [9, 10] авторы приводят алгоритм поиска адаптивной траектории утяжеления, который основывается на степени изменения величины определителя матрицы Якоби при различных изменениях мощности в узлах расчетной модели ЭС. В предлагаемом алгоритме при изменении активной мощности в каком-либо узле ЭС реактивная мощность остается неизменной. Однако, как следует из1, изменение реактивной мощности необходимо учитывать при определении величины предельного по САУ перетока активной мощности в КС. Также в рассматриваемом алгоритме не учитывается изменение мощности генерации при изменении мощности нагрузки, т.е. используется несбалансированная траектория утяжеления. Использование такой траектории утяжеления может приводить к изменению баланса мощности балансирующего узла, частных производных матрицы Якоби в балансирующем и соседних с ним узлах, и в конечном итоге - к неверному значению изменения определителя матрицы Якоби. Помимо этого, при таком подходе предел по САУ может достигаться не в исследуемом КС, а в других сечениях расчетной модели.
В работах [1-3, 8-10] предлагаемые алгоритмы преимущественно апробируются на тестовых ЭС. В расчетную модель реальной ЭС входят тысячи узлов и ветвей, в связи с чем размерность матрицы Якоби для схемы реальной ЭС будет большой, поэтому возникает проблема анализа определителя матрицы Якоби для ЭС большой размерности.
Проблема поиска адаптивной траектории утяжеления связана в том числе с внедрением и функционированием в Единой энергетической системе (ЕЭС) России системы мониторинга запасов устойчивости (СМЗУ), которая предназначена для определения в режиме реального времени величин МДП и АДП в КС [11]. На данном этапе развития технологии в СМЗУ приме-
няется алгоритм пересчета нормировочных коэффициентов вектора изменения режима (ВИР) [12]. Данный алгоритм позволяет на каждом цикле расчета изменять параметры утяжеления, т.е. приращения мощностей в заранее определенных узлах утяжеления, учитывая текущие значения мощностей нагрузки и генерации, а также технологические ограничения генераторов электрических станций. Использование адаптивной траектории утяжеления сделало бы возможным не только определение параметров узлов утяжеления на каждом расчетном цикле СМЗУ, но и адаптивный выбор узлов утяжеления с учетом текущей схем-но-режимной ситуации.
Таким образом, проблема поиска адаптивной траектории утяжеления с определением опорных узлов для схем большой размерности до сих пор остается актуальной.
В статье ставится задача алгоритмизации процесса поиска оптимальной траектории утяжеления в текущем режиме ЭС, которая будет являться адаптивной к сложившейся схемно-режимной ситуации.
Целью исследования является поиск адаптивной траектории утяжеления для текущей схемно-режимной ситуации ЭС, соответствующей наименьшему предельному по САУ перетоку активной мощности.
2. МАТЕРИАЛЫ И МЕТОДЫ ИССЛЕДОВАНИЯ
Классическим подходом к анализу САУ является известный метод малых колебаний. Корни характеристического уравнения позволяют судить об устойчивом или неустойчивом состоянии ЭС. Появление хотя бы одного корня с положительной вещественной частью будет свидетельствовать об апериодическом нарушении устойчивости и смене знака свободного члена характеристического уравнения. При наличии в расчетной схеме балансирующего узла, рассматриваемого как шины бесконечной мощности, задании для генераторов иг = const, использовании PU-модели генератора, применении в расчетах устой-
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):Б82-Б9Б
чивости и установившихся режимов одних и тех же статических характеристик нагрузки, значения свободного члена характеристического уравнения и определителя матрицы Якоби совпадают между собой, что дает возможность использовать определитель матрицы Якоби в качестве критерия САУ.
Матрица Якоби представляет собой квадратную матрицу размерностью
(2« "(«баз + «Узд )"«баз ) [13]. Структура матрицы Якоби и размерность ее блоков приведена на рис. 1.
На рис. 1 использованы обозначения: п - общее число узлов в расчетной модели ЭС, «баз - число базисно-балансирующих узлов, «Гзд - число PU-
узлов (узлы, моделирующие генераторы электростанций, а также плавно регулируемые средства компенсации реактивной мощности, например, синхронный компенсатор, статический тиристорный компенсатор, управляемый шунтирующий реактор).
При расчете определителя матрицы большой размерности использовалось LU -разложение матрицы Якоби с учетом слабой заполненности матрицы.
П - Пбааз П - «баз- «Узд
S
I
'"ч £ s
<
M / 1
4 1 dP/dô ' dP/dU
s 1
<- _ +-----
! с 1
s 1
1 iS dQ/dô 1 dQ/dU
s 1
1 1
s
I
S <N
^- -
2П - («баз + Пу-д) - «баз
Рис. 1. Структура матрицы Якоби Fig. 1. Jacobian matrix structure
В специализированных математических пакетах имеется функция sparse, которая позволяет хранить слабозаполненные матрицы с помощью покоординатного хра-
нения, когда хранятся только ненулевые элементы матрицы и соответствующие строчные и столбцовые индексы [14].
LU - разложение матрицы А имеет следующий вид [15]:
A = LU,
(1)
где A - исходная матрица; L - нижняя унитреугольная матрица (на главной диагонали расположены единицы, все элементы, расположенные выше главной диагонали, равны нулю); U - верхняя треугольная матрица (все элементы, расположенные ниже главной диагонали, равны нулю).
Получив LU - разложение матрицы, можно вычислить ее определитель как произведение диагональных элементов матрицы U , т.е.
det ( A) = det (LU) = det (L ) det (U ) = = det (U ) = un ...unn.
Блок-схема алгоритма поиска адаптивной траектории утяжеления в части выбора нагрузочных узлов, изменение мощности в которых наиболее сильно влияет на определитель матрицы Якоби, приведена на рис. 2 с учетом сведений, представленных в работах [9, 10].
Алгоритм поиска адаптивной траектории утяжеления состоит в следующем:
1. Рассчитывается исходное значение определителя матрицы Якоби detJ для рассматриваемой схемно-режимной ситуации в следующей последовательности:
1.1 выгружается массив значений производных матрицы Якоби (дР/дб, dP/dU, dQ/d5, dQ/dU) из ПК RastrWin в пакет Microsoft Excel;
1.2 выполняется макрос на языке VBA, результатом которого являются три массива, реализующие покоординатное хранение матрицы Якоби;
1.3 выполняется LU - разложение матрицы Якоби и определяется значение ее определителя путем перемножения диагональных элементов верхней треугольной матрицы U .
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
Начало
1. Расчет исходного значения определителя матрицы Якоби det J для рассматриваемой схемно-режимной ситуации
>
2. Увеличение мощности нагрузки в i - м рассматриваемом узле энергорайона, принимающего мощность по КС, на величину ДPШi+jДQШi
1 г
3. Увеличение загрузки по активной мощности на одном из генераторов в части ЭС, передающей мощность по КС, на величину ДРг
г
4. Расчет нового значения определителя матрицы Якоби det J
г
5. Определение значения отношения (det J - det Jt )/АРш
нет 6. Рассмотрены все нагрузочные
узлы энергорайош?-
да
7. Сортировка нагрузочных узлов по убыванию значения отношения (det J - det J )/АРн
8. Выбор нагрузочных узлов, у которых значение отношения J - det Ji )/ДРш максимально (метод Аг-средних, метод иерархической агломеративной кластеризации, критерии _Стьюдента ^ и Фишера)_
Конец
Рис. 2. Блок-схема алгоритма поиска адаптивной траектории утяжеления Fig. 2. Flow diagram of the search algorithm of an adaptive heavy-loaded trajectory
2. В узле части ЭС, принимающей мощность по КС, увеличивается значение активной мощности нагрузки на величину ДРН:
P
ДР-~ =-P^P. баз , (2)
где ДРН1 - приращение активной мощности нагрузки в /-м узле; РН!. - исходная активная мощность нагрузки /-го узла; Рнбаз - базовое
значение активной мощности нагрузки (вып
бирается произвольно); ^р - суммарная
г=1
активная мощность нагрузки всех рассматриваемых нагрузочных узлов; п - количество рассматриваемых узлов.
Пропорционально значению ДРН (с
сохранением коэффициента мощности) увеличивается значение реактивной мощности нагрузки на величину ДQн:
(3)
где - приращение реактивной мощности нагрузки в /-м узле; ДРН1 - приращение активной мощности нагрузки в /-м узле; - коэффициент мощности.
3. На генераторе в части ЭС, передающей мощность по КС, увеличивается загрузка по активной мощности на величину ДРг:
ДР = ДР, + Р,
(4)
У Р
/ , ш
где ДР - приращение активной мощности генерации на одном из генераторов в передающей части ЭС; ДРШ. - приращение активной мощности нагрузки в /-м узле; Ршт -
величина потерь активной мощности в электросетевых элементах (принимается равной 1-2 МВт).
4. Рассчитывается новое значение йвг^ после увеличения нагрузки в /-м узле
на величину ДРн1. и увеличения загрузки одного из генераторов в передающей части ЭС на величину ДРГ.
5. Определяется значение отношения
i=1
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):Б82-Б9Б
A det J det J - det J,
AP
AP
(5)
где Айви - изменение значения определителя матрицы Якоби; АРШ - приращение активной мощности нагрузки в ш узле; detJ - исходное значение определителя матрицы Якоби; detJ - значение определителя матрицы Якоби после увеличения нагрузки в ш узле на величину АРШ. и увеличения загрузки одного из генераторов в передающей части ЭС на величину Ар.
6. Процесс повторяется в той же последовательности для других исследуемых узлов.
7. Рассматриваемые нагрузочные узлы сортируются по убыванию значения
А detJ
отношения
AP
8. Выбираются узлы, для которых
А detJ
значение отношения -является макАР
симальным.
Нагрузочные узлы разного класса напряжения рассматриваются отдельно, а узлы с активной мощностью нагрузки равной нулю, и узлы, в которых реактивная мощность нагрузки больше или равна активной, исключаются из рассмотрения, т.к. увеличение мощности в таких узлах приведет к высоким значениям отношения АdetJ/АРШ, выпадающим из общей выборки узлов.
Для выбора нагрузочных узлов, значения отношений АdetJ/АРШ которых являются максимальными, необходимо провести границу между двумя выборками, взяв при этом значение отношения АdetJ/ АРШ, выступаемого в качестве граничного: все нагрузочные узлы с отношением АdetJ/ Ар, большим, чем граничное, попадут в первую выборку, а со значением меньшим, чем граничное, - во вторую. В связи с этим возникает задача определения значения граничного отношения, от величины которого будет зависеть,
какое количество нагрузочных узлов попадет в первую выборку, сколько будут являться сенсорными и станут опорными узлами.
Для разбиения исходной выборки в качестве начального приближения были использованы методы кластерного анализа с указанием требуемого числа кластеров (два кластера): метод ^средних (k-means) и метод иерархической агломеративной кластеризации с использованием метода Уорда в качестве метрики расстояния между кластерами.
Метод ^средних является простым и распространенным методом кластеризации, используемым в совокупности с другими математическими методами для решения различных технологических задач в ЭС [16-19]. Суть его заключается в минимизации суммарного квадратичного отклонения точек кластера от центроида этого кластера:
min
11 К7 ^
i=1 x«EÄ
(6)
где к - количество кластеров; х(') - j-е наблюдение; - -й кластер; д - центроид кластера ^.
Метод ^средних по умолчанию использует Евклидово расстояние в ^мерном пространстве в качестве метрики расстояния между точками кластера:
p( Х,У H Iх - У II=
n 2
Ж-УP), (7)
p=i
x
где р(х,у) - Евклидово расстояние; хр, ур - рассматриваемые точки кластера; п -
размерность пространства кластера.
Одним из недостатков алгоритма ^ средних является необходимость задания требуемого числа кластеров. Однако в нашем случае оно известно заранее: первый кластер - сенсорные нагрузочные узлы, которые войдут в адаптивную траекто-
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
рию утяжеления и будут опорными, второй кластер - все остальные нагрузочные узлы. Качество проведенной кластеризации во многом зависит от выбранных начальных центроидов. Для отыскания центроидов кластеров использован метод k-meaпs++ [20, 21].
Метод иерархической агломератив-ной кластеризации состоит в рассмотрении исходных точек наблюдений в качестве отдельных кластеров. Затем эти кластеры объединяются в более крупные, а те в свою очередь в еще более крупные, и т.д. Таким образом, происходит иерархическое формирование кластеров. В качестве метрики расстояния между кластерами использован метод Уорда, который на каждом шаге кластеризации объединяет два кластера так, что это приводит к минимальному увеличению дисперсии [22]:
Д=Е(х -х )2"
I 2 _ 2 (8)
-Е(х-а) -Е(х-ъ) >
х. е Л х. е В
где Д - метрика Уорда; х - '-е наблюдение; Л - первый объединяемый кластер; В - второй объединяемый кластер; х -
центр объединенного кластера; а - центр кластера А; ъ - центр кластера В.
Для подтверждения правильности результатов кластеризации использованы критерии Стьюдента и Фишера, которые позволяют проверить гипотезу об однородности двух выборок путем проверки равенства средних значений и дисперсий в двух выборках [23].
3. РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ И ИХ ОБСУЖДЕНИЕ
Для подтверждения работоспособности алгоритм поиска адаптивной траектории утяжеления был апробирован на КС № 1. Схема ЭС с КС № 1 приведена на рис. 3. При проведении расчетов рассматривалась схема ремонта устройства продольной компенсации (УПК) 500 кВ (УПК 500 кВ зашунтировано).
КС № 1 состоит из сетевых элементов:
- ВЛ 500 кВ ПС 500 кВ № 3 - ПС 500 кВ № 6;
- ВЛ 500 кВ ГЭС № 2 - ПС 500 кВ № 5 - первая цепь;
- ВЛ 500 кВ ГЭС № 2 - ПС 500 кВ № 5 - вторая цепь.
ГЭС №1
и-
ГЭС №2
_ Рассматриваемый энергорайон
УПК 500 кВ
У
I ПС 500 кВ №5
ПС 110 кВ №5
кВ пс 110 Пело
кВ №9 кВ №10
ТЭС №5 ТЭС №6
Рис. 3. Схема энергосистемы с контролируемым сечением № 1 Fig. 3. Diagram of a power system with a controlled section no. 1
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):Б82-Б9Б
Сенсорные нагрузочные узлы Sensor load nodes
^асс номинального Общее количество ^личество
напряжения, кВ нагрузочных узлов сенсорных узлов
6-10 6 3
110 114 14
220 8 6
Часть ЭС, обведенная на схеме пунктирной линией, принимает мощность по КС № 1, при этом нагрузочные узлы сгруппированы по классам номинального напряжения: 6-10 кВ, 110 кВ, 220 кВ. Мощность генерации изменялась на генераторах ГЭС № 2. В таблице приведено общее количество исследуемых узлов и количество тех узлов, изменение мощности в которых наиболее сильно влияет на изменение определителя матрицы Якоби.
Например, для нагрузочного узла № 1738 (1 СШ 110 кВ ПС 110 кВ № 14):
1. Определим исходное значение определителя матрицы. Значение определителя составляет detJ = 2,695947 -1СГ229 по пунктам 1.1-1.3 алгоритма. Отметим, что малое значение показателя степени обусловлено выполнением деления каждого диагонального элемента верхней треугольной матрицы U на 10000 для исключения переполнения переменной типа double, представляющей собой значение определителя матрицы Якоби. Поскольку далее в расчетах используется не само значение определителя матрицы Якоби, а его изменение, то эта операция не влияет на качество расчетов.
2. Увеличим значение активной мощности нагрузки узла № 1738 на величину APH по формуле (2):
P
\р * н. 1738 р
^ н.1738 n 1 н.баз
У p
/ , Hi
7,825 145G
• 31,825 = G, 172 МВт.
Пропорционально значению АР (с
сохранением коэффициента мощности) увеличиваем значение реактивной (р)
мощности нагрузки на величину АQШ по (3):
А£ш.1738 =АРШ.1738 љР=
= 0,172 • 0,72 = 0,124 МВ • Ар.
3. На одном из генераторов ГЭС № 2 увеличим загрузку по активной мощности на величину Ар,, формула (4):
Ар = АР + P =
АР г.ГЭС№2 АР н.1738 + р пот
= G, 172 + 1 = 1,172 МВт.
4. Определим новое значение опре-
делителя detJ17з8 после увеличения нагрузки в узле на величину АРн1738 и увеличения загрузки одного из генераторов ГЭС № 2 на величину Аргэс№2. Значение определителя составляет
с1гих1Ъ% =2,577738 • 10 229.
5. Вычислим значение отношения АdetJ / АР,
н.1738
по формуле (5):
А det J Q229 _ det J det J1738 jq229 _
АР ~ АР =
АР н.1738 АР н.1738
_( 2,695947 - 2,577738) • 1 G" = G ,172
-229
• 1G 229 = G, 687.
Значение
отношения
AdetJ/Арл738 домножено на 1G229 для
удобства дальнейших расчетов.
Аналогичным образом алгоритм применялся к остальным нагрузочным узлам. Для разделения полученной выборки узлов на два кластера по методам ^ средних и иерархической агломеративной кластеризации использованы формулы (6)-(8). Методы выделили одни и те же нагру-
590
ISSN 1814-3520
i=1
зочные узлы в первый кластер - кластер сенсорных узлов. С использованием критериев Стьюдента и Фишера было доказано, что полученные две выборки в результате процесса кластеризации действительно являются статистически различными.
На рис. 4 и 5 приведены численные значения отношений Др для вы-
борок узлов с классами напряжений 220 кВ и 110 кВ, соответственно.
На рис. 4 и 5 темно-серым цветом выделены значения отношений
ДdetJ/ДРн. узлов, которые попали в первый кластер, являются сенсорными и должны быть использованы в адаптивной траектории утяжеления как опорные узлы, светло-серым цветом - значения отношений ДdetJ/ДРн. нагрузочных узлов, которые попали во второй кластер. На рис. 4 приведены не все номера нагрузочных узлов 110 кВ в связи с их большим количеством.
2,5
<U
о
^ 1,5
<
х 1 5 1
к ч <и и
0,5
х
<и £
к
СП
X
к
к о ч . <и
Ц 5а
11 £ ^ ас <3
СП
0,6 0,5 0,4 0,3 0,2 0,1 0
1767 1693 1742 1744 1694 1750 1751 1746 Номера узлов в расчетной модели ЭС
Рис. 4. Выборка узлов напряжением 220 кВ Fig. 4. Sample of 220 kV nodes
,-н о oo i/ч oo m oo h ю h t о ю h h h oo h о
0
oo
7
4 0 0 2
oo чо vo m о (N
0 6 5 1 0
^T i/o m oo
0 4
(N oo
69
4 0
m m vo (N On О О ю ю oo ю h in
Номера узлов в расчетной модели ЭС
Рис. 5. Выборка узлов напряжением 110 кВ Fig. 5. Sample of 110 kV nodes
ISSN 1814-3Б20
Б91
2
0
При утяжелении режима по найденной траектории предел по САУ составил 2089 МВт. Предел по САУ с использованием траектории утяжеления, включающей в себя нагрузочные узлы, в которых увеличение мощности прогнозируется в соответствии с СиПР, составил 2072 МВт. Предельный переток с использованием адаптивной траектории утяжеления на 27 МВт превышает предельный переток, полученный с использованием траектории в соответствии с СиПР. Разница между значениями предельных перетоков не превышает амплитуду нерегулярных колебаний в 100 МВт. Разница в значениях предельных перетоков объясняется различием исходных расчетных моделей, которые использованы при утяжелении режима.
4. ЗАКЛЮЧЕНИЕ
На основе выполненного исследования можно сделать вывод о корректности созданного алгоритма в части поиска адап-
1. Крюков А.В. Предельные режимы электроэнергетических систем. Иркутск: Изд-во ИрГУПС, 2012. 236 с.
2. Паздерин А.В., Чусовитин П.В., Шабалин Г.С., Юферев С.В. Определение запасов устойчивости и управляющих воздействий для обеспечения статической устойчивости в задаче противоаварийного управления на основе Обобщенного метода Ньютона // Современные направления развития систем релейной защиты и автоматики энергосистем: сб. науч. тр. V Междунар. науч.-техн. конф. (г. Сочи, 1-5 июня 2015 г.). Сочи, 2015. Т. 2. С. 50-57.
3. Shabalin G., Pazderin A., Bannykh P., Balakh E. Voltage stability analysis using quadratic objective function taking into account equality constraints // International Conference on the Science of Electrical Engineering (Eilat, 16-18 November 2016). New Jersey: IEEE, 2017. P. 1-5. https://doi.org/10.1109/ICSEE.2016.7806166
4. Аржанников С.Г., Вторушин А.С., Захаркин О.В., Ландман А.К., Петров А.Э., Попова Е.Ю. Алгоритмическое обеспечение ПТК верхнего уровня ЦСПА ОЭС Сибири и перспективы развития // Известия НТЦ Единой энергетической системы. 2013. № 1. С. 91-98.
5. Аксаева Е.С., Глазунова А.М. Оценивание максимально допустимого перетока контролируемого сечения электроэнергетической системы // Известия
тивной траектории утяжеления с последующим ее применением при определении предельных по САУ перетоков активной мощности в КС. Для получения более достоверных результатов и накопления статистических данных предложенный алгоритм планируется апробировать на других ремонтных схемах для КС № 1, а также для других КС, в том числе и для тех, где в качестве изменяемого режимного параметра используется мощность генерации в передающей и приемной частях ЭС. Кроме того, при удовлетворении требований к быстродействию алгоритм может быть апробирован на каждом цикле расчета в СМЗУ для получения адаптивной траектории утяжеления.
В статье не были исследованы смежные КС по отношению к КС № 1, не был рассмотрен вопрос о возможном нарушении устойчивости в смежных КС. Но в настоящее время алгоритм поиска адаптивной траектории утяжеления совершенствуется для решения этих задач.
шй список
НТЦ Единой энергетической системы. 2018. № 1. С. 23-34.
6. Глазунова А.М. Применение метода оценивания состояния для вычисления предельных режимов электроэнергетической системы // Электричество. 2013. № 10. С. 23-29.
7. Glazunova A.M., Aksaeva E.S. Available transfer capability determination on the basis of a trade-off approach // Manchester PowerTech (Manchester, 18-22 June 2017). New Jersey: IEEE, 2017. P. 1-6. https://doi.org/10.1109/PTC.2017.7981014
8. Войтов О.Н., Воропай Н.И., Гамм А.З., Голуб И.И., Ефимов Д.Н. Анализ неоднородностей электроэнергетических систем. Новосибирск: Наука, Сиб. изд. фирма РАН, 1999. 256 с.
9. Ефремова И.Ю., Глускин И.З. Разработка структурной схемы адаптивного пускового органа автоматики предотвращения нарушения устойчивости // Вестник Ивановского государственного энергетического университета. 2017. № 6. С. 15-24.
10. Ефремова И.Ю., Глускин И.З. Адаптивная настройка пускового органа противоаварийной автоматики для транзитов с промежуточными отборами // Электричество. 2017. № 2. С. 13-17. http://dx.doi.org/10.24160/0013-5380-2017-2-13-17
11. Томалев А.А. Система мониторинга запасов устойчивости (СМЗУ). Опыт эксплуатации и перспективные направления развития СМЗУ в ОЭС
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
Сибири // Электроэнергетика глазами молодежи: сб. науч. тр. юбилейной X Междунар. науч.-техн. конф. (г. Иркутск, 16-20 сентября 2019 г.). Иркутск: Изд-во ИРНИТУ, 2019. Т. 2. С. 83-86.
12. Александров А.С., Максименко Д.М., Михайленко А.Ф., Неуймин В.Г. Развитие системы мониторинга запасов устойчивости с контролем динамической устойчивости для учета действия противоаварийной автоматики // Известия НТЦ Единой энергетической системы. 2017. № 1. С. 64-72.
13. Bonini Neto A., Amancio Alves D. Singularities Analysis of the Jacobian Matrix Modified in the Continuation Power Flow: Mathematical Modeling // IEEE Latin America Transactions. 2016. Vol. 14. Issue 12. P. 4750-4756. https://doi.org/10.1109/TLA.2016.7817006
14. Shah M. Sparse Matrix Sparse Vector Multiplication - A Novel Approach // 44th International Conference on Parallel Processing Workshops (Beijing, 1-4 September 2015). New Jersey: IEEE, 2015. P. 67-73. https://doi.org/10.1109/ICPPW.2015.18
15. Razik L., Schumacher L., Monti A., Guironnet A., Bureau G. A comparative analysis of LU decomposition methods for power system simulations // Milan PowerTech (Milan, 23-27 June 2019). New Jersey: IEEE, 2019. P. 1-6. https://doi.org/10.1109/PTC.2019.8810616
16. Lami B., Bhattacharya K. Clustering Technique Applied to Nodal Reliability Indices for Optimal Planning of Energy Resources // IEEE Transactions on Power Systems. 2016. Vol. 31. Issue 6. P. 4679-4690. https://doi.org/10.1109/TPWRS.2015.2507061
17. Kim Do-In, White A., Shin Yong-June. PMU-Based Event Localization Technique for Wide-Area Power System // IEEE Transactions on Power Systems. 2018. Vol. 33. Issue 6. P. 5875-5883.
https://doi.org/10.1109/TPWRS.2018.2824851
18. Grigoras G., Neagu B.-C., Scarlatache F., Ciobanu R.C. Identification of pilot nodes for secondary voltage control using K-means clustering algorithm // 26th International Symposium on Industrial Electronics (Edinburgh, 19-21 June 2017). New Jersey: IEEE, 2017. P. 106-110. https://doi.org/10.1109/ISIE.2017.8001231
19. Yasrebi Nayini S.E., Geravand S., Maroosi A. A novel threshold-based clustering method to solve k-means weaknesses // International Conference on Energy, Communication, Data Analytics and Soft Computing (Chennai, 1-2 August 2017). New Jersey: IEEE, 2018. P. 47-52. https://doi.org/10.1109/ICECDS.2017.8389496
20. Olukanmi P.O., Twala B. K-means-sharp: Modified centroid update for outlier-robust k-means clustering // Pattern Recognition Association of South Africa and Robotics and Mechatronics (Bloemfontein, 30 November - 1 December 2017). New Jersey: IEEE, 2018. P. 14-19. https://doi.org/10.1109/RoboMech.2017.8261116
21. Niu Kun, Gao Zhipeng, Jiao Haizhen, Deng Nanjie. K-means+: A developed clustering algorithm for big data // 4th International Conference on Cloud Computing and Intelligence Systems (Beijing, 17-19 August 2016). New Jersey: IEEE, 2016. P. 141 -144. https://doi.org/10.1109/CCIS.2016.7790241
22. Puggina Bianchesi N.M., Romao E.L., Lopes M.F.B.P., Balestrassi P.P., De Paiva A. P. A Design of Experiments Comparative Study on Clustering Methods // IEEE Access. 2019. Vol. 7. P. 167726-167738. https://doi.org/10.1109/ACCESS.2019.2953528
23. Pruim R. Foundations and Applications of Statistics. An Introduction Using R. Second Edition. American Mathematical Society, 2018. 820 p.
References
1. Kryukov AV. Limit Load Operation of Electric Power Systems. Irkutsk: Irkutsk State Transport University; 2012, 236 p. (In Russ.)
2. Pazderin AV, Chusovitin PV, Shabalin GS, Yuferev SV. Determination of Stability Margins and Control Actions to Ensure Steady-State Stability in the Emergency Control Problem using Newton's Method. In: Sov-remennye napravleniya razvitiya sistem releinoi zash-chity i avtomatiki energosistem: Trudy V mezhdunarod-noy nauuchno-tekhnicheskoy konferentsii = Actual Development Trends of Power System Relay Protection and Automation: Proceedings of 5th International Scientific and Technical Conference. 01-04 June 2015, Sochi. Sochi; 2015, vol. 2, p. 50-57. (In Russ.)
3. Shabalin G, Pazderin A, Bannykh P, Balakh E. Voltage Stability Analysis using Quadratic Objective Function Taking into Account Equality Constraints. In: International Conference on the Science of Electrical Engineering. 16-18 November 2016, Eilat, New Jersey: IEEE; 2017, p. 1-5. https://doi.org/10.1109/ICSEE.2016.7806166
4. Arzhanikov SG, Vtorushin AS, Zakharkin OV, Land-
man AK, Petrov AE, Popova EYu. Algorithmic Support of the Upper Level Device of the Sibirean Centralized Emergency Control System and its Development Prospects. Izvestiya NTC Edinoj Enegreticheskoj Sistemi = STC of Unified Power System Proceedings. 2013;1:91 — 98. (In Russ.)
5. Aksaeva ES, Glazunova AM. Estimation of Total Transfer Capability for a Controlled Cutset of Electric Power System. Izvestiya NTC Edinoj Enegreticheskoj Sistemi = STC of Unified Power System Proceedings. 2018;1:23-34. (In Russ.)
6. Glazunova AM. Using a State Estimation Method for Calculating Limit Load Operation Conditions of an Electric Power System. Elektrichestvo. 2013;10:23-29. (In Russ.)
7. Glazunova AM, Aksaeva ES. Available Transfer Capability Determination on the Basis of a Trade-Off Approach. In: Manchester PowerTech. 18-22 June 2017, Manchester. New Jersey: IEEE; 2017, p. 1-6. https://doi.org/10.1109/PTC.2017.7981014
8. Voitov ON, Voropai NI, Gamm AZ, Golub II, Efimov DN. Analysis of Electrical Power System Heterogenei-
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
ties. Novosibirsk: Nauka. The Siberian Publishing Company of the Russian Academy of Sciences; 1999, 256 p. (In Russ.)
9. Efremova IU, Glouskin IZ. Development of an Adaptive Fault Detector Structural Scheme of Automatic Stability Control System. Vestnik Ivanovskogo Gosudar-stvennogo Energeticheskogo Universiteta = Vestnik of Ivanovo State Power Engineering University. 2017;6:15-24. (In Russ.)
10. Efremova IYu, Gluskin IZ. Adaptive Tuning of the Emergency Control System for Transmission Systems with Intermediate Power Takeoff Points. Elektrichestvo. 2017;1:13-17. (In Russ.) http://dx.doi.org/10.24160/0013-5380-2017-2-13-17
11. Tomalev AA. System Stability Reserves Monitoring (SMZU). Using Experience and Development Prospects of SMZU in UES of Siberia. In: Elektroenergetika glazami molodezhi: trudy yubileinoi X mezhdunarodnoy nauuchno-tekhnicheskoy konferentsii = Power Industry through the Eyes of Young People: Proceedings of Anniversary X International Scientific and Technical Conference. 16-20 September 2019, Irkutsk. Irkutsk: Irkutsk National Research Technical University; 2019, vol. 2, p. 83-86. (In Russ.)
12. Alexandrov AS, Maximenko DM, Mikhaylenko AF, Neuymin VG. Evolution of Stability Margin Monitoring System with Control of Dynamic Stability for Taking into Account Emergency Automatic Actions. Izvestiya NTC Edinoj Enegreticheskoj Sistemi = STC of Unified Power System Proceedings. 2017;1:64-72. (In Russ.)
13. Bonini Neto A, Amancio Alves D. Singularities Analysis of the Jacobian Matrix Modified in the Continuation Power Flow: Mathematical Modeling. IEEE Latin America Transactions. 2016;14(12):4750-4756. https://doi.org/10.1109/TLA.2016.7817006
14. Shah M. Sparse Matrix Sparse Vector Multiplication - A Novel Approach. In: 44th International Conference on Parallel Processing Workshops. 1-4 September 2015, Beijing. New Jersey: IEEE; 2015, p. 67-73. https://doi.org/10.1109/ICPPW.2015.18
15. Razik L, Schumacher L, Monti A, Guironnet A, Bureau G. A Comparative Analysis of LU Decomposition Methods for Power System Simulations. In: Milan PowerTech. 23-27 June 2019, Milan. New Jersey: IEEE; 2019, p. 1-6.
Критерии авторства
Бацева Н.Л., Сухоруков В.А. заявляют о равном участии в получении и оформлении научных результатов и в равной мере несут ответственность за плагиат.
Конфликт интересов
Авторы заявляют об отсутствии конфликта интересов.
Все авторы прочитали и одобрили окончательный вариант рукописи.
https://doi.org/10.1109/PTC.2019.8810616
16. Lami B, Bhattacharya K. Clustering Technique Applied to Nodal Reliability Indices for Optimal Planning of Energy Resources. IEEE Transactions on Power Systems. 2016;31 (6):4679—4690. https://doi.org/10.1109/TPWRS.2015.2507061
17. Kim Do-In, White A, Shin Yong-June. PMU-Based Event Localization Technique for Wide-Area Power System. IEEE Transactions on Power Systems. 2018;33(6):5875—5883. https://doi.org/10.1109/TPWRS.2018.2824851
18. Grigoras G, Neagu B-C, Scarlatache F, Ciobanu RC. Identification of Pilot Nodes for Secondary Voltage Control using K-Means Clustering Algorithm. In: 26th International Symposium on Industrial Electronics. 1921 June 2017, Edinburgh. New Jersey: IEEE; 2017, p. 106-110. https://doi.org/10.1109/ISIE.2017.8001231
19. Yasrebi Nayini SE, Geravand S, Maroosi A. A Novel Threshold-Based Clustering Method to Solve K-Means Weaknesses. In: International Conference on Energy, Communication, Data Analytics and Soft Computing. 12 August 2017, Chennai. New Jersey: IEEE; 2018, p. 47-52. https://doi.org/10.1109/ICECDS.2017.8389496
20. Olukanmi PO, Twala B. K-Means-Sharp: Modified Centroid Update for Outlier-Robust K-Means Clustering. In: Pattern Recognition Association of South Africa and Robotics and Mechatronics. 30 November - 1 December 2017, Bloemfontein. New Jersey: IEEE; 2018, p. 14-19.
https://doi.org/10.1109/RoboMech.2017.8261116
21. Niu Kun, Gao Zhipeng, Jiao Haizhen, Deng Nanjie. K-Means+: A Developed Clustering Algorithm for Big Data. In: 4th International Conference on Cloud Computing and Intelligence Systems. 17-19 August 2016, Beijing. New Jersey: IEEE; 2016, p. 141-144. https://doi.org/10.1109/CCIS.2016.7790241
22. Puggina Bianchesi NM, Româo EL, Lopes MFBP, Balestrassi PP, De Paiva AP. A Design of Experiments Comparative Study on Clustering Methods. IEEE Access. 2019;7:167726-167738. https://doi.org/10.1109/ACCESS.2019.2953528
23. Pruim R. Foundations and Applications of Statistics. An Introduction Using R. Second Edition. American Mathematical Society; 2018, 820 p.
Authorship criteria
Batseva N.L., Sukhorukov V.A. declare equal participation in obtaining and formalization of scientific results and bear equal responsibility for plagiarism.
Conflict of interests
The authors declare that there is no conflict of interests regarding the publication of this article.
The final manuscript has been read and approved by all the co-authors.
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):Б82-Б9Б
СВЕДЕНИЯ ОБ АВТОРАХ
Бацева Наталья Ленмировна,
кандидат технических наук, доцент, доцент отделения электроэнергетики и электротехники Инженерной школы энергетики, Национальный исследовательский Томский политехнический университет, 634050, г. Томск, пр. Ленина, 30, Россия; Н e-mail: batsevan@tpu.ru
Сухоруков Василий Александрович,
магистрант,
Национальный исследовательский Томский политехнический университет, 634050, г. Томск, пр. Ленина, 30, Россия; e-mail: vas60005596@yandex.ru
INFORMATION ABOUT THE AUTHORS
Natalia L. Batseva,
Cand. Sci. (Eng.),
Associate Professor of the Department of Electric Power and Electrical Engineering, School of Energy and Power Engineering, National Research Tomsk Polytechnic University, 30, Lenin Ave, Tomsk 634050, Russia; H e-mail: batsevan@tpu.ru
Vasiliy A. Sukhorukov,
Master Degree Student,
National Research Tomsk Polytechnic University, 30, Lenin Ave, Tomsk 634050, Russia; e-mail: vas60005596@yandex.ru
ВЕСТНИК ИРКУТСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА 2020;24(3):582-595
