Моделирование больших систем на основе методов декомпозиции и компактной

обработки разреженных матриц

В.Н. Гридин,
дир., д.т.н., проф., info@ditc.ras.ru,
ЦИТП РАН , г. Одинцово,
В.И.Анисимов,
проф., д.т.н., проф., vianisimov@inbox.ru,
СПбГЭТУ, г. Санкт-Петербург

Рассматривается методика моделирования больших систем на базе математического описания системы в виде блочно-диагональной окаймленной структуры. Показывается, что реализация такого подхода обеспечивает существенное увеличение эффективности программного обеспечения по затратам требуемой оперативной памяти и по скорости решения задачи моделирования. Для дальнейшего повышения эффективности вычислительного процесса предлагается использовать компактную форму описания разреженных блочных матриц отдельных подсистем, на которые в результате декомпозиции разделяется исходная система. Устанавливается методика построения вычислительного процесса для компактной формы описания разреженных блочных матриц на основе связанных списков и на основе метода индексно-адресных матриц.

The technique of modeling large systems based on the mathematical description of the system in the form of block-diagonal-lined structure. It is shown that the implementation of such an approach provides a significant increase in the efficiency of the software required cost of memory and the speed of solving the problem of modeling. To further increase the efficiency of computing process are encouraged to use compact form block matrices sparse description of the individual subsystems, which are divided as a result of decomposition of the original system. Sets the method of construction of the computational process to form a compact description of the block sparse matrices based on linked lists and on the basis of index-addressable arrays.

Если моделируемая система имеет слабосвязанную структуру, то наиболее эффективным способом организации вычислительных процессов является декомпозиция исходной системы на ряд составляющих подсистем с использованием блочно-диагональной окаймленной структуры. При этом существенно повышается производительность процессов, реализующих моделирование системы, что позволяет более эффективно использовать вычислительные ресурсы. Для дальнейшего повышения эффективности вычислительного процесса целесообразно использовать компактную форму описания разреженных блочных матриц отдельных подсистем, на которые в результате декомпозиции разделяется исходная система. Для этой цели могут быть использованы различные методы сжатия данных, каждый из которых имеет свои достоинства и недостатки. При этом наибольший интерес представляют методы, основанные на связанных списках, и метод индексно-адресных матриц, поскольку эти методы предоставляют возможность не только хранить исходное описание в компактной форме, но и позволяют в случае необходимости изменять формат начального описания в процессе обработки данных. Такая необходимость появляется при решении уравнений вследствие неизбежного появления новых ненулевых элементов в процессе этого решения. Все другие известные методы используют фиксированный формат, что не позволяет непосредственно изменять этот формат в процессе обработки сжатых данных на последующих шагах вычислительного процесса.

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

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

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

Для каждого узла схемы можно ввести узловую переменную, определяющую состояние этого узла и из всех введенных узловых переменных v_p можно образовать вектор узловых потенциалов . В каждой подсистеме можно выделить внутренние переменные x_i и составить для каждой подсистемы вектор . При этом, общий вектор всех переменных системы будет иметь вид , где вектор будет включать в себя узловые переменные узлов, а векторы - переменные k-й подсистемы.

С учетом принятой структуры вектора матричное уравнение всей системы в расширенном базисе узловых потенциалов может быть записано в виде, содержащем окаймленную блочную матрицу, в диагональных блоках которой описываются отдельные подсистемы, на которые разделяется моделируемая на основе декомпозиционного подхода система.

W11					W10	X1		S1
	…				…	…		…
		Wkk			Wk0	Xk		Sk
			…		…	…	+	…	= 0
				Wmm	Wm0	Xm		Sm
W01	…	W0k	…	W0m	W00	X0		S0

Рис.1. Окаймленная блочно-диагональная структура описания системы

Система уравнений для k-й подсистемы согласно структуре k-й блочной строки матрицы системы имеет вид

(1)

Отсюда можно получить выражения для вектора внутренних переменных отдельных подсистем моделируемой схемы

,

Введем обозначения

(2)

Тогда выражение для вектора внутренних переменных моделируемой системы будет иметь следующий вид

, (3)

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

Подставляя сюда вектор внутренних переменных подсистем , определяемый выражением (2), получим

Обозначим

(4)

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

Отсюда получим выражение для расчета узловых переменных узлов разделения

(5)

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

В соответствии с математическим описанием задачи в виде окаймленных матриц, вычислительный процесс решения задачи разделяется на следующие этапы:

1. Для каждой подсистемы на основе введенной информации формируются массивы параметров блочных матриц моделируемой системы , , , , ,

2. Выполняется алгоритм Гаусса – Жордана для каждой k-й подсистемы, в результате чего блочная матрица преобразуется в единичную диагональную матрицу, а блочные матрицы и преобразуются соответственно в блочные матрицы и .

3. Вычисляются матрицы , и на основе выражения (5) проводится расчет вектора узловых переменных.

4. Последовательно для каждой k-й подсистемы выполняется на основании выражения (3) расчет векторов внутренних переменных отдельных подсистем.

Эффективность процедуры моделирования систем на основе декомпозиции может быть оценена коэффициентом экономии оперативной памяти и степенью повышения быстродействия.

Если V_D - объем памяти, необходимой для расчета с использованием декомпозиции, V - объем памяти, требуемый для расчета обычными методами и система разделена на m подсистем, число переменных в каждой подсистеме равно n_k, а число переменных связи составляет n, то коэффициент экономии оперативной памяти α=V/V_D подхода основе декомпозиции может быть определен выражением

Если система разделена на 10 одинаковых подсистем (m=10) и число переменных в каждой подсистеме равно числу переменных связи n_k = n, то получим значение коэффициента экономии оперативной памяти α=6.

При оценке степени повышения быстродействия на основе декомпозиции следует отметить, что время решения системы уравнения определяется числом мультипликативных операций и пропорционально величине nдля метода Гаусса и n_/2 для метода Гаусса-Жордано, где n число решаемых совместных уравнений..

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

,

где b=1/3 или b=1/2

Если моделируемая система разбита на 10 одинаковых подсистем и число переменных в каждой подсистеме равно числу переменных связи, т.е. n_k=n₀, то получим значение степени повышения быстродействия при решении задачи моделирования β = 30.

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

К списковым схемам относится метод Кнута, для реализации которого необходимо создать следующие компактные массивы:

WZ – для значений ненулевых элементов исходной матрицы, которые располагаются в этом массиве по относительному адресу a, имеющему значение порядкового номера элемента в массиве WZ;

WI – для номеров строк ненулевых элементов ;

WJ – для номеров столбцов ненулевых элементов ;

NR – для хранения относительного адреса a следующего ненулевого элемента строки;

NC – для хранения относительного адреса a следующего ненулевого элемента столбца;

ER – для хранения относительного адреса a начального элемента очередной строки;

EC – для хранения относительного адреса a начального элемента очередного столбца.

Достоинством метода Кнута является возможность записи элементов массива WZ в любом порядке, и как следствие, возможность введения в описание дополнительных ненулевых элементов, которые неизбежно появляются при решении уравнений. Следует также отметить, что классический метод Кнута позволяет сканировать исходную матрицу, как по строкам, так и по столбцам, что является также достоинством метода. Вместе с тем некоторым недостатком метода является значительное число массивов, необходимых для хранения информации о ненулевых элементов исходной матрицы.

Если учесть, что массивы WZ, WI, WJ, NR, NC имеют длину m, определяемую числом ненулевых элементов, а массивы ER, EC имеют длину n, определяемую порядком исходной матрицы, то эффективность использования памяти для данных T типа double, может быть определена для метода Кнута выражением:

Из полученного выражения видно, что с уменьшением эффективность метода неограниченно возрастает.

Практическая реализация методов Кнута связана с необходимостью полной перекодировки всех массивов при ведении в описание дополнительных элементов, что связано с большими техническими трудностями в процессе построения программного обеспечения и снижением производительности его дальнейшего функционирования. Возникающие проблемы могут быть устранены путем перехода к объектно-ориентированному подходу описания схемы Кнута на основе существующей в современных языках программирования C# и Java методики работы с коллекциями объектов класса ArrayList предназначенный для поддержки динамических массивов, размеры которых могут изменяться в процессе работы. Класс ArrayList реализует ряд интерфейсов и помимо свойств и методов, определенных в интерфейсах, имеет также и собственные свойства и методы.

При этом для решения задачи построения компактного описания необходимо создание объектов rows и columns класса ArrayList. Каждый элемент объекта rows соответствует определенной строке матрицы и ссылается либо на ненулевой элемент в данной строке, либо указывает в null, если в строке элементов нет, а каждый элемент объекта columns соответствует определенному столбцу матрицы и ссылается либо на ненулевой элемент в данном столбце, либо указывает в null, если в столбце отсутствуют элементы.

Описание элементов в списке задается в классе Element, члены которого val, i, j описывают значение элемента w_ij исходной матрицы и значения, определяющие его координаты i, j, а ссылки nextI, nextJ указывают на следующий элемент в строке и на следующий элемент в столбце соответственно.

class Element

{

public double val;

public int i;

public int j;

public Element nextI;

public Element nextJ;

public Element(double val, int i, int j)

{

this.val = val;

this.i = i;

this.j = j;

}

}

Компактное описание всей разреженной матрицы выполняется в классе Matrix, структура которого имеет следующей вид:

public class Matrix

{ Element el;

ArrayList rows = new ArrayList();

ArrayList columns = new ArrayList();

public Matrix(){}

//Методы включения в список новых элементов

//Методы для доступа к элементам списка

}

Реализация методов включения в список новых элементов и методов для доступа к элементам списка базируется на использовании стандартных методов, определенных в интерфейсах класса ArrayList.

Общая структура, отображающая связь списков с объектами rows и columns, приведена на рис.1.

Рис.1. Связь списков с объектами rows и columns

Эффективность списковой схемы при использовании класса коллекций по сравнении с классической схемой Кнута несколько снижается, так как для каждого объекта класса Element помимо значений val, i, j необходимо также хранить значения cсылок nextI и nextJ. При этом эффективность полной схемы определяется выражением

Метод индексно-адресных матриц основан на введении некоторой целочисленной матрицы A, которая в точности повторяет структуру исходной разреженной матрицы, входящей в уравнение системы моделирования. В качестве элементов индексно-адресная матрица содержит порядковый номер a ненулевых элементов, которые перечисляются в некотором массиве WZ. При этом если порядковый номер ненулевого элемента исходной матрицы равен a, то в индексно-адресную матрицу вводится значение A(i,j)=a.

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

При реализации метода индексно-адресных матриц достаточно провести сканирование индексно-адресной матрицы и выбрать из нее очередной ненулевой элемент. В соответствии с ее численным значением, затем необходимо выбрать соответствующий порядковый элемент из массива WZ. Эффективность метода индексно-адресных матриц в случае, когда для хранения индексно-адресной матрицы A используется тип данных длиной два байта, а для хранения каждого значащего элемента исходной матрицы, используются данные длиной 8 байтов, может быть определена выражением:

=

Из приведенного соотношения видно, что коэффициент экономии памяти для данного метода увеличивается с уменьшением коэффициента разреженности.

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

=

При практической реализации задачи построения программного обеспечения для работы с индексно-адресными матрицами целесообразно применить для ее решения существующую в современных языках программирования C# и Java технологию работы с коллекциями, и использовать класс ArrayList, предназначенный для поддержки динамических массивов, размеры которых могут изменяться в процессе работы. Класс ArrayList реализует интерфейсы ICollection, IList, IEnumerable, ICloneable и помимо свойств и методов, определенных в интерфейсах, имеет также и собственные свойства и методы. Для этого класса необходимо создать два объекта – для параметров схемы (например, WZ) и для задающих источников (например, SZ), ссылки на которые удобно объявить в специально созданном классе глобальных переменных (например, GV). В классе глобальных переменных следует также объявить ссылку типа short на индексно-адресный массив A. Создание самих объектов выполняется в произвольной функции после ввода значений всех переменных для размерности схемы.

При построении программного обеспечения для работы с индексно-адресной матрицей необходимо решать две основные задачи. Первая задача связана с формированием компактного описания системы, вторая – с приведением матрицы к диагональной единичной матрице на основе алгоритма Гаусса-Жордано.

Целью формирования компактного описания системы является занесение в массивы WZ и SZ некоторых значений, определяемых параметрами компонентов и структурой схемных связей. Так, если в компактный массив WZ необходимо занести некоторое значение z, которое в полном описании разреженной комплексной матрицы должно было бы заносится в i–ю строку и j–й столбец этой матрицы, то при занесении в компактный массив WZ такая задача решается путем вызова функции set(i, j, z), которая имеет описание

void set(int i, int j, Complex z)

{ int a = GV.A[i, j];

if (a != 0) GV.WZ[a] = (Complex)GV.WZ[a] + z;

else

{ GV.WZ.Add(z);

GV.A[i, j] = (short)(GV.WZ.Count - 1); }

}

Аналогичным образом формируется и компактный массив SZ.

При приведении матрицы к диагональной единичной матрице на основе алгоритма Гаусса-Жордано при использовании полного описания разреженной матрицы w необходимо выполнять k-й шаг преобразования элементов матрицы на основе соотношения

w[i, j] = w[i, j] - w[i, k] *w[k, j]/w[k, k];

При использовании компактного описания такая задача решается путем вызова функции step(i, j), которая имеет описание

void step(int i, int j)

{ if (GV.A[i, j] == 0)

{ GV.WZ.Add(new Complex(0, 0));

GV.A[i, j] = (short)(GV.WZ.Count - 1); }

GV.WZ[GV.A[i, j]] = (Complex)GV.WZ[GV.A[i, j]] -

(Complex)GV.WZ[GV.A[i, k]]*

(Complex)GV.WZ[GV.A[k, j]]/ (Complex)GV.WZ[GV.A[k, k]];

}

В результате выполнения этой функции в массиве WZ будет выполнен k-й шаг виртуального преобразования к диагональной единичной матрице. Аналогичным образом решается задача преобразования массива SZ.

После виртуального приведения на основании компактного описания системы диагональных блочных матриц к единичным матрицам, нетрудно выполнить расчет моделируемой системы в соответствии с выражениями (2),(3),(4) и (5), при этом здесь нет необходимости перехода к компактной форме описания задачи, поскольку все используемые в этих выражениях матрицы не характеризуются сильной разреженностью.

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