CC BY
43
УДК 519.683.8
РАСЧЕТ УСТАНОВИВШЕГОСЯ РЕЖИМА ЭЛЕКТРИЧЕСКОЙ СЕТИ
СРЕДСТВАМИ GPU
Ф.Н. ЯСИНСКИЙ, Н.Б. ИЛЬИЧЕВ, А.И. КУЛЕШОВ, Д.П ХАРИТОНОВ ФГБОУВПО «Ивановский государственный энергетический университет им. В. И. Ленина»,
Иваново, Российская Федерация E-mail: [email protected]
Авторское резюме
Состояние вопроса: На данный момент задача получения установившегося режима является решенной и результат получается за конечное время. Однако при больших объемах сети и массовых расчетах установившегося режима необходимо сократить время вычисления.
Материалы и методы: Для расчета установившегося режима используется метод узловых потенциалов. При решении системы нелинейных уравнений методом Ньютона выполняется линеаризация уравнений посредством определения матрицы Якоби и решается система линеаризованных уравнений при помощи GPU. Результаты: Рассматривается алгоритм расчета установившегося режима электросети. Описывается общая структура видеокарты. Приводится анализ распараллеливания приведенного алгоритма. Выводы: Показано, что можно получить ускорение расчета установившегося режима электросети.
Ключевые слова: установившийся режим, GPU, параллельный алгоритм.
CALCULATION OF STEADY MODE OF ELECTRICAL NETWORK WITH MEANS OF GPU
F.N. YASINSKIY, N.B. IL'ICHEV, A.I. KULESHOV, D.P. KHARITONOV Ivanovo State Power Engineering University, Ivanovo, Russian Federation E-mail: [email protected]
Abstract
Background: At present the problem of getting the steady mode is solved, and the result can be received for the short period of time. However, for big volumes of networks and mass calculation of steady mode it is necessary to decrease the calculation time.
Materials and methods: The nodal-voltage method is used for calculation of steady mode. The authors carry out the linearization of equations by means of finding the matrix of Jacobi while solving the system of nonlineal equations with the Newton method. The linearization of equations helps to solve linearized equations with GPU.
Results: The article considers the calculation algorithm of steady mode of electrical network. The description of the graphical card general structure is given. The analysis of the given algorithm paralleling is carried out. Conclusions: It is shown that it is possible to get calculation acceleration of steady mode of electrical network.
Key words: steady mode, GPU, parallel algorithm.
Расчеты установившихся режимов в электроэнергетических системах являются наиболее массовыми расчетами. Усложнение расчетных моделей, связанное с более детальным представлением сети, значительно увеличивает размерность задачи и затраты времени на получение решения. В то же время современная вычислительная техника имеет неиспользованный в настоящее время потенциал производительности вычислительных систем, связанный с применением многоядерных процессоров и наличием матричных процессоров в составе видеокарт. То, что данный ресурс пока не используется, обусловлено тем, что распараллеливание процесса расчета установившегося режима связано с разработкой специальных вычислительных алгоритмов. Известно, что расчет установившегося режима, как правило, состоит в решении системы нелинейных уравнений узловых напряжений, которые могут быть записаны в форме баланса мощностей или в форме ба-
ланса токов. При этом матрица коэффициентов левой части уравнения является слабозапол-ненной. Таким образом, речь идет о распараллеливании процесса решения системы нелинейных уравнений.
Архитектура GPU и особенности работы с устройством. Для расчета установившегося режима (УР) в качестве входных данных используется матрица собственных и взаимных проводимостей. Эта матрица является разреженной (слабозаполненной), что накладывает определенные ограничения на выбор математических методов, применяемых для расчета УР. При составлении матрицы собственных и взаимных проводимостей учитываются различные типы ветвей: трансформаторы, линии, генераторы и пр.
Для начала рассмотрим оборудование, с которым придется работать, учтем особенности его работы, которые могут пригодиться нам для эффективного расчета УР.
Тенденция последних лет к увеличению числа полупроводниковых элементов на устройстве, которые можно было бы использовать в вычислениях, привела к созданию GPGPU (General-Purpose Graphics Processing Units) -техники использования графического процессора видеокарты для выполнения расчётов в приложениях для общих вычислений, которые обычно проводит центральный процессор. На рис. 1 представлены схематичные изображения архитектуры центрального и графического процессоров, демонстрирующие разницу этих двух устройств.
Рис. 1. Схематичное изображение архитектуры процессоров: а - центрального; б - графического
Графический процессор обладает большим количеством ядер, но следует отметить, что ядра GPU предназначены только для вычислений и имеют меньшее количество инструкций по сравнению с CPU. Соответственно, исходя из различной организации архитектур устройств логично предположить, что наличие меньшего числа ядер с большим количеством инструкций рационально для использования параллелизма по задачам (MIMT - Multiple Instructions Multiple Threads), а использование устройства с большим количеством параллельно работающих ядер оптимально для распараллеливания по данным (SIMT - Single Instruction Multiple Threads).
Для понимания особенностей реализации алгоритмов для GPU необходимо иметь представление об архитектуре устройства с точки зрения логической компоновки вычислительных элементов. На рис. 2 представлена архитектура GPU в терминах технологии CUDA (альтернатива CUDA - OpenCL - имеет очень похожую архитектуру, только используется другая терминология, но сути программирования под GPU с использованием OpenCL это не меняет).
Пото< р. 1 I ____ Поим p. q
Рис. 2. Архитектура GPU CUDA
Также стоит отметить наличие различных типов памяти в GPU-устройстве (рис. 3).
Рис. 3. Модель памяти GPU
Различные типы памяти имеют различные скорости доступа. Так, например, самая медленная память - глобальная (global), далее идут локальная (local) и частная (private). Но и объемы у различных типов памяти разные. Так, глобальная память - это вся оперативная память компьютера, как правило, достигает объема в несколько гигабайт. А локальная и частная память намного меньше - как правило, несколько десятков килобайт.
Алгоритм расчета УР. Теперь рассмотрим основные этапы расчета установившегося режима:
1) получение матрицы узловых прово-димостей на вход (процесс формирования этой матрицы можно найти в [1]);
2) задание начальных приближений напряжений узлов (искомых переменных) (в качестве первого приближения берутся номинальные напряжения узлов, которые в дальнейшем корректируются на последующих итерациях выбранного математического метода);
3) определение величины небалансов мощностей;
4) формирование матрицы Якоби;
5) решение системы линейных уравнений;
6) определение новых значений напряжений и фазовых углов напряжений;
7) проверка условий окончания расчета и возврат на п.3, если это необходимо.
Остановимся поподробнее на некоторых из перечисленных выше пунктах.
Так, уравнения небалансов мощностей имеют следующий вид:
Ф,(и,8) = и2 д,, + £ и, и] х
М ] (1)
х [д,у 003(8 - 8]) + Ьу 3т(8, -8]) ] - Р = 0; т
у,(и,8) = -и,2Ь„ +У и,и, х
и ] (2)
х [д(] з1п(8,- - 8]) - Ьу 003(8, - 8]) ] - О,- = 0,
где к - число узлов в сети без учета балансирующего узла и узлов с заданным модулем напряжения.
Формирование матрицы Якоби производится по следующим формулам:
—Р т
др = X и ; иу Ь 003(8, - 8у) - д ] з!п(8, - 8у) ]; (3)
др = -и, иу [Ьу 003(8, - 8у) - д у з!п(8, - 8у) ]; (4) —Р т
-и = 2Ц д,, + Хиу [ду 003(8 - 8у) + Ьу 3'п(8 - 8у)]; ди, у=1
(5)
—р
-и = и, [ду 003(8 -8у) + Ьу з1п(8 -8у) ]; (6)
—О т
-О- = X и, иу [ду 003(8, - 8у) + Ьу 31п(8, - 8у) ]; (7)
= -X и, иу [ду 003(8, - 8у) + Ьу 31п(8, - 8у)]; (8)
—О,
—8
у у=1
—О т
■О = 2и ,Ь„ -Хиу [Ьу 003(5,-5у) - ду 3^(8,-8 у)]; —и, у=1
(9)
—о
—и
= -X и у [Ьу 003(8, - 8 у) - д у 31п(8, - 8 у) ] . ( 10)
у у=1
—Р —Р
Объединяя элементы вида—'-; —'- в
—8, —8у
—Р
—р —Р
в мат-
матрицу—, элементы вида ,
—8 —и, —иу
—Р
рицу— и аналогичным образом формируя
—и
—О —О
матрицы — и — , уравнение баланса мощ-
—8 —и
ностей в блочно-матричной форме будет выглядеть следующим образом:
—V ]дик = _дРк -1. '^Ц.]д8к + ^ = -Д Ок-1.
(11)
(12)
Следующим шагом является решение системы линейных уравнений, полученных в
результате выполнения предыдущего шага. Решив полученную систему уравнений, получаем новые значения приращения напряжения ди и фазовых углов Д8. После этого можно получить очередное приближение этих двух параметров, воспользовавшись следующими формулами:
ик=ик-+дик-1; 8к =8к-1 + д8к-1.
Получив новые значения, проверяем на сходимость (проверка на сходимость включает в себя больше, чем просто проверку на достижение значений небаланса желаемой точности. Более подробно этот шаг описан в [1, с. 70]) и в случае, если желаемая точность не достигнута, переходим к следующей итерации, повторив вышеупомянутые действия, начиная с шага пересчета небаланса мощностей.
На рис. 4 приводится схема сети, для которой выполняется расчет установившегося режима. Подробное описание формирования исходных данных для этой сети и пошаговый алгоритм расчета для однопроцессорного варианта можно найти в [1].
Рис. 4. Примерная схема для расчета установившегося режима
Распараллеливание алгоритма расчета УР. Выше был приведен алгоритм выполнения расчета УР электроэнергетической системы, представляющий собой алгоритм решения системы нелинейных уравнений. Согласно этому алгоритму, линеаризация уравнений выполняется посредством получения матрицы Якоби. Таким образом, для решения системы нелинейных уравнений необходимо решить систему линейных уравнений.
В качестве алгоритма решения системы нелинейных уравнений будем использовать метод Ньютона [7, 8].
Ниже рассмотрим, как можно распараллелить расчет матрицы Якоби и подробно остановимся на решении системы линейных уравнений.
Согласно формулам (3)-(10), можно выделить 8 значений, которые необходимо получить. Причем каждое из этих значений, во-первых, не зависит от других значений, во-вторых, если оперировать матричной терминологией, не зависит от значений по столбцам, «накапливая» данные по строкам. Таким образом, само собой напрашивается распараллеливание данных по строкам, предоставив каж-
дому процессу (нити, потоку - в различных русскоязычных источниках используется разная терминология для термина thread) считать свою строку, разграничивая пересчет диагональных и недиагональных элементов. То есть каждый процесс будет считать одну строку, во время которой будет пересчитано 8 значений. Посчитаем теоретическое ускорение, которое можно получить на данном этапе. Для простоты вычислений примем, что не используются структуры данных для хранения разреженной матрицы (т.е. хранятся и нулевые элементы). Тогда, если в строке содержится N значений, то пересчет одной строки матрицы Якоби составит 8N. А пересчет всей матрицы Якоби при однопроцессорной обработке составит 8-N-N. При использовании N процессов получаем, что каждый процесс считает 8 N значений. Теоретически получаемое ускорение по количеству выполняемых операций будет рассчитываться как
N.
8 • N2 8 • N
Это довольно грубый расчет, не учитывающий время на передачу между различными типами памяти (см. рис. 3.) и различное время обращения к памяти на CPU и GPU (в зависимости от типа памяти). Целью подсчета ускорения по операциям было показать наличие ускорения как такового на данном этапе.
Однако вся выгода от ускорения на этапе получения матрицы Якоби может быть «съедена» неправильным выбором метода расчета системы линейных уравнений и частыми передачами данных с ядра видеокарты в оперативную память (об этом будет сказано ниже).
Ключевым моментом, таким образом, становится выбор расчета системы линейных уравнений. Существует множество способов решения системы линейных уравнений, которые можно подразделить на 2 больших класса -прямые и итерационные математические методы. К прямым методам относятся:
• метод Гаусса;
• метод Жордана-Гаусса;
• метод Крамера;
• разложение Холецкого;
• метод Томаса.
К итерационным методам относятся:
• метод Якоби (метод простой итерации);
• метод Гаусса - Зейделя;
• метод релаксации;
• многосеточный метод;
• метод Монтанте;
• метод Абрамова;
• метод обобщенных минимальных невязок;
• метод бисопряженных градиентов;
• стабилизированный метод бисопряженных градиентов;
• квадратичным метод сопряженных градиентов;
• метод квази-минимальных невязок.
Обоснование использования того или
иного метода представляет собой отдельную задачу, на решении которой мы останавливаться не будем. Описание некоторых этих методов можно найти в [3, 4, 5].
На данный момент в системе ЕпегдуС8 (работа выполняется для модификации расчетного блока программного комплекса) используется метод Гаусса, дающий точное решение за N шагов. Однако реализация этого метода на видеокарте является нетривиальной задачей, так как в процессе решения системы данным методом матрица коэффициентов раскладывается на верхнюю треугольную и нижнюю треугольную матрицы (Ш-разложение, подробное описание этого метода см. в [4, 6]). При Ш-разложении имеет место зависимость по данным, то есть расчет следующего значения разложенной матрицы невозможен до тех пор, пока не будут посчитаны предыдущие значения. Следовательно, толку от распараллеливания разложения не будет, так как каждый процесс будет ждать пересчета других значений, превратив параллельный расчет Ш-разложения в последовательный. Одним из вариантов выхода из данной ситуации является использование особого алгоритма расчета Ш-разложения, позволяющего параллельно проводить эту операцию. Информация о параллельном Ш-разложении представлена в [6].
Другим вариантом выхода из этой ситуации является использование итерационного метода. Однако тут необходимо учитывать следующее. При использовании разреженных матриц используются особые структуры для их хранения, что ведет к тому, что, например, операция транспонирования матрицы представляет собой нетривиальную и довольно трудоемкую задачу. Другой момент - некоторые итерационные математические методы требуют предварительного расчета собственных чисел матрицы, что ведет к дополнительным затратам машинного времени. И наконец, итерационные методы могут быть нестабильны . Стабилизация того или иного метода также является большой областью исследований в численных методах и зачастую носит очень конкретный характер, варьируясь от задачи к задаче.
Ниже приведем два алгоритма итерационных методов: метода простых итераций (как самого простого и понятного) и метода бисоп-ряженных градиентов стабилизированного (БЮС81аЬ - метод, по которому в русскоязычной литературе мало информации). Недостатком этих методов является то, что они предназначены лишь для положительно определенных матриц и не всегда сходятся.
Основная формула для метода простых итераций (метода Якоби):
¡-1 m
X = ^ -^а]Х] -^аЦХ] , (13)
]=1 i=+1
и Ь аа
где Ь =-*-; а] =-].
аа ап
Согласно (13), данный метод очень хорошо подходит под парадигму Б1МТ. Каждый процесс считает свою строку, собирая в самом конце невязку в общую переменную. Процесс идет до тех пор, пока отношение суммарной невязки к количеству задействованных процессов не достигнет необходимой точности.
Данный метод прост в реализации, прост в понимании, по нему довольно много информации в русскоязычной литературе, однако он довольно требователен к исходным данным и долго сходится.
Алгоритм метода ВЮСБ1аЬ следующий:
к = -1,
еыбрэтъ -Ро " г0?
посчитать Го = — Ах q .
u_l = = 0, p_t = L
пс"0то|жтъ пока 11 11| не достаточно мало
k-k+i.
i>k ('';■■ 'У). Рк = -pk/pk-i'
Чк = 1Ъ - ihUk-1- <Ъ ~ -4(11 ■
"к = fk - PkUk-l,
7 к = (<-k-fk)- at = pkhk-
■*'jt+l = Я* + ''(.■+! =Tk- OfcCfc, rk+i =f к -akATuk
Данный метод, являющийся модификацией метода сопряженных градиентов (CG -Conjugate Gradient), отличается повышенной стабильностью, однако все равно не всегда сходится. Способы повышения стабильности данного метода можно найти в [7]. В основе данного алгоритма лежит одна из базовых векторных операций - saxpy: saxpy (a, x, y) = ax + y,
где a - число; x, y- векторы.
Распараллеливание данной операции очень просто и не требует особых трудозатрат. Однако, ввиду того что метод BiCGStab требует синхронизации во время расчета для получения промежуточных значений, это замедляет скорость его расчета.
Если сравнивать методы Якоби и BiCGStab, то второй метод сходится за меньшее количество итераций, однако выгода от его использования становится очевидна лишь при больших размерах матрицы, так как при малых размерах матрицы метод Якоби, не требующий синхронизации значений во время выполнения, оказывается быстрее.
Стоит обратить внимание на метод итераций Чебышева [2], не требующий межпроцессорной синхронизации во время выполнения. Однако недостатком данного метода является необходимость знать спектр матрицы (максимальное и минимальное собственные числа) перед началом расчета.
Заключение
Алгоритм расчета установившегося режима электрической сети и предложенные подходы к распараллеливанию этого расчета позволяют сделать вывод, что выбор численного метода остается за пользователем и осуществляется в зависимости от используемых структур для хранения данных и возможностей параллельной обработки данных на используемой машине. Рекомендации, касающиеся выбора математического метода с учетом особенностей работы с вычислительным блоком на видеокарте, даны с учетом того, что выбор видеокарты для расчетов обусловлен наличием большого количества процессоров (технология легких ядер) на борту.
Список литературы
1. Кулешов А.И., Прахин Б.Я. Расчет и анализ установившихся режимов электроэнергетических систем на персональных компьютерах: учеб. пособие / Иван. гос. энерг. ун-т. - Иваново, 2001. - 171 с.
2. Чадов С.Н. Некоторые вопросы численного моделирования динамических систем / ГОУВПО «Ивановский государственный энергетический университет им. В.И. Ленина». - Иваново, 2010. - 120 с.
3. Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. - Новосибирск: Изд-во НГТУ, 2000. - 70 с.
4. Голуб Дж., Ван Лоун Ч. Матричные вычисления: пер. с англ. - М.: Мир, 1999. - 548 с.
5. BiCGstab(L) for linear equations involving unsym-metric matrices with complex spectrum, gerard l.g. sleijpen, diederik r. fokkema, Electronic Transactions on Numerical Analysis. - September 1993. - Vol. 1. - Р. 11-32.
6. Семушин И.В. Численные методы алгебры: учебное пособие для вузов. -Ульяновск: УлГУ, 2008. 178 с.
7. Radolph E. Bank, Tony F. Chan. An analysis of the composite step biconjugate gradient method // Numerische mathematic. - Vol. 66. - № 1. - Р. 295-319.
8. http://ru.wikipedia.org/wiki/Метод_Ньютона
9. http://en.wikipedia.org/wiki/Newton's_method
References
1. Kuleshov, A.I., Prakhin, B.Ya. Raschet i analiz ustanovivshikhsya rezhimov elektroenergeticheskikh sistem na personal'nykh komp'yuterakh [Calculation and Analysis of Steady Modes of Electric Power Systems with using of PCs], Ivanovo, 2001, 171 p.
2. Chadov, S.N. Nekotorye voprosy chislennogo modeli-rovaniya dinamicheskikh sistem [Specific Issues of Numerical Simulation in Dynamic Systems], Ivanovo, 2010, 120 p.
3. Balandin, M.Yu., Shurina, E.P. Metody resheniya SLAU bol'shoy razmernosti [Solving Methods of SLAU with great Dimention], Novosibirsk: Izd-vo NGTU, 2000, 70 p.
4. Golub, J., Van Loan, Ch. Matrichnye vychisleniya [Matrix Computations], Moscow: Mir, 1999, 548 p.
5. BiCGstab(L) for linear equations involving unsymmet-ric matrices with complex spectrum, gerard l.g. sleijpen, diederik r. fokkema, Electronic Transactions on Numerical Analysis, vol. 1, pp. 11-32, September 1993.
6. Semushin, I.V. Chislennye metody algebry [Computational Algebra Methods], Ul'yanovsk: UlGU, 2008, 178 р.
7. Radolph, E. Bank, Tony, F. Chan. An analysis of the composite step biconjugate gradient method. Numerische mathematic, vol. 66, 1, pp. 295-319.
8. http://ru.wikipedia.org/wi й/Метод_Ньюто на
9. http://en.wikipedia.org/wiki/Newton's_method
Ясинский Федор Николаевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
доктор физико-математических наук, профессор кафедры высокопроизводительных вычислительных систем,
телефон (4932) 26-98-29.
Ильичев Николай Борисович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
кандидат технических наук, доцент кафедры электрических станций, подстанций и диагностики электрооборудования, e-mail: [email protected]
Кулешов Анатолий Иванович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
кандидат технических наук, доцент кафедры электрических станций, подстанций и диагностики электрооборудования, e-mail: [email protected]
Харитонов Денис Петрович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
аспирант кафедры высокопроизводительных вычислительных систем,
телефон (4932) 26-98-29,
e-mail: [email protected]
