CC BY
15
ПОСТРОЕНИЕ И АНАЛИЗ КВАДРАТИЧНОЙ АППРОКСИМАЦИИ НЕЛИНЕЙНОЙ МОДЕЛИ ДВУХРАЙОННОЙ ЭЛЕКТРОЭНЕРГЕТИЧЕСКОЙ СИСТЕМЫ В ПРОСТРАНСТВЕ СОСТОЯНИЙ1
Кутяков Е. Ю.2
(ФГБУН Институт проблем управления им. В.А. Трапезникова РАН, Москва)
Представлены результаты математического моделирования двухрайонной четырёхгенераторной электроэнергетической системы (ЭЭС) в пространстве состояний с учётом членов второго порядка малости. Разработано программное обеспечение для построения квадратичных аппроксимаций указанной энергосистемы. Исходная модель ЭЭС основана на нелинейной модели синхронного генератора 6-го порядка, которая объединяет электромеханические уравнения движения ротора, динамику электродвижущей силы (ЭДС) обмотки возбуждения, а также динамику трёх демпферных обмоток. Автоматические регуляторы возбуждения и системные стабилизаторы в модели не учитываются. Для построения квадратичной аппроксимации используется метод линеаризации Карлемана, который широко применяется при представлении моделей нелинейных динамических систем в билинейном виде. Приводится подробный анализ спектра матрицы динамики полученной квадратичной аппроксимации исследуемой ЭЭС. Изложены результаты численных экспериментов, в которых продемонстрировано, что по сравнению с традиционно использующейся линеаризованной моделью разработанная квадратичная аппроксимация даёт качественно новую информацию о переходных процессах исходной ЭЭС, возникающих как при отклонении отдельных переменных состояния, так и при возмущении через входы системы. Полученная модель может использоваться для аналитического исследования нелинейного поведения ЭЭС вблизи её рабочей точки, в частности, для анализа межрайонных колебаний и параметрических резонансов. Отметим, что применённый в данной работе метод построения квадратичной аппроксимации (линеаризация Карлемана) может быть распространён и на ЭЭС более крупных размеров, а разработанное программное обеспечение допускает модернизацию на случай учёта членов третьего и более высоких порядков малости.
Ключевые слова: нелинейный анализ ЭЭС, нелинейная модель ЭЭС, линеаризация Карлемана, межрайонные колебания, нелинейные модальные взаимодействия, ряд Тейлора.
1 Работа выполнена при финансовой поддержке РНФ (Проект №19-19-00673).
2 Евгений Юрьевич Кутяков, н.с. (evgeniykutyakov@gmail.com).
Управление большими системами. Выпуск 100 1. Введение
Хорошо известно, что крупные ЭЭС демонстрируют существенно нелинейное поведение при отклонениях от установившегося режима, одним из проявлений которого являются резо-нансы различной природы, а также межрайонные колебания -низкочастотные электромеханические колебания, как правило слабо демпфированные и возникающие между удалёнными друг от друга группами когерентных генераторов [1, 28]. Классические средства анализа устойчивости ЭЭС «в малом», такие как модальный анализ [19], используют линеаризованные модели и не всегда пригодны для исследований указанных явлений, поскольку последние связаны с нелинейными модальными взаимодействиями и проявляются при больших возмущениях, ввиду чего их анализ требует более точных моделей, которые бы учитывали не только линейные члены, но и члены более высокого порядка малости [23].
В промышленной практике нелинейные модели ЭЭС, которые достаточно точно описывают реальные переходные процессы, обычно используют при исследовании переходной устойчивости с применением численных методов интегрирования. В ходе этих исследований выполняется построение динамического отклика системы во временной области при воздействии на неё больших возмущений, которые вносятся в систему путём реализации специальных сценариев, имитирующих некоторые аварии и нештатные ситуации, наиболее часто встречающиеся в практике эксплуатации ЭЭС. Набор таких сценарных аварий немногочислен и обычно ограничивается короткими замыканиями между фазами и на землю, отключением линий электропередач, крупных генерирующих узлов или потребителей [13]. Подобные исследования, очевидно, не гарантируют, что система будет демонстрировать требуемое качество переходных процессов при возмущениях, отличных от сценарных. Эта особенность численных методов исследования ЭЭС может расцениваться как недостаток, особенно в сочетании с ограниченностью наборов реализуемых сценариев.
Вместе с тем существуют и в последнее время получают активное развитие различные аналитические методы исследования, позволяющие учесть некоторые нелинейности энергосистем. Так, в работе [18] авторы применили спектральный анализ Гильберта для оценки нелинейных модальных взаимодействий и анализа межрайонных колебаний, в [24] применён метод мод Купмана для задачи идентификации когерентных генераторов, в работе [5] используется метод разложения по динамическим модам для идентификации опасных переходных процессов в крупных ЭЭС. Перечисленные методы объединяет то, что все они используют наборы данных, полученные либо в ходе измерений на реальных ЭЭС, либо путём интегрирования систем нелинейных дифференциально-алгебраических уравнений, описывающих динамику энергосистемы.
В части аналитических методов, основанных на моделях ЭЭС, следует отметить метод нормальных форм Пуанкаре как получивший наибольшее распространение при решении различных теоретических и практических задач. Хотя этот метод известен достаточно давно, примером первого его использования в сфере электроэнергетики, по всей видимости, является работа [27], в которой нелинейная энергосистема, представленная с помощью ряда нелинейных преобразований в нормальной форме, сопоставляется с исходной нелинейной моделью по переходным процессам при больших возмущениях. Кроме задач анализа [27] и прогнозирования [25] межрайонных колебаний, метод нормальных форм показал хорошие результаты и в задачах настройки регуляторов систем возбуждения генераторов [17], оптимального расположения системных стабилизаторов [16], анализа устойчивости систем с возобновляемыми источниками энергии [2], исследования резонансов [6] и др.
Родственным к методу нормальных форм Пуанкаре является метод линеаризации Карлемана, который заключается в расширении вектора переменных состояния исходной нелинейной модели за счёт включения в этот вектор произведений его первоначальных элементов в различных комбинациях. Достаточно подробное введение в этот метод приведено в [22] и в [26]. В последней работе, кроме прочего, показана эквивалентность
линеаризации Карлемана методу нормальных форм Пуанкаре для случая автономной динамической системы.
Основное отличие метода линеаризации Карлемана от всех перечисленных ранее заключается в том, что конечный результат применения этого метода имеет вид классической динамической системы в пространстве состояний, что несомненно является преимуществом, особенно при решении задач управления энергосистемами. Кроме этого, данный метод позволяет учитывать не только нелинейности между переменными состояния, но и нелинейности между переменными состояния и входами системы (управлениями), то есть применим к классу билинейных динамических систем, которые активно используются в электроэнергетике [8-10].
Несмотря на хорошо описанную теорию и широкое применение линеаризации Карлемана в различных задачах [11, 12, 20, 21], готовые аппроксимационные модели ЭЭС либо программное обеспечение (ПО) для их получения в свободном доступе найти проблематично. Немногочисленные работы [3, 4], в которых использовались аппроксимированные с применением указанного метода модели энергосистем, не содержат подробного описания самих аппроксимаций. Вместе с тем, наличие последних (или соответствующего ПО) могло бы быть достаточно удобным в теоретических исследованиях, особенно на фоне сложившейся специфической практики в сфере электроэнергетики. Она состоит в том, что в научных исследованиях вместо реальных ЭЭС традиционно используются специально разработанные тестовые энергосистемы. На сегодняшний день уже существуют готовые серии таких систем, наиболее распространённой из которых является библиотека тестовых ЭЭС IEEE. Для компьютерной обработки этих систем применяется специализированное ПО, которое обычно позволяет выполнять построение нелинейных моделей энергосистем и их линейных приближений, а также построение динамического отклика нелинейной модели и проведение различных аналитических исследований линейной модели (например, модальный анализ). На работу с квадратичными аппроксимациями данное ПО не рассчитано. В этой связи создание библиотеки квадратичных ап-
проксимаций общеиспользуемых тестовых ЭЭС, а также ПО для построения квадратичных аппроксимаций в условиях отсутствия аналогов таких библиотек и ПО можно считать актуальной задачей.
В связи с этим ряд предыдущих работ автора настоящей статьи был посвящён построению квадратичных и кубических аппроксимаций простейшей одномашинной модели ЭЭС [7, 14], а в [15] впервые была рассмотрена многомашинная двухрайон-ная система, для нелинейных уравнений динамики которой получено разложение в ряд Тейлора до квадратичных членов.
Настоящая же работа является продолжением указанной серии публикаций. В этой статье для нелинейной модели двух-районной четырёхгенераторной ЭЭС представлено описание техники перехода от многочлена Тейлора, полученного в работе [15], к пространству состояний; дано описание программного обеспечения, предназначенного для автоматизированного получения квадратичных аппроксимаций рассматриваемой ЭЭС при различных начальных параметрах генераторов и сети; представлены результаты верификации полученной квадратичной аппроксимации.
Дальнейшее изложение статьи организовано следующим образом: во втором разделе представлено краткое изложение техники линеаризации Карлемана; третий раздел содержит подробное описание нелинейной математической модели двух-районной ЭЭС; в четвёртом разделе изложена техника построения квадратичной аппроксимации исследуемой двухрайонной ЭЭС, а также краткое описание программного обеспечения, реализующего данную технику; пятый раздел посвящён верификации и анализу полученной квадратичной аппроксимации; шестой раздел содержит основные выводы проведённого исследования.
2. Линеаризация Карлемана
Определение 1. Кронекеровской степенью порядка к вектора г = [г 22 ... г„]т, г 6 К, будем называть вектор 2Ш, определяемый следующим выражением:
к-\ раз
где ® - произведение Кронекера.
Рассмотрим динамическую систему следующего вида:
(1) ^=/(*)+ т
где х е Кп - вектор переменных состояния; и е К™ - вектор входов (управлений) системы; / () = [/1 /2 ... /]т, g () = [gl g2 ... gn]т -произвольные, вообще говоря, нелинейные вектор-функции векторных аргументов, аналитические на всей области определения.
Пусть существует пара точек {хо, ио} такая, что / (хо) + g(uо) = 0, т.е. система (1) в рассматриваемой паре точек находится в равновесии. Тогда динамику отклонения точки х от равновесной точки хо можно записать через отклонения функций /(х) и g(u) от соответствующих равновесных значений / (хо) и g(uо):
(2) ф-Х0) = Лх = / (х)_ / х )+ § (м )_ § (Мо ) = д/ ^ )+ ^ (Мо )
т
Поведение отклонений Л/(хо) и Лg(uо) в некоторых окрестностях точек хо, ио может быть аппроксимировано многочленами Тейлора порядка N и М соответственно:
(3) Гы (Дх) = ±А4Х01 Дх",
к=1 к!
(4) Ом (Ди^^Ди
к=1 к!
где Ак(хо), Вк(ио) - матрицы частных производных порядка к вектор-функций /( ) и g ( ) в точках хо и ио соответственно (для краткости в дальнейшем аргументы будем опускать, пользуясь обозначениями вида Ак, Вк). В частном случае при к = 1 эти матрицы представляют собой матрицы Якоби соответствующих функций, а при к = 2 - формируются из матриц Гессе скалярных функций/, gi следующим образом:
®к
(5)
А =
±в\ ®(е'пИ ( / ))
г=1
±еП ®((/2))
/=1
52 =
®(етн (
1=1
®(е;я(Я2))
/=1
Ь'т ®(е'тИ (Яп ))
где ек - вектор-строка из к элементов, состоящий из нулей с единицей на позиции I; Н£), - матрицы Гессе скалярных функций £ и gi соответственно.
Далее будем предполагать, что М = 1 (т.е. функция входа линеаризуется) и N 6 {1, 2, ...}. Тогда аппроксимация системы (2) до степени N по х и до степени 1 по и может быть представлена в виде
N А
(6) Аг = ^ (Ах-) + Ц (Ам) = £ — Л*8* + ВгМ.
к=1 к\
Согласно определению 1, использование кронекеровских степеней в выражении (6) будет обеспечивать возникновение в правой части последнего новых переменных состояния, не содержащихся в исходном векторе Ах. Уравнения динамики для этих новых переменных состояния могут быть получены с использованием правила дифференцирования произведения функций:
(7)
й (Ах®* )
сИ
= Дхи
>Ах + Дх"
) Ах.
Тогда выражение (6) в новом пространстве состояний может быть представлено следующим билинейным уравнением [22]:
(8) М = ААх + ЫАхАи + ВАи,
где Ат =
Аг' Аг'
®1
®2
Аг®
А =
А А • Аы
о ^2,1 ' А
о О • • А
о о •
' О О ■■■ О ' В'
N = В2\ О ■■■ О , В = О
О '•• вм О О
О - нулевая матрица.
Матрицы с двойной индексацией, образующие отдельные
блоки матриц А и N , вычисляются по формулам
з
(9) ^=2]/®--® 4
к=1
(10) ВЛ=±1<8>-
позиция к
' позиция А:
где 7 - единичная матрица.
3. Описание нелинейной модели двухрайонной ЭЭС
Электроэнергетическая сеть обычно описывается системой дифференциально-алгебраических уравнений, общий вид которых для к-го динамического устройства может быть представлен следующим образом:
(11)
где хк, Ук и 1к - соответственно вектор переменных состояния, вектор напряжений на клеммах и вектор впрыскиваемых в систему токов к-го динамического устройства; / ( ), g ( ) - нелинейные функции.
Аналитическое представление статических устройств ЭЭС, очевидно, будет состоять только из алгебраических уравнений, связывающих впрыскиваемые токи с напряжениями на клеммах.
В данной работе рассматривается тестовая двухрайонная ЭЭС, общая схема которой представлена на рис. 1.
Рис. 1. Общая схема исследуемой двухрайонной тестовой ЭЭС [13]
Эта система широко известна в сфере электроэнергетики как «схема Кундура» и активно используется в теоретических и прикладных исследованиях, связанных с идентификацией межрайонных колебаний и их надлежащем гашении. Подробное описание линейной части этой системы для различной степени детализации генераторов приведено, например, в [13].
Динамические устройства в рассматриваемой ЭЭС представлены в виде четырёх синхронных генераторов, подключенных в узлах №1-4, а статические устройства - в виде двух нагрузок, в узлах №7 и 9. Генераторы вместе с нагрузками образуют две локальные области протяжённостью несколько десятков километров (первая область образована генераторами в узлах №1-2 и нагрузкой в узле №7, а вторая область - оставшимися двумя генераторами и нагрузкой), которые соединены между собой парой линий (между узлами №7 и №9) протяжённостью 220 км.
3.1. МОДЕЛЬ ГЕНЕРАТОРА
Количество и структура дифференциальных уравнений, описывающих динамику генератора, определяются исходя из целей моделирования и требуемой точности модели. Самые простые, применяемые на практике модели имеют второй порядок и описывают только медленные процессы, связанные с электромеханическими колебаниями ротора генератора. Более сложные модели позволяют учитывать относительно быстрые
процессы изменения потокосцеплений демпферных обмоток и обмотки возбуждения, а также влияние на эти процессы автоматического регулятора и системного стабилизатора. Количество переменных состояния в таких моделях может быть от трёх до девяти.
В данной работе используется относительно простая, но в то же время достаточно подробная модель генератора шестого порядка, которая учитывает динамику потокосцепления обмотки возбуждения, а также динамику потокосцеплений трёх демпферных обмоток. В качестве переменных состояния выбраны: отклонение частоты вращения ротора генератора Аюг, угол поворота ротора генератора 3, потокосцепление обмотки возбуждения уа, а также потокосцепления демпферных обмоток Ум у^ и ^2q. Входами системы являются механический крутящий момент Tm на приводном валу генератора и напряжение обмотки возбуждения Efd. В указанном пространстве состояний уравнения динамики (первое уравнение в (11)) для рассматриваемой модели имеют вид (12)
(12)
а д®г 1 (
йг — 1 2 Иу
йд аЛаг ,
йг
_ аоКА
йг Ей
й¥и ао Ки
йг Ей
ао Кц
йг
й¥2Ч _ а0К2Ч
йг ц 2Ч
Ааг - Тт ),
-¥аа ) + -Т— Е
т
айы
(у ),
-(у -ущ);
где ^ - коэффициент демпфирования; H - постоянная инерции ротора; юо - базовая частота сети в рад/с; Rfd - сопротивление обмотки возбуждения; Lfd - индуктивность обмотки возбужде-
ния; Ladu - ненасыщенное значение взаимной индуктивности между обмоткой статора и обмоткой ротора на d-оси; Ru, Rq R2q - сопротивления демпферных обмоток; L1d, L\q, L2q - индуктивности демпферных обмоток; tyad, tyaq - взаимные потокосцеп-ления между обмотками ротора и статора, численные значения которых вычисляются по выражениям (13)-(14):
Í „г \
Wfd , Wn
(13) Wad = L'd.
(14) w = L"
v / r aq aqs
-ь
-i„
Lfd L1a
Aq
W2 q
L
\
2 q /
где id и iq - токи обмоток статора, вычисляемые по выражениям:
L„
(15)
K
R,
L"
ads
f ^
Wfd , Wid
V Lfd L1d
- Vj sin (ó)- vR cos(ó)
K
f
I"
aqs
L1q
W2q
R,
(16)
i =—-q K
v iq (
L
+ vR sin (ó) - v cos (ó)
'2q
L"
ads
WfL + Wk.
Lid
V Lfd
\ \ - V sin (ó) - vñ cos (ó)
+-
L
K
f
L"
aqs
К
W2 q
L
\
\
+ VR sin (ó)- V cos (ó)
2 q /
ин-
где ^ = + Ц; ья, = Цц, +1;; ^ = ЬйЬя, + , -
дуктивность утечки статора; Д - сопротивление статора; Уд, V/ - соответственно действительная и мнимая часть сетевого напряжения Ук на клеммах к-го генератора; и Ц - насыщенные значения взаимных индуктивностей между обмотками статора и обмотками ротора на й- и ^-осях соответственно, которые вычисляются по следующим выражениям:
(17) Ц* - 1 " 1
1
1 1
■ +-+ -
L" =■
aqs
1
Lads L
fd
L
1 1
■ +-+ -
1d
L L L
aqs 1q
"2q
id =
Уравнения токов (второе уравнение в (11)) для рассматриваемой модели генератора имеют вид
(18) |1к = 8'П (^) + С0§ ' 111 = 1Ч зт С0$ (£);
где 1е и II - соответственно действительная и мнимая часть тока 1к, впрыскиваемого к-м генератором в сеть.
3.2. МОДЕЛЬ НАГРУЗКИ
Нагрузки в узлах №7 и №9 в общем случае могут описываться различными линейными и нелинейными функциями. В данной работе предполагается, что токи, потребляемые к-й нагрузкой из сети, нелинейно зависят от напряжений следующим образом [13, с. 796-797]:
( , , \т ( , , У
(19)
Т = Р
±К(к) Ь0 (к)
Т = Р
±1 (к) Ь 0 (к)
V
V
0 (к)
+ О
2 + О10 (к)
К
V
V
V
0 (к)
V
0 (к)
Т(к) _п
2 ОЩк)
V
К
V
0(к)
^ (к) \К\2
УЧ к) .
I |2 ;
V
где 1к(к) и 1дк) - соответственно действительная и мнимая часть тока 1к, потребляемого к-й нагрузкой из сети; уад и Увд - соответственно действительная и мнимая часть напряжения Ук в узле подключения к-й нагрузки в сеть; Рщк) и Qu)(k) - соответственно установившиеся значения активной и реактивной мощности к-й нагрузки; Уо(к) - установившееся значение напряжения в узле подключения к-й нагрузки в сеть. Параметры т и п предполагаются равными соответственно 1 и 2 [13, с.814].
4. Построение квадратичной аппроксимации
4.1. РАЗЛОЖЕНИЕ В РЯД ТЕЙЛОРА
Нелинейная модель рассматриваемой четырёхмашинной ЭЭС состоит из четырёх блоков уравнений вида (12), четырёх
V
блоков уравнений вида (18), и двух блоков уравнений (19). Разложение этой системы в ряд Тейлора в равновесной точке до членов 2-го порядка для уравнений генераторов и до членов 1 -го порядка для уравнений нагрузок можно записать в векторно-матричном виде следующим образом:
' Ах = Ав Ах + ВвАГ + 42)£Ах®2 + В^МАУ®2 + (20) ] + Н0Р(АУ®Ах) + В1Аи;
А1 = С0 Ах + Ц А V + С£2) £Ах®2 + ^Р (А V ® Ах).
В приведённой системе использованы следующие обозначения:
- Дх[24х1] - вектор отклонений переменных состояния ЭЭС от равновесного значения. Формируется из векторов переменных состояния отдельных генераторов. Общая структура вектора и его отдельных составляющих:
Ах =
хо(0
хб(1) хо( 2) хо(3) хо(4)
_А®(0 А5(г) А^ (г) (г) А^Цг) А^2
- Д/[12х1] - вектор отклонений узловых токов ЭЭС от равновесного значения. Каждому к-му активному узлу системы соответствует два компонента вектора Д/: действительная Тед и мнимая /дк) составляющая. Неактивные узлы 5, 6, 8, 10 и 11, которые не впрыскивают и не потребляют ток, исключаются из уравнений сети с помощью соответствующих преобразований, т.е. исходная система на рис.1 замещается эквивалентной ше-стиузловой системой. Общая структура вектора имеет вид
нагрузочные узлы
А/ = АА1(2) А1(3) А1(4) А1(1) А1(9)
генераторные узлы
- ДК[12х1] - вектор отклонений узловых напряжений, структура которого аналогична структуре вектора Д/;
- Дм[8х1] - вектор отклонений входных сигналов (управлений) ЭЭС от равновесных значений. Формируется из векторов
управлений отдельных генераторов. Общая структура вектора и его отдельных составляющих:
Ап =
и
П„,Л =
0(1) АТ
0(2) и0(3) и0( 4)
АЕ
УЫ (г
О)
- Ло [24x24], 5д[24х12], Л%[24 х144], В(2)о[24х24], ^[24x48], В\ [24x24] -
блочные матрицы коэффициентов Тейлора при разложении в ряд уравнений динамики генераторов вида (12). Общая структура данных матриц (за исключением матрицы В [24x24]), структура отдельных блоков, а также выражения для вычисления их ненулевых элементов представлены в работе [15]. Матрица В1 [24x8] имеет следующую структуру:
А = ^
в(1)
в(2) в
=
1
2И
0 0
Куы (г
ь
в
0 0 0
000
- с
о [12x24], Оо[12х12].
"'аЫи(г)
С(2)о[12 xl44], 5,D[l2x48] — блочные матрицы коэффициентов Тейлора при разложении в ряд уравнений токов вида (18)—(19). Общая структура данных матриц, структура отдельных блоков, а также выражения для вычисления их ненулевых элементов представлены в работе [15];
- Е [144x576], .^[24x144], ^[48x288] - блочные матрицы специального вида, состоящие из нулей и единиц и имеющие следующую структуру:
Е,.
Е =
(
ЕЕ
Е
м =
р =
м„ч м^ м^ м,
[(1)
(4).
'[16x96]
а
[16x48]
а
рр
р1) р2)
[8х96]
рр
р3) р 4)_
а
[8x48]
а
[48x96]
Т
Т
Т
г )[36х144]
ек ,/[6х 6]
ы,
Р
)[4х 24]
т
'к ,1[2х 2]
)[12х 48]
ек ,1 = Рк ,1 =
к=1...6
1=1.. .24 к=1.. .2
Ы...12
Аг=1...2
/=1..
Рк ,1[6х 6]
\хб] при/ = 4(^-1)+7, 0[6х6] при остальных I;
\<2] пРи 1 = 6 (к " 1) + ^
0[2х2] при остальных I;
где 1 - единичная матрица; О - нулевая матрица; 1 = 1.. .4 - номер генератора; к - номер строки; I - номер столбца.
4.2. ПЕРЕХОД К ПРОСТРАНСТВУ СОСТОЯНИЙ
Для представления системы (20) в пространстве состояний необходимо исключить из её первого уравнения вектор узловых напряжений и выполнить группировку по векторам Дх и Ах02.
Учитывая, что Д/ = УДУ (где У - матрица проводимости сети), из второго уравнения системы (20) получим:
(21) АV = (У -Ц0 -80Р(1[12х12] 0 Ах))-1 (СцАх + сЦ2ЕАх02).
Поскольку выражение (21) нелинейно по Дх и по Ах02, то при его использовании в первом уравнении системы (20) группировка по указанным векторам невозможна. Поэтому, применяя к выражению (21) (а также к выражениям для АV02 и для ДУ 0 Дх, которые тоже используются в системе (20)) аппроксимацию многочленом Тейлора второго порядка, получим следующую систему в пространстве состояний:
(22) Ах = ДАх + ^Ах:®2 + ВгАи ,
Д =Ап+ ВпВ1т + В^МВ21т + НвРНПп,
4 = БвВщг + АЦ} Е + ВЦ ыв2 ^ + И0РИщг
?(2 ь
где Вйи, В2йи, Инп - матрицы коэффициентов Тейлора линейных членов при аппроксимации многочленом Тейлора второго порядка выражения (21), выражения для АУ02 и выражения для АУ 0 Ах соответственно; В^г, В^дг, ИЩг - матрицы коэффициентов Тейлора квадратичных членов при аппроксимации многочленом Тейлора второго порядка выражения (21), выражения для АУ02 и выражения для АУ 0 Ах соответственно.
Используя процедуру линеаризации Карлемана, описанную в разделе 2 настоящей работы, запишем уравнение (22) в билинейном виде:
(23) Ах = ААх +
к=1
где Ах =
-В Ап.
Ах Г А А2' Г О О' Г в 1
, А = 1 , NN. = В ( к ) _ В21 , в = 1
_Ах02_ ? О А21 _ ? ^ ' к О О
-й стол-
Ах = А 01+10 А, в£) = в(к)01+10 в(к), В1(к) - к-й
бец матрицы В1.
Выражение (23) является искомой квадратичной аппроксимацией исследуемой нелинейной модели ЭЭС в пространстве состояний. Выражение (22) без второго слагаемого является линейной аппроксимацией.
4.3. ОПИСАНИЕ ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ ДЛЯ ПОСТРОЕНИЯ КВАДРАТИЧНОЙ АППРОКСИМАЦИИ
В рамках данного исследования разработано программное обеспечение, позволяющее автоматизировано вычислять матрицы Л, Ык и В билинейного уравнения (23), а также матриц Л1 и В1 линейной части уравнения (22) для различных параметров исследуемой двухрайонной ЭЭС. Все необходимые для этого функции и скрипты реализованы в среде МайаЬ и располагаются в открытом доступе на ресурсе ОйНиЬ [29]. Пакет программ имеет следующий состав и функционал:
1) арргох.ш - базовый скрипт, при выполнении которого запускается процесс вычисления матриц Л, Ык и В билинейного
уравнения (23), а также матриц А1 и В1 линейной части уравнения (22);
2) id4mac.m - функция, вызываемая скриптом арргох.ш и содержащая параметры генераторов, линий и нагрузок для рассматриваемой двухрайонной ЭЭС. Данная функция поддерживает ручное изменение пользователем значений параметров ЭЭС. Установка желаемых параметров системы должна быть выполнена до запуска скрипта арргох.т;
3) fun_diagmat.m - функция, вызываемая скриптом арргох.т и осуществляющая вычисление матрицы В1, а также матриц с нижним индексом В из системы (20);
4) tгp_simu.m - скрипт построения графиков переходных процессов линейной (уравнение вида (22) без второго слагаемого) и квадратичной (уравнение (23)) аппроксимации ЭЭС. Возмущающее воздействие реализуется либо через отклонение переменных состояния от равновесных значений, либо через отклонение компонент вектора Дм.
4.4. УСТРАНЕНИЕ ЛИНЕЙНОЙ ЗАВИСИМОСТИ В МАТРИЦАХ ДИНАМИКИ И ИЗБЫТО ЧНОСТИ В ВЕКТОРЕ ПЕРЕМЕННЫХ СОСТОЯНИЯ
Вследствие использования произведения Кронекера при формировании вектора Ах , парные произведения разноимённых переменных состояния будут дублироваться, т.е. входить в указанный вектор дважды. С целью устранения такой избыточности и уменьшения размерности системы (23) дублирующиеся переменные состояния вектора Ах, а также соответствующие им строки и столбцы матриц А, Ык и В подлежат удалению. Данная процедура реализована в виде функции Шп_^исе, содержащейся в скрипте арргох^.
Кроме этого необходимо отметить, что в силу различных особенностей моделирования многомашинных энергосистем спектр матрицы динамики линеаризованной модели всегда содержит одно или два близких к нулю собственных значения, действительные части которых могут быть в том числе и положительными. Одно из этих собственных значений связано с от-
сутствием однозначности в абсолютной величине угла поворота ротора. Эта неоднозначность проявляется в том, что изменение угла ротора всех генераторов модели на одинаковую величину не оказывает влияния на устойчивость модели, что, однако, не соответствует положению вещей в реальных ЭЭС. Второе близкое к нулю собственное значение возникает в тех моделях, в которых крутящие моменты генераторов не зависят от отклонений угловых скоростей вращения роторов (те случаи, когда в модели не учтены регуляторы скорости и коэффициенты демпфирования Кэ для всех генераторов равны нулю).
Спектр матрицы динамики квадратичной аппроксимации естественным образом «унаследует» эти два собственных значения, а также их всевозможные попарные суммы. Эти собственные значения могут оказывать негативное влияние на дальнейший анализ системы, а в некоторых случаях их наличие в спектре является крайне нежелательным.
Для устранения собственного значения, связанного с неоднозначностью абсолютной величины угла поворота ротора, необходимо для одного из генераторов, называемого референс-ным (индекс К), положить [13]
Л '
(т.е. предположить, что для референсного генератора отклонение угла поворота ротора от равновесного значения не меняется), а для остальных генераторов принять
^^ = 0о(Д0г(О-Д0г(д)), 7 = 1,..., и;
Устранение второго близкого к нулю собственного значения выполняется аналогичным образом с тем лишь отличием, что указанные выше операции применяются в отношении уравнений угловых ускорений роторов генераторов (первого уравнения системы (12)).
В разработанном пакете программ реализация описанной выше процедуры представлена в виде функции &п_НМер, содержащейся в скрипте арргох.ш.
5. Верификация
В первой части данного раздела представлены результаты сравнения переходных процессов линейной и квадратичной аппроксимаций исследуемой ЭЭС при различных видах и величинах возмущений, а во второй части приведён анализ спектров матриц динамики указанных аппроксимаций.
Исходные параметры генераторов, нагрузок и линий исследуемой ЭЭС приняты в соответствии с [13, с. 813-814] и представлены в таблице 1. Численные значения переменных состояния и иных величин, характеризующих рабочую точку аппроксимации, представлены в таблице 2.
Таблица 1. Исходные параметры ЭЭС
Параметры генераторов в системе относительных единиц 900 МВА, 20 кВ
Номер генератора по рис. 1 Параметр
X, X* X*
1, 2, 3, 4 1,8 1,7 0,2 0,3 0,55 0,25
X' К Т 0 Т Тч 0 ГТ1ГГ Т* 0 гггГГ Тч 0
1, 2, 3, 4 0,25 0,0025 8,0 с 0,4 с 0,03 с 0,05
А ■^Ъаг и ка Н
1, 2 0,015 9,6 0,9 0 6,5
3, 4 6,175
Параметры генераторных узлов
Параметр Номер узла по рис.1
1 2 3 4
Активная мощность 700 МВт 700 МВт 719 МВт 700 МВт
Реактивная мощность 185 МВАр 235 МВАр 176 МВАр 202 МВАр
Напряжение на клеммах 1,03 Z 20,2° 1,01 Z 10,5° 1,03 Z -6,8° 1,01 Z -17,0°
Таблица 1 (продолжение)
Параметры нагрузочных узлов
Параметр Номер узла по рис.1
7 8
Активная мощность, МВт 967 1767
Реакт. индуктивная мощность, МВАр 100 100
Реакт. емкостная мощность, МВАр 200 350
Параметры линий в системе относительных единиц
100 МВА, 230 кВ
Параметр Значение
Удельное активное сопротивление 0,0001 о.е./км
Удельное реактивное сопротивление 0,001 о.е./км
Удельная проводимость параллельного 0,00175 о.е./км
емкостного сопротивления
Параметры трансформаторов в системе относительных единиц, 900 МВА, 20/230 кВ
Реактивное сопротивление трансформатора 0,15
Таблица 2. Значения величин, определяющих рабочую точку
Значения переменных состояния
Переменная состояния Номер генератора
1 2 3 4
Люг, о.е. 0 0 0 0
3, рад 1,0794 0,8990 0,6300 0,4463
№, о.е. 0,1340 0,1382 0,1346 0,1346
уы, о.е. 1,7843 1,7814 1,7619 1,7398
Уи, о.е. -0,6165 -0,5969 -0,6300 -0,6122
У2«, о.е. -0,6165 -0,5969 -0,6300 -0,6122
Значения входных величин, о.е.
Тт 0,7778 0,7778 0,7985 0,7778
2,0097 2,0736 2,0195 2,0185
Отметим, что линейная аппроксимация модели ЭЭС, описанной в третьем разделе, имеет размерность 24, а квадратичная аппроксимация после устранения избыточности в векторе переменных состояния имеет размерность 324.
5.1. АНАЛИЗ ПЕРЕХОДНЫХ ПРОЦЕССОВ ЛИНЕЙНОЙ И КВАДРАТИЧНОЙ АППРОКСИМАЦИИ
На рис. 2А представлены переходные процессы для угла поворота ротора генератора №4 линейной (сплошная линия) и квадратичной (штриховая линия) аппроксимаций при отклонении угла поворота ротора генератора №2 на величину Д&32 = 0,01 рад. Как следует из рисунка, при такой величине возмущения оба графика совпадают.
При увеличении возмущения в 100 раз, т.е. при Дёо2 = 1 рад, график переходного процесса линейной аппроксимации изменяется только по масштабу и пропорционально величине возмущения (что и отражает свойство линейности), в то время как график переходного процесса квадратичной аппроксимации изменяется и по масштабу, и по форме - см. рис. 2Б.
Рис. 2. Переходные процессы угла поворота ротора генератора №4 при отклонении угла поворота ротора генератора №2 А) на величину 0,01 рад; Б) на величину 1 рад
Аналогичное поведение исследуемых аппроксимаций наблюдается и при рассмотрении переходных процессов других переменных состояния. Например, на рис. 3 представлены переходные процессы потокосцепления демпферной обмотки генератора №4 (уы).
На рис. 4 показаны динамические отклики рассматриваемых аппроксимаций по углу поворота ротора генератора №4 при возмущении системы через один из её входов, а именно, через отклонение величины Д£де2) напряжения обмотки возбуждения генератора №2.
■ Ю4 (А)
1 -1-1-1-1-1-1-
О 2 4 6 Й 10 12 14
Время, с
-0.04 V__I_=_I_I_I_I_I_:
0 2 4 6 В 10 12 14
Время, с
Рис. 3. Переходные процессы потокосцепления демпферной обмотки генератора №4 при отклонении угла поворота ротора генератора №2 А) на величину 0,01 рад; Б) на величину 1 рад
При этом реализовывалось два варианта возмущающего воздействия: в первом, менее «нагруженном» варианте величина отклонения составляла Д£де2) = 0,5 о.е. на интервале времени ^ = 1...1,1 с - см. рис. 4А; во втором, более «нагруженном»
варианте величина отклонения составляла Д£де2) = 1,5 о.е. на интервале времени t = 1.. .2 c - см. рис. 4Б.
Рис. 4. Переходные процессы угла поворота ротора генератора №4 при отклонении напряжения обмотки возбуждения генератора №2 А) на величину 0,5 о.е. на интервале времени 1.. .1,1 с; Б) на величину 1,5 о.е. на интервале времени 1.. .2 с
По форме графиков можно заметить, что тенденция усиления различий между переходными процессами линейной и квадратичной аппроксимаций по мере возрастания величины возмущающего воздействия сохраняется и при данном способе внесения возмущения.
5.2. АНАЛИЗ СПЕКТРА ЛИНЕЙНОЙ И КВАДРАТИЧНОЙ АППРОКСИМАЦИИ
Линейная аппроксимация исследуемой ЭЭС имеет порядок 24. Все собственные значения (моды) её матрицы динамики представлены в таблице 3, где также указаны их частоты и коэффициенты демпфирования, вычисляемые по следующему выражению:
ЯеЯ
к
damp
у!( ЯеД)2 +(1т А)2 где X - собственное значение.
Таблица 3. Собственные значения матрицы динамики линейной аппроксимации
Номер моды Значение Частота, Гц Коэф. демпфирования
1, 2 0 + 0i 0 -
3 -0,0981 + 0,0000i 0 1
4,5 -0,1112 ± 3,4417i 0,5478 0,0323
6 -0,1131 + 0,0000i 0 1
7 -0,2609 + 0,0000i 0 1
8 -0,2733 + 0,0000i 0 1
9,10 -0,4919 ± 6,8336i 1,0876 0,0718
11,12 -0,5059 ± 7,0238i 1,1179 0,0718
13 -3,3974 + 0,0000i 0 1
14 -4,1233 + 0,0000i 0 1
15 -5,2838 + 0,0000i 0 1
16 -5,2989 + 0,0000i 0 1
17 -30,9588 + 0,0000i 0 1
18 -32,4025 + 0,0000i 0 1
19 -34,0501 + 0,0000i 0 1
20 -35,5215 + 0,0000i 0 1
21,22 -37,8645 ± 0,1486i 0,0236 ~ 1
23,24 -37,9892 ± 0,0613i 0,0098 ~ 1
Как следует из таблицы, спектр линейной аппроксимации содержит пять пар колебательных мод, при этом только три
из них являются слабо демпфированными (т.е. с коэффициентом демпфирования <0,1). Это пары мод №4-5, №9-10 и №11-12. Первая из указанных пар является межрайонной модой, поскольку её частота лежит в промежутке от 0,1 Гц до 1 Гц, а две другие пары мод (№9-10 и №11-12) являются межмашинными. Оставшиеся две пары колебательных мод №21-22 и №23-24 имеют частоты менее 0,1 Гц, и поэтому не относятся ни к межмашинным ни к межрайонным колебаниям.
Квадратичная аппроксимация исследуемой ЭЭС после устранения избыточности в векторе переменных состояния имеет размерность 324. Расположение собственных значений матрицы динамики этой аппроксимации на комплексной плоскости в разных масштабах представлено на рис. 5-6, где также отмечены моды матрицы динамики линейной аппроксимации.
Рис. 5. Расположение собственных значений матриц динамики линейной и квадратичной аппроксимации на комплексной плоскости
Рис. 6. Расположение собственных значений матриц динамики линейной и квадратичной аппроксимации с действительными частями от 0 до -6
Спектр квадратичной аппроксимации содержит 100 пар колебательных мод, при этом 95 из них являются новыми по отношению к спектру линейной аппроксимации. Из этих 95 пар новых мод слабо демпфированными являются 12 пар. Их характеристики представлены в таблице 4, а расположение на комплексной плоскости отражено на рис. 7. В последнем столбце таблицы 4 также указана связь с модами линейной аппроксимации из таблицы 3.
Как следует из приведённой таблицы, учёт квадратичных членов приводит к появлению в спектре двух новых межрайонных слабо демпфированных колебаний (пары мод 158, 159 и 160, 161 в таблице 4), и десяти межмашинных слабо демпфированных колебаний (остальные моды таблицы 4). Новые межрайонные колебания имеют частоты около 0,5 Гц с коэффициентами демпфирования около 0,06 и получаются из спектра линейной
аппроксимации как сумма межрайонной моды №4, 5 и апериодических мод №3 и №6 таблицы 3.
Таблица 4. Слабо демпфированные собственные значения матрицы динамики квадратичной аппроксимации_
№ мод Значение Частота, Гц ^Ратр Номера мод по таблице 3
109, 110 -0,2225 ± 6,88331 1,0955 0,0323 4+4, 5+5
45, 46 -0,6032 ± 10,27531 1,6354 0,0586 4+9, 5+10
43, 44 -0,6171 ± 10,46551 1,6656 0,0589 4+11, 5+12
158, 159 -0,2093 ± 3,44171 0,5478 0,0607 4+3, 5+3
160, 161 -0,2243 ± 3,44171 0,5478 0,0650 4+6, 5+6
41, 42 -0,9838 ± 13,66731 2,1752 0,0718 9+9, 10+10
39, 40 -0,9978 ± 13,85751 2,2055 0,0718 9+11, 10+12
37, 38 -1,0118 ± 14,04761 2,2358 0,0718 11+11,12+12
112, 113 -0,6040 ± 7,02381 1,1179 0,0857 11+3, 12+3
128, 129 -0,5900 ± 6,83361 1,0876 0,0860 9+3, 10+3
120, 121 -0,6190 ± 7,02381 1,1179 0,0878 11+6, 12+6
126, 127 -0,6050 ± 6,83361 1,0876 0,0882 9+6,10+6
15
["Ж1 Моды 37. 39. 41 Мода 43 N 1 1
О Линейная аппроксимация X Квадратичная аппроксимация
I ^ Моды 112,1 126, 12В 'Л од а 45 Мода 11 X | % — Мода 20, Мода 109 X 9
Моды 161 Мода 4 3. 15В
О -'-'-'-'-'-
-1 -0.0 -0.6 -0.4 -0.2 о
Действительная часть
Рис. 7. Расположение собственных значений матриц динамики линейной и квадратичной аппроксимации слабо демпфированных мод с положительными мнимыми частями
Частоты новых межмашинных слабо демпфированных колебаний находятся в диапазоне примерно от 1,1 до 2,2 Гц с коэффициентами демпфирования около 0,06...0,09. Эти колебания получаются путём попарного суммирования мод линейной аппроксимации, как это указано на соответствующих строках в пятом столбце таблицы 4.
Кроме этого, в спектре квадратичной аппроксимации присутствуют ещё 4 пары межрайонных и 4 пары межмашинных мод, которые хоть и не являются слабо демпфированными, однако близки к таковым, и также требуют внимания. Их характеристики приведены в таблице 5.
Остальные моды квадратичной аппроксимации имеют коэффициент демпфирования более 0,48, т.е. достаточно хорошо гасятся, и в данной работе не анализировались.
Таблица 5. Собственные значения матрицы динамики квадра-
тичной аппроксимации с ксЪшр = 0,1.0,2
№ мод Значение Частота, Гц ^kdamp Номера мод по Табл. 3
175, 176 -0,3722 ± 3,44171 0,5478 0,1075 4+7, 5+7
122, 123 -0,7668 ± 7,02381 1,1179 0,1085 11+7, 12+7
130, 131 -0,7528 ± 6,83361 1,0876 0,1095 9+7, 10+7
124, 125 -0,7792 ± 7,02381 1,1179 0,1103 11+8, 12+8
181, 182 -0,3846 ± 3,44171 0,5478 0,1110 4+8, 5+8
132, 133 -0,7652 ± 6,83361 1,0876 0,1113 9+8, 10+8
152, 153 -0,6171 ± 3,58211 0,5701 0,1698 5+11, 4+12
173, 174 -0,6032 ± 3,39201 0,5398 0,1751 5+9, 4+10
В заключение отметим, что объединение спектра линейной аппроксимации и всевозможных попарных сумм её отдельных собственных значений в точности даёт спектр матрицы динамики квадратичной аппроксимации, что является одним из признаков корректности построения последней, и, следовательно, выполненного в разделе 4 моделирования.
6. Выводы
Основным результатом данного исследования является разработанная квадратичная аппроксимация нелинейной модели четырёхгенераторной ЭЭС. Наряду с подробным математическим описанием этой аппроксимации, которое содержится в настоящей статье, было также разработано и представлено в открытом доступе специальное программное обеспечение, осуществляющее построение линейного и квадратичного приближений рассмотренной энергосистемы в пространстве состояний.
Численно было показано, что учёт квадратичных членов приводит к существенному пополнению слабо демпфированной части спектра матрицы динамики системы. А именно, при трёх колебательных модах в спектре линейной аппроксимации, спектр квадратичного приближения содержит уже 12 новых слабо демпфированных мод, две из которых являются межрайонными. Такое изменение в составе мод даёт качественно новую информацию о процессах в системе, в первую очередь о нелинейных модальных взаимодействиях, что наглядно демонстрируют результаты сравнения переходных процессов линейной и квадратичной аппроксимаций, которые также представлены в данной работе.
Полученная аппроксимация может применяться в теоретических и практических исследованиях, связанных с нелинейным анализом устойчивости ЭЭС, в первую очередь в исследовании межрайонных колебаний. По сравнению с повсеместно использующимся для этих целей методом нормальных форм Пуанкаре, использованный в данной работе метод линеаризации Карлема-на позволяет учитывать нелинейности в том числе и в управлении, что формирует новизну результатов данной работы. Кроме этого, разработанная квадратичная модель и пакет программ сами по себе являются новыми, поскольку ни в отечественной ни в зарубежной литературе аналогов этим результатам нет.
Несмотря на вполне законченный результат, данная работа требует развития в части уточнения исходной нелинейной модели генераторов, а именно, их дооснащения регуляторами
напряжения и системными стабилизаторами. Это существенно
расширит область применения модели не только на задачи анализа, но и на задачи синтеза регуляторов, что является предметом дальнейших исследований.
Литература
1. ABED E.H., VARAIYA P.P. Nonlinear oscillations in power systems // Electric Power Energy System. - 1984. - Vol. 6, No. 1. - P. 37-43.
2. AMANO H., YOKOYAMA A. Rotor angle stability analysis using normal form method with high penetrations of renewable energy sources-energy index for multi-swing stability // Proc. of the 20th Power Systems Computation Conference (PSCC-2018), Dublin. - 2018. -P. 1-6.
3. ARROYO J., BAROCIO E., BETANCOURT R., MESSINA A.R. A
bilinear analysis technique for detection and quantification of nonlinear modal interaction in power systems // Proc. of the IEEE Power Engineering Society General Meeting, Montreal, Canada. - 2006. - P. 1-10.
4. ARROYO J., BETANCOURT R., MESSINA A.R., BAROCIO E. Development of bilinear power system representations for small signal stability analysis // Electric Power Systems Research. - 2007. -Vol. 77. - P. 1239-1248.
5. BAROCIO E., PAL B.C., THORNHILL N.F., MESSINA A.R. A dynamic mode decomposition framework for global power system oscillation analysis // IEEE Trans. on Power Systems. - 2015. - Vol. 30, No. 6. - P. 2902-2912.
6. BETANCOURT R.J., BAROCIO E., ARROYO J., MESSINA A.R. A real normal form approach to the study of resonant power systems // IEEE Trans. on Power Systems. - 2006. - Vol. 21, No. 1. - P. 431-432.
7. DUSHIN S.V., ABRAMENKOV A.N., KUTYAKOV E.J., IS-KAKOV A.B., SALNIKOV A.M. Developing a weakly nonlinear power system model using the Carleman bilinearization procedure // Proc. of the IEEE 2nd Int. Conf. on Control Systems, Mathematical Modeling, Automation and Energy Efficiency (SUMMA2020). - Lipetsk: IEEE, 2020. - P. 1-5.
8. GIL-GONZÁLEZ W., MONTOYA O D. Passivity-based pi control of a SMES system to support power in electrical grids: A bilinear approach // Journal of Energy Storage. - 2018. - Vol. 18. - P. 459-466.
9. HE F. Nonlinear adaptive control in the design of power system stabilisers: PhD thesis - Adelaide, 1991. - P. 349.
10. HE F., GIBBARD M.J. Design of an adaptive bilinear power system stabilizer // Automatica. - 1997. - Vol. 33, No. 4. - P. 663-668.
11. HERNÁNDEZ-ORTEGA M.A., MESSINA A.R. Recursive linearization of higher-order for power system models // IEEE Trans. on Power Systems. - 2021. - Vol. 36, No. 2. - P. 1206-1216.
12. HERNÁNDEZ-ORTEGA M.A., MESSINA A.R. Recursive linearization of Carleman-based nonlinear power system models // Electric Power Systems Research. - 2022. - Vol. 209. - P. 1-10.
13. KUNDUR P. Power Systems Stability and Control. - New York: McGraw-Hill, Inc., 1994. - 1176 p.
14. KUTYAKOV E.Y., DUSHIN S.V., ABRAMENKOV A.N., IS-KAKOV A.B., SALNIKOV A.M. Quadratic approximation of nonlinear models of the synchronous machine using the bilinear representation // Proc. of the IEEE 13th Int. Conf. «Management of Large-Scale System Development» (MLSD). - Moscow: IEEE, 2020. - P. 1-5.
15. KUTYAKOV E.Y., DUSHIN S.V., ISKAKOV A.B., ABRAMENKOV A.N. Taylor series decomposition of nonlinear model of two-area four-generator power system // Proc. of the IEEE 3nd Int. Conf. on Control Systems, Mathematical Modeling, Automation and Energy Efficiency (SUMMA2021). - Lipetsk: IEEE, 2021. - P. 1-6.
16. LIU S., MESSINA A.R., VITTAL V. A normal form analysis approach to siting power system stabilizers (PSSs) and assessing power system nonlinear behavior // IEEE Trans. on Power Systems. - 2006. -Vol. 21, No. 4. - P. 1755-1762.
17. LOMEI H., ASSILI M., SUTANTO D., MUTTAQI K.M. A new approach to reduce the nonlinear characteristics of a stressed power system by using the normal form technique in the control design of the excitation system // IEEE Trans. on Industry Applications. - 2017. -Vol. 53, No. 1. - P. 492-500.
18. MESSINA A.R., VITTAL V. Nonlinear, non-stationary analysis of inter-area oscillations via Hilbert spectral analysis // IEEE Trans. on Power Systems. - 2006. - Vol. 21, No. 3. - P. 1234-1241.
19. PEREZ-ARRIAGA I.J., VERGHESE G.C., SCHWEPPE F.C. Selective modal analysis with applications to electric power systems, Part I: Heuristic introduction // IEEE Trans. on Power Apparatus and Systems. - 1982. -Vol. 101, No. 9. - P. 3117-3125.
20. RAUH A., MINISINI J., ASCHEMANN H. Carleman linearization for control and for state and disturbance estimation of nonlinear dynamical processes // IFAC Proceedings Volumes. - 2009. - Vol. 42, No. 13. - P. 455-460.
21. ROTELLA F., ZAMBETTAKIS I., DAUPHIN-TANGUY G., BORNE P. Modelling and simulation of non-linear systems based on infinite bilinear realizations // IFAC Proceedings Volumes. - 1987. -Vol. 20, No. 12. - P. 61-66.
22. RUGH W.J. Nonlinear System Theory. The Volterra/Wiener Approach. - Baltimore: The Johns Hopkins University Press, 1981. - 344 p.
23. SANCHEZ-GASCA J.J., VITTAL V., GIBBARD M.J., MESSINA A.R., VOWLES D.J., LIU S., ANNAKKAGE U.D. Inclusion of higher order terms for small-signal (modal) analysis: committee report-task force on assessing the need to include higher order terms for small-signal (modal) analysis // IEEE Trans. on Power Systems. - 2005. - Vol.20, No. 4. - P. 1886-1904.
24. SUSUKI Y., MEZIC' I. Nonlinear Koopman modes and coherency identification of coupled swing dynamics // IEEE Trans. on Power Systems. - 2011. - Vol. 26, No. 4. - P. 1894-1904.
25. THAPAR J., VITTAL V., KLIEMANN W., FOUAD A.A. Application of the normal form of vector fields to predict interarea separation in power systems // IEEE Trans. on Power Systems. - 1997. - Vol. 12, No. 2. - P. 844-850.
26. TSILIGIANNIS C.A., LYBERATOS G. Normal forms, resonance and bifurcation analysis via the Carleman linearization // Journal of Mathematical Analysis and Applications. - 1989. - Vol. 139, No. 1. -P. 123-138.
27. VITTAL V., KLIEMANN W., STARRETT S.K., FOUAD A.A. Analysis of stressed power systems using normal forms // Proc. of the IEEE Int. Symposium on Circuits and Systems, San Diego, CA, USA. -1992. - P. 2553-2556.
28. WINKELMAN J.R., CHOW J.H., BOWLER B.C., AVRAMOVIC' B., KOKOTOVIC' P.V. An analysis of interarea dynamics of multi-machine systems // IEEE Trans. on Power Apparatus and Systems. - 1981. -Vol. PAS-100, No. 2. - P. 754-763.
29. https://github.com/e-kut/Bilinear-4machine-EPS-model (дата обращения: 15.09.2022).
DEVELOPING AND ANALYSIS OF A QUADRATIC AP-PROXIMATION OF A NONLINEAR TWO-AREA POWER SYSTEM MODEL IN THE STATE SPACE
Evgeniy Kutyakov, V.A. Trapeznikov Institute of Control Sciences of RAS, Moscow, research associate (evgeniykutyakov@gmail.com).
Abstract: We present the results of mathematical modeling of a two-area four-generator electric power system (EPS) in the state space, with accounting the second order terms. Software for construction of quadratic approximations of the EPS was developed. The model is based on a 6-th order nonlinear synchronous generator model that integrates the electromechanical equations of the rotor motion, the excitation winding electromotive force (EMMF) dynamics, and the dynamics of the three damping windings. Automatic voltage regulators and system stabilizers are not considered. For construction of the quadratic approximation, the Carleman linearization method is used, which is widely used for representation of nonlinear dynamic systems models in bilinear form. A detailed analysis of the dynamics matrix spectrum of the obtained quadratic approximation is presented. The results of numerical experiments are presented to demonstrate that in comparison with the traditionally used linearized model the developed quadratic approximation gives qualitatively new information about the transient processes of the initial EPS model, arising both at the deviation of some state variables and at the disturbance through system inputs. The resulting model can be used for analytical studies of the nonlinear behavior of the EPS near its operating point, in particular, for the analysis of interarea oscillations and parametric resonances. It should be noted that the method for constructing a quadratic approximation (Carleman linearization) applied in this work can be extended to EPS of larger sizes, and the developed software allows upgrading for taking into account terms of the third and higher orders of smallness.
Keywords: nonlinear analysis of EPS, nonlinear model of EPS, Carleman linearization, inter-area oscillations, nonlinear modal interactions, Taylor series.
УДК 004.942 ББК 31.27
DOI: 10.25728/ubs.2022.100.6
Статья представлена к публикации членом редакционной коллегии П.В. Пакшиным.
Поступила в редакцию 04.10.2022.
Опубликована 30.11.2022.
