Подводка

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

При моделировании возникает такое понятие как шины бесконечной мощности (ШБМ). Это шины, напряжение на которых не зависит от тока - в терминах ТОЭ это идеальный источник ЭДС. Эти шины являются опорной точкой и для многих задач, предположение о неизменности напряжения в какой-либо точке не является большой ошибкой и упрощает анализ.

От ШБМ идут линии к другим подстанциям, далее к нагрузке. Так и строится сеть. Инженер может взять модель линии и соединить её с ШБМ. Получится такое соединение (ёмкость с источником ЭДС):

Однако, такое соединение не рекомендуется и даже описано в Modeling Best Practices от Mathworks. Там прямо указано, что такое соединение порождает дифференциально-алгебраические системы уравнений индекса 2 и может приводить к проблемам численного моделирования.

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

Дифференциально-алгебраические системы уравнений (ДАУ)

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

Здесь вектор функция y входит в уравнения со своей производной (это могут быть напряжения на конденсаторах или токи в катушках). А вектор функция z производной не имеет (это остальные токи и напряжения, которые не связаны с запасом энергии). В этой системе из второго уравнения хотя бы в теории возможно выразить z, подставить в первое уравнение и получить ОДУ.

Для схемы подключения конденсатора к источнику ЭДС через резистор уравнения запишутся как:

Здесь u - дифференциальная переменная (напряжение на ёмкости), i - алгебраическая (ток в контуре), e(t) - известный источник ЭДС (функция, таблица и т.д.). Это дифференциально-алгебраическая система уравнений индекса 1.
Если же резистор занулить, то мы получим более сложную систему вида:

где в алгебраическом уравнении нет z. Теоретически можно было бы найти z из первого уравнения, после подстановки туда как y так и der(y). Это совсем другой порядок вычислений.
Для нашего примера уравнения выглядят довольно тривиально, тем не менее их структура именно такая, как описано выше.

Формулы для решения

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

Решим эту СЛАУ
equation (5).png

Выглядит довольно логично. При уменьшении шага интегрирования h значение напряжения, практически перестает меняться u_m+1≈u_m, а ток вычисляется через падение напряжения на резисторе.

Если же занулить резистор (и тут первые намеки на проблемы) то:
equation (6).png

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

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

Начальные условия (скрытые ограничения)

Еще немного вводных.

Решение системы
equation (7).png

требует задания начальных условий (значений в нулевой момент времени) для неизвестных функций u(t), i(t). При решении систем ОДУ эти начальные условия могут быть любыми (как начальные значения интеграторов 1/s. В системах ДАУ появляются нюансы. Например, в приведенной выше системе, хорошо бы, что бы начальные значения тока и напряжения удовлетворяли алгебраическому ограничению (второму уравнению), т.е. были с ним согласованы. Мы ожидаем, что u₀, i₀ выбраны такими, что

В системах индекса 1, мы видим явное ограничение на выбор начальных условий.

В системах индекса 2 ситуация усложняется, тем, что в них кроме явных ограничений, есть еще и скрытые. Например в нашей второй системе:

можно явно указать, что u~0~ должно быть равно e(0), но условия на ток нет. Если же продифференцировать второе уравнение, то можно получить, что согласованным является значение тока: equation (10).png

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

Симуляции

Будем решать обе системы уравнений несколькими способами, чтобы увидеть основные моменты. Все вычисления специально сделаны в single точности, чтобы более явно видеть проблемы. Для жестких задач подохдят методы дифференцирования назад (BDF). Будем использовать метод второго порядка. Также электротехничекие задачи часто решают неявным методом трапеций (тоже второго порядка). Применим и его.

#ПЛАНИРОВАЛОСЬ, ЧТО ВЕСЬ КОД БУДЕТ СКРЫТ МАСКАМИ. ОН НЕ ВАЖЕН, ВАЖНЫ ТОЛЬКО ГРАФИКИ И ТЕКСТ

# @markdown Код функции, производящей численное решение
using LinearAlgebra

function voltage_source(Um, phi, t)
    Float32.(Um * cos.(2*pi*50*t .+ phi));
end

function der_voltage_source(Um, phi, t)
    Float32.(-Um *2*pi*50* sin.(2*pi*50*t .+ phi));
end

# Функция для решения задачи разными методами на разных сетках
function solveDAE(N, r, Tmax, u0, params, method, fun)
    Nr = N * r #количество интервалов на сгущенной сетке
    h = Float32(Tmax / Nr) #шаг интегрирования
    C, R  = params

    current = zeros(Float32, Nr + 1)
    voltage = zeros(Float32, Nr + 1)
    time = zeros(Float32, Nr + 1)

    voltage[1], current[1] = u0

    for k in 2: Nr + 1
        t = (k-1) * h
        t_prev = (k-2) * h
        time[k] = t;

        if method == "ImplicitEuler"
            A = [Float32(1) -h/C;
                Float32(1) R];

            b = [voltage[k-1];
                fun(t)];

        elseif method == "trapz"
            A = [Float32(1) -h/C/2;
                Float32(1) R];

            b = [voltage[k-1] + h/2/C*current[k-1];
                fun(t) + fun(t_prev) - voltage[k-1] - current[k-1] * R];

        elseif method == "BDF2"
            if k == 2
                A = [Float32(1) -h/C;
                    Float32(1) R];

                b = [voltage[k-1];
                    fun(t)];
            else
                A = [Float32(1) -2/3*h/C;
                    Float32(1) R];

                b = [4/3*voltage[k-1] - voltage[k-2]/3;
                    fun(t)];
                    
            end
        end

        solve = A \ b
        voltage[k] = solve[1]
        current[k] = solve[2]
    end

    return voltage[1: r: end], current[1: r: end], time[1: r: end]
end;

# Начальные значения
# @markdown Начальные значения
const Um = Float32(100);
const phi = Float32(pi / 4);
const N = Int32(200);
const Tmax = Float32(0.2);
const R = Float32(1.2);
const C = Float32(1e-2);

const e0 = voltage_source(Um, phi, 0);
const u0 = Float32(0)
const i0 = (e0 - u0) / R;

Индекс 1. Согласованные начальные условия.

Метод трапеций (неявный) очень популярен при электротехнических расчетах. Он не подавляет колебания и поэтому подходит для моделирования схем с осцилляторами. А сети переменного тока и есть схемы с осцилляторами.

vol, cur, tm = solveDAE(N, 1, Tmax, (u0, i0), (C, R), "trapz", t->voltage_source(Um, phi, t));
gr()
plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Согласованные начальные условия. \n Индекс 1. Метод трапеций.",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
# @markdown Решение и построение графиков для метода трапеций
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

Выводы по графику

Шаг в 1 мс хорошо передает колебания;
Нет артефактов и особенностей;
Согласованные начальные условия позволяют графикам начинаться корректно (без изломов).

Посмотрим как изменяется погрешность при сгущении сетки.

# @markdown Функция для расчета погрешности по мере сгущения сетки
function сompute_error(stadium, N, fun)
    err = zeros(Float32, stadium-1);
    steps = zeros(Float32, stadium-1);
    local prevCur
    
    for s in 0: stadium-1
        vol, cur, tm = fun(2^s);

        if s > 0
            err[s] = norm(cur - prevCur);
            steps[s] = 2^s * N;
        end

        prevCur = cur;
    end

    return err, steps
end;

# @markdown Построение графика уменьшения абсолютной погрешности для метода трапеций
err, steps = сompute_error(10, N, x->solveDAE(N, x, Tmax, (u0, i0), (C, R), "trapz", t->voltage_source(Um, phi, t)))
p1 = plot(steps, err,
    title = "Абсолютная погрешность в зависимости \n от количества интервалов",
    label = "trapz",
    xaxis = :log,
    yaxis = :log,
    lw=3,
    marker = :circle,
    markersize = 5,
    # xlims = (1e2, 1e8),
    # ylims = (1e-10, 1e4),
    legend = :bottomleft,
    xlabel = "Количество интервалов",
    ylabel = "Норма разности двух \n последующих решений")

Комментарий к графику

Видно, что метод работает как метод 2-ого порядка вплоть до 5000 шагов на интервал моделирования (наклон в прямолинейной части равен минус 2). Далее точность не увеличивается из-за ошибок округления.

Посмотрим на еще один метод второго порядка основанный на формулах дифференцирования назад. Он в отличие от метода трапеций уже не только A но и L устойчивый.

# @markdown Добавляем на грфик зависимость для метода BDF2
err, steps = сompute_error(10, N, x->solveDAE(N, x, Tmax, (u0, i0), (C, R), "BDF2", t->voltage_source(Um, phi, t)))
plot!(p1, steps, err,
    label = "BDF2",
    lw=3,
    marker = :circle,
    markersize = 5)

Метод BDF2 также работает как метод второго порядка примерно с теми же величинами шагов.

Индекс 2. Согласованные начальные условия

# @markdown Зануляем резистор и пересчитываем начальные условия
const R = Float32(0);
const u0 = e0;
const i0 = C * der_voltage_source(Um, phi, 0)

err, steps = сompute_error(12, N, x->solveDAE(N, x, Tmax, (u0, i0), (C, R), "trapz", t->voltage_source(Um, phi, t)))
p2 = plot(steps, err,
    title = "Абсолютная погрешность в зависимости от \n количества интервалов. Индекс 2.",
    label = "trapz_index2",
    xaxis = :log,
    yaxis = :log,
    lw=3,
    marker = :circle,
    markersize = 5,
    legend = :bottomleft,
    xlabel = "Количество интервалов",
    ylabel = "Норма разности двух \n последующих решений")
    
err, steps = сompute_error(12, N, x->solveDAE(N, x, Tmax, (u0, i0), (C, R), "BDF2", t->voltage_source(Um, phi, t)))
plot!(p2, steps, err,
    label = "BDF2_index2",
    lw=3,
    marker = :circle,
    markersize = 5)

Комментарий к графику

Метод трапеций перестал работать. Мы сгущаем сетку, а результат всё хуже и хуже;
Метод BDF2 сначала (пока шаг интегрирования был не очень маленький) работал корректно и сгущение сетки улучшало результат. Но потом погрешность начала увеличиваться (из-за наличия шага интегрирования в знаменателе численных формул).

vol, cur, tm = solveDAE(N, 128, Tmax, (u0, i0), (C, R), "trapz", t->voltage_source(Um, phi, t));

# @markdown График одного из некорректных решений методом трапеций
plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Некорректное решение ДАУ \n индекса 2 методом трапеций",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

vol, cur, tm = solveDAE(N, 32768, Tmax, (u0, i0), (C, R), "BDF2", t->voltage_source(Um, phi, t));

# @markdown График одного из решений методом BDF2 (при слишком малом шаге)
plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Проявление численных ошибок \n метода BDF2",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

Комментарии к графикам

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

Индекс 1. Несогласованные условия

Сделаем начальные условия нулевыми, хотя ни для индекса 1 ни для индекса 2 это неверно для нулевого момента времени.

# @markdown Построение графика уменьшения абсолютной погрешности для рассматриваемых методов для задачи индекса 1 и несогласованных начальных условиях
const R = Float32(1.2);
err, steps = сompute_error(12, N, x->solveDAE(N, x, Tmax, (0, 0), (C, R), "trapz", t->voltage_source(Um, phi, t)))
p3 = plot(steps, err,
    title = "Абсолютная погрешность в зависимости от \n количества интервалов при несогласованных условиях",
    label = "trapz_index1",
    xaxis = :log,
    yaxis = :log,
    lw=3,
    marker = :circle,
    markersize = 5,
    legend = :bottomleft,
    xlabel = "Количество интервалов",
    ylabel = "Норма разности двух \n последующих решений")
    
err, steps = сompute_error(12, N, x->solveDAE(N, x, Tmax, (0, 0), (C, R), "BDF2", t->voltage_source(Um, phi, t)))
plot!(p3, steps, err,
    label = "BDF2_index1",
    lw=3,
    marker = :circle,
    markersize = 5)

Комментарий к графику

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

vol, cur, tm = solveDAE(N, 1, Tmax, (0, 0), (C, R), "trapz", t->voltage_source(Um, phi, t));

# @markdown график для метода трапеций и несогласованных начальных условиях
plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Незатухающие колебания при несогласованных \n условиях методом трапеций (индекс 1)",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

Индекс 2. Несогласованные условия.

# @markdown Зануляем резистор
const R = Float32(0);

# @markdown График для несогласованных условий, метод трапеций, индекс 2
vol, cur, tm = solveDAE(N, 1, Tmax, (0, 0), (C, R), "trapz", t->voltage_source(Um, phi, t));

plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Нарастающие колебания при \n несогласованных условиях методом трапеций (индекс 2)",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

# @markdown График выбросов. Индекс 2. BDF2. Несогласованные условия.
vol, cur, tm = solveDAE(N*8, 1, Tmax, (0, 0), (C, R), "BDF2", t->voltage_source(Um, phi, t));

plot(tm, [voltage_source(Um, phi, tm), vol, cur],
title = "Выбросы при несогласованных условиях. \n Индекс 2. Метод BDF2.",
lw = 3,
color = ["#0072BD" "#D95319" "#EDB120"],
legend = :bottomright,
xlabel = "Время, с",
ylabel = "Напряжение, В; Ток, А",
labels = ["Источник" "Напряжение на ёмкости" "Ток"])

Комментарии к графикам

Метод трапеций на задачах индекса 1 выдавал колебания, которые не уменьшали амплитуду при уменьшении шага. На задачах индекса 2, эти колебания уже могут увеличиваться приводя к в корне неверным результатам;
Метод BDF2 производит при несогласованных условиях выброс(ы) в начале симуляции. Этот выброс !увеличивается! при уменьшении шага. Это может стать поблемой для решателей с переменным шагом, которые могут захотеть уточнить первое значение, уменьшат шаг и вместо уточнения, получат увеличеине выброса. Такой алгоритм может "споткнуться" с самого старта.

Подводя итог

Задачи индекса 2 гораздо чувствительнее к согласованности начальных условий;
Согласовать эти условия сложно из-за скрытых ограничений на значения переменных или их производных;
В таких задачах появляется численное дифференцирование, которое являетя более неустойчивой опрецией, чем интегрирование. Это означает, что шаг мельчить ради достижения необходимой точности может не получиться;
Решатели с переменным шагом могут застыть на первом шаге из-за выбросов, которые не убрать уменьшением шага.

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

Замечание в документации Engee еще раз об этом предостерегает:

Модель Lumped parameter pi-section имеет параллельный конденсатор на обоих концах. Это означает, что ее нельзя подключать непосредственно к идеальному источнику напряжения, то есть источнику без внутреннего сопротивления.

Компания Mathworks предупреждает

Подводка

Дифференциально-алгебраические системы уравнений (ДАУ)

Формулы для решения

Начальные условия (скрытые ограничения)

Симуляции

Индекс 1. Согласованные начальные условия.

Выводы по графику

Комментарий к графику

Индекс 2. Согласованные начальные условия

Комментарий к графику

Комментарии к графикам

Индекс 1. Несогласованные условия

Комментарий к графику

Индекс 2. Несогласованные условия.

Комментарии к графикам

Подводя итог

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник