ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2022, том 58, № 2, с.280-288

ХРОНИКА

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ

ДИНАМИКИ И УПРАВЛЕНИЯ ПРИ МОСКОВСКОМ

ГОСУДАРСТВЕННОМ УНИВЕРСИТЕТЕ

им. М.В. ЛОМОНОСОВА∗)

Ниже публикуются краткие аннотации докладов, состоявшихся в осеннем семестре 2021 г.

(предыдущее сообщение о работе семинара дано в журнале “Дифференц. уравнения”. 2021.

Т. 57. № 2; за дополнительной информацией обращаться по адресу: nds@cs.msu.su)∗∗).

DOI: 10.31857/S0374064122020145

А.С. Фурсов, Ю.М. Мосолова (Москва, Россия, ВМК МГУ). “Перспективы использо-

вания нейрорегуляторов для стабилизации переключаемых систем” (13.09.2021).

Рассматривается непрерывная скалярная переключаемая линейная система

x = A_σx + b_σu, σ ∈ S_τ,

(1)

где σ : R₊ → I = {1, . . . , m} - кусочно-постоянная функция (переключающий сигнал) с ко-

нечным числом разрывов (переключений) на любом конечном промежутке; S_τ - множество

всех тех переключающих сигналов σ, для которых время между каждыми двумя соседни-

ми переключениями не меньше τ (τ > 0); I - множество индексов, нумерующих режимы

функционирования системы (1); A_σ = A ◦ σ - композиция отображения A: I → {A₁, . . . , A_m}

(A_i ∈ R^n×n) и переключающего сигнала σ, а b_σ = b◦σ - аналогичная композиция для отобра-

жения b: I → {b₁, . . . , b_m} (b_i ∈ Rⁿ); пары матриц (A_i, b_i), i = 1, m, определяют режимы

функционирования системы (1); x ∈ Rⁿ - вектор состояния, u ∈ R¹ - управляющий вход.

Требуется стабилизировать систему (1) в некоторой окрестности D(0) нуля, т.е. постро-

ить обратную связь u = u(x) такую, что для каждого σ ∈ S_τ и любого x₀ ∈ D(0) норма

решения x(t) с начальным вектором x(0) = x₀ системы (1) стремится к нулю при t → ∞. При

этом вектор-функция x(t) определяется как решение линейной дифференциальной системы

с кусочно-постоянными коэффициентами

x = Aσ(t)x + bσ(t)u, x(0) = x₀, σ ∈ S_τ.

Один из возможных подходов к решению поставленной задачи заключается в построении

переключаемого регулятора

u = -k_σx,

(2)

каждый режим которого является стабилизирующей обратной связью u = -k_ix для, соот-

ветственно, i-го режима

x = A_ix + b_iu переключаемой системы (1), т.е. обеспечивает асимп-

тотическую устойчивость замкнутой системы

x = (A_i - b_ik_i)x. Такую обратную связь всегда

можно построить в предположении полной управляемости пары (A_i, b_i) [1, с. 249].

Основная проблема данного подхода состоит в том, что для реализации указанного ре-

гулятора необходимо знать моменты переключения режимов и номера активных режимов в

каждый момент времени (т.е. функцию σ ) - для того чтобы обеспечить синхронность пере-

ключений регулятора и отображения A_σ. Однако, как правило, такая информация о пере-

ключаемой системе в режиме реального времени не доступна.

∗) Семинар основан академиками РАН С.В. Емельяновым и С.К. Коровиным.

∗∗) Составитель хроники А.В. Ильин.

280

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

281

Возможное решение, позволяющее использовать данный регулятор, заключается в постро-

ении “наблюдателя”, который определял бы номер активного режима для текущего момента

времени. В качестве такого “наблюдателя” предлагается использовать нейронную сеть [2, с. 55].

Пару - переключаемый регулятор и нейронная сеть - далее будем называть нейрорегуля-

тором. Перенумеруем (от 1 до m²) все возможные режимы

x = (Ai - bikj )x (i, j = 1, m)

системы (1), замкнутой переключаемым регулятором (2). Основная работа нейронной сети

состоит в следующем: выбирается дискретная последовательность моментов времени εl , l =

= 0, 1, 2 . . . , где ε - достаточно малая положительная константа (ε < τ), и далее по каждой

паре измерений вектора состояния (x(εl), x(ε(l + 1))) нейронная сеть должна определять но-

мер текущего режима замкнутой переключаемой системы в момент t = ε(l + 1). Используя

данную информацию, осуществляется переключение векторов k_i стабилизирующей обратной

связи в переключаемом регуляторе.

Можно выделить несколько основных шагов настройки нейрорегулятора:

1) выбор структуры нейросети;

2) выбор величины ε c учётом времени задержки τ;

3) выбор регуляторов u = -k_ix, обеспечивающих нужное быстродействие с учётом пара-

метра τ;

4) построение обучающей выборки для нейросети с учётом окрестности D(0), в которой

осуществляется стабилизация переключаемой системы;

5) обучение нейросети на парах (x(0), x(ε));

6) проверка работоспособности нейросети по определению текущего активного режима;

7) моделирование переключаемой системы (1), замкнутой нейрорегулятором.

Данный алгоритм прошёл апробацию на примере переключаемой системы вида (1) вто-

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

изложенного выше алгоритма стабилизации системы (1).

Работа выполнена при финансовой поддержке Российского научного фонда (проект 22-21-

00162).

Литература. 1. Андреев Ю.Н. Управление конечномерными линейными объектами. М.,

1976. 2. Хайкин С. Нейронные сети. М., 2006.

Н.А. Магницкий (Россия, Москва, ФИЦ ИУ РАН). “Вывод дифференциальных уравне-

ний Максвелла из уравнений классической механики” (27.09.2021).

Считается, что система линейных дифференциальных уравнений электромагнитного поля

в вакууме - система уравнений Максвелла - “не вытекает из каких-либо более общих теоре-

тических положений, но является обобщённой записью наблюдавшихся на опыте фактов” [1,

c. 41]. Покажем, что это не так и что система уравнений Максвелла является следствием систе-

мы уравнений вакуума (эфира), предложенной автором в [2-4] и описывающей эфирную среду

в трёхмерном евклидовом пространстве двумя нелинейными уравнениями, которые следуют

из уравнений классической механики Ньютона

∂ρ

dρ u

∂ρu

+ div (ρu) = 0,

+ (u · ∇)(ρu) = 0,

(1)

∂t

в них векторы u(t, r) и r(t) являются соответственно скоростью и координатами распростра-

нения локальных возмущений плотности ρ(t, r) среды (эфира) в каждый момент времени t.

Определения векторов напряжённости электрического поля E и магнитной индукции B

даны в работах автора [2-4] и выражаются формулами

B = c∇ × (ρu), E = (u · ∇)(ρu),

(2)

где c - скорость распространения возмущений в эфире (скорость света).

Три уравнения Максвелла из четырёх получаются достаточно просто из системы уравнений

эфира (1). Применим к выражениям (2) оператор дивергенции ∇ · . Получим два уравнения

из системы уравнений Максвелла:

∇·B=0

∇ · E = 4πσ,

(3)

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

282

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

здесь σ имеет смысл плотности заряда, определяемой возмущениями плотности эфира

4πσ = ∇ · (|u|∇(ρ|u|) - u × (∇ × (ρu))).

Применим оператор ротора c∇× ко второму уравнению системы (1). Придём к первому

уравнению из системы уравнений Максвелла

∂B

+ c∇ × E = 0,

(4)

∂t

являющемуся обобщением закона индукции Фарадея.

Наиболее сложным является вывод второго уравнения из системы уравнений Максвелла,

содержащего плотность тока j и плотность тока смещения. Применим ко второму уравнению

системы (1) оператор производной вдоль кривой (u · ∇). Получим

(

)

∂E

∂u

+ (u · ∇)((u · ∇)(ρu)) -

· ∇ (ρu) = 0.

∂t

Используя далее вытекающее из системы (1) уравнение ∂u/∂t = (u·∇)u-(∇·u)u и известные

правила действия с оператором набла ∇, придём к уравнению

)

∂E

⁽ |u|²

-∇×

B + 4πj = 0,

(5)

∂t

в котором j имеет смысл вектора плотности электрического тока и

4πj = (b · ∇)u + u(∇ · b) - b(∇ · u) - ∇ × (u × (∇(u · a) - (a · ∇)u - a × (∇ × u))) +

+ ∇ × (u(u · (∇ × a))) - (((∇ · u)u - (u · ∇)u) · ∇)a, a = ρu, b = (u · ∇)ρu.

В случае |u| ≈ c и экспериментально определяемых σ и j уравнения (3)-(5) переходят в

классическую линейную систему уравнений Максвелла

1 ∂B

1 ∂E

4πj

∇×E=-

∇×B=

∇ · B = 0,

∇ · E = 4πσ.

(6)

c ∂t

Таким образом, систему уравнений (3)-(5) можно интерпретировать как систему обобщён-

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

эфира (1) инвариантны относительно преобразования Галилея. Одной из причин потери та-

кой инвариантности в классических линейных уравнениях Максвелла является линеаризация

обобщённых нелинейных уравнений Максвелла (3)-(5) при |u| ≈ c. Другой причиной являет-

ся применение операторов дифференцирования при переходе к векторам E и B. Последнее

приводит также к потере продольной компоненты в винтовой электромагнитной волне эфира

(фотоне)

u(t, r) = u₀ sin(ωt - ky)i_x + ci_y + u₀ cos(ωt - ky)i_z, ω = kc,

(7)

распространяющейся со скоростью света c в направлении y в эфире постоянной плотности,

где u₀- амплитуда поперечной компоненты скорости, а i_x, i_y, i_z - единичные базисные век-

торы декартовой системы координат. Нетрудно проверить, что вектор (7) является решением

системы уравнений (1), записанной в декартовой системе координат. Вычисляя теперь векторы

индукции магнитного поля и напряжённости электрического поля по формулам (2), получаем

B = c∇ × (ρu) = u₀ρω(sin(ωt - ky)i_x + cos(ωt - ky)i_z),

(8)

E = (u · ∇)(ρu) = u₀ρω(-cos(ωt - ky)i_x + sin(ωt - ky)i_z).

(9)

Векторы (8) и (9) ортогональны между собой, имеют равные амплитуды |E| = |B| = u₀ρω

и лежат в плоскости (x, z), ортогональной направлению движения эфирной волны (фотона).

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

283

Нетрудно проверить, что векторы (8) и (9) являются также решениями классической линейной

системы уравнений Максвелла (6) в случае отсутствия зарядов и токов. Следовательно, в этом

частном случае линейная система уравнений Максвелла (6) даёт верный результат в виде

плоской монохроматической циркулярно поляризованной электромагнитной волны.

В другом частном случае, когда решением системы уравнений (1) является продольная

волна возмущений плотности эфира

u(t, r) = ci_y, ρ(t, r) = ρ₀(1 + α sin(ωt - ky)), ω = kc, α ≪ 1,

(10)

вычисленные по формулам (2) векторы индукции магнитного поля и напряжённости электри-

ческого поля имеют вид

B = c∇ × (ρu) = 0, E = (u · ∇)(ρu) = -cωρ₀αcos(ωt - ky)i_y.

Следовательно, этот случай не описывается линейной системой уравнений Максвелла. Вполне

вероятно, что продольная волна возмущений плотности эфира (10) является волной ради-

антного электричества Николы Тесла. Вывод из системы (1) других известных физических

законов и уравнений можно найти в работах автора [5, 6].

Литература. 1. Левич В.Г. Курс теоретической физики. Т. 1. М., 1969. 2. Магницкий Н.А.

Математическая теория физического вакуума. М., 2010. 3. Магницкий Н.А. К электродина-

мике физического вакуума // Сложные системы. 2011. Т. 1. № 1. С. 83-91. 4. Magnitskii N.A.

Mathematical theory of physical vacuum // Comm. Nonlin. Sci. and Numer. Simul. 2011. V. 16.

P. 2438-2444. 5. Magnitskii N.A. Theory of compressible oscillating ether // Results in Phys. 2019.

V. 12. P. 1436-1445. 6. Магницкий Н.А. Теория сжимаемого осциллирующего эфира. М., 2021.

Е.И. Атамась, А.И. Роговский (Москва, Россия, ВМК МГУ). “О свойствах обобщения

относительного порядка для систем с запаздыванием” (11.10.2021).

Рассматриваем линейную стационарную управляемую систему с соизмеримыми запазды-

ваниями. Такая система имеет следующий вид:

∑

x(t) =

A_ix(t - iτ) +

B_iu(t - iτ), y(t) =

C_ix(t - iτ).

(1)

i=0

Здесь x(t) ∈ Rⁿ, y(t), u(t) ∈ R^l, матрицы A_i, B_i, C_i, i = 1, l, постоянны и имеют соответ-

ствующие размеры, τ > 0. Будем рассматривать квадратные системы, т.е. системы, размер-

ности входа и выхода которых одинаковы. Далее введём оператор запаздывания и, заменив

его символом Δ, получим систему в алгебраической форме [1]:

x = A(Δ)x + B(Δ)u, y = C(Δ)x,

(2)

где A(Δ), B(Δ), C(Δ) - полиномиальные матрицы размеров n×n, n×l, l×n соответствен-

но. Систему (2) будем называть также системой (A, B, C).

Определение 1. Вектор r = (r₁,... ,r_l)^т ∈ N^l называется вектором относительного

порядка системы (2), если выполнены следующие условия:

1) C_i∗(Δ)Ari-¹(Δ)B(Δ) = 0∗) и если r_i > 1, то C_i∗(Δ)Aj-1(Δ)B(Δ) = 0, где через C_i∗(Δ)

обозначена i-я строка матрицы C(Δ);

2) определитель матрицы

⎛

⎞

C1∗(Δ)Ar1-¹(Δ)B(Δ)

H(Δ) =^⎝

⎠

C_l∗(Δ)Arl-¹(Δ)B(Δ)

является ненулевым полиномом.

∗) т.е. указанная строка отлична от строки, целиком состоящей из нулевых полиномов.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

284

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

Если, кроме этого, определитель det H(Δ) обратим, т.е. является ненулевым числом, то

относительный порядок называется чистым.

Понятие относительного порядка имеет фундаментальное значение в математической тео-

рии управления, в частности, при решении задач наблюдения и обращения [2]. К сожалению,

не для всех систем относительный порядок определён, поэтому приходится рассматривать

обобщения этого понятия. В докладе приводятся некоторые результаты о свойствах таких

обобщений для систем с запаздыванием, близкие к результатам работы [3].

Определение 2. Вектор r = (r₁,... ,r_l)^т ∈ N^l называется вектором неполного относи-

тельного порядка (НОП ) системы (2), если для него выполнено условие 1) определения 1.

Определение 3. Вектор r = (r₁,... ,r_l)^т ∈ N^l называется вектором главного неполного

относительного порядка (ГНОП), если r - вектор НОП и для любых попарно различных

индексов i₁, i₂, . . . , i_k таких, что ri1 = ri2 = . . . = riq , полином

G[Ci1∗Ari1 -¹B, Ci₂∗Ari2 -¹B, . . . , Ciq Ariq -¹B]

ненулевой, где G[v₁, v₂, . . . , v_p] = det(v_ivтj)ki,j=1 - определитель матрицы Грама.

Определение 4. Систему {A,B,C} назовём слабо приводимой, если при любой унимо-

дулярной матрице T для системы {A, B, T C} определён вектор НОП.

Утверждение 1. Система {A, B, C} является слабо приводимой тогда и только тогда,

когда матрица

V (Δ) = [CB, CAB, . . . , CAn-1B]

имеет полный ранг l при всех Δ, за исключением, быть может, конечного числа значений.

Утверждение 2. Пусть система {A, B, C} является слабо приводимой. Тогда найдётся

такая унимодулярная матрица T, что система {A, B, T C} имеет вектор ГНОП.

Утверждение 3. Пусть r - вектор ГНОП слабо приводимой системы {A, B, C}. Тог-

да при любой унимодулярной матрице T для системы {A,B,TC} определён вектор НОП

r, причём выполнено покомпонентное неравенство ord (r) ≥ ord (r), где ord (r) обозначает

вектор с теми же компонентами, что и у r, но упорядоченными по невозрастанию.

Работа выполнена при финансовой поддержке гранта Президента Российской Федерации

для молодых учёных-кандидатов наук (проект МК-4905.2021.1.1).

Литература. 1. Morse A.S. Ring models for delay-differential systems // IFAC Proc. Volumes.

1974. V. 7. № 1. P. 439-445. 2. Ильин А.В., Атамась Е.И., Фомичев В.В. О приведении систем

с запаздыванием к форме с выделением нулевой динамики // Докл. РАН. 2018. Т. 480. № 1.

С. 11-15. 3. Фомичев В.В., Краев А.В., Роговский А.И. Обобщение понятия относительного

порядка и его свойства // Дифференц. уравнения. 2016. Т. 52. № 8. С. 1099-1108.

В.В. Фомичев, Н.И. Денисова (Москва, Россия, ВМК МГУ, ИПМ им. М.В. Келдыша).

“Наблюдатели для квазистационарных систем при внешних возмущениях” (15.11.2021).

Хорошо известны подходы к построению наблюдателей линейных стационарных систем [1].

Как правило, для построения таких наблюдателей требуется достаточная гладкость парамет-

ров системы. Большой интерес вызывает случай линейных нестационарных систем. С возника-

ющими проблемами можно бороться с помощью, например, адаптивных наблюдателей [2]. Для

некоторых линейных нестационарных систем можно построить даже экспоненциально устой-

чивый наблюдатель [3]. Значительно интереснее обстоит дело в задаче о синтезе наблюдателей

для систем с возмущениями при непрерывности или однократной дифференцируемости пара-

метров системы. Если для стационарного случая известны различные методы решения этой

задачи (при различных предположениях о системе), то в нестационарном случае продвижений

в её решении очень мало.

В докладе предпринята попытка перенести некоторые результаты, полученные для линей-

ных стационарных систем с возмущением, на нестационарный случай.

Итак, рассматриваем систему

x = A(t)x + D(t)f, y = C(t)x,

(1)

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

285

где x ∈ Rⁿ, u ∈ R^k, f ∈ R^m, y ∈ R^l, известные матрицы A(t), C(t), D(t) - гладкие по t и

имеют соответствующий порядок, возмущение f мажорируется константой.

Требуется по информации об известном входе u(t) и известном выходе y(t) построить

асимптотическую оценку x(t) вектора состояния x(t).

Далее предполагаем, что (1) является гипервыходной системой, т.е. системой, для которой

l > m (число известных выходов больше числа неизвестных входов). Предполагаем также,

что rank (CD) = m, т.е. матрица CD имеет полный ранг. Кроме того, считаем, что последнее

условие выполняется “равномерно по выходам”, т.е. что выход можно разделить:

(

)

(

)

′

(

)

y^′

det

C^′(t)D(t)

для всех t.

C^′′ x⁼

y^′′

Применяя сначала стационарное преобразование координат с произвольной матрицей F

такой, что

(

)

(

)

x^′

T =

det T (t) = 0 для всех t, T x =

⇔x=T^′x^′ +T^′′y^′,

C^′

y^′

а затем делая нестационарную замену x^′ = x^′ - F D(C^′D)-1y^′, получаем систему

x^′ = A₁₁(t)x^′ + A₁₂(t)y^′,

y′ = A21(t)x′ + A22(t)y′ + A24f.

Рассмотрим “запасные” выходы, которые после преобразований принимают вид

y^′′ = C^′′T^′x^′ + C^′′(T^′FD(C^′D)-1 + T^′′)y^′,

откуда получаем итоговую систему без явного влияния возмущения с известным выходом, для

которой можно пробовать строить наблюдатель:

x^′ = A₁₁(t)x^′ + A₁₂(t)y^′,

Cx^′,

(2)

где

C =C^′′T^′.

Для иллюстрации предложенного подхода рассмотрим квазистационарный случай системы

третьего порядка, когда n = 3, C и D - постоянные матрицы. Тогда

⎛

⎞

⎛

⎞

(

)

a₁₁(t) a₁₂(t) a₁₃(t)

A(t) =^⎝a₂₁(t) a₂₂(t) a₂₃(t)⎠, D = ⎝d2⎠, C =

(3)

a₃₁(t) a₃₂(t) a₃₃(t)

d₃

Лемма 1. Система (1), (3) наблюдаема тогда и только тогда, когда a²¹³(t) + a²²³(t) = 0.

В силу вида исходной матрицы C можно выполнить невырожденное стационарное преоб-

разование

⎛

⎞

(

)

(

)

x^′

Tx=

1⎠x=

C^′ x=⎝0

y^′

В итоге получим систему (2) с матрицами

(

)

-d₂a₁₂(t) + a₂₂(t)

-d₂a₁₃(t) + a₂₃(t)

A₁₁(t) =

−d₃a₁₂(t) + a₃₂(t)

-d₃a₁₃(t) + a₃₃(t)

(

)

a₂₁(t) - d₂a₁₁(t) - d²²a₁₂(t) + d₂a₂₂(t) - d₂d₃a₁₃(t) + d₃a₂₃(t)

A₁₂(t) =

a₃₁(t) - d₃a₁₁(t) - d₂d₃a₁₂(t) + d₂a₃₂(t) - d²³a₁₃(t) + d₃a₃₃(t)

C = (1, 0).

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

286

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

Для построения наблюдателя необходима наблюдаемость пары {A₁₁(t)

C }.

Лемма 2. Пара {A₁₁(t)

C} наблюдаема, когда a₂₃(t) ≡ d₂a₁₃(t).

Замечание. В случае ненаблюдаемости исходной системы (3) пара {A₁₁(t)

C} также не

будет наблюдаемой.

Итак, имеем наблюдаемую систему

x^′ = A₁₁(t)x^′ + A₁₂(t)y^′,

Cx^′.

Строим наблюдатель

x=A11(t)x + A₁₂(t)y^′ - L

Cx^′ - y), где L = (l₁,l₂)^т.

Для матрицы

(

)

ã₁₁(t) - l₁(t)

ã₁₂(t)

A_L(t) = A₁₁(t) -

C =

ã₂₁(t) - l₂(t)

ã₂₂(t)

рассмотрим функцию Ляпунова V = e^тP e,

ė = A_Le с матрицей P = diag[h₁,h₂], h₁,h₂ > 0.

Теорема. Задача наблюдения может быть решена при ã₂₂(t) = a₃₃(t) - d₃a₁₃(t) > 0.

В этом случае l₁ = l₁(t) = ã₁₁(t) + δ, δ > 0, а функция l₂ выбирается достаточно малой.

Работа выполнена при финансовой поддержке Российского научного фонда (проект 22-21-

00288).

Литература. 1. Fomichev V.V., Vysotskii A.O. Algorithm for designing a cascade asymptotic

observer for a system of maximal relative order // Differ. Equat. 2019. V. 67. № 3. P. 553-560.

2. Куок Дат Во, Бобцов А.А. Адаптивный наблюдатель переменных состояния линейных не-

стационарных систем с параметрами, заданными не точно // Автоматика и телемеханика.

2020. № 12. С. 100-110. 3. Зубер И.Е. Синтез экспоненциально устойчивого наблюдателя для

линейных нестационарных систем с одним выходом // Автоматика и телемеханика. 1995. № 5.

С. 42-49.

Д.В. Злобин (Москва, Россия, ВМК МГУ). “Анализ эффективности автоматически сге-

нерированного кода для вычисления матриц Якоби и Гессе в задачах нелинейного програм-

мирования” (13.12.2021).

Пусть нелинейный динамический объект задан в виде системы обыкновенных дифферен-

циальных уравнений

x(t) = f(x(t), u(t)),

x(0) = x₀, t ∈ [0, T ],

а на его фазовый вектор наложены нелинейные ограничения

GL ≤ G(x(t), u(t), T ) ≤ GU ,

где f(x, u), x(t), x₀ ∈ Rⁿ; u(t) ∈ R^m; G(x, u, T ), G_L, G_U ∈ R^r; T ∈ R₊.

Для этого объекта требуется решить задачу оптимального управления (ОУ) с нелинейным

функционалом качества:

min F (x(t), u(t), T ),

u(t)∈U

где F (x, u, T ) ∈ R, а U - множество допустимых управлений.

Один из способов решения такой задачи состоит в переходе от неё к задаче нелинейного

программирования (НЛП) при помощи дискретизации по времени фазового вектора, управ-

ления, уравнений динамики, ограничений и функционала качества [1, с. 91-216]. Полученная

задача НЛП решается одним из численных методов. Если ограничения и функционал качества

задачи НЛП дважды непрерывно дифференцируемы, то можно использовать решатель IpOpt

[2], являющийся открытой эффективной реализацией метода внутренней точки.

Для применения решателя IpOpt пользователь должен предоставить реализацию на языке

C++ функций для вычисления функционала качества, ограничений, градиента функциона-

ла качества, ненулевых элементов матрицы Якоби для ограничений и ненулевых элементов

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

287

матрицы Гессе для лагранжиана задачи. Эти функции могут быть написаны вручную или сге-

нерированы при помощи специального программного обеспечения на основе высокоуровневого

описания задачи НЛП.

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

полученных при помощи пакетов sympy2ipopt [3, 4] и CasADi [5]. Пакет CasADi, основанный на

идеях автоматического дифференцирования, ориентирован на общий случай задач НЛП. В от-

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

структуры задач НЛП, полученных при дискретизации задач ОУ. Это делает результирующий

код более компактным и потенциально более быстрым.

Для проведения тестирования были выбраны четыре тестовые задачи ОУ из набора COPS3

[6]: “robot”, “rocket”, “steering”, “glider”. Применялся компилятор gcc 9.3.0, вычисления проводи-

лись на компьютере с процессором Intel ® Core^TM i7-9750H.

Результаты экспериментов приведены в табл. 1 и 2. В первом случае сборка осуществлялась

с использованием флага оптимизации “-O3”, а во втором - флага “-Ofast”.

Таблица 1. Время (в мкс), затраченное на вычисление автоматически сгенерированных функций,

сборка осуществлялась с флагом “-O3”

sympy2ipopt

CasADi

Функции

robot rocket steering glider robot rocket steering glider

Функционал качества

Ограничения

118

144

146

Градиент

Матрица Якоби

430

364

370

144

140

385

Матрица Гессе

1050

223

2359

404

183

770

Таблица 2. Время (в мкс), затраченное на вычисление автоматически сгенерированных функций,

сборка осуществлялась с флагом “-Ofast”

sympy2ipopt

CasADi

Функции

robot rocket steering glider robot rocket steering glider

Функционал качества

Ограничения

117

Градиент

Матрица Якоби

265

136

127

344

Матрица Гессе

297

177

711

По сравнению с фагом “-O3” флаг “-Ofast” включает дополнительные оптимизации, кото-

рые не соответствуют спецификациям IEEE и ISO для математических функций. В частности,

не устанавливается переменная “errno” после выполнения математических операций, которые

занимают одну инструкцию, предполагается, что в процессе выполнения не будет ошибок де-

ления на нуль, переполнения и т.п.

Анализ результатов позволяет сделать следующие выводы. Во-первых, использование па-

кета sympy2ipopt даёт выигрыш в скорости, за исключением вычисления элементов матриц

Якоби и Гессе при сборке с флагом “-O3” для задач “robot” и “glider”. Во-вторых, использование

флага компилятора “-Ofast” существенно увеличивает скорость выполнения кода для вычис-

ления элементов матриц Якоби и Гессе, сгенерированного при помощи пакета sympy2ipopt.

Степень влияния флага “-Ofast” объясняется способом представления математических вы-

ражений, который использует программный пакет. Пакет CasADi основан на представлении

выражений в виде направленного графа, а пакет sympy2ipopt использует средства библиотеки

SymPy [7] (символьные вычисления на языке Python), в которой выражения представляются в

виде дерева. Из-за этого многократно вычисляются одни и те же подвыражения, что приводит

к большим накладным расходам на проверки, если не указан флаг компилятора “-Ofast”.

Полученные результаты показывают, что изменение представления математических выра-

жений в пакете sympy2ipopt позволит добиться высокой скорости выполнения сгенерирован-

ного кода без использования “агрессивной” оптимизации “-Ofast”.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022

288

О СЕМИНАРЕ ПО ПРОБЛЕМАМ НЕЛИНЕЙНОЙ ДИНАМИКИ

Литература. 1. Betts J.T. Practical Methods for Optimal Control and Estimation Using

Nonlinear Programming. Philadelphia, 2010. 2. Wächter A., Biegler A. On the implementation of a

primal-dual interior point filter line search algorithm for large-scale nonlinear programming // Math.

Programming. 2006. V. 106. P. 25-57. 3. Злобин Д.В. Использование разреженной структуры

матриц Якоби и Гессе для ускорения численного решения задач оптимального планирования

траекторий // Интеллектуальные системы. Теория и приложения. 2021. T. 25. № 4. C. 285-288.

4. GitHub. symp2ipopt. 2021. Режим доступа: https://github.com/zlobin-d/sympy2ipopt (Дата

обращения: 13.12.2021). 5. Andersson J.A.E., Gillis J., Horn G., Rawlings J.B., Diehl M. CasADi -

a software framework for nonlinear optimization and optimal control // Math. Progr. Comput. 2019.

V. 11. P. 1-36. 6. Dolan E.D., Moré J.J., Munson T.S. Benchmarking Optimization Software with

COPS 3.0 // Technical Memorandum. ANL/MCS-TM-273. Argonne National Laboratory, 2004.

7. Meurer A., Smith Ch.P., Paprocki M.,Čert´ik O. SymPy: symbolic computing in Python // Peer

J. Comp. Sci. 2017. V. 3. P. e103.

Р.Р. Бегишев, А.В. Ильин, А.И. Роговский (Москва, Россия, ВМК МГУ). “Достаточ-

ные условия существования периодических решений у нелинейной системы третьего порядка”

(20.12.2021).

Рассматривается нелинейная система третьего порядка

1 = z2,

2 = -z1 + μ(-γz2 + (Δ - z1)f(z3)),

3 = -βz3 + α(cos(ζ arctg(z1)) + 1),

(1)

где числа α, β, γ, Δ, ζ положительны, а f(z) - возрастающая, непрерывно дифференци-

руемая и ограниченная на R функция. Требуется найти параметры системы, при которых у

неё существует нестационарное периодическое решение.

Отметим, что вопрос о существовании периодических решений у системы (1) с конкретной

функцией f(z₃) исследовался в работах [1-4]. В настоящем докладе результаты работы [4]

обобщаются на случай произвольной функции f(z₃).

Используя результаты, приведённые в монографии [5, гл. XIV], доказана следующая

Теорема. Пусть функция H(p), p ∈ R, задана равенством

∫2π

H(p) = γπp + (Δ - p cos(s))f(ω(s, p)) sin(s) ds,

где

(

∫

2π

∫

)

ω(s, p) = e^-βs

e^βq cos(ζ arctg (p cos(q)))dq + α e^βq cos(ζ arctg (p cos(q)))dq

e2πβ -1

и числа α, β, γ, Δ, ζ таковы, что найдётся положительное p, для которого выполняют-

ся соотношения H(p) = 0 и dH(p)/dp = 0. Тогда найдётся такое μ > 0, что система (1)

имеет отличное от константы периодическое решение при всех μ ∈ [0,μ].

Литература. 1. Fomichev V.V., Iline A.V., Rogovskii A.I., Todorov G.D., Sofronov Ya.P.

Search for periodic regimes in an energy-harvester model by simulation // Comput. Math. and

Model. 2020. V. 31. № 1. P. 293-307. 2. Фурсов А.С., Митрев Р.П., Крылов П.А., Тодоров Т.С.

О существовании периодического режима в одной нелинейной системе // Дифференц. уравне-

ния. 2021. Т. 57. № 8. С. 1104-1115. 3. Фурсов А.С., Митрев Р.П., Крылов П.А., Тодоров Т.С. О

существовании колебательных режимов в одной нелинейной системе с гистерезисами // Диф-

ференц. уравнения. 2020. Т. 56. № 8. С. 1103-1121. 4. Фомичев В.В., Ильин А.В., Рогов-

ский А.И., Бегишев Р.Р., Митрев Р.П., Тодоров Т.С. О существовании периодических решений

у нелинейной системы третьего порядка // Дифференц. уравнения. 2021. Т. 57. № 10. С. 1367-

1383. 5. Коддингтон Э.А., Левинсон Н. Теория обыкновенных дифференциальных уравнений.

М., 1958.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58

№2

2022