Переходные процессы в линейных цепях
andreyilinskiyУравнения для переходных процессов в линейной электрической цепи
В документации Engee есть пример, где описывается применение пакета DifferentialEquations.jl для решения задачи по переходным процессам для электрической схемы с двумя индуктивностями и двумя конденсаторами.
В данном скрипте я бы хотел еще раз обратиться к этой же задаче и решить ее с более общих позиций.
Исходный пример
В исходном примере задача решена методом переменных состояния. В качестве переменных состояния выбраны токи через индуктивности , и напряжения на конденсаторах , . Алгебраические переменные , , выражены через переменные состояния. Первоначальная дифференциально-алгебраическая система уравнений тем самым сведена к дифференциальной и использована функция ODEProblem для постановки задачи по решению системы обыкновенных дифференциальных уравнений.
Проблема может быть в том, что необходимо выразить алгебраические переменные через переменные состояния. Для небольшой линейной схемы это можно сделать вручную. Но математически процесс выражения алгебраических переменных, через известные с прошлого шага переменные состояния равносилен решению СЛАУ, сложность решения которой быстро увеличивается с увеличением количества уравнений. Если же система нелинейная то аналитически выразить алгебраические переменные может быть невозможно.
Составим систему DAEs
Вместо уравнений из примера (2 узловых, 3 контурных, два для связи тока и напряжения на кондесаторах) напишем уравнения, которые проще сформировать автоматически, анализируя топологию схемы.
Введём напряжения относительно земли для каждого узла (в том числе в котором соединены 2 ветви - см. рисунок). Запишем 2 узловых уравнения, закон Ома для каждого элемента схемы и уравнения конденсаторов:
\begin{equation*}
\begin{cases}
i_1-i_2-i_3=0\
i_3-i_4-i_5=0\
u_t-u_1=i_1R_1\
u_1-u_2=L_1i^\prime_1\
u_2=i_2R_2\
u_2-u_3=U_{C1}\
C_1 U^\prime_{C1}=i_3\
u_3-u_4=i_4 R_3\
u_4=L_2i^\prime_4\
u_3-u_5=U_{C2}\
C_2 U^\prime_{C2}=i_5\
u_5=i_5 R_4
\end{cases}
\end{equation*}
Мы получили дифференциально-алгебраическую систему уравнений индекса 1 (semi-explicit DAE of index 1) с пятью неизвестными токами и семью неизвестными напряжениями. В [1] такие системы охарактеризованы как:
Systems of the form ... are arranged automatically by tools of computer aided design (CAD).
Будем решать эту систему напрямую, без дополнительных преобразований.
Использование записи в виде Mass-Matrix Differential-Algebraic Equations
Пакет DifferentialEquations.jl предлагает для таких систем использовать либо запись с матрицей масс, либо ставить задачу как неявную и использовать DAEProblem вместо ODEProblem. Для сокращения длины скрипта воспользуемся только первым способом.
Введем искомую вектор функцию:
и запишем систему в виде
Нулевые строки соответствуют алгебраическим уравнениям, ненулевые - дифференциальным.
Далее я буду переиспользовать код из оригинального примера.
using DifferentialEquations
# Ввод параметров электрической цепи
R1 = 10.0
R2 = 15.0
R3 = 20.0
R4 = 10.0
f = 50.0
C1 = 10e-6
C2 = 30e-6
L1 = 15e-3
L2 = 20e-3
# Определение дифференциально-алгебраической системы
function circuitDAEsemiexplicit(dy, y, p, t)
R1, R2, R3, R4 = p
i1, i2, i3, i4, i5, u1, u2, u3, u4, u5, uc1, uc2 = y
U = 100 * sin(2 * π * 50 * t)
dy[1] = u1 - u2
dy[2] = i1 - i2 - i3
dy[3] = i3 - i4 - i5
dy[4] = u4
dy[5] = U - u1 - i1*R1
dy[6] = u2 - i2*R2
dy[7] = u2 - u3 - uc1
dy[8] = u3 - u5 - uc2
dy[9] = u3 - u4 -i4*R3
dy[10] = u5 - i5*R4
dy[11] = i3
dy[12] = i5
return nothing
end
# mass matrix
M = Diagonal([L1, 0, 0, L2, 0, 0, 0, 0, 0, 0, C1, C2])
# Начальные условия для y
y0 = zeros(12)
# Диапазон времени для расчета, пусть будет два периода промышленной частоты 50 Гц
tspan = (0.0, 0.04)
# Задание параметров
paramsDAE = (R1, R2, R3, R4)
# Создание объекта ODEProblem с передачей параметров
fn = ODEFunction(circuitDAEsemiexplicit, mass_matrix = M);
probDAEsemiexplicit = ODEProblem(fn, y0, tspan, paramsDAE);
Для DAEs предпочтительнее использовать неявные схемы интегрирования. Поэтому в Engee разный набор решателей для физических (acasual) и "обычных" (casual) моделей. Здесь Tsit5() даст ошибку, так как не подходит для систем, записанных с mass matrix.
# Решение дифференциальных уравнений
solDAE = solve(probDAEsemiexplicit, Rosenbrock23());
# Выбираем переменные для сравнения
i1_sol = [solDAE.u[k][1] for k in 1:length(solDAE.t)];
i4_sol = [solDAE.u[k][4] for k in 1:length(solDAE.t)];
uc1_sol = [solDAE.u[k][11] for k in 1:length(solDAE.t)];
uc2_sol = [solDAE.u[k][12] for k in 1:length(solDAE.t)];
Добавим код из примера для сравнения результатов.
# Определение дифференциальной системы
function circuit!(du, u, p, t)
U, R1, R2, R3, R4, C1, C2, L1, L2 = p
Uc1, Uc2, i1, i4 = u
# Алгебраические уравнения связи
i2 = (u[1] + u[2] + (u[3] - u[4]) * R4) / (R2 + R4)
i3 = u[3] - i2
i5 = i3 - u[4]
# Дифференциальные уравнения
du[1] = i3 / C1
du[2] = i5 / C2
du[3] = 1 / L1 * (U(t) - i2 * R2 - u[3] * R1)
du[4] = 1 / L2 * (u[2] + i5 * R4 - u[4] * R3)
end
# Начальные условия для u
u0 = [0.0, 0.0 , 0.0, 0.0]
# Задание параметров
U(t) = 100 * sin(2 * π * f * t)
params = (U, R1, R2, R3, R4, C1, C2, L1, L2)
# Создание объекта ODEProblem с передачей параметров
prob = ODEProblem(circuit!, u0, tspan, params)
# Решение дифференциальных уравнений
sol = solve(prob, Tsit5());
# Построение графиков
plot(solDAE.t, [uc1_sol uc2_sol i1_sol i4_sol], label = ["Uc1_mass" "Uc2_mass" "i1_mass" "i4_mass"], layout = 4, lw=3)
plot!(sol, label = ["Uc1_ode" "Uc2_ode" "i1_ode" "i4_ode"], layout = 4, ls=:dash, lw=3)
Вывод
- По графиками мы получили визуальное совпадение. Конечно, полного совпадения не может быть так как, по сути, разные начальные системы уравнений решались разными численными методами;
- Способ с mass matrix при решении semi-explicit DAE of index 1 позволяет сразу решать автоматически сгенерированную систему уравнений без исключения алгебраических переменных.
Литература
[1] Numerical Analysis and Simulation of Ordinary Differential Equations. Roland Pulch