Программная реализация редуцированной вычислительной схемы численного интегрирования системы линейных дифференциальных уравнений на основе методов Эйлера

ЖАДОВ А.Д., ВОСТРИКОВ А.В., БОРИСОВ Н.И., ТУМКОВСКИЙ С.Р.

Программная реализация редуцированной вычислительной схемы
численного интегрирования системы линейных дифференциальных уравнений на основе методов Эйлера

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

Одними из важнейших проблем длительной эксплуатации космических аппаратов (КА) являются электризация и связанные с этим процессом электростатические разряды (ЭСР). В 30% случаев потеря КА происходит из-за электризации. ЭСР порождают токи на поверхности КА, которые наводят в бортовой кабельной сети наводки величиной до 10 В. Помехи такой величины способны привести к отказам бортовой радиоэлектронной аппаратуры. Для сведения: к минимуму возможных сбоев электроники КА необходимо смоделировать картину растекания тока по поверхности КА. Далее проводится расчет наводок в кабелях и даются рекомендации по их экранированию или оптимальному маршруту прокладки. С этой целью в МИЭМ были разработаны структурная электрофизическая модель (СЭМ) КА [1] и программное обеспечение ПО «Satellite-MIEM» для ее реализации.

ПО «Satellite-MIEM» синтезирует структурную электрофизическую модель (СЭМ) на базе заданной полигональной 3D-модели КА (см. рис. 1). Полигональная модель состоит из совокупности элементарных фигур – треугольников или прямоугольников, преобразуемых средствами программы в поверхностную сетку: совокупность связанных узлов. Сетка является равномерной, то есть шаги по всем координатам равны между собой. Ячейка сетки представляется прямоугольником или треугольником. Количество узлов СЭМ крупногабаритного КА будет равняться (1…2)х105. Каждая связь (ветвь) представляется в виде элементов электрической цепи (R, L, C), в целом образующих эквивалентную электрическую схему поверхности КА. При этом номиналы элементов одного типа одинаковы.

Рисунок 1. Преобразование полигональной модели в СЭМ

С помощью ПО «Satellite-MIEM», зная место электрического разряда, можно получить картину растекания тока. Для расчета переходных токов используется наиболее производительная программа расчета электрических схем LTSpice, однако для анализа схемы из 150 000 узлов требуется 83 часа (эксперимент проводился на ЭВМ с двуядерным процессором с тактовой частотой 1,8 ГГц на каждом ядре, объем оперативной памяти равняется 2 Гб) [2], что неприемлемо для предприятий космической отрасли. Используемые в ПО численные методы хорошо показывают себя только при расчете небольших моделей [3]. Поэтому возникает необходимость в создании новых методов расчета больших электрических схем. Ранее в [4-6] был предложен новый подход редукции модели линейной схемы, основанный на исключении из модели подвекторов, содержащих в себе фазовые переменные, значения которых не превышают 1-2% от величины приложенного к месту разряда тока. В данной работе авторы разрабатывают программу для реализации редуцированной вычислительной схемы и проверяют теоретические оценки эффективности вычислительной схемы на практике.

1. Построение редуцированной вычислительной схемы

Модель ЭЭС КА может быть сформирована в расширенном однородном координатном базисе (РОКБ) и записана в виде системы линейных обыкновенных дифференциальных уравнений [4]:

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

Отметим специфику модели схемы.

1) Матрица С является диагональной и вырожденной. Диагональ состоит из двух групп коэффициентов и коэффициенты внутри каждой группы одинаковы.

2) Матрица G – невырожденная, симметричная и разреженная матрица.

3) Вектор  содержит только один или два ненулевых коэффициента вида , где i может быть любым числом от 1 до n.

Необходимо с минимальными временными затратами выполнить решение системы уравнений в момент времени t*, т.е. вычислить числовой вектор .

2. Редукция вычислительной схемы методов Эйлера

Для решения задачи (1) формируется модель схемы в РОКБ следующим образом:

, (2)

Отметим, что матрицы – симметричные, подматрицы C11, L22, G33 – диагональные, при этом:

•  – подвектор напряжений в узлах, к которым подсоединены конденсаторы;
•  – подвектор токов, проходящих через индуктивные элементы;
•  – подвектор напряжений в узлах, к которым не подсоединены конденсаторы.

Модель (2) записывается в виде:

, (3)

, (4)

, (5)

Исключаем из модели подвектор . Из (5) получим:

, (6)

Подстановка  из (6) в (3) и (4) приводит к следующим выражениям:

(7)

 (8)

или

(9),

где r – сопротивление резисторов в схеме.

Модель может быть записана в виде системы дифференциальных уравнений, заданной в явной форме:

 , (10)

где 

Будем считать, что исходная задача имеет следующий вид:

, . (11)

При использовании численных методов процесс решения (11) заключается в вычислении вектора решения в моменты времени , т.е. , ,…, .

Рассмотрим одношаговый метод Эйлера для решения задачи (11). Явный метод Эйлера решения системы уравнений с учетом конкретного вида системы уравнений (11):

. (12)

Рассмотрим неявный метод Эйлера. С учетом конкретного вида системы уравнений (11) получим следующую формулу метода Эйлера:

 . (13)

Для построения новой вычислительной схемы необходимо преобразовать задачи (12) и (13) в некоторую редуцированную модель (макромодель) КА, количество уравнений в которой будет равно

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

 , (14)

 , (15)

где , .

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

  (16)

или

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

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

• На первом этапе строить редуцированную вычислительную схему.
• На втором этапе проводить расчет картины растекания токов по поверхности КА с помощью вычислительной схемы.

Иерархическую структуру разработанной программы покажем на рис. 2.

Рисунок 2. Иерархическая структура разработанной программы

Разработанное ПО для получения требуемого результата выполняет ряд этапов:

1. Конвертирование ЭЭС КА из программы LTSpice в удобный для дальнейших преобразований матричный вид.
2. Преобразование матриц для ускорения скорости расчета по редуцированной вычислительной схеме.
3. Расчет ЭЭС с помощью вычислительной схемы.
4. Трансляция файла полученных результатов и построение модели КА с учетом картины растекания токов.

Рассмотрим алгоритм конвертирования файла LTSpice в удобный для дальнейших преобразований матричный вид и результат его выполнения (см. рис. 3):

1. Объявление и обнуление переменных и массивов, выделение памяти.
2. Открытие файлов ввода данных и вывода результата (вывод ошибки при невозможности открытия).
3. Построчное чтение файла и запись его в символьный массив.
4. Поиск в файле величины и узла включения источника возмущающего воздействия.
5. Проход по всему массиву, нахождение в нем C-переменных и запись их номиналов в массив С в порядке возрастания по диагонали, подсчет количества С, фиксирование новых позиций для С-переменных, поставить источник возмущающего воздействия в узел включения первого из элементов С.
6. Проход по всему массиву, нахождение в нем L-переменных и запись их номиналов с отрицательным знаком в массив С в порядке возрастания по диагонали, после позиций элементов С.
7. Проход по всему массиву, нахождение в нем R-переменных и узлов их включения. Если узел связан с фиксированной позицией С, заменить его на номер этой позиции, иначе присвоить номер: сумма всех элементов С + сумма всех элементов L + порядок вставки. Порядок вставки увеличить на 1.
8. Заполнение массива G: считываем номиналы текущих R и по данным из шестого шага узлам включения ставим их на соответствующие позиции в матрице симметрично – G[узел1][узел 2]=-1/номинал; G[узел2][узел1]=-1/номинал, а также заполнение диагонали соответствующих узлов (прибавление к элементу диагонали величины, обратной номиналу).

Рисунок 3. Пример результата конвертирования файла LTSpice в удобный для дальнейших преобразований матричный вид согласно (2)

Преобразования матрицы А (16) отличаются от предложенных в [5], где была получена треугольная матрица с окаймлением с помощью метода определяющих величин. Предложенный алгоритм позволяет получить более выгодную с точки зрения дальнейших вычислений треугольную матрицу, но посредством арифметических операций со строками матрицы, вследствие чего может быть потеряна точность вычислений. Алгоритм и пример результата (см. рис. 4) приведены ниже:

1. Редукция матриц С и R.
2. Получение матрицы А путем перемножения матрицы G и C^-1.
3. Поиск в матрице А строк и столбцов, содержащих наибольшее количество ненулевых элементов.
4. Нулевые строки поместить в низ матрицы.
5. Остальные строки расположить в порядке возрастания количества элементов.
6. Расположить в порядке возрастания столбцы слева направо в матрице А так, чтобы максимальные по модулю элементы попадали на главную диагональ.
7. Нулевые элементы на диагонали исключить посредством сложения строки с нулевым элементом на диагонали и строки с максимальным элементом по модулю в этом же столбце.
8. Исключить ненулевые элементы над главной диагональю посредством вычитания текущей строки от нижестоящей строки, содержащей элемент диагонали в этом же столбце, и умножения на коэффициент, равный отношению текущего элемента и найденного ниже (идти по матрице необходимо с нижнего правого угла).

Рисунок 4. Результат преобразования матрицы

На современном компьютере (ЭВМ с двуядерным процессором с тактовой частотой 1,8 ГГц на каждом ядре, объем оперативной памяти равняется 2 Гб) проведены исследования временных затрат на анализ модели ЭЭС, содержащей до 3000 обыкновенных дифференциальных уравнений. Время анализа неявным методом Эйлера, явным методом Эйлера и по предложенной новой вычислительной схеме с учетом времени построения самой вычислительной схемы показано на рис. 5. Для модели, состоящей из 100 уравнений, неявному методу Эйлера потребовалось 0,16 с для анализа, расчет по новой вычислительной схеме при размерности подматрицы (10х10) занял 0,012 с, что в 13 раз быстрее (таблица 1). Для модели из 3000 уравнений расчет неявным методом Эйлера потребовал 28397 с для анализа, вычисления по новой схеме при размерности подматрицы (10х10) затратили 387,5 с, что в 75 раз быстрее. В случае больших систем использование редуцированной схемы оказывается выгоднее на два-три порядка других неявных численных методов. При этом результат расчета в выбранной локальной области традиционными численными методами и по редуцированной вычислительной схеме совпадает.

Рисунок 5. Время анализа модели ЭЭС разными методами

Таблица 1. Результат тестирования численных методов на ПК

Существенно отличается время анализа явным методом Эйлера и по предложенной новой вычислительной схеме без учета времени построения самой вычислительной схемы (такой вариант используется при многократном использовании макромодели). Для модели из 3000 уравнений явному методу Эйлера потребуется 51,3 с для анализа, вычисления по новой схеме при размерности подматрицы (10х10) потребуют 0,005 с, что в 10 260 раз быстрее, в случае подматрицы (100х100) вычисления по новой схеме потребуют 2,9 с, что в 17,7 раза быстрее. На практике подматрица будет содержать менее 1% коэффициентов от их общего количества модели, поэтому выгода может достигать четырех порядков по сравнению даже с явными численными методами.

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

1. Новиков Л.С., Бабкин Г.В., Морозов Е.П., Колосов С.А., Крупников К.К., Милев В.Н., Саенко В.С. Комплексная методология определения параметров электростатической зарядки, электрических полей и пробоев на космических аппаратах в условиях их радиационной электризации. Руководство для конструкторов. – ЦНИИМАШ, Королев 1995.
2. Востриков А.В., Абрамешин А.Е.Тестирование коммерческого программного обеспечения для моделирования и анализа эквивалентных электрических схем космических аппаратов. // «Технологии электромагнитной совместимости», №1, 2012 г. – С. 25-28.
3. Azmy S. Ackleh, Edward James Allen, R. Baker Kearfott, Padmanabhan Seshaiyer. //Classical and Modern Numerical Analysis: Theory, Methods and Practice/ 2009 by Chapman and Hall/CRC, 628 p.
4. Востриков А.В., Абрамешин А.Е., Борисов Н.И. Расчет наводок в бортовой кабельной сети космических аппаратов с помощью макромоделирования на основе методов Эйлера. //«Технологии электромагнитной совместимости», 2012. №1, 2012 г. – С. 19-24.
5. Востриков А.В., Абрамешин А.Е. Вычислительная схема ускоренного метода расчета наводок в бортовой кабельной сети космических аппаратов. //«Технологии электромагнитной совместимости», №3, 2012 г. – С. 22-28.
6. Востриков А.В., Борисов Н.И. Разработка эффективного метода анализа эквивалентных электрических схем космических аппаратов, использующего редукцию и разреженность матриц моделей. //В кн.: Новые информационные технологии в автоматизированных системах: материалы пятнадцатого научно-практического семинара/Отв. ред.: С.Р. Тумковский. – М.: Московский государственный институт электроники и математики, 2012. – С. 113-120.

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

Software implementation of the reduced computational schemes for the numerical integration of the system of linear differential equations based on the methods of Euler

Jadov A.D., Vostrikov A.V., Borisov N.I, Tumkovskiy S.R.

The paper considers the implementation of the basic algorithms for constructing the reduced computational schemes for the numerical integration of the system of linear differential equations by Euler's method. Shows the time and construction costs reduced by calculation scheme as compared with conventional numerical methods. Payment scheme for computing large mathematical models occurs at 2-3 orders of magnitude faster used in specialized software techniques.

Keywords: algorithm, the reduction of the system of differential equations, matrix transformation

Комментарии могут оставлять только зарегистрированные пользователи