Уравнения для переходных процессов в линейной электрической цепи

В документации Engee есть пример, где описывается применение пакета DifferentialEquations.jl для решения задачи по переходным процессам для электрической схемы с двумя индуктивностями и двумя конденсаторами.

В данном скрипте я бы хотел еще раз обратиться к этой же задаче и решить ее с более общих позиций.

Исходный пример

В исходном примере задача решена методом переменных состояния. В качестве переменных состояния выбраны токи через индуктивности $i_1$, $i_4$ и напряжения на конденсаторах $U_{C1}$, $U_{C2}$. Алгебраические переменные $i_2$, $i_3$, $i_5$ выражены через переменные состояния. Первоначальная дифференциально-алгебраическая система уравнений тем самым сведена к дифференциальной и использована функция 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. Для сокращения длины скрипта воспользуемся только первым способом.

Введем искомую вектор функцию: $$\mathbf{y}=\begin{pmatrix} i_1\\ i_2\\ i_3\\ i_4\\ i_5\\ u_1\\ u_2\\ u_3\\ u_4\\ u_5\\ U_{C1}\\ U_{C2}\\ \end{pmatrix}$$ и запишем систему в виде $$ \mathbf{My^\prime}=f(\mathbf{y}) $$ $$ \begin{pmatrix} L_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & L_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & C_1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & C_2\\ \end{pmatrix} \begin{pmatrix} i^\prime_1\\ i^\prime_2\\ i^\prime_3\\ i^\prime_4\\ i^\prime_5\\ u^\prime_1\\ u^\prime_2\\ u^\prime_3\\ u^\prime_4\\ u^\prime_5\\ U^\prime_{C1}\\ U^\prime_{C2}\\ \end{pmatrix}=\begin{pmatrix} u_1-u_2\\ i_1-i_2-i_3\\ i_3-i_4-i_5\\ u_4\\ u(t)-u_1-i_1R_1\\ u_2-i_2R_2\\ u_2-u_3-U_{C1}\\ u_3-u_5-U_{C2}\\ u_3-u_4-i_4 R_3\\ u_5-i_5 R_4\\ i_3\\ i_5\\ \end{pmatrix} $$ Нулевые строки соответствуют алгебраическим уравнениям, ненулевые - дифференциальным.

Далее я буду переиспользовать код из оригинального примера.

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)

Вывод

1. По графиками мы получили визуальное совпадение. Конечно, полного совпадения не может быть так как, по сути, разные начальные системы уравнений решались разными численными методами;
2. Способ с mass matrix при решении semi-explicit DAE of index 1 позволяет сразу решать автоматически сгенерированную систему уравнений без исключения алгебраических переменных.

Литература

[1] Numerical Analysis and Simulation of Ordinary Differential Equations. Roland Pulch