ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2023, том 59, № 6, с.814-827

ЧИСЛЕННЫЕ МЕТОДЫ

УДК 519.622.2+681.5.015

АЛГОРИТМ ПОДВИЖНОГО ОКНА

ДЛЯ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ

ДИНАМИЧЕСКИХ СИСТЕМ С ПРЯМОУГОЛЬНЫМИ

И ЭЛЛИПСОИДНЫМИ ОБЛАСТЯМИ

НЕОПРЕДЕЛЁННОСТИ ПАРАМЕТРОВ

Решена задача параметрической идентификации динамических систем с прямоугольными

и эллипсоидными областями неопределённости параметров для случая, когда эксперимен-

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

систем в каждый момент времени является параметрическим множеством. Построена це-

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

степень отклонения параметрических множеств состояний от экспериментальных интер-

вальных оценок. Для минимизации целевой функции разработан алгоритм подвижного

окна, относящийся к градиентным методам. В его основе лежит алгоритм адаптивной ин-

терполяции, позволяющий в рамках заданной области неопределённости параметров (ок-

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

Продемонстрирована эффективность и работоспособность предлагаемого алгоритма.

DOI: 10.31857/S0374064123060110, EDN: FIELWZ

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

разброса в значениях [1, 2]. В основном это связано с погрешностями проводимых измерений

или с внешними факторами, влияние которых нельзя исключить. Обычно эксперимент повто-

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

значности в данных. Наиболее простым описанием неопределённости является интервальная

оценка возможных значений.

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

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

таким образом, чтобы модель наилучшим образом воспроизводила эксперимент. Традиционно

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

кривой от экспериментальных точек.

Применение интервального аппарата [3-6] в задачах параметрической идентификации свя-

зано с предположением, что параметры модели могут быть интервальными. Преимущество в

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

ресующие величины, в отличие от классических моделей, которые аппроксимируют их. В рабо-

те [7] представлен итеративный интервальный метод для прогнозирования границ параметров

заданной модели в условиях неопределённости измеренных данных. В исследовании [8] пред-

лагается структурно-параметрическая идентификация линейных, либо допускающих линеа-

ризацию, динамических объектов по интервальным исходным данным. Статья [9] посвящена

развитию методов параметрической идентификации на основе нейронных сетей, применяе-

мых для решения обратных задач теплофизики в условиях интервальной неопределённости

параметров.

Ранее авторами был разработан интерполяционный подход к решению задачи интерваль-

ной параметрической идентификации для случая точечных экспериментальных данных [10].

В данной работе выполнено обобщение подхода на случай интервальных экспериментальных

данных и различных форм области неопределённости параметров динамической системы; рас-

смотрены динамические системы с прямоугольными и эллипсоидными областями неопределён-

ности параметров. Будет показано, что данные системы можно представить как системы с ин-

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

814

АЛГОРИТМ ПОДВИЖНОГО ОКНА

815

является параметрическим множеством. Задача параметрической идентификации заключает-

ся в нахождении такой области неопределённости параметров, в которой параметрические мно-

жества состояний будут содержать соответствующие экспериментальные интервальные оценки

или минимизировать отклонение от них.

Как и при обычной параметрической идентификации, здесь составляется целевая функция

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

модельных множеств от экспериментальных данных. Для минимизации целевой функции раз-

работан алгоритм подвижного окна, относящийся к градиентным методам. В его основе лежит

алгоритм адаптивной интерполяции [11, 12], позволяющий в рамках заданной области неопре-

делённости параметров (окна) в явном виде получать параметрические множества состояний

динамической системы. Ключевое преимущество алгоритма подвижного окна заключается в

том, что он стремится получить область неопределённости параметров как можно меньшего

объёма, что является важным на практике.

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

прямоугольными и эллипсоидными областями неопределённости параметров. Наиболее часто

математические модели динамических систем задаются в виде систем обыкновенных диффе-

ренциальных уравнений (ОДУ), поэтому рассматриваются именно они. В п. 2 формулируется

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

даны в виде интервалов, и в общем виде строится целевая функция. В следующем пункте

анализируется целевая функция и доказывается утверждение о том, что при определённых

свойствах параметрического множества состояний вместо экспериментальных интервальных

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

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

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

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

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

1. Прямая задача. Рассматривается задача Коши для системы ОДУ, состоящей из n

уравнений с m-мерной областью неопределённости Ω в начальных условиях:

dy_i(t)

= f_i(y₁(t),y₂(t),... ,y_n(t)), i = 1,n, t ∈ [t₀,t_N],

(y₁(t₀), y₂(t₀), . . . , y_m(t₀))^т ∈ Ω, y_i(t₀) = y0i, i = m + 1, n,

(1)

где t₀ - начальный момент времени; t_N - конечный момент времени; вектор-функция f =

= (f₁, f₂, . . . , f_n)^т удовлетворяет всем условиям, обеспечивающим единственность и существо-

вание решения при всех (y₁(t₀), y₂(t₀), . . . , y_m(t₀))^т ∈ Ω. Любую неавтономную систему ОДУ

с параметрами можно преобразовать к виду (1) путём добавления фиктивных уравнений.

В качестве области неопределённости Ω рассматривается или прямоугольный параллеле-

пипед, или эллипсоид, которые могут быть сориентированы определённым образом в прост-

ранстве. Область Ω задаётся в параметрическом виде

{

}

∏

Ω = z = (z₁,z₂,...,z_m)^т : z = z^c +

Rϕi,j Sg(ξ), ξ ∈ E

(2)

i=1 j=i+1

где zc = (zc1, zc2, . . . , zcm)т - центр области; Rϕi,j - матрица поворота на угол ϕi,j отно-

сительно плоскости, образованной осями i и j; S = diag {s1,s2,... ,sm} - диагональная

матрица масштабирования; E = {ξ = (ξ₁, ξ₂, . . . , ξ_m)^т : ξ_i ∈ [0, 1]} - единичный куб; g(ξ) =

= (g₁(ξ), g₂(ξ), . . . , g_m(ξ))^т - отображение единичного куба в куб с центром в начале координат

и длиной стороны 2 (тогда g_i(ξ) = 2ξ_i -1, i = 1, m) или в единичный шар. В последнем случае

в качестве g(ξ) удобно использовать переход из сферической системы координат в декартову

систему координат [13]:

∏

g₁(ξ) = ξ₁ cos(2πξ₂)

sin(πξ_j ), g₂(ξ) = ξ₁ sin(2πξ₂)

sin(πξ_j),

j=3

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

№6

2023

816

МОРОЗОВ, РЕВИЗНИКОВ

∏

g_i(ξ) = -ξ₁ cos(πξ_i)

sin(πξ_j ), i = 3, m.

j=i+1

Систему (1) можно записать в виде системы ОДУ с интервальными параметрами

dy_i(t)

= f_i(y₁(t),y₂(t),... ,y_n(t)), i = 1,n, t ∈ [t₀,t_N],

y_i(t₀) = z_i(z^c,s,ϕ,ξ), i = 1,m,

y_i(t₀) = y0i, i = m + 1,n,

(3)

∏_m ∏_m

где ξ ∈ [0, 1]^m; z(z^c, s, ϕ, ξ) = (z₁(z^c, s, ϕ, ξ), . . . , z_m(z^c, s, ϕ, ξ))^т = z^c +

Rϕi,j S ×

i=1

j=i+1

×g(ξ) - преобразование единичного куба в область определённой формы; z^c =(zc1,zc2,... ,z^cm)^т,

s = (s₁,s₂,...,s_m)^т, ϕ = (ϕ1,2,...,ϕ_i,j,...,ϕm-1,m), i < j, - параметры, полностью задающие

область неопределённости. Для сокращения записи обозначим наборы параметров z^c, s и ϕ

через θ = (θ₁, θ₂, . . . , θnp )^т, n_p = 2m + (m - 1)m/2.

Решение системы (3) в каждый момент времени t_k является параметрическим множеством

Y k = {y^k(ξ₁,ξ₂,...,ξ_m) : ξ_i ∈ [0,1], i = 1,m}.

(4)

Отметим, что на практике часто используется внешняя интервальная оценка множества (4),

которую можно рассматривать как одномерные проекции Y^k на каждое фазовое измерение:

W^k = [wk1,wk1] × [wk2,wk2] × ... × [w^kn,w^kn], w^ki = min

[y^ki(ξ)],

ξ∈[0,1]^m,

w^ki = max

[y^ki(ξ)], i = 1, n.

ξ∈[0,1]^m

Обратим внимание, что Y^k в общем случае зависит от параметров θ :

Y k(θ) = {y^k(z(θ,ξ)) : ξ ∈ [0,1]^m}.

В дальнейшем эти параметры могут опускаться.

Цель алгоритма адаптивной интерполяции заключается в построении для каждого момента

времени t_k вектор-функции P^k(ξ), интерполирующей y^k(z(θ, ξ)), ξ ∈ [0, 1]^m. В начальный

момент времени P⁰ определяется явным образом:

P⁰(ξ) = (z₁(θ,ξ),z₂(θ,ξ),... ,z_m(θ,ξ),y0m+1,... ,y0n)^т.

Построение Pk+1 по P^k сводится к интерполяции неявной функции

F (ξ₁, ξ₂, . . . , ξ_m) = (y₁(tk+1), y₂(tk+1), . . . , y_n(tk+1))^т,

заданной в виде системы ОДУ

dy_i(t)

= f_i(y₁(t),y₂(t),... ,y_n(t)), t ∈ [t_k,tk+1],

y_i(t_k) = P^ki(ξ₁,ξ₂,... ,ξ_m), i = 1,n.

(5)

Интерполяционный полином строится по некоторому набору узлов, образующих интерпо-

ляционную сетку. При построении Pk+1 сначала выполняется перенос решений, находящихся

в узлах сетки, с k-го слоя на (k + 1)-й слой. Далее в зависимости от погрешности интерполя-

ции происходит адаптация сетки путём добавления или удаления узлов. Интерполяционный

полином P может быть любым, только необходимо, чтобы была возможность контролиро-

вать погрешность интерполяции. В приведённых далее численных примерах использовалась

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

817

интерполяция на адаптивных разреженных сетках с нелинейным базисом [14-16] (подробное

описание см. в работах [17-19]).

Задание области неопределённости Ω в параметрическом виде (2) необходимо для удобства

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

сетка всегда строится над единичным кубом.

Отметим, что так как P^k(ξ) является явной функцией, то можно аналитически вычислить

производные решения по ξ, а также по z и θ. Для определения dy^ki(z(θ, ξ))/dz_j , j = 1, m,

i = 1,n, необходимо решить n систем линейных алгебраических уравнений (СЛАУ) с оди-

наковой матрицей коэффициентов, которая является матрицей Якоби преобразования z(θ, ξ)

относительно ξ :

∑

dP^k(ξ)

dy^ki(z(θ, ξ)) dz_l(θ, ξ)

j = 1,m, i = 1,n,

(6)

dξ_j

dz_l

dξ_j

l=1

а производные по θ уже вычисляются через dy^ki(z(θ, ξ))/dz_j :

∑

dy^k(z(θ, ξ))

dy^ki(z(θ, ξ)) dz_l(θ, ξ)

j = 1,n_p, i = 1,n.

(7)

dθ_j

dz_l

dθ_j

l=1

Выражения для производных являются явными функциями, что позволяет выполнять анализ

исходной системы (3) в области неопределённости без дополнительного интегрирования систем

ОДУ для вычисления градиента.

2. Интервальная задача параметрической идентификации. Пусть известны N экс-

периментальных оценок фазовых переменных для различных моментов времени t_k, k = 1, N :

Y k = [ŷk1, ŷk1] × [ŷk2, ŷk2] × ... × [ŷ^kn, ŷ^kn].

С геометрической точки зрения множеств

Y k является n-мерным прямоугольным парал-

лелепипедом (или, по-другому, брусом). Интервальная задача параметрической идентифика-

ции заключается в определении таких параметров θ области неопределённости Ω начальных

условий системы (3), чтобы

Y k ⊂ Y ^k, k = 1,N, или чтобы степень непересечения множеств

Y k и Y ^k, k = 1,N, была минимальна. Если экспериментальная информация о диапазонах

значений известна не для всех фазовых переменных, то при поиске неизвестных параметров

такие фазовые переменные не учитываются.

Сформулируем задачу минимизации отклонения модельного решения от эксперименталь-

ных данных. Минимизируется следующая целевая функция:

∑

J (θ) =

ρ(Y^k(θ)

Y k),

k=1

где ρ(Y^k

Y k) характеризует степень непересечения множеств

Y k и Y ^k, и если

Y k ⊂ Y ^k, то

ρ(Y^k

Yk) = 0. В частности, ρ(Y^k

Y k) может соответствовать квадрату расстояния между

множеством Y^k и самой удалённой от него точкой множества

ρ(Y^k

Y k) = max ρ_sp(Y ^k, ŷ^k)

(8)

ŷk

Y^k

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

ŷk₁

ŷk₂

ŷkn

∫

ρ(Y^k

Yk) =

ρ_sp(Y^k, ŷ^k)dŷ^kndŷkn-1 ... dŷk1,

(9)

ŷk

ŷkn

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

№6

2023

818

МОРОЗОВ, РЕВИЗНИКОВ

где ŷk = (ŷk1, ŷk2, . . . , ŷkn)т, ρsp(Yk, ŷk) = min

∥y^k - ŷ^k∥² - квадрат расстояния между мно-

y^k∈Y^k

жеством Y^k и точкой ŷk. Традиционно в качестве нормы ∥ · ∥ используется взвешенная

евклидова норма

∑

∥a∥² = p_ia2i,

i=1

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

изложения будем полагать p_i = 1). С учётом этого ρ_sp можно записать следующим образом:

∑

ρ_sp(Y^k(θ), ŷ^k) = min

[y^ki(z(θ, ξ)) - ŷ^ki]².

(10)

ξ∈[0,1]^m

i=1

Отметим, что в данной постановке задачи не требуется, чтобы получающаяся область

неопределённости Ω была минимальной. Тем не менее если очевидным образом Ω можно

уменьшить, то эта процедура будет выполняться.

3. Построение целевой функции. Рассмотрим вопрос вычисления ρ(Y^k(θ)

Y k). В слу-

чае с (8) необходимо решать максиминную задачу:

∑

ρ(Y^k(θ)

Y k) = max

min

[y^ki(z(θ, ξ)) - ŷ^ki]².

(11)

ŷk

Y^k

ξ∈[0,1]^m

i=1

Решение подобных задач требует применения специального математического аппарата [20, 21].

Однако за счёт того, что Y^k и

Y k могут обладать определёнными свойствами, задачу (11)

можно упростить.

Теорема. Если множество Y^k является выпуклым, то максимум по ŷ^k в задаче (11)

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

Yk:

Y k) = {(ŷk1 + (ŷk1 - ŷk1)j₁, ŷk2 + (ŷk2 - ŷk2)j₂,... , ŷ^kn + (ŷ^kn - ŷ^kn)j_n) : j_i ∈ {0,1}}.

(12)

Доказательство. От противного предположим, что максимум достигается в точке a ∈

∈

Yk \V

Yk). Возьмём две любые диаметрально противоположные точки b

Yk и c

относительно точки a: a = (b + c)/2, в которых максимум не достигается. Отметим частные

случаи: если одна из точек b или c, также как и a, является точкой максимума, то дальней-

шие рассуждения будут аналогичными, а если абсолютно во всех диаметрально противопо-

ложных точках достигается максимум, то он будет достигаться и в граничных точках V

Y k).

Обозначим через b^′ и c^′ точки множества Y^k, ближайшие к соответствующим точкам b

и c:

b^′ = arg min ∥y^k - b∥², c^′ = arg min ∥y^k - c∥²,

y^k∈Y^k

а через a^′ их линейную комбинацию a^′ = (b^′ + c^′)/2, которая тоже принадлежит Y^k по

свойству выпуклости. Так как a - точка максимума, то ∥b^′ - b∥² < ρ_sp(Y^k, a) и ∥c^′ - c∥² <

< ρ_sp(Y ^k,a); кроме того, по определению расстояния между множеством и точкой как мини-

мального имеем дополнительное неравенство ∥a^′ - a∥² ≥ ρ_sp(Y^k, a). Таким образом,

∥b^′ - b∥ < ∥a^′ - a∥ и

∥c^′ - c∥ < ∥a^′ - a∥.

Подставим выражения для a и a^′ :



₁



₁



∥b^′ - b∥ <^

(b^′ + c^′) -

(b + c)

∥c^′ - c∥ <

b^′ + c^′) -

(b + c)

^,

^2(

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

819

далее перегруппируем слагаемые

∥b^′ - b∥ <

∥(b^′ - b) + (c^′ - c)∥,

∥c^′ - c∥ <

∥(b^′ - b) + (c^′ - c)∥

и сложим неравенства

∥b^′ - b∥ + ∥c^′ - c∥ < ∥(b^′ - b) + (c^′ - c)∥.

В результате получаем противоречие с неравенством треугольника. Следовательно, в окрест-

ности точки максимума a не может существовать диаметрально противоположных точек b

и c, в которых максимум не достигается. Такое может быть только в случае, если a ∈ V

Y k).

Теорема доказана.

Количество точек в множестве V

Y k) конечно и равно 2ⁿ, т.е. если множество Y ^k выпук-

лое, то вычисление (11) сводится к решению 2ⁿ задач минимизации. Требование выпуклости

Y k не всегда может выполняться, однако если множество Ω изначально является эллипсо-

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

преобразования

y₁ = x₁(x²¹ + x²²), y₂ = x₂(x²¹ + x²²)

для квадрата и круга с центром в начале координат. Образ квадрата не будет являться вы-

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

В случае с (9) вычисление ρ(Y^k

Y k) также будет сводиться к решению определённого

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

честве ρ(Y^k

Y k) для практических задач предлагается использовать некоторый компромисс

между (8) и (9):

∑

ρ(Y^k(θ)

Yk)=

min

[y^ki(z(θ, ξ)) - ŷ^ki]².

(13)

ξ∈[0,1]

ŷk ∈V

Y^k)

i=1

С учётом (12) и (13), а также замены y^k(z(θ,ξ)) на P^k(ξ), целевая функция примет вид

∑

J (θ) =

min

[P^ki(ξ) - ŷ^kj,i]²,

(14)

ξ∈[0,1]

k=1 j∈{0,1}ⁿ

i=1

где j = (j₁, j₂, . . . , j_n) - мультииндекс, состоящий из 0 и 1,

ŷkj = (ŷkj,1, ŷkj,2, . . . , ŷ^kj,n)т - j-я

вершина прямоугольного параллелепипеда

Yk,

ŷkj,i = (1 - j_i)ŷ^ki - j_i ŷ^ki, i = 1, n.

4. Алгоритм подвижного окна. Рассмотрим вопрос минимизации целевой функции

(14). Для вычисления её значения необходимо решить набор задач минимизации по ξ для

явных функций:

∑

J^kj(z(θ,ξ)) =

[P^ki(ξ) - ŷ^kj,i]², j ∈ {0, 1}ⁿ, k = 1, N .

(15)

i=1

Здесь удобно применить градиентные методы [22, 23], так как аналитически возможно постро-

ить выражение для градиента:

dJ^kj(z(θ, ξ))

∑

dP^ki(ξ)

[P^ki(ξ) - ŷ^kj,i]

p = 1,m.

dξ_p

i=1

Обозначим через ξ^kj точку минимума (15):

ξ^kj = arg minJkj(z(θ,ξ)).

ξ∈[0,1]^m

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

№6

2023

8^∗

820

МОРОЗОВ, РЕВИЗНИКОВ

С учётом соотношений (6) и (7) можно также построить градиент для основной целевой функ-

ции (14):

∑

∑ n∑

∑

dy^ki(z(θ, ξ^kj)) dz_l(θ, ξ^kj)

dJ(θ)

[P^ki(ξ^kj) - ŷ^kj,i]

p = 1,...,n_p,

dθ

dz_l

dθ_p

k=1 j∈{0,1}ⁿ i=1

l=1

где dy^ki(z(θ, ξ^kj))/dz_l, l = 1, m, находятся из решения СЛАУ

∑

dy^ki(z(θ, ξ^kj)) dz_l(θ, ξ^kj)

dP^ki(ξ^kj)

p = 1,m.

dz_l

dξ_p

dξ

l=1

Однако, как было отмечено ранее при постановке задачи, размер получающейся области Ω

никак не учитывается в целевой функции. Применение классического метода градиентного

спуска, как правило, будет приводить к постоянному увеличению Ω на каждой итерации опти-

мизации, что уменьшает практическую ценность найденного таким образом решения. В связи

с этим предлагается алгоритм подвижного окна, который улучшает данную ситуацию и может

рассматриваться как градиентный метод.

Решение задачи интервальной параметрической идентификации можно рассматривать как

поиск для каждой точки ŷ^kj прообраза в пространстве начальных условий системы (1) и на-

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

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

мы ОДУ будет проходить через точку ŷ^kj или минимизировать отклонение от неё. Необходимо

дополнительно отметить, что ŷ^kj обычно содержит информацию не обо всех фазовых перемен-

ных, поэтому нельзя проинтегрировать систему ОДУ по времени в обратную сторону, чтобы

найти прообразы.

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

зависимости решения системы ОДУ от начальных условий в некоторой области (в окне). В

рамках окна для каждой точки ŷ^kj определяется прообраз, а также выполняется оценка, ку-

да этот прообраз должен переместиться далее, чтобы ещё более уменьшить соответствующее

отклонение. Следующее окно строится исходя так, чтобы оно включало в себя найденные про-

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

когда окно перестаёт перемещаться.

Выполним формальное описание алгоритма. Пусть Ω⁽⁰⁾ - начальное окно, которое задаёт-

ся с помощью параметров zc(0), s⁽⁰⁾ и ϕ⁽⁰⁾; h - параметр расширения окна; ε - параметр

остановки; s^min = (s^min1, s^min2, . . . , sminm), s^max = (s^max1, s^max2, . . . , smaxm) - минимальный и макси-

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

1. С помощью алгоритма адаптивной интерполяции выполняется интегрирование систе-

мы (3) с параметрами zc(it), s(it), ϕ(it), (θ(it)) и строятся соответствующие вектор-функ-

ции P^k(ξ).

2. Для каждой точки ŷ^kj определяется прообраз в пространстве начальных условий z^kj =

= z(θ(it),ξ^kj), где ξ^kj = arg minJkj(z(θ(it),ξ)).

ξ∈[0,1]^m

3. Для нахождения направления перемещения прообраза z^kj вычисляются производные

dJ^kj(z(θ(it), ξ^kj))/dz_l, l = 1, m, с помощью решения СЛАУ

∑

dJ^kj(z(θ(it), ξ^kj)) dz_l(θ(it), ξ^kj)

∑

dP^ki(ξ^kj)

[P^ki(ξ^kj) - ŷ^kj,i]

p = 1,m.

dξ_p

l=1

i=1

4. Определяются точки zn^kj, в направлении которых могут переместиться прообразы:

{

∥Δz^kj∥-1, если

∥Δz^kj∥ > ε,

zn^kj = z^kj + hΔz^k

иначе,

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

821

где

(

dJ^kj(z(θ(it), ξ^kj))

dJ^kj(z(θ(it), ξ^kj)))т

Δz^kj =

,...,-

dz₁

dz₂

dz_m

5. Определяются параметры zc(it+1), s^next и ϕ(it+1) окна, которое содержит все точки z^kj

и zn^kj, j ∈ {0,1}ⁿ, k = 1,N.

6. Применяются ограничения на размер окна:

s(it+1) = min(s^max,max(s^min,s^next)),

где min, max - поэлементные операции. Получается новое окно Ω(it+1).

7. Алгоритм останавливается, когда Ω(it+1) = Ω(it).

Отметим несколько важных практических моментов в алгоритме.

1. Если для нахождения ξ^kj применяются итерационные методы, то в качестве стартового

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

шился на предыдущем шаге: ξ(it+1)start = z-1(θ(it+1), zk(it)j), где z-1(θ, z) - обратная функция

к z(θ,ξ).

2. Окна, полученные на двух подряд идущих итерациях, имеют не пустое пересечение:

Ω(it+1)

^⋂Ω(it) = ∅. Для уменьшения вычислительных затрат при построении Pk(it+1), вместо

интегрирования системы ОДУ (5) в точках, которые попадают в пересечение, можно исполь-

зовать Pk(it).

3. Модуль антиградиента ∥Δz^kj∥ отличен от нуля только для точек, находящихся на гра-

нице окна, при этом антиградиент всегда будет направлен от окна, перпендикулярно границе.

4. Параметр s^min отвечает за то, чтобы окно не становилось вырожденным, а параметр

s^max, в большей степени, требуется для контроля вычислительных ресурсов в рамках одного

шага алгоритма.

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

затратами. Поэтому алгоритм можно применять следующим образом: сначала ищется окно

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

эллипсоидной формы.

6. Если решение задачи arg minJkj(z(θ(it),ξ)) не единственное, то в качестве ξkj имеет

ξ∈[0,1]^m

смысл брать точку, находящуюся ближе к центру области: (0.5, 0.5, . . . , 0.5)^т, что позволит

дополнительно уменьшить размер окна. Однако такая ситуация приводит к некоторой неод-

нозначности при построении Ω, и в этом случае критерием остановки алгоритма может вы-

ступать ограничение на модули градиентов: ∥Δz^kj∥ ≤ ε, j ∈ {0, 1}ⁿ, k = 1, N .

Основной интерес с позиции размера получающегося окна Ω представляют результаты

п. 5.

5. Построение окна. Рассмотрим задачу построения эллипсоида или бруса, содержащего

заданное множество точек. Известны постановки такой задачи с условием минимальности объ-

ёма соответствующей фигуры: задача Сильвестра [24] (найти наименьший по объёму шар), за-

дача о минимальном эллипсоиде, задача о минимальном прямоугольном параллелепипеде [25]

и т.д. Алгоритмы их решения (см. [26, 27]) зачастую обладают высокой вычислительной слож-

ностью и носят итеративный характер.

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

ный на статистическом анализе исходного множества точек, который даёт приемлемые резуль-

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

Пусть исходное множество состоит из N точек z^k = (zk1, zk2, . . . , z^km)^т, k = 1, N . Построим

выборочную матрицу ковариации

∑

Σ=

(z^k - μ)(z^k - μ)^т,

N-1

k=1

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

№6

2023

822

МОРОЗОВ, РЕВИЗНИКОВ

∑_N

где μ = N-1

z^k - выборочное математическое ожидание.

k=1

Обозначим через U = (u₁, u₂, . . . , u_m) матрицу, составленную из нормированных соб-

ственных векторов матрицы ковариации Σ. Для определения углов ϕ выполняется последо-

вательное обнуление недиагональных элементов матрицы U-1 = U^т :

V =U-1,

for i = 1, m:

for j = (i + 1), . . . , m:

ϕi,j = arctg (vi,j/vi,i),

V =VRϕi,j,

где Rϕi,j - матрица поворота.

Масштабные коэффициенты s и центр области z^c вычисляются по следующим формулам:

(d^max - d^min), z^c = μ +

U (d^min + d^max),

где d^min = min[U-1(zk - μ)] и dmax = max[U-1(zk - μ)] - предельные отклонения центри-

рованных точек исходного множества от начала координат в базисе U; min и max - поэле-

ментные операции.

В случае с эллипсоидом масштабные коэффициенты дополнительно умножаются на кор-

ректирующее значение η, равное максимальному среди всех точек расстоянию от центра коор-

динат в базисе U с учётом уже найденных s и z^c :

η = max∥diag {s}-1U-1(z^k - z^c)∥,

где diag {s} - диагональная матрица, ∥ · ∥ - евклидова норма.

Данный подход можно использовать на промежуточных шагах алгоритма подвижного ок-

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

нимальности объёма. Это позволит уменьшить вычислительные затраты.

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

ный прообраз в пространстве начальных условий и параметров, то верно утверждение о том,

что решение, полученное с помощью алгоритма подвижного окна, будет наилучшим с точки

зрения объёма полученной области неопределённости. Кроме того, алгоритм даёт некоторые

вероятностные оценки неопределённых величин.

6. Результаты. Приведём результаты решения нескольких задач параметрической иден-

тификации динамических систем с помощью представленного алгоритма подвижного окна.

Вместо экспериментальных интервальных оценок фазовых переменных использовались ква-

зиэкспериментальные: задавались параметры области неопределённости z^c, s, ϕ и решалась

прямая задача; в каждый заданный момент времени из полученного параметрического мно-

жества решений Y^k брался случайный прямоугольный параллелепипе

Y k. Для всех вершин

Y k, k = 1,N, вычислялись прообразы в пространстве начальных условий и определялись эта-

лонные (исходные) значения параметров z^c, s и ϕ. Во всех примерах значение параметра

остановки ε = 10-10.

Рассмотрим систему ОДУ с двумя неизвестными параметрами, которая соответствует мо-

дели Лотки-Вольтерры:

u^′ = αu - 2uv, v^′ = -v + βuv, u(0) = 1, v(0) = 3,

(

)

(

)

(

)(

)

α₀

cos ϕ

- sin ϕ

s_αg₁(ξ₁,ξ₂)

ξ₁ ∈ [0,1], ξ₂ ∈ [0,1], t ∈ [0,2.5],

(16)

β₀

sin ϕ cos ϕ

s_βg₂(ξ₁,ξ₂)

где g₁, g₂ определяют преобразование единичного квадрата или в квадрат с длиной стороны 2

и центром в начале координат: g₁(ξ₁, ξ₂) = 2ξ₁ - 1, g₂(ξ₁, ξ₂) = 2ξ₂ - 1, или в единичный круг:

g₁(ξ₁,ξ₂) = ξ₁ cos(2πξ₂), g₂(ξ₁,ξ₂) = ξ₁ sin(2πξ₂). Параметры области неопределённости: α₀,

β₀, s_α, s_β, ϕ.

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

823

На рис. 1 показаны сгенерированные квази-

экспериментальные интервальные оценки фазо-

вых переменных

Y k, k = 1,20, в равномерно

расположенных точках на интервале интегриро-

вания: t_k = k/8. Исходные значения параметров

в случае с прямоугольным параллелепипедом:

α₀ = 2.00, β₀ = 0.98, s_α = 0.21, s_β = 0.39,

ϕ = -0.67, в случае с эллипсоидом: α₀ = 2.00,

β₀ = 0.98, s_α = 0.28, s_β = 0.52, ϕ = -0.67.

Заметим, что площадь соответствующего прямо-

угольника (S = 4s_αs_β = 0.33) меньше, чем пло-

щадь эллипса (S = πs_αs_β = 0.46), это связано с

тем, что при генерации изначально использова-

лась прямоугольная область Ω.

Выполняется решение задачи параметриче-

ской идентификации. Шаг расширения окна

h = 0.2; начальные значения параметров: α⁽⁰⁾⁰ =

(0)

Рис. 1. Квазиэкспериментальные интервальные

= 2.50, β₀

= 0.50, sα0) = 0.05, s(0)β = 0.05,

оценки фазовых переменных системы (16) в раз-

личные моменты времени.

ϕ(0) = 0.00; ограничения на размер окна s^min и

s^max не задаются. На рис. 2 и 3 проиллюстри-

рован процесс работы алгоритма в случаях пря-

моугольной и эллипсоидной областей Ω.

Рис. 2. Иллюстрация решения задачи параметрической идентификации системы (16) в случае прямоугольной

области Ω.

В верхних четырёх строках на рис. 2 и 3 показаны квазиэкспериментальные оценки фа-

зовых переменных

Y k (темно-серый цвет) для k = 1,7,14,20 и соответствующие модельные

множества Y^k (светло-серый цвет) на всех итерациях алгоритма подвижного окна. Серые ли-

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

№6

2023

824

МОРОЗОВ, РЕВИЗНИКОВ

нии соединяют каждую вершин

Y k с ближайшей точкой в Y ^k. Сумма квадратов длин этих

линий по всем k фактически определяет значение целевой функции (14).

В нижней строке на рис. 2 и 3 демонстрируется движение окна Ω(it) (закрашенная фигура)

в пространстве параметров α × β. Исходное множество показано незакрашенной фигурой

(точное решение). Серые линии соединяют точки z^kj и zn^kj. Каждое окно покрывает линии,

полученные на предыдущем шаге. Со 2-й итерации линии начинают вырождаться, это говорит

о том, что для соответствующей экспериментальной точки найдены значения α и β, при

которых отклонение от модельного решения будет минимальным (в данном примере равно

нулю). На 4-й итерации все z^kj стали равны zn^kj, и после 5-й итерации окно более не будет

перемещаться.

Рис. 3. Иллюстрация решения задачи параметрической идентификации системы (16) в случае эллипсоидной

области Ω.

Найденные значения α₀, β₀, s_α, s_β, ϕ на последней итерации алгоритма совпадают с ис-

ходными значениями. Модельные множества полностью содержат в себе интервальные оценки

фазовых переменных. Отметим, что полученные множества Y^k в случае прямоугольной фор-

мы окна не являются выпуклыми (левая сторона множеств имеет вогнутость), таким образом,

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

Далее рассмотрим задачу с трёхмерной областью неопределённости:

u^′ = αu -

uv +

u², v^′ = -2v +

uv +

v²,

⎛

⎞

⎛

⎞

⎛

⎞

u(0)

u₀

s_ug₁(ξ₁,ξ₂,ξ₃)

⎝v(0)⎠= ⎝v0⎠+Rϕ

Rϕ1,3 Rϕ2,3

⎝svg2(ξ1,ξ2,ξ3)⎠,t∈[0,4],

(17)

1,2

α₀

s_αg₃(ξ₁,ξ₂,ξ₃)

где ξ₁ ∈ [0, 1], ξ₂ ∈ [0, 1], ξ₃ ∈ [0, 1]; Rϕ1,2 , Rϕ1,3 , Rϕ2,3 - 3 × 3-матрицы поворота; g1,

g₂, g3 задают преобразование единичного куба в куб с длиной стороны 2 и центром в начале

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

825

координат или в шар единичного радиуса. Параметры области неопределённости: u₀, v₀, α₀,

α₀, s_v, s_α, ϕ1,2, ϕ1,3, ϕ2,3.

На рис. 4 тёмно-серым цветом показаны сгенерированные квазиэкспериментальные интер-

вальные оценки фазовых переменных

Y k, k = 1,20, в равномерно расположенных точках

tk = k/5 на интервале интегрирования. При генерации использовалась прямоугольная об-

ласть неопределённости со следующими значениями параметров: u₀ = 4.5, v₀ = 3.0, α₀ = 4.0,

s_u = 0.1, s_v = 0.1, s_α = 0.1, ϕ1,2 = -π/6, ϕ1,3 = π/3, ϕ2,3 = π/4.

Рис. 4. Квазиэкспериментальные интервальные оценки фазовых пе-

ременных системы (17) в различные моменты времени (темно-серый

цвет) и полученные множества Y^k на последней итерации алгоритма

в случае области неопределённости в форме прямоугольного паралле-

лепипеда (светло-серый цвет).

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

движного окна. Шаг расширения h = 0.1; начальные значения параметров: u⁽⁰⁾⁰ = 3.5, v⁽⁰⁾⁰ =

= 4.0, α⁽⁰⁾⁰ = 3.5, su0) = 0.1, sv0) = 0.1, sα0) = 0.1, ϕ(0)1,2 = 0, ϕ(0)1,3 = 0, ϕ(0)2,3 = 0; ограничения

на размер окна: smaxu = 0.2, smaxv = 0.2, smaxα = 0.2. Алгоритм завершил работу за 20 ите-

раций. На рис. 5 показаны получающиеся окна в процессе работы алгоритма для каждой 4-й

итерации.

Рис. 5. Иллюстрация перемещения окна в процессе работы алгоритма

при решении задачи параметрической идентификации системы (17).

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

№6

2023

826

МОРОЗОВ, РЕВИЗНИКОВ

Получены значения для прямоугольной области неопределённости: u⁽²⁰⁾⁰ = 4.55, v⁽²⁰⁾⁰ =

= 2.98, α⁽²⁰⁾⁰ = 3.94, su20) = 0.11, sv20) = 0.11, sα20) = 0.07, ϕ(20)1,2 = 0.30, ϕ(20)1,3 = 0.29,

ϕ(20)2,3 = -0.37. Они не совпадают со значениями, использующимися при генерации.

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

точке на фазовой плоскости соответствует множество точек в пространстве u(0) × v(0) × α.

В связи с этим решение задачи параметрической идентификации является не единственным.

На рис. 4 показана фазовая плоскость на последней итерации алгоритма: все модельные мно-

жества Y^k (светло-серый цвет) полностью содержат в себе квазиэкспериментальные интер-

вальные оценки

Y k (темно-серый цвет). Множества Y ^k являются проекцией трёхмерного

деформирующегося прямоугольного параллелепипеда. Каждая точка на такой проекции бу-

дет соответствовать кривой в пространстве u(0) × v(0) × α.

На последней итерации алгоритма была построена также область неопределённости в фор-

ме эллипсоида: u⁽²⁰⁾⁰ = 4.55, v⁽²⁰⁾⁰ = 2.98, α⁽²⁰⁾⁰ = 3.94, su20) = 0.25, sv20) = 0.15, sα20) = 0.1,

ϕ(20)1,2 = 0.30, ϕ(20)1,3 = 0.29, ϕ(20)2,3 = -0.37. В этом случае все множества Y^k будут являться

выпуклыми, в отличие от множеств на рис. 4.

Представленный в работе алгоритм успешно справился с данными задачами: во всех слу-

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

шение системы ОДУ полностью содержит интервальные оценки фазовых переменных.

Заключение. В работе рассмотрена задача параметрической идентификации динамиче-

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

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

чи сводится к минимизации целевой функции, которая характеризует суммарное отклонение

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

Сформулировано и доказано утверждение о том, что если модельные множества являются вы-

пуклыми, то вместо экспериментальных интервальных оценок можно рассматривать только

крайние точки этих оценок. Для эффективной минимизации целевой функции разработан ал-

горитм подвижного окна, относящийся к градиентным методам. В его основе лежит алгоритм

адаптивной интерполяции, который позволяет в рамках определённого окна в пространст-

ве неизвестных параметров для каждого момента времени построить полином, интерполи-

рующий зависимость состояния динамической системы от параметров. Выполнена апробация

предложенного алгоритма на нескольких задачах, содержащих разное количество неизвестных

величин. Для всех задач успешно получены области значений неизвестных параметров, при

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Шенк Х. Теория инженерного эксперимента. М., 1972.

2. Martyshov M.N., Emelyanov A.V., Demin V.A. et al. Multifilamentary character of anticorrelated

capacitive and resistive switching in memristive structures based on (Co-Fe-B)x(LiNbO3)100-x nano-

composite // Phys. Rev. Appl. 2020. V. 14. № 3. P. 034016.

3. Moore R. Interval Analysis. Englewood Cliffs, 1966.

4. Moore R.E., Kearfott R.B., Cloud M.J. Introduction to Interval Analysis. Philadelphia, 2009.

5. Шарый С.П. Конечномерный интервальный анализ. Новосибирск, 2019.

6. Добронец Б.С. Интервальная математика. Красноярск, 2007.

7. Xiao N., Fedele F., Muhanna R.L. Inverse problems under uncertainties-an interval solution for the beam

finite element // 11th Intern. Conf. on Structural Safety & Reliability. New York, 2013. P 1-8.

8. Петрикевич Я.И. Структурно-параметрическая идентификация динамических объектов по интер-

вальным исходным данным: дис

канд. техн. наук. М., 2006.

9. Дилигенская А.Н., Самокиш А.В. Параметрическая идентификация в обратных задачах теплопро-

водности в условиях интервальной неопределённости на основе нейронных сетей // Вестн. Самар-

ского гос. техн. ун-та. 2020. Т. 28. № 4 (68). С. 6-18.

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

№6

2023

АЛГОРИТМ ПОДВИЖНОГО ОКНА

827

10. Морозов А.Ю., Ревизников Д.Л. Интервальный подход к решению задач параметрической иденти-

фикации динамических систем // Дифференц. уравнения. 2022. Т. 58. № 7. С. 962-976.

11. Морозов А.Ю., Ревизников Д.Л. Алгоритм адаптивной интерполяции на основе kd-дерева для чис-

ленного интегрирования систем обыкновенных дифференциальных уравнений с интервальными

начальными условиями // Дифференц. уравнения. 2018. Т. 54. № 7. С. 963-974.

12. Морозов А.Ю., Ревизников Д.Л., Гидаспов В.Ю. Алгоритм адаптивной интерполяции на основе

kd-дерева для решения задач химической кинетики с интервальными параметрами // Мат. моде-

лирование. 2018. Т. 30. № 12. С. 129-144.

13. Морозов А.Ю. Интерполяционный подход в задачах моделирования динамических систем с эллип-

соидными оценками параметров // Тр. МАИ. 2022. № 124. С. 1-24.

14. Смоляк С.А. Квадратурные и интерполяционные формулы на тензорных произведениях некоторых

классов функций // Докл. АН СССР. 1963. Т. 148. № 5. С. 1042-1045.

15. Bungatrz H-J., Griebel M. Sparse grids // Acta Numerica. 2004. V. 13. № 1. P. 147-269.

16. Gerstner T., Griebel M. Sparse grids // Encyclopedia of Quantitative Finance / Ed. R. Cont. New York,

2010.

17. Morozov A.Yu., Zhuravlev A.A., Reviznikov D.L. Sparse grid adaptive interpolation in problems of

modeling dynamic systems with interval parameters // Mathematics. 2021. V. 9. P. 298.

18. Морозов А.Ю., Ревизников Д.Л. Алгоритм адаптивной интерполяции на разреженных сетках для

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

неопределённостями // Дифференц. уравнения. 2021. Т. 57. № 7. С. 976-987.

19. Морозов А.Ю. Параллельный алгоритм адаптивной интерполяции на основе разреженных сеток для

моделирования динамических систем с интервальными параметрами // Программная инженерия.

2021. Т. 12. № 8. С. 395-403.

20. Демьянов В.Ф., Малоземов В.Н. Введение в минимакс. М., 1972.

21. Евтушенко Ю.Г. Некоторые локальные свойства минимаксных задач // Журн. вычислит. матема-

тики и мат. физики. 1974. Т. 14. № 3. С. 669-679.

22. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. М., 1985.

23. Пантелеев А.В., Летова Т.А. Методы оптимизации в примерах и задачах. М., 2005.

24. Sylvester J. J. A question in the geometry of situation // Quarterly J. of Math. 1857. V. 1. P. 79.

25. Васильев Н.С. О численном решении экстремальных задач построения эллипсоидов и параллеле-

пипедов // Журн. вычислит. математики и мат. физики. 1987. Т. 27. № 3. С. 340-348.

26. Шор Н.З., Стеценко С.И. Алгоритм последовательного сжатия пространства для построения опи-

санного эллипсоида минимального объёма // Исследование методов решения экстремальных задач.

Киев, 1990. С. 25-29.

27. Khachiyan L.G. Rounding of polytopes in the real number model of computation // Math. of Operations

Research. 1996. V. 21. № 2. P. 307-320.

Федеральный исследовательский центр

Поступила в редакцию 20.02.2023 г.

“Информатика и управление” РАН, г. Москва,

После доработки 20.02.2023 г.

Московский авиационный институт

Принята к публикации 18.04.2023 г.

(национальный исследовательский университет)

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

№6

2023