Локальные решатели — Ускоряем физическую модель в три раза

Введение

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

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

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

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

Адаптер питания

Рассмотрим достаточно простую модель адаптера питания, выполняющего преобразование переменного тока в постоянный (пульсирующий). Подключим его к источнику синусодиальной ЭДС с действующим значением 220 В и частотой 50 Гц. Ток с источника поступает на понижающий трансформатор, представленный упрощённой моделью в виде резистивно-индуктивного элемента и идеального трансформатора за ним. Получившееся напряжение, уже с меньшей амплитудой, попадает на диодный мост. Диоды открываются попеременно парами, а результирующий ток сглаживается параллельно подключенным фильтрующим конденсатором. К выводам конденсатора подключается обобщённая резистивная нагрузка величиной 700 Ом. Не забываем добавить "Electrical Reference" по обе стороны от трансформатора.

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

Заменим "Inductor" на три одинаковые подсистемы "50 RL". Они состоят из пяти одинаковых подсистем, каждая из которых включает 10 "Inductor". Так как они соединены последовательно, индуктивности каждого из 150 "Inductor" должны быть равны 1/150 от требуемой эквивалентной индуктивности исходного "Inductor". То же касается и последовательной резистивной составляющей.

N_L=150;      # Количество Inductor
L=1.7e-3/N_L; # Индуктивность одного Inductor, Гн
r_L=0.17/N_L; # Последовательное сопротивление одного Inductor, Ом

Аналогично поступаем и с "Capacitor", но их соединяем параллельно в подсистемах "50 C".

N_C=150;    # Количество Capacitor
C=1e-2/N_C; # Ёмкость одного Capacitor, Ф

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

Обратите внимание на наличие паразитных параллельной проводимости и последовательного сопротивления в "Inductor" и "Capacitor" соответственно. Если их сделать нулевыми, никакое число дополнительных "Inductor" и "Capacitor" не будет приводить к появлению новых степеней свободы. Действительно, последовательное соединение идеальных катушек индуктивности означает, что через все из них протекает один и тот же ток. Это не позволяет нам выбирать произвольные значения тока через каждую катушку независимо. То же самое происходит при параллельном соединенении идеальных конденсаторов. Математически мы получаем систему дифференциально-алгебраических уравнений индекса 2, а затем индекс успешно понижается ядром Engee. В итоге мы получаем простейшую систему алгебраических уравнений, большинство из которых классифицируются как выходные. Это очень полезная функциональность вычислительного ядра, которая значительно ускоряет расчёт и повышает его стабильность, но в нашей постановке чисто математического эксперимента будет только мешать.

Чтобы исключить хотя бы часть факторов, влияющих на производительность модели, которые могут только усложнить анализ эффекта от использования нескольких локальных решателей, установим дискретную синхронизацию, аналитический расчёт матрицы Якоби, быструю версию неявного метода Эйлера с периодом дискретизации 0.0001:

В качестве глобального решателя можно явно выбрать "discrete" с тем же шагом 0.0001, так как мы знаем, что вне физической сети нет непрерывных состояний:

Но даже если этого не сделать, Engee самостоятельно заменит выбранный пользователем непрерывный решатель на дискретный для ускорения работы модели. Результаты расчёта при этом останутся теми же самыми.

Запустим модель:

example_path = @__DIR__; # Получаем абсолютный путь директории, содержащей текущий скрипт
cd(example_path); # Переходим в директорию примера

model_name = "local_solvers_partitioned_vs_whole_one.engee"
engee.open(model_name); # Открываем модель
model = engee.load(model_name); # Загружаем модель для быстрого доступа

one_data = engee.run(model); # Выполняем прогон и получаем результаты расчёта

# Вспомогательная функция для получения времени старта и симуляции модели
function get_model_times(model) 
    prepare_time, compile_time, init_time, sim_time = 0.0, 0.0, 0.0, 0.0

    logs = engee.get_logs(model)

    prepare_time  = parse(Float64, (match(r"\d+\.\d+", logs[1][:content])).match)
    compile_time  = parse(Float64, (match(r"\d+\.\d+", logs[2][:content])).match)
    init_time     = parse(Float64, (match(r"\d+\.\d+", logs[3][:content])).match)
    sim_time      = parse(Float64, (match(r"\d+\.\d+", logs[4][:content])).match)

    return prepare_time + compile_time + init_time, sim_time
end

one_start_time, one_simulation_time = get_model_times(model);

Посмотрим на результаты расчёта:

plot(eachcol(collect(one_data["Voltage Sensor.V"]))..., ylabel="Напряжение нагрузки, В", xlabel="Время, с", titlefont=font(10), xguidefont=font(10), leg=false)

При холодном старте спустя где-то 0.2 с на нагрузке обеспечивается практически постоянное напряжение 12 В. Оно колеблется с двойной частотой относительно заданной во внешнем источнике, что обусловлено поочередной отработкой диодов.

one_start_time

39.2778

секунд занял запуск модели на расчёт, а

one_simulation_time

32.596

секунд занял сам расчёт. Попробуем ускорить модель.

Разбиение на две физических сети

Разобьём модель на две практически равные части сразу после диодного моста.

Поставим в левой сети датчик тока, который будет передавать сигнал через единичное усилительное звено на идеальный источник тока в правой сети. А в правой сети будем с уже имеющегося датчика напряжения снимать актуальное значение и передавать его в идеальный источник ЭДС в левой сети. Обратите внимание на "Unit Delay", необходимый для разрыва получившейся алгебраической петли по напряжению.

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

Engee требует, чтобы в каждой электрической сети был хотя бы один собственный экземпляр "Electrical Reference". В данном случае его можно добавить в любом месте.

Обязательно добавляем копию "Solver Configuration" в правую сеть. Без него модель не запустится, так как каждая физическая сеть обязательно должна иметь собственный единственный экземпляр "Solver Configuration". Теперь обе физических сети считаются независимыми локальными решателями.

Запустим модифицированную модель:

model_name = "local_solvers_partitioned_vs_whole_two.engee"
engee.open(model_name);
model = engee.load(model_name);
two_data = engee.run(model);
two_start_time, two_simulation_time = get_model_times(model);

two_start_time

23.279899999999998

секунд занял запуск модели с двумя локальными решателями, а

two_simulation_time

9.653

секунд продолжался сам расчёт. Модифицированная модель запустилась быстрее в

one_start_time/two_start_time

1.687197969063441

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

one_simulation_time/two_simulation_time

3.3767740598777576

раза! Наибольший вклад на этом этапе обычно вносится алгоритмом декомпозиции матрицы Якоби, который в нашем случае характеризуется временной сложностью , где — количество уравнений в системе. Отсюда, максимальное ускорение этой части симуляции за счёт решения двух вдвое меньших задач составляет четыре раза.

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

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

p1 = plot(eachcol(collect(one_data["Voltage Sensor.V"]))..., ylabel="Напряжение нагрузки, В", label="Одна сеть")
plot!(p1, eachcol(collect(two_data["Voltage Sensor.V"]))..., label="Две сети")

abs_error = abs.(two_data["Voltage Sensor.V"].value .- one_data["Voltage Sensor.V"].value)
p2 = plot(one_data["Voltage Sensor.V"].time, abs_error, ylabel="Абсолютная разница, В", leg=false)

rel_error = abs_error ./ abs.(one_data["Voltage Sensor.V"].value)
p3 = plot(one_data["Voltage Sensor.V"].time, rel_error .* 100, xlabel="Время, с", ylabel="Относительная разница, %", leg=false)
plot(p1, p2, p3, layout=(3,1), leg=false, size=(1400,800), titlefont=font(10), xguidefont=font(10))

Без масштабирования графика и тем более без накладывания графиков друг на друга разницу практически невозможно заметить. Она носит колебательный характер, в первые миллисекунды составляет немногим более 1%, но дальше очень быстро падает и почти на всём интервале моделирования колеблется около 0.01%. Откуда она взялась? Мы изменили математическую модель внесением задержки в падение напряжения правой части схемы, когда добавили единичное запаздывание для разрыва алгебраической петли.

Заключение

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

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

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

Локальные решатели — Ускоряем физическую модель в три раза