CC BY
14ISSN 2782-2001 Системы анализа и обработки данных том 85, № 1, 2022, с. 25-40
http://journals.nstu.ru/vestnik Analysis and data processing systems Vol. 85, No. 1, 2022, pp. 25-40
ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ
INFORMATICS, COMPPUTER ENGINEERING AND MANAGEMENT
УДК 681.51:519.6 Б01: 10.17212/2782-2001-2022-1-25-40
Идентификация квадратичного ядра уравнения Вольтерра для моделирования нелинейных динамических систем*
Ю.Е. ВОСКОБОЙНИКОВ1'2'", В.А. БОЕВА1-4
1 630008, РФ, г. Новосибирск, ул. Ленинградская, 113, Новосибирский государственный архитектурно-строительный университет
2 630073, РФ, г. Новосибирск, пр. Карла Маркса, 20, Новосибирский государственный технический университет
а voscob@mail.ru ь v.boyeva@sibstrin.ru
В последние два десятилетия для описания динамики стационарных нелинейных систем используются интегральные модели «вход-выход», в которых ядрами являются члены ряда Вольтерра. Наиболее часто используются линейный член (импульсная переходная функция зависит от одной переменной) и квадратичный член, зависящий от двух переменных. Для выделения в выходном сигнале идентифицируемой системы двух его составляющих - выхода линейной «подсистемы» и выхода «квадратичной» подсистемы - проводят активный эксперимент, в котором на вход системы подается специальная комбинация прямоугольных импульсов. После выделения выхода «квадратичной» подсистемы идентификация квадратичного члена ряда Вольтерра сводится к решению двумерного интегрального уравнения первого рода. В литературе приводятся формулы обращения, в которых функция квадратичного ядра получается в результате арифметических операций с производными второго порядка от выходного сигнала. Дифференцирование функций является некорректно поставленной задачей, когда малые погрешности задания функции (шумы измерения) вызывают большие ошибки в производных (особенно в производных второго порядка). В работе предлагается для устойчивого вычисления производных использовать сглаживающие кубические сплайны. Для вычисления смешанной производной второго порядка строится сплайн с двумя переменными - сглаживающий бикубический сплайн. Основной проблемой, возникающей на практике при обработке данных реального эксперимента, является выбор параметра сглаживания, от величины которого зависит ошибка сглаживания зашумленных данных. Как правило, в эксперименте не известна величина дисперсии шума измерения. Поэтому в работе предлагается для выбора параметра сглаживания в построенных сплайнах (особенно в бикубическом) использовать алгоритм, основанный на методе ¿-кривой, где не требуется задание дисперсии шума измерения. Предлагаемый алгоритм идентификации имеет высокую вычислительную эффективность. Выполненный вычислительный
*
Статья получена 14 июля 2021 г. Исследование выполнено при финансовой поддержке РНФ в рамках научного проекта №22-21-00409.
эксперимент показал маленькую методическую ошибку (порядка 1 %) и хорошую устойчивость к шумам измерений выходных сигналов идентифицируемой системы. Для уменьшения случайной составляющей ошибки идентификации предлагается использовать постобработку локально-пространственным комбинированным фильтром.
Ключевые слова: интегральные модели нелинейных систем, ядра Вольтерра, формулы обращения, дифференцирование зашумленных данных, сглаживающий кубический сплайн, сглаживающий бикубический сплайн, выбор параметра сглаживания, локально-пространственный комбинированный фильтр
Интегральные представления нелинейных динамических систем в терминах «вход-выход» нашли широкое применение при их моделировании (в частности, интегральные модели Гаммерштейна, Винера - Гаммерштейна [1]). Для идентификации динамических характеристик нелинейной системы в условиях априорной неопределенности о внутренней структуре в последние два десятилетия активно применяется математический аппарат интегростепенных рядов Вольтерра (подробные обзоры см. в [2, 3]). Особенно часто используется ряд, ограниченный квадратичным членом ряда Вольтерра. В этом случае связь между входным сигналом ф(т) моделируемой стационарной системы и ее выходным сигналом /(?) можно представить следующим интегральным соотношением [4]:
f(t) = J -x)dx + JJ *2(т, " s)dxds, t e[0,T ] , (1)
где ki(x), &2(т,s) - линейное и квадратичное ядра Вольтерра соответственно. Нужно построить метод и алгоритм идентификации этих ядер по измеренным значениям входного и выходного сигналов идентифицируемой системы.
Очевидно, что сигнал f (t) есть сумма двух неизвестных сигналов: fi(t) и fi(t) - первое и второе слагаемые в (1), которые можно интерпретировать как выходные сигналы линейной и квадратичной «подсистем». Следовательно, на первом этапе идентификации необходимо из f (t) выделить сигналы fi (t), f 2 (t) по отдельности. Для произвольного входного сигнала такая задача неразрешима. Поэтому обратимся к простой, но достаточно эффективной методике, предложенной в работах [5, 6], которая подразумевает проведение активной идентификации, когда на вход системы подаются сигналы заданной формы и заданной амплитуды.
Обозначим через фн (t) входной сигнал, являющийся функцией Хэви-сайда:
1. ВВЕДЕНИЕ И ПОСТАНОВКА ЗАДАЧИ
t
11
0
00
Определим следующий набор входных сигналов:
фа (t) = афн (t), Ф а (t) = -афн (t), 0 < t < T ,
(2)
Фа V) = а(Фя (0-Фя ^ - V)), 0 < V < t < Т ;
(3)
ф-а (t,V) = -а(фн ^) - Фя ^ - V)), 0 < V < t < Т.
При подаче на вход системы сигналов (2) имеем следующие выходные сигналы:
г г г
/^) = а{х + а2 Ц£2(х, s)dхds, 0 < t < Т;
0 00 t 11
(4)
f (-a)(t) = -a jkx(x)dx + a2 j jk2(x,s)dxds, 0 < t < T.
0 00
При подаче входных сигналов (3) имеем
t t t f (a)(t, v) = a j kj(x)dx + a2 j j k2(x, s)dxds, 0 <v <t <T;
t-v t-v t-v
t t t
(5)
f( a)(t, v) = -a j kj(x)dx + a2 j j k2(x, s)dxds, 0 <v <t <T.
t-v t-v t-v
Осуществляя операцию вычитания в (4) и операцию сложения в (5), получаем следующие полезные соотношения:
fk,(x)dx = М) = f(a)(t)-f 1 -al(t);
0 2a
j j k2(x, s)dxds = f2(t, v) = f (t V) + f~ (t, v), 0 < v < t < T.
t-vt-v 2a
(6)
(7)
Из выражения (6) следует формула обращения
) = /№), t е [0,Т], (8)
а из (7) - формула обращения [4, 6]
/2» ^ V) + /^2 ' V) ^ - V, t) =-2-' 0 < V < t < Т, (9)
где используются следующие обозначения:
2 2
/1(0 = "/), /2*' а, V) V), /2v2 " (t, V) = ^г />(/, V). (10)
dt дtдv 2v dV2
К сожалению, реализация на практике полученных формул обращения сталкивается с принципиальной трудностью - некорректностью операции
дифференцирования [7]. Одно из проявлений некорректности заключается в больших ошибках вычисления производной даже при очень малых погрешностях задания дифференцируемой функции. Отметим, что операция сложения в (7) погрешностей регистрации двух функций приводит к увеличению дисперсии суммарной погрешности задания функции /V). Таким образом, устойчивое дифференцирование зашумленных данных становится актуальной задачей для реализации формулы (9) на практике.
В работе [8] был построен устойчивый алгоритм идентификации на основе формулы (8), где для устойчивого вычисления первой производной использовался сглаживающий кубический сплайн (СКС) с дефектом единица с выбором параметра сглаживания из условия минимума среднеквадратической ошибки сглаживания. В случае идентификации квадратичного ядра ^(т, использование сглаживающих сплайнов существенно усложняется. Во-первых, для вычисления смешанной производной второго порядка /2^" V) нужно строить уже сглаживающий бикубический сплайн (СБС), являющийся функцией двух переменных I, V . Во-вторых, краевые условия задаются уже не в двух крайних точках интервала построения СКС, а на четырех прямых, являющихся границами прямоугольной области построения СБС. В-третьих, из-за разной «гладкости» функции /V) по разным переменным необходимо выбрать уже два параметра сглаживания из условия минимума ошибки сглаживания.
Эти затруднения обусловили следующие основные задачи, которые не получили своего решения в соответствующих научных публикациях и которые решаются в настоящей работе:
• разработка алгоритма, позволяющего строить сглаживающий бикубический сплайн с большим числом комбинаций краевых условий на разных границах области построения сплайна;
• раздельный выбор двух параметров сглаживания (по каждой переменной сплайна) из условия минимума ошибки сглаживания на основе метода Х-кривой;
• исследование устойчивости полученного решения к случайным ошибкам измерения функции /2 (^, V);
• применение алгоритмов постобработки результатов идентификации для повышения ее точности.
2. АЛГОРИТМ ИДЕНТИФИКАЦИИ КВАДРАТИЧНОГО ЯДРА
Обратимся к формуле (7) и предположим, что значения двумерных функций / (а)(^ V), / (_а)(^, V) измеряются при значениях аргументов , / = 1...; Vj, у = 1...NV , т. е. в узлах прямоугольной сетки
, Vj |, } = 1. Nt, у = 1. NV . Заметим, что узлы и ту могут иметь неодинаковый и неравный шаг. Тогда значения функции /2 V) также определяются в узлах этой прямоугольной сетки, и для учета возможных погрешностей
(шумов) измерений принимается следующая модель для зашумленных значений , Vj ):
, Vj) = /2 , Vj) + Ц|,у, I = 1. Nt, у = 1. N, (11)
где ц. у - случайный шум измерения с нулевым средним и дисперсией стЦ (равноточные измерения). Требуется по исходным данным | /2(^, V у )| вычислить значения производных /^,V), / 2" ^,V) в заданных узлах.
Для устойчивого вычисления этих производных обратимся к сглаживающим кубическим сплайнам [9], широко используемым при обработке экспериментальных данных [10, 11]. Предположим, что на некотором интервале [^1, Р2] задано ^ узлов V = Vl < V2 <... < VNv = и в этих узлах измерены значения функции (сигнала) /(V):
/у = /V ) + Ц у, у = 1. Nv , (12)
где цу - случайный шум измерений с нулевым средним и дисперсией Стц (равноточные измерения). Сглаживающий кубический сплайн SNV а (V) с дефектом единица на каждом отрезке , Vj+1) можно представить кубическим многочленом вида [10]
SNV,а(V) = ау + Ьу (V - Vj) + Су (V - Vj )2 + ёу (V - Vj )3 , (13)
причем функция SNV а (V) должна быть дважды непрерывно дифференцируема на всем интервале [V, ^2 ]. Заметим, что в отличие от интерполяционного сплайна (проходящего через точки (уу, /у)) сглаживающий кубический сплайн SNV а (V) в общем случае не проходит через эти точки, а проходит более «плавно» в некоторой окрестности этих точек, зависящей от величины параметра сглаживания а , обеспечивая тем самым сглаживание (фильтрацию) шума измерений.
Для однозначного вычисления коэффициентов сплайна а у, Ьу, Су, ёу задают краевые условия в узлах V, VNv . Наиболее часто используются следующие условия [9, 11]:
• условия на нулевые вторые производные сплайна (естественные краевые условия):
S;vV,аЫ = 0; S'N,vа^) = 0; (14)
• условия на первые производные сплайна:
SNV ,а ) = ; SNV ,а (Уы, ) = , (15)
а также комбинация этих условий (например, слева условие (15), справа - (14)). Показано [9], что СКС, построенный при этих условиях, доставляет минимум функционалу:
^ (5) = а { |5»|2 ^ +Х Р-1(//-S(vJ■))2, (16)
VI М
где pJ - весовые множители, отражающие точность J-го измерения fj (в случае равноточных измерений задаются одинаковыми).
Для вычисления коэффициентов сплайна (при заданном параметре сглаживания) составляется система линейных алгебраических уравнений с пяти-диагональной матрицей относительно некоторого вектора (как правило, это значения второй производной сплайна в узлах {vJ}), через который затем
находятся все коэффициента сплайна (подробнее см. [9, 11]).
Параметр сглаживания а «управляет» гладкостью сплайна, и ошибка сглаживания (как и ошибка дифференцирования) существенно зависит от величины этого параметра [11, 12]. Существует значение параметра (назовем его оптимальным), для которого ошибка сглаживания (в принятой норме) минимальна [12]. Временно предположим, что приемлемое (с точки зрения минимума ошибки сглаживания) значение параметра сглаживания может быть найдено (выбор параметра рассматривается в следующем разделе).
Замечание 1. Из вида интеграла (7) следует, что функция /2(^ V) принимает ненулевые значения для аргументов, удовлетворяющих условию V < t. Для отрицательных значений V, t функция равна нулю в силу условия технической реализуемости системы, т. е.
V) = 0, если V < 0, / < 0.
Для уменьшения влияния разрыва первого рода при значениях V = t предлагается дополнить значения функции /V) для V > t по правилу
/2 ^, t + Лv) = /2 t) + (/2 (^ ^ -/2 ^ t - Л^, Лv > 0 .
*
Дополненную таким образом функцию будем обозначать как /2 ^, V) .
Первоначально остановимся на алгоритме вычисления значений производной /2 2 V) , который можно представить следующими шагами.
Шаг 1. Задаются краевые условия, комбинация которых в крайних точках
V!, VNv интервала построения определяется исходя из имеющейся априорной
*
информации о функции /2 ^, V) (при отсутствии такой достоверной информации следует обратиться к естественным условиям (14)).
Шаг 2. Для каждого г = 1... формируется набор данных
V, Н? = /2*(г,VJ), У = 1..N) ,
выбирается параметр сглаживания а1(|) и строится СКС S1(|) ( ) (V), по ко-
NV ,а1(|)
торому затем вычисляется первая производная /2У (^, Vj) = = ^1((|) 1(|) (V)!= Ь1(р (оценка производной /2У (^,V,-)), где Ь1(;г) - ко-
—V NV ,а1; V=Vj ■>
эффициент сплайна S1(|) (■)(v) в представлении (13).
Nv ,а1(|)
Шаг 3. Для каждого i = 1. Nt формируется набор данных
V, = ^(ti,Vj), у = 1,...,N),
выбирается параметр сглаживания а2(|) и строится СКС S2(|) ..)(v), первая
Nv ,а 2(|)
производная которого является оценкой /^ (^,V ) =—S2(|) (|) (V)! = =
у —V Nv ,а2( ) =^7
= Ь2(|) для второй производной , Vj), где Ь2(|) - коэффициент сплайна
S 2(|) (|) (V) в представлении (13).
Nv ,а 2(|)
Таким образом, вычислены оценки , Vj) второй производной
,) в узлах vj для ti, | = 1,...,Nt.
Перейдем к построению бикубического сглаживающего сплайна (следуя методике работы [14]) для вычисления смешанной производной (ti,Vj)
в виде следующих шагов.
Шаг 7. Для каждого у = 1. N формируется набор данных (фиксируется значение Vj):
{ti, А0-)=/^, Vj), |=1. N.},
выбирается параметр сглаживания а3( у) и по исходным данным (фиксируется
значение V;) строится СКС S3(( Л (t), по которому затем вычисляется
7 Nt ,а3( у)
первая производная /2Г(^,v,■) = —S3(у (.) (0| = Ь3((оценка производ-
— Nt ,а3(-') ''='г
ной /2Г(ti,v,■)), где Ь3(у) - коэффициент сплайна S3(у ( в представле-
7 Nt ,а3( ■/)
нии (13).
Шаг 2. Для каждого / = 1. N формируется набор данных (фиксируется значение ti):
V, / 4/) = /2г (ti V), ] = 1.. N
выбирается параметр сглаживания а4(г) и строится СКС 54(г) (.)(v), первая
Nv ,а 4(г)
производная которого является оценкой (¡г,V ) = —54(г) ( ) (V)! = =
—V Nv ,а4(г) lV=VJ
= М^Р для смешанной производной /2'^,Vj), где Ь4(.) - коэффициент
сплайна 5 4(г) () (V) в представлении (13).
Nv ,а 4(г)
Шаг 1 повторяется для Vу, у = 1... Nv, а шаг 2 - для , г = 1... Nt. После вычисления оценок /^ , VJ■), /2'^ (¡у , VJ■) по формуле (9) находится оценка ^(¡г - Vу, ) для значений Vу < .
Замечание 2. Формула обращения (9) определяет значение квадратичного ядра к2^,V) для аргументов 0 < V < t < Т, т. е. для значений аргумента V< t. Прямая V = t является осью симметрии ядра £2^,V) (следствие одномерного входного сигнала), и поэтому для определения значений ядра для V = t + Лv > t, где Лv > 0, предлагается симметричное дополнение значений ядра по формуле £2(¡, t + Лv) = £2^ + Дv,t).
Замечание 3. Так как построение СКС по переменной V требует примерно С0рег • Nv , арифметических операций [8, с. 345], где С0рег « 30 то предлагаемый алгоритм вычисления производных потребует порядка с0рег^N1
операций. Следовательно, предложенные алгоритмы вычисления производных имеют высокую вычислительную эффективность даже при большой размерности сетки , Vу).
3. ВЫБОР ПАРАМЕТРА СГЛАЖИВАНИЯ
Ранее было отмечено, что величина параметра сглаживания а (см. (16)) существенно влияет на величину ошибки сглаживания СКС. Если бы дисперсия ст^ шума измерений (см. (12)) была достоверно известна (хотя бы с точностью 5.. .8 %), то алгоритм выбора, построенный на основе проверки критерия оптимальности линейного алгоритма фильтрации, позволял бы с приемлемой точностью 5.8 % оценить значения оптимального параметра сглаживания, минимизирующего величину среднеквадратической ошибки сглаживания (см. [11, с. 60-67; 12]). Очевидно, что ситуация, когда дисперсия шума неизвестна, наиболее характерна при решении практических задач идентификации. Поэтому для выбора параметра в этом случае обратимся к методу ¿-кривой, который рассматривается в зарубежных публикациях (например, в [14, 15]) для выбора параметра регуляризации в алгоритмах решения линейных некорректных задач. В работе [16] была сделана следующая модификация метода ¿-кривой для выбора параметра сглаживания. Кратко о сути этого алгоритма выбора.
Введем в рассмотрение следующие функционалы (см. (16)):
Щ . . VNv
2
Р (а) =1 Р-1(/у - Sи,а ^у))2, у (а) = { (V) —V . (17)
М v1
Тогда ¿-кривой (форма которой напоминает начертание латинской буквы Ь) называется параметрическая кривая с координатами (р (а), у (а)). Можно показать, что кривизна Ь-кривой определяется следующей формулой:
кь(а) = 2 р'(а) у''(а)-рЧа) у'(а)3 , а8)
(р '(а) )2 + (у '(а) )2
2
где рр (а) = 1п р (а), у (а) = 1п у (а). Для эффективного вычисления значения функционала у(а) предлагается следующая формула:
У(а) = 'I:1 (4с2 к + 12С| к2 + 12—2к3 ), i=1
где к| = ti+l -ti, i = 1...п -1, с|, - коэффициенты СКС в представлении (13), вычисленные при заданном параметре а . В качестве параметра сглаживания будем принимать величину аь , для которой кривизна кь (а) принимает максимальное значение.
Для вычисления значения аь можно использовать численные методы одномерной минимизации. При этом для нахождения кривизны по формуле (18) нужны значения производных функционалов р (а), у (а) в заданной точке а . Производные можно вычислить либо с использованием интерполяционных сплайнов р(а), у (а) (подробнее см. в [16]), либо через конечные разности соответствующего порядка.
В работе [16] выполнен обширный вычислительный эксперимент для ответа на вопрос: велик ли проигрыш по ошибке сглаживания при использовании аь вместо оптимального а0рГ (который можно определить только в вычислительном эксперименте)? Эксперимент проводился с функциями, являющимися «типовыми» выходными сигналами динамической системы при подаче на вход ступенчатых сигналов. Анализ результатов эксперимента показал, что алгоритм выбора параметра сглаживания на основе метода Ь-кривой позволяет достаточно хорошо оценить оптимальное значение параметра сглаживания. Увеличение ошибки сглаживания при использовании параметра а ь в среднем не превышает 5...15 % по сравнению с а0рГ, вычисление которого на практике невозможно.
Таким образом, для вычисления параметров сглаживания а1(|), а2(|),
а3(|), а4(|) предлагается использовать алгоритм выбора параметра сглаживания на основе метода Ь-кривой.
4. ЧИСЛЕННЫЕ ИССЛЕДОВАНИЯ АЛГОРИТМА ИДЕНТИФИКАЦИИ
В качестве тестового квадратичного ядра (т, я) была взята функция, используемая для описания динамики одного типа теплообменников [6]. Поверхность этой функции показана на рис. 1, а, а изолинии - на рис. 1, б. Граница временного интервала Т = 1, количество узлов N = 80, = 80.
Рис. 1. Тестовое квадратичное ядро Fig. 1. The test quadratic kernel
Первоначально определим методическую ошибку алгоритма идентификации. Для этого вычислялись значения интеграла (7) в узлах ^, i = 1... Nt, Vj, j = 1...Nv, которые интерпретировались как точные значения функции fi(tj,Vj). Эти данные, представленные в виде матрицы F размером 80 х 80 с элементами Fi j = f (ti,vj), являлись исходными для предлагаемого
алгоритма идентификации. Так как эти исходные данные принимались как точные, то вместо СКС строились интерполяционные кубические сплайны (включая бикубический) с краевыми условиями (14). По этим сплайнам вычислялись оценки производных f2'v (ti, Vj), (ti, Vj), а затем по формуле (9)
строилась оценка для квадратичного ядра (см. замечание 2). На рис. 2 показаны изолинии этой оценки, имеющей относительную ошибку идентификации
K2 " K2 -
- = 0,011, где K2, K2 - матрицы, составленные из значений
5 К =
11К 2||
точного ядра к2^г-, V у) и его оценок к2^г-, Уу ) соответственно, Ц - евклидова
норма матрицы. Примерно такая же ошибка наблюдалась при других размерах сеток по I, V. Поэтому можно сделать вывод, что предложенный алгоритм идентификации имеет низкую методическую ошибку.
б
а
Рассмотрим влияние шумов измерения функции f 2_(t, v) на точность идентификации. Для этого все элементы «точной» матрицы F искажались нор-
||f - Fil
мально распределенным шумом с относительным уровнем 5р = 11 11,
F
где F - матрица с «зашумленными» элементами. Сформированная таким образом матрица F использовалась в качестве исходных данных для ранее описанного алгоритма идентификации. Задавалось три уровня шума: 0.02, 0.04, 0.06. Параметр сглаживания на всех шагах вычисления производных выбирался описанным выше методом L-кривой. Основная составляющая ошибки идентификации имела случайный характер, и поэтому для ее уменьшения предлагается ввести этап постобработки, который заключается в фильтрации матрицы
K2 (результат работы алгоритма идентификации) локально-пространственным комбинированным двумерным фильтром [17], который эффективно удаляет артефакты вейвлет-фильтрации [18]. Результатом работы этого фильтра
Рис. 2. Оценка ядра к2(т, s), построенная по точным данным
Fig. 2. The Kernel estimation k2(%,s) on accurate data
^кф
Рис. 3. Оценка ядра k2(t, s) после постобработки
Fig. 3. The Kernel estimation k2(x, s) after post-processing
является матрица 1С ф , элементы которой определялись следующим соотношением:
K,2 = aver (K А :
VI, j 2
i\,i 2eA\
¿> мф {> мф
киЪ — Kj\, j 2
<д
F
(19)
т. е. усреднялись только те значения медианного фильтра ^гмф2, которые попадали в прямоугольную апертуру А фильтра с размерами К1 х К2 и
удовлетворяли еще дополнительному условию
Kмф — кмф Ki1,i 2 Kj1, j 2
< д
K :
что предотвращало сглаживание контрастных деталей изображения при правильном выборе порога А к • Результат медианной фильтрации Kj^j 2 определялся выражением
j 2 = med (*2ЙЙ ) • (20)
Операция med ((i 2 ) означала вычисление медианы для элементов мат-
i\,i 2eB\ 12'
рицы 1^2, равных значениям оценки /¿(tj, vj) квадратичного ядра и попавших в прямоугольную апертуру B фильтра с размерами Ly х L2 . Применение в данном вычислительном эксперименте комбинированного фильтра с размерами 7 х 7 (med), 7 х 7 (Aver) и А к = 3 на этапе постобработки позволило в 2... 2,5 раза уменьшить величину относительной ошибки идентификации. От-
K 2 - Kкф
носительная ошибка идентификации 5 К =-т-т-для измерения уровня
||К 2|,
шума 0,02 была равна 0,042 (без фильтрации 0,082), для уровня шума 0,04 относительная ошибка равна 0,068 (без фильтрации 0,133), для уровня шума 0,06 относительная ошибка равна 0,098 (без фильтрации 0,186). Эти результаты позволяют сделать вывод о целесообразности фильтрации результатов идентификации. На рис. 3 показаны изолинии оценки ¿2кф (и , V-) для уровня шума 0,04.
ЗАКЛЮЧЕНИЕ
Непараметрическая идентификация ядер Вольтерра для моделирования нелинейных систем является весьма сложной для практики задачей. При использовании в активном эксперименте известных формул обращения выходного сигнала идентифицируемой системы возникают принципиальные трудности, обусловленные некорректностью задачи дифференцирования и неустойчивостью получаемых решений к шумам измерений. Для устойчивого вычисления производных второго порядка в работе строятся сглаживающие сплайны одной и двух переменных, что обеспечивает высокую вычислительную эффективность и низкую методическую ошибку (примерно один процент) идентификации по точным данным. Для минимизации ошибки сглаживания экспериментальных данных сплайнами предложен алгоритм выбора параметров сглаживания сплайнов, основанный на методе ¿-кривой. Для уменьшения случайной ошибки идентификации в работе используется постобработка результатов идентификации, заключающаяся в фильтрации результатов идентификации локально-пространственным комбинированным фильтром. Выполненный вычислительный эксперимент показал приемлемую для практики точность идентификации квадратичного ядра Вольтерра предлагаемым подходом «алгоритм обращения + алгоритм фильтрации».
СПИСОК ЛИТЕРАТУРЫ
1. Каминскас В. Идентификация динамических систем по дискретным наблюдениям. -Вильнюс: Мокслас, 1985. - 152 с.
2. Giannakis G.B., Serpedin E. A bibliography on nonlinear system identification // Signal Processing. - 2001. - Vol. 81, no. 3. - P. 533-580.
3. Volterra-series-based nonlinear system modeling and its engineering applications: a state-of-the-art review / C.M. Cheng, Z.K. Peng, W.M. Zhang, G. Meng // Mechanical Systems and Signal Processing. - 2017. - Vol. 87. - P. 430-364.
4. Солодуша С.В. Квадратичные и кубичные полиномы Вольтерра: идентификация и приложение // Вестник СПбГУ. Прикладная математика. Информатика. Процессы управления. -2018. - Т. 14, вып. 2. - С. 131-144. - DOI: 10.21638/11702/spbu10.2018.205.
5. Solodusha S. V. New classes of Volterra integral equations of the first kind related to the modeling of the wind turbine dynamics // 15th International Conference on stability and oscillations of nonlinear control systems (Pyatnitskiy's Conference). - Moscow, 2020. - P. 35-39. -DOI: 10.1109/STAB49150.2020.9140662.
6. Солодуша С.В. Методика построения интегральных моделей динамических систем: алгоритмы и приложения в энергетике: дис. ... д-ра техн. наук: 05.13.18. - Иркутск, 2018. - 353 с.
7. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. - М.: Наука, 1986. -
285 с.
8. Воскобойников Ю.Е., Боева В.А. Алгоритмы непараметрической идентификации сложных технических систем // Научный вестник НГТУ. - 2020. - № 4 (80). - С. 47-64. -DOI: 10.17212/1814-1196-2020-4-47-64.
9. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. - М.: Наука, 1980. - 345 с.
10. Wang Y. Smoothing splines methods and applications. - A Chapman & Hall, 2011. - 347 p. -(Monographs on Statistics and Applied Probability; vol. 121).
11. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. - Новосибирск: Наука, 1984. - 238 с.
12. Voskoboynikov Yu£., Boeva V.A. Synthesis of smoothing cubic spline in non-parametric identification technical systems' algorithm // IOP Conference Series: Materials Science and Engineering. - 2020. - Vol. 953. - P. 012035. - DOI: 10.1088/1757-899X/953/1/012035.
13. Воскобойников Ю.Е., Боева В.А. Устойчивый алгоритм вычисления смешанных производных в задачах непараметрической идентификации нелинейных систем // Современные наукоемкие технологии. - 2021. - № 4. - С. 25-29. - DOI: 10.17513/snt.38610.
14. RezghiM., Hosseini S.M. A new variant of L-curve for Tikhonov regularization // Journal of Computational and Applied Mathematics. - 2012. - Vol. 231 (2). - Р. 914-924.
15. Cultrera A., Callegaro L. A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems // IOP SciNotes. - 2020. - Vol. 1 (2). - Р. 25004.
16. Воскобойников Ю.Е., Боева В.А. Метод L-кривой для оценивания оптимального параметра сглаживающего кубического сплайна // Международный научно-исследовательский журнал. - 2021. - № 11 (113), ч. 1. - С. 6-13. - DOI: 10.23670/IRJ.2021.113.11.003.
17. Воскобойников Ю.Е. Артефакты вейвлет-фильтрации изображений и их устранение // Автометрия. - 2020. - Т. 56, № 6. - С. 3-10. - DOI: 10.15372/AUT20200601.
18. Meyer Y. Wavelets: algorithms and applications. - Philadelphia: SIAM, 1993. - 133 p.
Воскобойников Юрий Евгеньевич, доктор физико-математических наук, заведующий кафедрой прикладной математики Новосибирского государственного архитектурно-строительного университета (Сибстрин), профессор кафедры автоматики Новосибирского государственного технического университета, заслуженный работник высшей школы РФ. Основное направление научных исследований - методы и алгоритмы решения некорректных задач интерпретации экспериментальных данных; методы и алгоритмы фильтрации сигналов и изображений; идентификация динамических систем. Является автором более
300 научных работ, в том числе пяти научных монографий по решению некорректно поставленных задач. E-mail: voscob@mail.ru
Боева Василиса Андреевна, аспирант Новосибирского государственного архитектурно-строительного университета (Сибстрин) по профилю «Системный анализ, управление и обработка информации», ассистент кафедры прикладной математики Новосибирского государственного архитектурно-строительного университета (Сибстрин). Основное направление научных исследований - методы и алгоритмы непараметрической идентификации реальных технических систем. Является автором 30 научных статей. E-mail: v.boyeva@sibstrin.ru
Voskoboinikov Yuri Evgenievich, graduate of the NSTU (NETI), Department of Automation, Doctor of Physics and Mathematics, Head of the Department of Applied Mathematics at the Novosibirsk State University of Architecture and Civil Engineering (Sibstrin), Professor at the Department of Automation of the Novosibirsk State Technical University, Honored Worker of the Higher School of the Russian Federation. The main direction of scientific research is methods and algorithms for solving ill-posed problems of interpreting experimental data; methods and algorithms for filtering signals and images, as well as identification of dynamic systems. He is the author of over 300 scientific papers including 5 scientific monographs on solving ill-posed problems. E-mail: voscob@mail.ru
Boeva Vasilisa A., postgraduate student, NSTU Department of System Analysis, Control and Data Interpretation; teaching assistant at the Department of Applied Mathematics, Novosibirsk State University of Architecture and Civil Engineering (Sibstrin). Her research interests are currently focused on non-parametric identification methods and algorithms for real engineering systems. She has 18 scientific articles. E-mail: v.boyeva@sibstrin.ru
DOI: 10.17212/2782-2001-2022-1-25-40
Identification of the quadratic kernel of the Volterra equation for modeling non-linear dynamic systems'
Yu.K. VOSKOBOINIKOV1-2-", V-А. BOEVA2b
1 Novosibirsk State Technical University, 20 K. Marx Prospekt, Novosibirsk, 630073, Russian Federation
2 Novosibirsk State University of Architecture and Civil Engineering (Sibstrin), 113 Leningradskaya Street, Novosibirsk, 630008, Russian Federation
a voscob@mail.ru b v.boyeva@sibstrin.ru
Abstract
In the last two decades, integral models have been used to describe the dynamics of stationary nonlinear systems in which the terms of the Volterra series are the kernels. The most commonly used are the linear term (the impulse transition function depends on one variable) and the quadratic term (depending on two variables). An active experiment in which a special combination of rectangular pulses is fed to the input of the system is carried out to select two of its components in the output signal of the identified system - the output of the linear "subsystem" and the output of the "quadratic" subsystem. After isolating the output of the "quadratic" subsystem, the identification of the quadratic term of the Volterra series is reduced to solving a two-dimensional integral equation of the first kind. In the literature, inversion formulas are given in which the quadratic kernel function is obtained as a result of arithmetic operations with second-order derivatives of the output signal. Differentiation of functions is an incorrectly posed problem, when small errors in the assignment of a function (measurement noise) cause large errors in
* Received 14 July 2021.
The reported study was funded by RFBR, project number 20-38-90041. The reported study was funded by RSF, project number 22-21-00409.
derivatives (especially in second-order derivatives). The paper proposes the use of smoothing cubic splines for stable calculation of derivatives. To calculate the mixed second-order derivative, a spline with two variables is built - a smoothing bicubic spline. The main problem that arises in practice when processing the data of a real experiment is the selection of a smoothing parameter, on the value of which a smoothing error of noisy data depends. As a rule, the value of the variance of the measurement noise is unknown in the experiment. Therefore, in this work, it is proposed to use an algorithm based on the L-curve method to select the smoothing parameter in the constructed splines (especially in the bicubic one), which does not require setting the variance of the measurement noise. The proposed identification algorithm has a high computational efficiency. The performed computational experiment showed a small methodical error (about 1 %) and good resistance to noise in measurements of the output signals of the identified system. To reduce the random component of the identification error, it is proposed to use post-processing with a local-spatial compound filter.
Keywords: integral models of non-linear systems, Volterra kernels, inversion formula, differentiation of noise contaminated data, smoothing cubic spline, smoothing bicubic spline, smoothing parameter selection, local-spatial compound filter
REFERENCES
1. Kaminskas V. Identifikatsiya dinamicheskikh sistem po diskretnym nablyudeniyam [Identification of dynamic systems based on discrete observations]. Vilnius, Moxlas Publ., 1985. 152 p. (In Russian).
2. Giannakis G.B., Serpedin E. A bibliography on nonlinear system identification. Signal Processing, 2001, vol. 81, no. 3, pp. 533-580.
3. Cheng C.M., Peng Z.K., Zhang W.M., Meng G. Volterra-series-based nonlinear system modeling and its engineering applications: a state-of-the-art review. Mechanical Systems and Signal Processing, 2017, vol. 87, pp. 430-364.
4. Solodusha S.V. Kvadratichnye i kubichnye polinomy Vol'terra: identifikatsiya i prilozhenie [Quadratic and cubic Volterra polynomials: identification and application]. Vestnik Sankt-Peterburg-skogo Universiteta. Prikladnaya matematika. Informatika. Protsessy upravleniya = Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2018, vol. 14, iss. 2, pp. 131-144. DOI: 10.21638/11701/spbu10.2018.205.
5. Solodusha S.V. New classes of Volterra integral equations of the first kind related to the modeling of the wind turbine dynamics. 15th International Conference on stability and oscillations of nonlinear control systems (Pyatnitskiy's Conference). Moscow, 2020, pp. 35-39. DOI: 10.1109/ STAB49150.2020.9140662.
6. Solodusha S.V. Metodikapostroeniya integral'nykh modelei dinamicheskikh sistem: algoritmy i prilozheniya v energetike. Diss. dokt. tekhn. nauk [Methods for constructing integral models of dynamic systems: algorithms and applications in energy. Dr. eng. sci. diss.]. Irkutsk, 2018. 353 p.
7. Tikhonov A.N., Arsenin V.Ya. Metody resheniya nekorrektnykh zadach [Methods for solving ill-posed problems]. Moscow, Nauka Publ., 1986. 285 p.
8. Voskoboynikov Yu.E., Boeva V.A. Algoritmy neparametricheskoi identifikatsii slozhnykh tekhnicheskikh sistem [Non-parametric identification algorithms for complex engineering systems]. Nauchnyi vestnik Novosibirskogo gosudarstvennogo tekhnicheskogo universiteta = Science bulletin of the Novosibirsk state technical university, 2020, no. 4 (80), pp. 47-64. DOI: 10.17212/1814-11962020-4-47-64.
9. Zav'yalov Yu.S., Kvasov B.I., Miroshnichenko V.L. Metody splain-funktsii [Methods of splines]. Moscow, Nauka Publ., 1980. 345 p.
10. Wang Y. Smoothing splines methods and applications. A Chapman & Hall, 2011. 347 p.
11. Voskoboynikov Yu.E., Preobrazhensky N.G., Sedelnikov A.I. Matematicheskaja obrabotka jeksperimenta v molekuljarnoj gazodinamike [Mathematical experiment proceeding in molecular gas dynamics]. Novosibirsk, Nauka Publ., 1984. 238 p.
12. Voskoboynikov Yu.E., Boeva V.A. Synthesis of smoothing cubic spline in non-parametric identification technical systems' algorithm. IOP Conference Series: Materials Science and Engineering, 2020, vol. 953, p. 012035. DOI: 10.1088/1757-899X/953/1/012035.
13. Voskoboynikov Yu.E., Boeva V.A. Ustoichivyi algoritm vychisleniya smeshannykh pro-izvodnykh v zadachakh neparametricheskoi identifikatsii nelineinykh sistem [Stable algorithm of mixed derivatives calculation for non-parametric identification problems in nonlinear systems]. Sov-remennye naukoemkie tekhnologii = Modern High Technologies, 2021, no. 4, pp. 25-29. DOI : 10.17513/snt.38610.
14. Rezghi M., Hosseini S.M. A new variant of L-curve for Tikhonov regularization. Journal of Computational and Applied Mathematics, 2012, vol. 231 (2), pp. 914-924.
15. Cultrera A., Callegaro L. A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems. IOP SciNotes, 2020, vol. 1 (2), p. 25004.
16. Voskoboynikov Yu.E., Boeva V.A. Metod L-krivoi dlya otsenivaniya optimal'nogo parametra sglazhivayushchego kubicheskogo splaina [L-curve method for evaluating the optimal parameter of a smoothing cubic spline]. Mezhdunarodnyi nauchno-issledovatel'skii zhurnal = International Research Journal, 2021, no. 11 (113), pt. 1, pp. 6-13. DOI: 10.23670/IRJ.2021.113.11.003.
17. Voskoboinikov Yu.E. Artefakty veivlet-fil'tratsii izobrazhenii i ikh ustranenie [Artefacts of wavelet filtration of images and their elimination]. Avtometriya = Optoelectronics, Instrumentation and Data Processing, 2020, vol. 56, no. 6. DOI: 10.15372/AUT20200601. (In Russian).
18. Meyer Y. Wavelets: algorithms and applications. Philadelphia, SIAM, 1993. 133 p.
Для цитирования:
Воскобойников Ю.Е., Боева В.А. Идентификация квадратичного ядра уравнения Вольтерра для моделирования нелинейных динамических систем // Системы анализа и обработки данных. -2022. - № 1 (85). - С. 25-40. - DOI: 10.17212/2782-2001-2022-1-25-40.
For citation:
Voskoboinikov Yu.E., Boeva V.A. Identifikatsiya kvadratichnogo yadra uravneniya Vol'terra dlya modelirovaniya nelineinykh dinamicheskikh sistem [Identification of the quadratic kernel of the Volterra equation for modeling non-linear dynamic systems]. Sistemy analiza i obrabotki dannykh = Analysis and Data Processing Systems, 2022, no. 1 (85), pp. 25-40. DOI: 10.17212/2782-2001-20221-25-40.
ISSN2782-2001, http://journals.nstu.ru/vestnik Analysis and data processing systems Vol. 85, No 1, 2022, pp. 25-40
