Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 504, № 1, стр. 56-59

Наиболее эффективным способом расчета динамических режимов больших электронных схем является их декомпозиция с последующим расчетом моделируемой схемы по частям [1–3]. Введем для каждого узла схемы узловое напряжение и из всех введенных узловых переменных ${{{v}}_{p}}$ образуем вектор узловых потенциалов X₀ = V₀ = ${{[...,{{{v}}_{p}},...]}^{t}}$. В каждой подсхеме выделим внутренние переменные x_i и образуем для каждой подсхемы вектор ${{X}_{k}} = {{[...,{{x}_{i}},...]}^{t}}$. Тогда результирующий вектор всех переменных схемы будет иметь вид X = = ${{[X_{1}^{t},...,X_{k}^{t},...,X_{m}^{t},X_{0}^{t}]}^{t}}$.

Уравнение всей схемы в расширенном базисе узловых потенциалов при этом может быть записано в виде

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

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

Обозначим

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

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

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

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

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

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

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

Построение математического описания для расчета динамических режимов на основе декомпозиции схемы выполняется на основе алгебро-дифференциального уравнения схемы, записанного в неявном виде $F(\dot {X},X,t) = 0$, которое формируется в расширенном базисе узловых потенциалов. Для формирования алгебро-дифференциальных уравнений необходимо включить все полюса емкостных компонентов схемы в число y-полюсов, а все полюса индуктивных компонентов – в число z-полюсов схемы, выделив эти элементы из структуры многополюсников. При этом уравнение моделируемой схемы можно записать в виде

Заменяя оператор s операцией дифференцирования, получим алгебро-дифференциальное уравнение моделируемой системы в неявной форме, которое для момента времени t = k + 1 имеет вид $F({{\dot {X}}^{{k + 1}}},{{X}^{{k + 1}}},{{t}^{{k + 1}}}) = 0$. Для исключения вектора производной ${{\dot {X}}^{{k + 1}}}$ следует воспользоваться формулами коррекции ${{\dot {X}}^{{k + 1}}} = a{{X}^{{k + 1}}} + b$ для многошаговых неявных методов высших порядков Адамса-Маултона или Гира [4, 5], что приводит к уравнению $F(a{{X}^{{k + 1}}} + b,{{X}^{{k + 1}}},{{t}^{{k + 1}}}) = 0$, которое можно решить методом Ньютона-Рафсона относительно вектора ${{X}^{{k + 1}}}$, поскольку значения коэффициентов a и b к моменту времени t = k + 1 уже определены на предыдущих шагах расчета. В качестве начального приближения ${{({{X}^{{k + 1}}})}^{i}} = {{({{X}^{{k + 1}}})}^{0}}$ при этом обычно принимается значение вектора на предыдущем шаге X^k.

Процесс исключения значений производных целесообразно выполнять непосредственно в процессе формирования дискретизированных компонентных уравнений $i_{C}^{{k + 1}} = Y_{C}^{k}u_{C}^{{k + 1}} + j_{C}^{k}$ и $u_{L}^{{k + 1}} = R_{L}^{k}i_{L}^{{k + 1}} + e_{L}^{k}$ для двухполюсников C и L, а также для частотно-зависимых управляемых источников [6].

В процессе обработки алгебро-дифференциального уравнения $F({{\dot {X}}^{{k + 1}}},{{X}^{{k + 1}}},{{t}^{{k + 1}}}) = 0$ необходимо решать методом Ньютона–Рафсона уравнение $F(a{{X}^{{k + 1}}} + b,{{X}^{{k + 1}}},{{t}^{{k + 1}}}) = 0$, которое в общем случае является нелинейным. Для решения этого уравнения целесообразно выполнить схемотехническую интерпретацию метода Ньютона–Рафсона, основанную на линеаризации нелинейных компонентов схемы в каждой точке итерационного процесса.

Так, если нелинейные компоненты в общем случае описываются уравнением P₂= F(Q₂), то разлагая это уравнение в ряд Тейлора в точке $Q_{2}^{{}} = Q_{2}^{i}$ и отбрасывая высшие производные, получим уравнение вида

где $W_{2}^{i} = {{\left. {\frac{{\partial {{P}_{2}}}}{{\partial {{Q}_{2}}}}} \right|}_{{Q_{2}^{i}}}}$, $S_{2}^{i} = P_{2}^{i} - W_{2}^{i}Q_{2}^{i}$ – линеаризованные в точке ${{Q}_{2}} = Q_{2}^{i}$, ${{P}_{2}} = P_{2}^{i}$ матрицы неавтономных и автономных параметров нелинейных компонентов.

Пусть нелинейная схема описывается уравнением

где W₁ и S₁ – матрица и задающий вектор линейной части схемы, S₂ = P₂ = F(Q₂) – вектор фиктивных задающих источников, отображающих нелинейные характеристики компонентов.

Тогда на основании уравнения линеаризованных компонентов (8) и уравнения нелинейной схемы (9) можно получить матричное уравнение моделируемой схемы

где W₁ и S₁ – матрица и задающий вектор линейной части схемы, $W_{2}^{i}$ и $S_{2}^{i}$ – матрица и задающий вектор линеаризованной в точке ${{Q}_{2}} = Q_{2}^{i}$ нелинейной части схемы.

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

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

1. Ввод начальных значений ${{t}^{j}} = {{t}^{0}}$, ${{(X_{k}^{j})}^{i}} = X_{k}^{0}$, ${{(X_{0}^{j})}^{i}} = X_{0}^{0}$.

2. ${{t}^{{j + 1}}} = {{t}^{j}} + h$.

3. Начало 1-го цикла k = 1, m.

4. Дискретизация частотно-зависимых компонентов в точке $Q_{s}^{j}$.

5. Ввод начального значения для линеаризации ${{(X_{k}^{{j + 1}})}^{i}} = {{X}_{{пр}}}$.

6. Линеаризация нелинейных компонентов в точке ${{(Q_{f}^{{j + 1}})}^{i}}$.

7. Формирование блочных матриц ${{(W_{{kk}}^{j})}^{i}},$ ${{(W_{{k0}}^{j})}^{i}},$ ${{(W_{{0k}}^{j})}^{i}},$ ${{(W_{{00}}^{j})}^{i}}$ и векторов ${{(S_{k}^{j})}^{i}},$ ${{(S_{0}^{j})}^{i}}$ уравнения (1).

8. Преобразование диагональных подматриц ${{(W_{{kk}}^{j})}^{i}}$ уравнения (1) к единичной матрице и расчет ${{(\bar {W}_{{k0}}^{j})}^{i}},$ ${{(\bar {S}_{k}^{j})}^{i}}$.

9. Вычисление поправок ${{({{\Delta }_{k}}W_{{00}}^{j})}^{i}}\, = \,{{(W_{{0k}}^{j})}^{i}}{{(\bar {W}_{{k0}}^{j})}^{i}},$ ${{({{\Delta }_{k}}S_{0}^{j})}^{i}} = {{(W_{{0k}}^{j})}^{i}}{{(\bar {S}_{k}^{j})}^{i}}$ и расчет новых значений блочной матрицы и задающего вектора для узлов связи

10. Окончание 1-го цикла k = 1, m.

11. Решение уравнения ${{(\bar {W}_{{00}}^{j})}^{i}}{{(X_{0}^{{j + 1}})}^{{i + 1}}} + {{(\bar {S}_{0}^{j})}^{i}}$ = 0 для переменных узлов связи.

12. Начало 2-го цикла k = 1, m.

13. Расчет вектора ${{(X_{k}^{{j + 1}})}^{{i + 1}}}$ = $ - {\kern 1pt} {{(\bar {W}_{{k0}}^{j})}^{i}}{{(X_{0}^{{j + 1}})}^{{i + 1}}}$ – – ${{(\bar {S}_{k}^{j})}^{i}}$ внутренних переменных подсхем.

14. Конец 2-го цикла k = 1, m.

15. Вычисление нормы $N = {\text{||}}{{(X_{k}^{{j + 1}})}^{{i + 1}}} - {{(X_{k}^{{j + 1}})}^{i}}{\text{||}}$.

16. Если N > ε, то переопределение векторов внутренних переменных подсхем ${{(X_{k}^{{j + 1}})}^{i}}\, = \,{{(X_{k}^{{j + 1}})}^{{i + 1}}}$ и вектора переменных узлов связи ${{(X_{0}^{{j + 1}})}^{i}}\, = \,{{(X_{0}^{{j + 1}})}^{{i + 1}}}$ и возврат к 5, иначе вывод векторов внутренних переменных подсхем ${{(X_{k}^{{j + 1}})}^{{i + 1}}}$ и переменных связи ${{(X_{0}^{{j + 1}})}^{{i + 1}}}$.

17. Если ${{t}^{{j + 1}}} < {{t}_{{{\text{оконч}}}}}$, то переход к 2, иначе конец расчета.

Реализация программного обеспечения для расчета динамического режима на основе организации автономных вычислительных процессов для отдельных подсхем оказывается целесообразной в случаях, когда отдельные подсхемы имеют существенно отличные свойства от других подсхем. Сущность автономного вычислительного процесса заключается в том, что для каждой подсхемы выбирается свой шаг расчета h_k и с этим шагом выполняется расчет вектора внутренних переменных подсхем$X_{k}^{{j + 1}}$. При этом выбираются точки синхронизации с шагом H > h_k, в которых расчет процессов в подсхемах прекращается и начинается расчет вектора переменных узлов связи $X_{0}^{{j + 1}}$ с последующим уточнением вектора внутренних переменных $X_{k}^{{j + 1}}$.

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

1. Ввод начальных значений ${{T}^{j}} = {{T}^{0}}$, $X_{k}^{j} = X_{k}^{0}$, $X_{0}^{j} = X_{0}^{0}$.

2. Начало 1-го цикла k = 1, m.

3. ${{t}^{j}} = {{T}^{j}}$.

4. ${{t}^{{j + 1}}} = {{t}^{j}} + {{h}_{k}}$.

5. Дискретизация частотно-зависимых компонентов в точке $Q_{s}^{j}$.

6. Формирование блочных матриц $W_{{kk}}^{j},$ $W_{{k0}}^{j},$ $W_{{0k}}^{j},$ $W_{{00}}^{j}$ и векторов $S_{k}^{j},$ $S_{0}^{j}$ подсхем уравнения (1).

7. Преобразование диагональных блочных матриц $W_{{kk}}^{j}$ к единичной матрице и расчет ${{(\bar {W}_{{k0}}^{j})}^{i}},$ ${{(\bar {S}_{k}^{j})}^{i}}$.

8. Расчет вектора внутренних переменных подсхем $X_{k}^{{j + 1}} = - \bar {W}_{{k0}}^{j}X_{0}^{j} - \bar {S}_{k}^{j}$.

9. Если ${{t}^{{j + 1}}} < = H,$ то ${{t}^{j}} = {{t}^{{j + 1}}}$, $X_{k}^{j} = X_{k}^{{j + 1}}$ и возврат к 4, иначе переход к 10.

10. Вычисление поправок ${{\Delta }_{k}}W_{{00}}^{j}\, = \,W_{{0k}}^{j}\bar {W}_{{k0}}^{j},$ ${{\Delta }_{k}}S_{0}^{j}$ = = $W_{{0k}}^{j}\bar {S}{{_{k}^{j}}^{i}}$ и расчет новых значений блочной матрицы и задающего вектора переменных узлов связи

11. Окончание 1-го цикла k = 1, m.

12. ${{t}^{{j + 1}}} = {{T}^{j}} + H$.

13. Решение уравнения $\bar {W}_{{00}}^{j}X_{0}^{{j + 1}} + \bar {S}_{0}^{j} = 0$ для переменных узлов связи.

14. Начало 2-го цикла k = 1, m.

15. Уточнение вектора внутренних переменных подсхем $X_{k}^{{j + 1}} = - \bar {W}_{{k0}}^{j}X_{0}^{{j + 1}} - \bar {S}_{k}^{j}$.

16. Окончание 2-го цикла k = 1, m.

17. Вывод ${{t}^{{j + 1}}}$, $X_{k}^{{j + 1}}$, $X_{0}^{{j + 1}}$.

18. Если ${{t}^{{j + 1}}} < {{t}_{{{\text{оконч}}}}}$, то переопределение переменных ${{T}^{j}} = {{t}^{{j + 1}}}$, $X_{k}^{j} = X_{k}^{{j + 1}}$, $X_{0}^{j} = X_{0}^{{j + 1}}$ и возврат к 2, иначе конец расчета.

Для нелинейных схем алгоритм дополняется шагами линеаризации компонентов.

Организация автономных вычислительных процессов дает значительный выигрыш при расчете схем, в которых можно выделить “быстрые” и “медленные” подсхемы, поскольку для расчета “быстрых” подсхем можно выбрать малые шаги расчета h_k . Пусть, например, схема содержит 2 подсхемы – “быструю”, которую можно рассчитывать с шагом h₁ = 10^–6 с, и “медленную”, которая рассчитывается с шагом H = h₂ = 1 с. Если t_оконч = 10 с и для расчета подсхем и переменных узлов связи на одном шаге необходимо выполнить N_во вычислительных операций, то в этом случае число вычислительных операций при автономном вычислительном процессе за время t_оконч = 10 с будет равно N_во1 = (10⁶ + 1 + 1)10N_во = = 10⁷N_во. При сквозном вычислительном процессе в этом случае имеем N_во2 = 3 × 10⁶ 10 N_во = 3 × 10⁷N_во. Следовательно, при использовании автономного итерационного процесса время решения будет в три раза меньше, чем при организации сквозного процесса вычислений.

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

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Работа выполняется в рамках темы № FFSM-2019-0001.