МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ ПАРАЛЛЕЛЬНЫХ РАСЧЕТОВ СЛОЖНЫХ ТЕХНИЧЕСКИХ СИСТЕМ НА ОСНОВЕ ПРИМЕНЕНИЯ ТЕНЗОРНОЙ ТЕОРИИ СЕТЕЙ

(г. Псков, Псковский государственный политехнический институт)

Одно из направлений реализации тензорной теории сетей – повышение эффективности вычислений. В статье рассматривается математическое обеспечение тензорной диакоптики для ее применения в параллельных расчетах. Основное внимание уделено математическим моделям тензорной теории сетей и методу расчета моделей по частям. Приводятся результаты тестовых расчетов для одной из наиболее доступных в настоящее время технологий – кластерной технологии.

Одно из направлений реализации тензорной теории сетей – повышение эффективности вычислений. В 1955 году была опубликована работа Г. Крона [1], в которой сообщалось об обращении матрицы 256х256 методом тензорной диакоптики на вычислительной машине CPC (Card Programme Calculator), построенной IBM. Так как компьютер СРС мог обращать матрицы с максимальным размером 16х16 (за 1час 10мин), то матрица была разделена на 16 частей. Общий расчет занял около 16 часов.

Развитие параллельных расчетов определяется появившимся в последнее время разнообразным техническим обеспечением для их проведения. В настоящее время наиболее известны кластерные решения, grid-расчеты, расчеты на многоядерных архитектурах (CPU и GPU), расчеты на кристаллах с программируемой архитектурой (ПЛИС), нейрокомпьютерах. Можно предположить и дальнейший рост архитектурного разнообразия параллельного компьютерного «железа».

НЕ нашли? Не то? Что вы ищете?

В статье рассматривается математическое обеспечение тензорной диакоптики [2] для ее применения в параллельных расчетах. Основное внимание уделено математическим моделям тензорной теории сетей и методу расчета моделей по частям. Приводятся результаты тестовых расчетов для одной из наиболее доступных в настоящее время технологий – кластерной технологии.

Математические модели. Основные составляющие математического моделирования включают [3] метод электрических аналогий для систем различной физической природы и построение резистивных сеток для разностной аппроксимации дифференциальных уравнений в частных производных и обыкновенных дифференциальных уравнений. Цель метода электроаналогий для систем различной физической природы – выделить два вида двойственных переменных, описывающих техническую систему: потоковые переменные и продольные переменные (в тензорной терминологии контрвариантные и ковариантные переменные). В электрической системе к ним относят ток i и напряжение U, их механическими аналогами могут быть сила и скорость (момент вращения и угловая скорость), тепловыми двойственными переменными будут тепловой поток и температура и т. п. Далее предполагается, что двойственные переменные связаны компонентными дифференциальными уравнениями (обыкновенными и в частных производных). В результате разностной аппроксимации дифференциальных уравнений, строятся алгебраические уравнения, которым ставится в соответствие резистивная сетка. Различные методы разностной аппроксимации могут быть выражены в форме уравнения обобщенной электрической ветви - одного из фундаментальных понятий тензорного анализа сетей. Одновременно методы расчета, алгоритмически «упаковываясь» в обобщенную ветвь, представляют базовый набор ветвей, наиболее употребительных в электротехнических расчетах: линейные и нелинейные RLC ветви, взаимные индуктивности и управляемые источники. На рис.1 приведена схема обобщенной ветви.

 

Рис. 1. Обобщенная ветвь U+e = z ( i + J )

Приведем пример построения обобщенной ветви для обыкновенного дифференциального уравнения, связывающего двойственные переменные для электрической емкости C: i=C ∙ dU/dt . После разностной аппроксимации по неявному методу Эйлера 1-го порядка получаем:

,

где индекс (n+1) соответствуют текущему шагу, n - предыдущему расчетному шагу. Емкость можно представить резистивной ветвью, соответствующей рис.1 с сопротивлением z=∆t/С, с последовательно включенным источником напряжения, величина которого равна напряжению на емкости в предыдущий момент времени и взятому с обратным знаком e=-Un и источником тока J=0.

Для нелинейных обыкновенных дифференциальных уравнений так же может быть выполнена алгоритмическая «адаптация» к обобщенной ветви. Приведем пример для метода простых итераций для емкостной нелинейной ветви. Ток в емкости описывается выражением: ic = dQ/dt, где Q − заряд на емкости. Заряд на емкости при численных расчетах можно представить через статическую емкость, зависящую от напряжения: Q=C(Uc) ∙ Uc. При численных расчетах производную некоторой функции можно представить приближенно как сумму этих функций, вычисленных в предыдущие моменты времени. Например, для метода обратного дифференцирования 2-го порядка для производной заряда в текущий, (n+1)-й момент времени:

где ∆t − расчетный шаг; Qn+1 – значение заряда при t=tn+1, Qn – значение заряда при t=tn, Qn-1 – значение заряда при t=tn−1. Коэффициенты α0=−1,5, α1=2, α2=−0,5. В результате получаем разностно-итерационное уравнение:

,

которое соответствует обобщенной ветви рис.1, если ввести обозначения:

, , J=0.

Так как в текущий (n+1)-й момент времени емкость неизвестна, то ее значение итерационно уточняется от начального значения, известного на предыдущем (к-1)-м шаге.

Более подробно различные методы построения разностно-итерационных резистивных ветвей изложены в [4].

Приведем пример построения резистивной сетки для уравнения в частных производных. Пусть распределение магнитной индукции В вдоль оси х описывается уравнением:

В данном случае уравнение описывает распределение магнитной индукции вдоль движущейся со скоростью V проводящей пластины с удельным сопротивлением ρ, толщиной а, расположенной в зазоре δ’ стального магнитопровода. Магнитное поле создается токовым слоем с заданной линейной плотностью J1.

Аппроксимация на основе центральных разностей с шагом h дает следующее выражение:

Для перехода к сеточной электрической модели разностное уравнение преобразуется к виду:

где приняты следующие обозначения:

B+=B(x+h, t); B=B(x, t); B-=B(x-h, t); J+=J(x+h, t); J-=J(x-h, t);

Разностному уравнению можно поставить в соответствие следующую электрическую схему:

Рис.2. Схемный аналог для центральных разностей.
 
 

Как видно на рис.2, резистивная сетка для решения приведенного уравнения в частных производных будет состоять из Т-образных подсхем, в которые входят 2 проводимости, одна емкость, один источник тока, управляемый напряжением и один независимый источник тока. Расчет полной схемы можно выполнить как во временной области для анализа переходных процессов, так и в частотной – для анализа установившихся режимов. Все ветви подсхемы рис.2 могут быть представлены в виде обобщенной ветви рис.1.

Таким образом, модель сложной технической системы можно представить в виде резистивной схемы замещения, состоящей из обобщенных ветвей. Весь набор обобщенных ветвей представляет собой так называемую элементарную схему, описываемую уравнением в векторно-матричной форме: U+e = z ( i + J ), где z – в общем случае несимметричная матрица сопротивлений ветвей. Вектор токов i можно трактовать как вектор, задающий систему координат пространства отдельных ветвей элементарной схемы. Любая схема, собранная из отдельных ветвей, образует новое пространство с контурными и узловыми координатами. Вектор тока i' схемы в этом пространстве задает систему координат, которая связана с координатами отдельных ветвей соединенной схемы:

i = C · i'

где С – матрица преобразований пространств. В данном случае матрица является топологической, в ее столбцах записывается информация о ветвях, входящих в контуры. В тензорном анализе сетей С – тензор преобразований, который может быть представлен в матричной форме. Матрица С в этом случае всегда квадратная. Так как столбцы матрицы С соответствуют координатам открытых контуров (узловых пар) и замкнутых контуров, то матрица называется контурной. В тензорном анализе сетей показано, что обратная контурная матрица С равна узловой транспонированной матрице А, т. е.:

А= (Сt)-1 или Сt ∙ А = 1.

Равенство единичной матрице произведения Сt ∙ А говорит об ортогональности узлового и контурного пространств. Узловая матрица А связывает вектор напряжений пространства элементарной схемы и вектор напряжений пространства собранной схемы:

U = A ∙ U'.

Уравнения напряжений собранной схемы можно сформировать автоматически, если известен контурный тензор преобразований С и компоненты уравнения элементарной схемы. Форма уравнений напряжений остается инвариантной по отношению к любым преобразованиям элементарной схемы и имеет вид:

U'+e' = z' ( i' + J'),

где z' = Сt ∙ z ∙С – матрица сопротивлений схемы в узловых и контурных координатах,

e' = Сt ∙ е – вектор источников напряжений в координатах соединенной схемы,

J' = Аt ∙ J – вектор источников тока в узловых и контурных координатах соединенной схемы.

Несмотря на одинаковый вид тензорного уравнения для всех собранных схем, в конкретной системе координат матрично-векторные уравнения будут разными, и будут соответствовать той или иной конкретной схеме.

Так как вид уравнения соединенной схемы такой же, как и вид уравнения элементарной схемы, то можно осуществить сборку из отдельных схем того или иного устройства по тем же правилам преобразований координат, которые были использованы при сборке схем из элементарных ветвей. Блочно-иерархический принцип описания технического устройства реализуется правилами преобразований тензорной теории сетей. На рис.3 приведены схемные преобразования, осуществляющие переход от элементарной схемы к обобщенному устройству и далее, к другим уровням.

Рис.3. Схемная интерпретация преобразования систем координат при блочно-иерархическом описании устройств.

Тензорные правила преобразования обратимы (матрицы преобразований квадратные, невырожденные), поэтому от уравнений устройств можно перейти к уравнениям отдельных схем. На рис.4 приведены связи уравнений схем и уравнений преобразований для разных уровней иерархии.

Обратим внимание на то, что декомпозиционные представления о технической системе изначально заложены в тензорной теории сетей: тензоры преобразований координатных систем могут нести дискретную топологическую информацию о соединении частей схемы в целую схему и далее связывать схемы в более сложные иерархии. Это отличает сетевые тензоры от «классических» тензоров: последние ассоциируются с непрерывными преобразованиями.

Рис.4 Связь тензорных уравнений для разных уровней иерархии устройства.

Приведем вид топологических тензоров преобразований для простого, но достаточно общего примера схемы, приведенной на рис.5а. Схема состоит из обобщенных ветвей. На рисунке приведена нумерация и обозначения ряда топологических объектов: цифры в кружочках соответствуют номерам узлов, ветви обозначены символом z и порядковым номером ветви. Направление ветви задано направлением источника напряжения ветви. На рис.5б приведен топологический граф, где утолщенными линиями показано дерево графа, пунктирными линиями показаны открытые контуры. Маршруты «токов» открытых контуров проходят по ветвям дерева в общий узел. Общему узлу присваивается нулевой номер, остальные номера узлов равны номерам координат открытых контуров. Закрытые контуры обозначены символами a, b, и c. Токи в закрытых контурах равны токам в ветвях хорд (на топологическом графе обозначены тонкими линиями).

б)
 

а)
 
 

Рис.5. Пример схемы из обобщенных ветвей (а) и ее топологический граф (б).

Если известен список соединений ветвей схемы (такие списки автоматически формируются в различных графических редакторах электрических схем), то наиболее просто формируется узловой тензор преобразований. Часть тензора, перечисляющая узловые координаты автоматически заполняется по столбцам. Столбец соответствует ветви схемы. Цифра 1 ставится в той строке, которая соответствует началу ветви. Цифра -1 ставится в строке, соответствующей концу ветви. Строка с нулевым узлом исключается из рассмотрения. Другая часть узлового тензора, которая соответствует контурным координатам, заполняется единицами на пересечении имени контура и имени хорды, ток в которой равен контурному току. Списки ветвей хорд и ветвей дерева должны быть предварительно найдены по тому или иному известному алгоритму. Наиболее простой способ построения контурного тензора преобразования – вычисление обратного узлового тензора преобразования. Сформированные топологические матрицы схем хранятся в упакованном виде в отдельных файлах вместе с информацией о ветвях схемы.

Ниже приведены контурный и узловой тензоры преобразований в матричной форме.

1

2

3

a

b

c

z1

1

z2

1

z3

-1

1

-1

C = z4

1

1

-1

z5

-1

-1

1

-1

-1

z6

1

z1

z2

z3

z4

z5

z6

1

1

-1

-1

2

1

-1

-1

3

-1

1

1

At = a

1

b

1

c

1

Компоненты уравнений напряжений схемы формируются автоматически либо непосредственно по приведенным выше матричным соотношениям, либо по индексной информации методами «поэлементного вклада», методами схемоанализа и т. п. Уравнение обобщенной схемы в блочном виде, с выделением блоков для контурных и узловых координат выглядит следующим образом:

o

k

Uo

+

eo

=

o

Zoo

Zok

∙(

0

+

Jo

)

0

ek

k

Zko

Zkk

ik

Jk

где индексами обозначены k - контурные координаты, о - узловые координаты. Из уравнений могут быть исключены либо контурные, либо узловые координаты. Для декомпозиционных расчетов при узловом способе разъединения схем подсхемы должны быть представлены в узловом базисе:

Uo + ĕo = Žoo ∙ Jo,

где ĕo = eo − Zok·(Zkk)-1·ek – вектор узловых источников напряжения,

Žoo = Zoo – Zok (Zkk)-1 Zko – приведенная матрица узловых сопротивлений

Для проведения декомпозиционных расчетов при контурном способе разбиения, из уравнений обобщенной схемы исключаются узловые координаты. Уравнение подсхем в контурном базисе:

ĕk = Zkk · ik,

где ĕk=ek − Zko·Jo − Zkk · Jk – вектор приведенных контурных источников напряжения.

Топологически одна и та же схема может быть представленной как чисто узловая (радиальная) и как чисто контурная. На рис.6 приведены оба варианта эквивалентных схем, соответствующих схеме на рис.5а.

Расчет моделей по частям. При расчете по частям схема замещения технического устройства разделяется на изолированные подсхемы удалением ветвей связей (это могут быть короткозамкнутые перемычки или двухполюсники). Получившиеся отдельные подсхемы приводятся к эквивалентным радиальным подсхемам (рис.7а). Основные вычислительные затраты при расчете параметров эквивалентных радиальных подсхем связаны с расчетом обратных матриц контурных сопротивлений (Zkk)-1. Ветви связи, объединяя подсхемы в граничных узлах, образуют цепь пересечений. Отдельные ветви цепи пересечений приведены на рис.7б. Эти ветви представляют собой аналог элементарной схемы, ветви которой затем объединяются в суммарную цепь рис.7в.

 

Рис.6. Эквивалентные схемы в узловом (а) и контурном (б) базисах.

Этапы расчета суммарной схемы можно представить одной формулой:

Ũo=Uo – Žoō∙ Cōpk ∙ (Сpкō ∙ Žōō ∙ Cōpк )-1 Cpkō ∙ Uō ,

где использованы следующие индексы: ō - координаты граничных узлов подсхем, pk - контурные координаты суммарной цепи, Cōpk – топологическая матрица, связывающая граничные узлы ō и контуры pk суммарной схемы. В приведенном соотношении для простоты изложения сопротивления связей представлены короткозамкнутыми перемычками.

Для расчета по приведенной формуле требуется вектор всех узловых потенциалов изолированных подсхем Uo, подвектор граничных потенциалов изолированных подсхем Uō, топологическая матрица Cōpk, матрица приведенных узловых сопротивлений граничных узлов подсхем Žōō, матрица приведенных узловых сопротивлений подсхем Žoō. Основные вычислительные затраты требуются при вычислении обратной матрицы контурных сопротивлений суммарной схемы (Zpкpк)-1.

 

Рис.7. Схемы замещения при расчете по частям.

а) – радиальные подсхемы, б) – цепь пересечений, в) – суммарная цепь

При декомпозиционных расчетах данные о подсхемах организованы в последовательные массивы. Специальные информационные массивы содержат сведения о повторяемости подсхем по структуре и по параметрам, месте их расположения, количестве элементов и т. д. На рис.8а приведен алгоритм параллельного расчета наиболее ресурсоемкой части вычислений – расчета матриц приведенных сопротивлений подсхем при их радиальном эквивалентировании. Этот этап расчета включает в себя также расчет обратных матриц контурных сопротивлений подсхем.

Вычисления реализованы по кластерной технологии в учебном компьютерном классе, в состав которого вошли: 10 компьютеров Pentium4, объединенные в радиальную локальную сеть Fast Ethernet (соединение узлов через сетевой коммутатор). Программное обеспечение реализовано на фортране-95 (в двух вариантах – с применением компиляторов g95 и Intel) с использованием библиотеки межпроцессорного обмена данными MPICH 2.1.0.7. Расчеты выполнялись для тестовой схемы, состоящей из набора одинаковых по структуре подсхем решетчатого вида. Число ветвей в подсхеме – 1024. Количество различающихся по сопротивлениям подсхем при расчетах варьировалось от 1 до 8. На рис.8б приведены сравнительные результаты расчета для программы, написанной для однопроцессорного компьютера с использованием последовательного алгоритма расчета по частям и для программы с тем же алгоритмом, но с распараллеливанием расчета матриц сопротивлений подсхем.

 

Рис.8. Параллельный алгоритм расчета (а) и оценка его эффективности (б).

Можно отметить, что увеличение количества подсхем в 4 раза, увеличило время счета последовательной программы в 4,6 раза, время счета параллельной программы увеличилось в 1,35 раза для компилятора g95. Для компилятора Intel эти результаты практически такие же, но абсолютное время расчета почти в 1,5 раза меньше.

ЛИТЕРАТУРА

Kron G. Inverting a 256x256 matrix. Engineering (London), v.178, March II, 1955. Исследование сложных систем по частям - диакоптика. – М.: Наука, 19с.

3.  Сохор модели и алгоритмы тензорного анализа сетей. Учебно-методическое пособие. - Псковск. гос. политехн. ин-т. - Псков: Издательство ППИ, 20с.

Ильин автоматизации схемотехнического проектирования. – М.: Энергия, 19с.