CC BY
148
Оригинальная статья / Original article УДК 621.311.001
http://dx.doi.org/10.21285/1814-3520-2018-1 -134-150
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЦИФРОВОЙ ДИФФЕРЕНЦИАЛЬНОЙ ЗАЩИТЫ ТРАНСФОРМАТОРА В СРЕДЕ MATLAB SIMULINK
© М.В. Андреев1, Н.Ю. Рубан2, А.А. Суворов3
Национальный исследовательский Томский политехнический университет, Российская Федерация, 634050, г. Томск, пр-кт Ленина, 30.
РЕЗЮМЕ. ЦЕЛЬ. В данной работе изучен вопрос влияния конкретных видов преобразователей и фильтров в составе измерительной части на функционирование совокупной математической модели цифровой дифференциальной защиты трансформатора. МЕТОДЫ. Для разработки математических моделей защиты использовался метод направленных графов, а для их формализации в виде программного кода - язык программирования Си. В качестве платформы для интеграции созданных моделей и кодов использовался программный комплекс MATLAB Simulink. РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ. Представленные результаты демонстрируют функционирование отдельных элементов математической модели цифровой дифференциальной защиты трансформатора и их влияние на работу защиты в целом. Выявлены некоторые особенности MATLAB Simulink, которые необходимо иметь в виду при реализации математических моделей в данном комплексе. ВЫВОДЫ. Показано, что для полного и достоверного моделирования цифровой релейной защиты необходимо не просто воспроизводить функции структурных элементов, а иметь представление, как минимум, об их видах и характеристиках. В противном случае адекватность математической модели будет низкой. Известная проблема численных методов для решения жестких систем дифференциальных уравнений имеет место даже для таких небольших по математическим меркам моделей, как математические модели релейной защиты.
Ключевые слова: цифровая релейная защита, математическое моделирование, MATLAB Simulink, S-Function Builder, численные методы, измерительная часть.
Формат цитирования: Андреев М.В., Рубан Н.Ю., Суворов А.А. Математическое моделирование цифровой дифференциальной защиты трансформатора в среде Matlab Simulink // Вестник Иркутского государственного технического университета. 2018. Т. 22. № 1. С. 134-150. DOI: 10.21285/1814-3520-2018-1-134-150
MATHEMATICAL SIMULATION OF NUMERICAL TRANSFORMER DIFFERENTIAL PROTECTION
IN MATLAB SIMULINK
M.V. Andreev, N.Yu. Ruban, A.A. Suvorov
National Research Tomsk Polytechnic University, 30 Lenina pr., Tomsk 634050, Russian Federation
ABSTRACT. PURPOSE. The paper studies the influence of specific types of converters and filters in the measuring part on the operation of the mathematical model of numerical transformer differential protection. METHODS. The method of oriented graphs is used for the development of mathematical models of protection, which are formalized in the form of programming codes using the C programming language. The MATLAB Simulink software complex is used as a platform for created models and codes integration. RESULTS AND THEIR DISCUSSION. The presented results demonstrate the operation of the specific elements of the mathematical model of the numerical transformer differential protection and their effect on the operation of protection as a whole. Some features of MATLAB Simulink have been revealed which must be taken into account when implementing mathematical models in this complex. CONCLUSIONS. It is shown that complete and reliable modeling of numerical relay protection requires rather to have an idea on structural element types and characteristics than to reproduce their functions. Otherwise, the adequacy of the mathematical model will be low. The known
Андреев Михаил Владимирович, кандидат технических наук, заведующий научно-исследовательской лабораторией «Моделирование электроэнергетических систем» кафедры электроэнергетических систем, e-mail: andreevmv@tpu.ru
Mikhail V. Andreev, Candidate of technical sciences, Head of the Scientific Research Laboratory of Electrical Power System Simulation, e-mail: andreevmv@tpu.ru
2Рубан Николай Юрьевич, кандидат технических наук, доцент кафедры электроэнергетических систем, e-mail: rubanny@tpu.ru
Nikolai Yu. Ruban, Candidate of technical sciences, Associate Professor of the Department of Electrical Power Systems, e-mail: rubanny@tpu.ru
3Суворов Алексей Александрович, ассистент кафедры электроэнергетических систем, e-mail: suvorovaa@tpu.ru Aleksei A. Suvorov, Assistant Professor of the Department of Electric Power Systems, e-mail: suvorovaa@tpu.ru
problem of the numerical methods for solving the stiff systems of differential equations is relevant even for such small (in mathematical terms) models as mathematical models of relay protection.
Keywords: numerical relay protection, mathematical simulation, MATLAB Simulink, S-Function Builder, numerical methods, measuring part
For citation: Andreev M.V., Ruban N.Yu., Suvorov A.A. Mathematical simulation of numerical transformer differential protection in MATLAB SIMULINK. Proceedings of Irkutsk State Technical University. 2018, vol. 22, no. 1, pp. 134-150. (In Russian) DOI: 10.21285/1814-3520-2018-1-134-150
Введение
Известно, что математическое моделирование является единственным общедоступным способом получения информации по всему спектру процессов в конкретном оборудовании, включая устройства релейной защиты (РЗ) и автоматики, и электроэнергетических системах (ЭЭС) в целом. С появлением и развитием различных средств моделирования множество задач, решение которых ранее было затруднено, становится актуальным. Одна из них - разработка новых методов и средств настройки устройств РЗ, учитывающих современное состояние и перспективы развития ЭЭС. Вывод о такой необходимости сформулирован на основе детального анализа существующих руководящих указаний по расчету уставок РЗ и применяемых для этого вычислительных комплексов (подробности отражены, например, в [1]).
Удовлетворительное решение обозначенной задачи невозможно без детального анализа работы устройств РЗ различных исполнений в конкретных режимах работы и изучения процессов изменения токов и напряжений в защищаемых объектах, измерительных трансформаторах тока (ИТТ) и напряжения (ИТН) и самих РЗ. Получаемая в результате этого информация позволяет формировать параметры адекватной реальным условиям функционирования настройки РЗ.
Обеспечить указанную возможность в совокупности с адекватным симулятором ЭЭС позволяют математические модели РЗ. В публикациях [1, 2] сформулирован подход к разработке детализированных математических моделей РЗ, основное отличие которого от существующих [3-8] заключается в формировании единого математического описания (в виде передаточной
функции (ПФ) и соответствующей системы дифференциальных уравнений) одновременно для всей совокупности элементов, входящих в состав конкретной РЗ: измерительных трансформаторов, промежуточных преобразователей, частотных фильтров, реле и т.д. При этом исключается ошибка при некорректном или неполном учете взаимосвязей конкретного функционального элемента РЗ, поскольку все связи учитываются естественным образом при моделировании схемы целиком.
Постановка проблемы. Очевидным является неизбежный переход в ЭЭС на цифровые системы защиты, в связи с чем логичным является решение указанной ранее проблемы применительно к цифровой РЗ (ЦРЗ). В качестве конкретного объекта исследований выбрана цифровая дифференциальная защита трансформаторов (ЦДЗТ), поскольку ЦДЗТ относится к разряду основных защит и представляет значительный интерес в плане изучения.
Известно, что наибольшие погрешности преобразования входного сигнала формируются в измерительной части (ИЧ) ЦДЗТ, а соответственно адекватность функционирования модели защиты во многом зависит от полноты и достоверности моделирования элементов ИЧ. В работе [2] представлены результаты исследования интегральных частотных характеристик математических моделей ИЧ для различных сочетаний элементов в ее составе: ИТТ, промежуточных трансформаторов тока (ПТТ), аналоговых фильтров низких частот (ФНЧ), показавшее их конкретное влияние на результаты моделирования процессов в ИЧ в целом. В рамках других исследований (например, [9]) продемонстрировано, что для решения даже таких несложных по ма-
тематическим меркам моделей не все численные методы применимы. Использование наиболее простых численных методов (явный метод Эйлера, усовершенствованный (неявный) метод Эйлера) приводило к существенным искажениям выходного сигнала ИЧ, что в свою очередь было связано с постоянным перемагничиванием ИТТ, а соответственно изменением на каждом шаге интегрирования коэффициентов в дифференциальных уравнениях, формирую-
щих модель ИЧ. Применение более совершенного метода (неявный метод Рунге -Кутта 4-ого порядка) позволило в значительной степени снизить указанные искажения.
Проведенные ранее исследования выполнены для отдельных ИЧ. Целью же данной работы являлось изучение функционирования совокупной математической модели ЦДЗТ с учетом выявленных ранее особенностей.
Материалы и методы исследования
Математические модели ЦДЗТ реализованы в программном комплексе MATLAB Simulink. Для решения дифференциальных уравнений, описывающих процессы в ИЧ, использовался неявный метод Рунге - Кутта 4-го порядка. Программы решения, написанные на языке программирования Си, интегрированы в блоки S-Function Builder. Модель ЦДЗТ составлена для случая трехобмоточного трансформатора.
В качестве входного использовался сигнал, отражающий квазиустановившийся (0-0,2 с), аварийный (0,2-0,8 с) и после-аварийный (0,8-1 с) режимы со следующими характеристиками:
1) для нормального и послеаварий-ного режимов: амплитуда - 0,6 кА (сторона ВН), 0,4 кА (сторона СН), 7,2 кА (сторона НН), частота - 50 Гц;
2) для аварийного режима: постоянная времени затухания переходного процесса - Ta = 0,1 с, ударный ток - 10 кА (сторона ВН), 10 кА (сторона СН), 15 кА (сторона НН), установившийся ток КЗ -1,5 кА (сторона ВН), 1 кА (сторона СН), 9 кА (сторона НН), состав гармоник относительной основной: 2-я - 10%; 3-я - 5%; 4-я - 3%; 5-я - 2%.
Математические модели ЦДЗТ созданы для следующих сочетаний элементов ИЧ:
- ИТТ-активный ПТТ-ФНЧ Баттер-ворта 3-го порядка;
- ИТТ-активный ПТТ-ФНЧ Чебыше-ва 3-го порядка;
- ИТТ-пассивный ПТТ-ФНЧ Баттер-ворта 3-го порядка;
- ИТТ-пассивный ПТТ-ФНЧ Чебы-шева 3-го порядка.
Основой для создания соответствующих ПФ и образующей их системы дифференциальных уравнений являются схема замещения и направленный граф [1, 2].
Структура и функции преобразовательной (ПЧ) и вычислительно-логической частей (ВЛЧ) ЦДЗТ полностью совпадают для всех сочетаний. Параметры ИТТ взяты из [10] (коэффициенты трансформации соответственно для ИТТ сторон ВН, СН и НН приняты 600/5, 1000/5, 12000/5 соответственно). Параметры ФНЧ для частоты среза 250 Гц [11] определены в соответствии с методиками, представленными в [12]. В качестве активного ПТТ выбран трансформатор тока LTS 15^Р [13]. Пассивный ПТТ отличается от активного отсутствием операционного усилителя на выходе, замененного резистором с сопротивлением, подобранным таким образом, чтобы обеспечить максимальную близость амплитудно-частотных характеристик активного и пассивного ПТТ.
Схема моделирования совокупности «ИТТ-ПТТ-ФНЧ» с учетом характеристики намагничивания представлена на рис. 1, а. Особенность данной схемы заключается в параллельном расчете двух математических моделей - ИТТ, нагрузкой которого является ПТТ и ФНЧ, и совокупности «ИТТ-ПТТ-ФНЧ». Первая модель используется для определения тока намагничивания и
соответственно сопротивления цепи намагничивания в соответствии с характеристикой. Алгоритм реализации характеристики намагничивания с учетом гистерезиса представлен на рис. 1, Ь. Расчет и на его основе осуществлялся по указанному алгоритму на каждом шаге интегрирования. Вторая модель является основной и имеет два входных сигнала - первичный ток /1 и сопротивление 2^, полученное из расчета первой модели. Такая схема необходима, поскольку априори выходным сигналом в любом узле схемы, рассчитываемым с помощью ПФ, является напряжение. Очевидно, что для получения вторичного тока ИТТ, который в свою очередь необходим для реализации алгоритма, представленного на рис. 1, Ь, нужно иметь ток и на входе, и на выходе. Это достигается делением выходного напряжения ИТТ на сопротивление
нагрузки, которое представляет собой ПФ «ПТТ-ФНЧ», составленной при условии тока на входе и напряжения на выходе.
Исходная характеристика намагничивания ИТТ в виде зависимости магнитной индукции от напряженности магнитного поля б=/(Н) для стали 3413, заданная таблично и изначально содержащая 90 точек, представлена на рис. 2. Поскольку математическая модель ИЧ основана на схеме замещения [1, 2, 9], то для учета нелинейности использовалась иная форма характеристики намагничивания в виде зависимости сопротивления цепи намагничивания от тока намагничивания 2у,=!(1у). С использованием метода кусочно-линейной аппроксимации число точек было увеличено до 46000 для обеспечения плавного изменения сопротивления намагничивания ИТТ.
b
Рис. 1. A - схема моделирования совокупности «ИТТ-ПТТ-ФНЧ» с учетом характеристики намагничивания; b - алгоритм реализации характеристики намагничивания Fig. 1. A - Simulation flow diagram of the measuring part "MCT -ICT - LFF" (measuring current transformer (MCT), intervening current transformer (ICT), low frequency filter (LFF) taking into account magnetization characteristic; b - algorithm of magnetization characteristic implementation
а
Рис. 2. Семейство петель гистерезиса Fig. 2. Family of hysteresis loops
Для учета в модели явления гистерезиса были сформированы три петли в результате анализа диапазона изменения тока намагничивания в соответствии с начальной характеристикой и условного деления его на три участка: нормальный режим работы (ряд 2 и ряд 3 на рис. 2), установившееся КЗ без насыщения ИТТ (ряд 4 и ряд 5 на рис. 2), насыщение ИТТ (ряд 6 и ряд 7 на рис. 2 - предельная петля гистерезиса).
Как известно, процесс аналого-цифрового преобразования (АЦП) состоит из двух подпроцессов: дискретизация по времени и квантование по уровню. Частота дискретизации (1200 Гц = 50*24) рассчитана исходя из значения базовой частоты (50 Гц) и числа выборок АЦП (принято 24) и удовлетворяет требованию теоремы Ко-тельникова.
Алгоритм ЦДЗТ работает с векторными величинами, для формирования которых необходима предварительная цифровая фильтрация с частотой среза, равной базовой частоте - 50 Гц. Для формирования самих векторов использовался метод, отраженный на рис. 3. Первостепенно определяется период T изменения преобразуемого сигнала путем контроля перехода сигнала через ноль сначала в точке 1, а затем в точке п: T = tn - и. Для каждой переменной величины рассчитывается среднеквадратичное значение Аск за период изменения сигнала. Для определения фазы вектора генерируется независимый синхросигнал с частотой 50 Гц. Фаза вектора представляет собой отношение времени от первого перехода через ноль в точке 1 до синхросигнала в точке 2, т.е. времени tA к времени полного периода Т.
Рис. 3. Формирование амплитуды и фазы вектора Fig. 3. Calculation of vector amplitude and phase
Результаты исследования и их обсуждение
На рис. 4 представлены входные сигналы ЦДЗТ для трех сторон трансформатора.
На рис. 5-8 представлены амплитуды и фазы векторов, полученные на основе выходных токов ИЧ для четырех указанных ранее сочетаний. На рисунках можно заметить, что расчет амплитуд и фаз начался не в 0 с, а примерно в 0,045 с, т.е. спустя два периода. По факту расчет векторов начался с первого перехода через 0 кривой тока, что произошло примерно в момент времени 0,005 с. Когда первый период сигнала закончился, сформировалось среднеквадратичное значение амплитуды вектора тока. Однако блок S-Function Builder, рассчитывающий амплитуду и фазу вектора, выдал это значение только в начале нового периода, что дало смещение 0,02 с. Для расчета фазы нужно иметь две метки времени, поэтому сдвиг происходит дважды по одной и той же причине - задержка блока S-Function Builder. Для синхронизации кривых фазы и амплитуды значение амплитуды вектора сдвинуто искусственно на 0,02 с вперед, что связано с отсутствием возможности при моделировании в MATLAB Sim-ulink сдвига сигнала назад во времени. Если данную синхронизацию не провести, то заведомо появится ошибка при расчете дифференциального и тормозного токов с использованием указанных векторов в реагирующем органе ЦДЗТ. Следует отметить, что при программной реализации не в MATLAb от этого эффекта можно избавиться, обратившись к конкретным членам массива данных, чего не получается сделать при работе с блоками S-Function Builder. На рис. 9 показано то, как формируется амплитуда и фаза вектора.
Неодновременность формирования амплитуд и фаз векторов токов с разных сторон силового трансформатора связана с неодинаковым сдвигом по фазе кривых токов после преобразования входного сигна-
ла в измерительной части.
Следует отметить, что корректировки сигналов не производилось, потому что входные сигналы были подобраны так, чтобы идеально компенсировать друг друга, и, соответственно, необходимость корректировки априори исключалась.
Проводилось вычисление: - дифференциального тока:
1d _\ 1A ВН + 1A СН + 1 A НН
- тормозного тока:
Is 1 1 A ВН 1 + 1 1ACE 1 + 1 1A НН \
Сформированные дифференциальный и тормозной токи представлены на рис. 10 и 11.
Из рис. 10, 11 видно, что различные варианты ИЧ совершенно по-разному влияют на результаты моделирования выходного сигнала. В предыдущих работах, например в [2], обращено внимание на существенное различие интегральных частотных характеристик ИЧ. Графики, представленные на рис. 10, 11, визуально подтверждают то, как это различие сказывается на формировании контрольных токов реагирующего органа ЦДЗТ. Рис. 12 показывает сформированные в итоге из этих токов аварийные характеристики /диф=/г (/торм) для каждого из указанных ранее сочетаний элементов ИЧ. Расчет характеристики срабатывания (ХС) ЦДЗТ сделан на основе руководящих указаний Siemens для SIPROTEC 7UT6 [14, 15]. Результаты сведены в таблицу.
Рис. 4. Сигналы на входе ЦДЗТ Fig. 4. Input signals of numerical transformer differential protection
Рис. 5. Амплитуды и фазы векторов на основе выходных токов ИЧ «ИТТ-активный ПТТ-ФНЧ Баттерворта 3-его порядка» Fig. 5. Amplitudes and phases of vectors based on output currents of the measuring part "MCT-active ICT-Butterworth 3rd-order LFF"
Рис. 6. Амплитуды векторов на основе выходных токов ИЧ «ИТТ-активный ПТТ-ФНЧ Чебышева 3-его порядка» Fig. 6. Amplitudes of vectors based on output currents of the measuring part "MCT- active ICT-Chebyshev 3rd-order LFF"
Рис. 7. Амплитуды векторов на основе выходных токов ИЧ «ИТТ-пассивный ПТТ-ФНЧ Баттерворта 3-его порядка» Fig. 7. Amplitudes of vectors based on output currents of the measuring part "MCT- passive ICT-Butterworth 3rd-order LFF"
m Энергетика
Ses Power Engineering
Рис. 8. Амплитуды векторов на основе выходных токов ИЧ «ИТТ-пассивный ПТТ-ФНЧ Чебышева 3-его порядка» Fig. 8. Amplitudes of vectors based on output currents of the measuring part "MCT- passive ICT-Chebyshev 3rd-order LFF"
I, I A / il llll Ш|,|, /
JjjJJjJj^r-iiJ llll / 1 ЛЛАЛЛЛЛЙЛЛЛЛЛЛЛЛЛА 1 = =r
1! |!ШЩШ llllllllr 14 u i 11 111 I 1 I I t, с
0.1 0.3 0.5 0.7 0.9
Рис. 9. Формирование амплитуды и фазы вектора в S-Function Builder ИЧ «ИТТ-активный ПТТ-ФНЧ Баттерворта 3-его порядка» Fig. 9. Formation of vector amplitude and phase in S-Function Builder of the measuring part "MCT-active ICT-Butterworth 3rd-order LFF"
Рис. 10. Дифференциальные токи для различных сочетаний ИЧ Fig. 10. Differential currents for different types of measuring part combination
Рис. 11. Тормозные токи для различных сочетаний ИЧ Fig. 11. Restraint currents for different combinations of the measuring part
«ИТТ-активный ПТТ-ФНЧ Чебышева Э-ега порядка» / «ИТТ-пассивный ГПТ-ФНЧ Чебышева 3-вго порядка» /
-МСТ- active ICT-Chebyshev 3rd-order LFF' 'MCT- passive ICT-Chebyshev 3rd-ordef LFP
Рис. 12. Аварийные характеристики, зафиксированные реагирующими органами ЦДЗТ
с названными ранее ИЧ Fig. 12. Fault characteristics measured by the tripping elements of numerical transformer differential
protection with the measuring part named earlier
Уставки ЦДЗТ Settings of numerical transformer differential protection
Параметр / Parameter Значение / Value
Ток срабатывания / Tripping current (IVCT), о.е. 0,2
Ток срабатывания дифференциальной отсечки / Tripping current of a differential cutoff (!уст ДО), о.е. 7,5
Основание первой ветви характеристики срабатывания / Base of the first branch of the tripping characteristic (Iхс1), о.е. 0
Наклон первой ветви характеристики срабатывания / Slope of the first branch of the tripping characteristic (Кхс1), о.е. 0,25
Основание второй ветви характеристики срабатывания / Base of the second branch of the tripping characteristic (Iхс2), о.е. 2,5
Наклон второй ветви характеристики срабатывания / Slope of the second branch of the tripping characteristic (Кхсг), о.е. 0,5
Начало зоны дополнительного торможения / Start of the additional stabilization zone (!доп), о.е. 5
Наклон ветки, определяющей зону дополнительного торможения / Slope of the branch determining the additional stabilization zone (Кдоп), о.е. 0,15
Доля второй гармоники / Percentage of the second harmonic (I2), % 15
Доля третьей (пятой) гармоник / Percentage of the third (fifth) harmonic (I3 (I5)), % 20
Вид ХС представлен на рис. 13. К завести на него несколько сигналов, поэто-сожалению, возможности блока XYplot из му зона дополнительного торможения ХС библиотеки MATLAB Simulink не позволяют вынесена отдельно.
b
Рис. 13. Характеристика срабатывания ЦДЗТ: а - основная характеристика; b - зона дополнительного торможения Fig. 13. Tripping characteristic of numerical transformer differential protection: a - main characteristic; b - additional stabilization zone
На рис. 14 представлены сигналы, сформированные реагирующими органами ЦДЗТ для каждого из указанных ранее сочетаний элементов ИЧ.
Как видно на рис. 14, не все реагирующие органы сформировали сигналы на срабатывание - не сработал вариант с ИЧ «ИТТ-активный ПТТ-ФНЧ Баттерворта 3-го порядка». Время выдачи команды на срабатывание у остальных сильно различается. Отчасти это связано с различием формируемых аварийных характеристик,
отчасти - с работой блока, реализующего блокировку по гармоническим составляющим, основанного на алгоритме быстрого преобразования Фурье, который наиболее часто применяется в ЦРЗ. Функционирование блока для случая ИЧ «ИТТ-активный ПТТ-ФНЧ Баттерворта 3-го порядка» отображено на рис. 15. Для сравнения на рис. 16 представлена работа аналогичного блока, составленного с использованием элемента Fourier из библиотеки MATLAB Simulink.
Рис. 14. Сигнал о срабатывании (логическая 1) / несрабатывании (логический 0) ЦДЗТ для различных сочетаний ИЧ Fig. 14. Tripping (logical 1) / non-tripping (logical 0) signal of numerical transformer differential protection
for different combinations of the measuring part
а
Рис. 15. Доли гармонических составляющих относительно основной гармоники (сверху) на основе мгновенных значений дифференциального тока (снизу) Fig. 15. Percentage of harmonic components relative to the fundamental harmonic (at the top) on the basis of the instantaneous values of the differential current (at the bottom)
Рис. 16. Доли гармонических составляющих, полученные с помощью блока Fourier,
относительно основной гармоники Fig. 16. Percentage of harmonic components obtained by means of the Fourier block relative to the fundamental harmonic
Как видно из рис. 15, доли второй, третьей и пятой гармонических составляющих, полученные относительно первой гармоники и рассчитанные на основе кривой дифференциального тока, значительно превышают уставки в 15, 20 и 20% соответственно. Это приводит к блокировке срабатывания ЦДЗТ. Аналогичная ситуация сло-
жилась и при использовании стандартного блока быстрого преобразования Фурье Fourier из библиотеки MATLAB Simulink. Получившаяся форма дифференциального тока объясняется следующими основными факторами: 1) сохранившимися искажениями формы выходного тока ИЧ вследствие численного интегрирования дифференци-
альных уравнений; 2) нелинейностью ИТТ в целом и существенным различием магнитных характеристик ИТТ с разных сторон силового трансформатора. Если второй фактор ожидаем и не является проблемой моделирования, то первый связан с вопросом применимости численных методов [16-18]. Авторами настоящей статьи предложен вариант решения данной проблемы, указанный, например, в [9], заключающийся в использовании метода неявного интегрирования аналоговым способом системы дифференциальных уравнений, образующих модель РЗ. В настоящее время этот вариант уже взят в разработку. Различие в работе созданного блока вычисления гармоник и стандартного блока МАНАБ 81т-и!тк становится понятным после рассмотрения последнего «под маской» (рис. 17).
В созданном блоке расчет гармоник осуществляется один раз в период, который определяется, как показано на рис. 3. Задержка в вычислении связана с особенностями МАНАБ 81ти!1пк, отмеченными ранее при описании блока формирования векторов. Данный метод также работает с первого перехода сигнала через ноль. В библиотечном блоке, как видно из рис. 17, сначала мгновенные значения входного сигнала перемножаются с искусственно
сформированными кривыми, изменяющимися по закону синуса (sin(wt)) и косинуса (cos(wt)), для формирования в дальнейшем амплитуды и фазы вектора (или действительной и мнимой частей в случае записи вектора в комплексной форме). Блок Mean осуществляет вычисление среднего арифметического за период, поэтому при использовании Fourier вначале априори будет появляться пауза, равная периоду, обратно пропорциональному частоте основной гармоники. Вычислив количество точек за заданный период (поскольку сигнал дискретный), далее по ходу изменения сигнала блок Mean осуществляет пересчет среднего арифметического исходя из полученного фиксированного числа точек. По мнению авторов настоящей статьи, этот подход не совсем корректен, поскольку при большой доле высокочастотных составляющих возможны частые переходы через ноль и, соответственно, отличие фактического периода изменения сигнала от рассчитанного на первом этапе. Написанный программный код, реализующий алгоритм быстрого преобразования Фурье, каждый раз заново рассчитывает период сигнала и определяет на нем доли гармонических составляющих, а значит, указанный недостаток исключен.
Рис. 17. Блок Fourier «под маской» Fig. 17. Fourier block "under mask"
Выводы
1. Ранее проведенные исследования [2, 9] позволили определить конкретное влияние отдельных элементов (ИТТ, ПТТ и ФНЧ) на результаты моделирования ИЧ в целом. Тем не менее требовало проверки предположение о том, что их влияние в зависимости от конкретного вида ПТТ и ФНЧ по-разному скажется на функционировании защиты в целом. Проведенные и представленные в данной статье исследования это доказали. Таким образом, для достаточно полного и достоверного моделирования ЦРЗ необходимо не просто воспроизводить функции структурных элементов, как это часто делается [5, 7, 8], а знать хотя бы их вид.
2. Проблема применимости численных методов [16-18] имеет место для решения таких небольших (по математическим меркам) математических моделей. Это подтверждается возникающими искажениями моделируемых сигналов, даже при использовании такого достаточно устойчивого метода, как неявный метод Рунге - Кутты 4-го порядка. Данное явление связано с частыми (на каждом шаге интегрирования) и иногда значительными (при КЗ) изменениями сопротивления цепи намагничивания ИТТ. Учет гистерезиса еще больше усугубляет ситуацию. Кусочно-линейная аппроксимация петель гистере-
зиса позволяет несколько сгладить искажения сигнала, но главным образом в установившемся режиме, когда ток намагничивания изменяется плавно. Предположительно эту проблему можно снять при неявном интегрировании аналоговым способом, что является задачей планируемых исследований.
3. Выявлены некоторые особенности MATLAB Simulink, которые необходимо учитывать при моделировании ЦРЗ: задержки, возникающие при использовании программируемых блоков; невозможность сдвига сигнала назад во времени (по факту входные и выходные сигналы S-Function Builder - это массивы данных, но обратиться к конкретному члену этого массива и задать его невозможно), что также создает видимость задержек в работе модели; несовершенство некоторых библиотечных элементов, которые необходимы для моделирования ЦРЗ. Тем не менее MATLAB Simulink позволил выполнить запланированные исследования и определить конкретные требования к математическим моделям ЦРЗ и средствам их реализации.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 16-38-60043 мол_а_дк.
Библиографический список
1. Андреев М.В., Рубан Н.Ю., Гордиенко И.С., Боровиков Ю.С., Гусев А.С., Сулайманов А.О. Всережим-ное математическое моделирование релейной защиты электроэнергетических систем. Томск: Изд-во Томского политехнического университета, 2016. 180 с.
2. Андреев М.В., Боровиков Ю.С., Гусев А.С., Рубан Н.Ю., Суворов А.А., Сулайманов А.О., Уфа Р.А. Влияние элементов измерительной части цифровой дифференциальной защиты трансформаторов на преобразование входного сигнала // Проблемы региональной энергетики. 2017. № 2 (34). С. 1-12.
3. Ершов Ю.А., Малеев А.В. Моделирование микропроцессорных релейных защит в среде MATLAB // Журнал Сибирского федерального университета. Серия: Техника и технологии. 2010. Т. 3. № 2.
С. 220-228.
4. Румянцев Ю.В. Комплексная модель для исследования функционирования цифровой дифференциальной защиты силового трансформатора // Энергетика. Известия высших учебных заведений и энергетических объединений СНГ. 2016. Т. 59. № 3. С. 203-224.
5. Ершов Ю.А., Киселев Д.Н. Исследование цифровых дифференциальных защит трансформаторов // Приоритетные научные направления: от теории к практике. 2016. № 34-1. С. 176-185.
6. Chen C.-S., Liu C.-W., Jiang J.-A. Application of combined adaptive Fourier filtering technique and fault detector to fast distance protection // IEEE Transactions on Power Delivery. 2006. Vol. 21. No. 2. P. 619-626.
7. Perez S.G.A., Sachdev M.S., Sidhu T.S. Modeling
relays for use in power system protection studies // Proc. Electrical and Computer Engineering Conf. 2005. P. 566-569.
8. Aktaibi A., Rahman M.A. MATLAB - A Fundamental Tool for Scientific Computing and Engineering Applications. Vol. 1. Chapter 10. DOI: http://dx.doi.org/10.5772/48624
9. Андреев М.В. Особенности математического моделирования цифровых устройств релейной защиты // Автоматизация и IT в энергетике. 2017. № 11 (100). С. 10-15.
10. Афанасьев В.В., Адоньев Н.М., Жалалис Л.В. и др. Трансформаторы тока. Л.: Энергия, 1980. 344 с.
11. Овчаренко Н.И. Цифровые аппаратные и программные элементы микропроцессорной релейной защиты и автоматики энергосистем. Библиотечка электротехника, приложение к журналу «Энергетик». Вып. 5-6 (89-90). М.: НТФ «Энергопрогресс», 2006. 120 с.
12. Johnson D., Johnson J., Moore H. A handbook of active filters, Englewood Cliffs. New Jersey, USA: Pren-
tice-Hall Inc. 1980. 128 p.
13. Овчаренко Н.И. Аналоговые элементы микропроцессорных комплексов релейной защиты и автоматики. Библиотечка электротехника, приложение к журналу «Энергетик». Вып. 9 (33). М.: НТФ «Энергопрогресс», 2001. 80 с.
14. Manual for Differential Protection 7UT6x V4.6, Siemens, 2014.
15. Instruction Manual for Numerical Differential Protection Relay for Transformers, Generators, Motors and Branches Points 7UT51 V3.0, Siemens, 2001.
16. Hall G., Watt J.M. Modern numerical methods for ordinary differential equations. Oxford, UK: Clarendon Press, 1976. 312 p.
17. Галанин М.П., Ходжаева С.Р. Методы решения жестких обыкновенных дифференциальных уравнений. Результаты тестовых расчетов // Препринты ИПМ им. М.В. Келдыша. 2013. № 98. 29 с.
18. Вержбицкий В.М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения). М.: Высш. шк. 2001. 382 с.
References
1. Andreev M.V., Ruban N.Yu., Gordienko I.S., Bo-rovikov Yu.S., Gusev A.S., Sulaimanov A.O. Vse-rezhimnoe matematicheskoe modelirovanie releinoi zashchity elektroenergeticheskikh sistem [All-regimes simulation of relay protection in electrical power systems]. Tomsk: Tomsk Polytechnic University Publ., 2016, 180 p. (In Russian)
2. Andreev M.V., Borovikov Yu.S., Gusev A.S., Ruban N.Yu., Suvorov A.A., Sulaimanov A.O., Ufa R.A. Impact of Measuring Part Elements of Transformer Differential Protection on Input Signal Processing. Problemy re-gional'noi energetiki [Problems of the regional energetics]. 2017, no. 2 (34), pp. 1-12. (In Russian)
3. Ershov Yu.A., Maleev A.V. Microprocessor relay protection modeling in MATLAB. Zhurnal Sibirskogo feder-al'nogo universiteta. Seriya: Tekhnika i tekhnologii [Journal of Siberian Federal University. Series "Engineering and Technologies"]. 2010, vol. 3, no. 2, pp. 220-228. (In Russian)
4. Rumyantsev Yu.V. A comprehensive model for the power transformer digital differential protection functioning research. Energetika. Izvestiya vysshikh uchebnykh zavedenii i energeticheskikh ob"edinenii SNG [Energetika. Proceedings of CIS higher education institutions and power engineering associations]. 2016, vol. 59, no. 3, pp. 203-224 (In Russian)
5. Ershov Yu.A., Kiselev D.N. Study of numerical transformer differential protection. Prioritetnye nauchnye napravleniya: ot teorii k praktike [Priority research areas: from theory to practice]. 2016, no. 34-1, pp. 176-185. (In Russian)
6. Chen C.-S., Liu C.-W., Jiang J.-A. Application of combined adaptive Fourier filtering technique and fault detector to fast distance protection. IEEE Transactions on power delivery. 2006, vol. 21, no. 2, pp. 619-626.
7. Perez S.G.A., Sachdev M.S., Sidhu T.S. Modeling relays for use in power system protection studies. Proc. Electrical and Computer Engineering Conf. 2005, pp. 566-569.
8. Aktaibi A., Rahman M.A. MATLAB - A Fundamental Tool for Scientific Computing and Engineering Applications. Vol. 1. Chapter 10. DOI: http://dx.doi.org/10.5772/48624
9. Andreev M.V. Features of mathematical modeling of digital relay protection devices. Avtomatizatsiya i IT v energetike [Automation and IT in power engineering]. 2017, no. 11 (100), pp. 10-15. (In Russian)
10. Afanas'ev V.V., Adon'ev N.M., Zhalalis L.V. i dr. Transformatory toka [Current transformers]. Leningrad: Energiya Publ., 1980, 344 p. (In Russian)
11. Ovcharenko N.I. Tsifrovye apparatnye i pro-grammnye elementy mikroprotsessornoi releinoi zash-chity i avtomatiki energosistem. [Digital hardware and software elements of power system microprocessor relay protection and automation]. Bibliotechka el-ektrotekhnika, prilozhenie k zhumalu "Energetik" [Electrical engineer library, supplement to the journal "Power Engineer"]. Issue 5-6 (89-90). Moscow: NTF «Ener-goprogress» Publ., 2006, 120 p. (In Russian)
12. Johnson D., Johnson J., Moore H. A handbook of active filters, Englewood Cliffs. New Jersey, USA: Prentice-Hall Inc. 1980, 128 p.
13. Ovcharenko N.I. Analogovye elementy mikro-protsessornykh kompleksov releinoi zashchity i avtomatiki [Analog elements of microprocessor complexes of relay protection and automation]. Bibliotechka elektrotekhnika, prilozhenie k zhurnalu "Energetik" [Electrical engineer library, supplement to the journal "Power Engineer"]. Issue 9 (33), Moscow: NTF "Ener-goprogress" Publ., 2001, 80 p. (In Russian)
14. Manual for Differential Protection 7UT6x V4.6, Siemens, 2014.
15. Instruction Manual for Numerical Differential Protection Relay for Transformers, Generators, Motors and Branches Points 7UT51 V3.0, Siemens, 2001.
16. Hall G., Watt J.M. Modern numerical methods for ordinary differential equations. Oxford, UK: Clarendon Press, 1976, 312 p.
17. Galanin M.P., Khodzhaeva S.R. Methods of solving
Критерии авторства
Авторы заявляют о равном участии в получении и оформлении научных результатов и в равной мере несут ответственность за плагиат.
Конфликт интересов
Авторы заявляют об отсутствии конфликта интересов.
Статья поступила 30.11.2017 г.
stiff ordinary differential equations. Results of test calculations. Preprinty IPM im. M.V. Keldysha. [Keldysh Institute preprints]. 2013, no. 98, 29 p. (In Russian) 18. Verzhbitskii V.M. Chislennye metody (matematich-eskii analiz i obyknovennye differentsial'nye uravneniya) [Numerical methods (mathematical analysis and ordinary differential equations)]. Moscow: Vysshaya shkola Publ., 2001, 382 p. (In Russian)
Authorship criteria
The authors declare equal participation in obtaining and formalization of scientific results and bear equal responsibility for plagiarism.
Conflict of interests
The authors declare that there is no conflict of interests regarding the publication of this article.
The article was received 30 November 2017
