УДК 681.51:519.6
Иследования эффективности использования сглаживающих кубических сплайнов в задачах непараметрической
идентификации
Ю.Е. Воскобойников1'2, В.А. Боева2 1 ФГБОУ ВПО НГТУ, 2ФГБОУ ВПО НГАСУ (Сибстрин), Новосибирск, Россия
Аннотация: Математические модели многих технических систем имеют вид интегрального уравнения Вольтерра I рода с разностным ядром. Для таких систем задача идентификации заключается в построении оценки для импульсной переходной функции системы по зарегистрированным значениям входного и выходного сигнала. Такая задача является некорректно поставленной, в которой могут нарушаться одно или несколько условий корректности. Для нахождения единственного и устойчивого решения используют различные (как детерминированные, так и статистические) методы регуляризации. В недавней работе авторов был предложен устойчивый алгоритм, суть которого заключается в переходе от решения вольтеровского уравнения I рода к решению уравнению II рода, что является уже корректной задачей. Вычисление первых производных входного и выходного сигналов, содержащихся в интегральном уравнении II рода, на основе аппарата сглаживающих кубических сплайнов позволило снизить негативное влияние шумов измерений на общую ошибку идентификации. Однако, применение сглаживающих сплайнов обусловило ряд вопросов, связанных с выбором параметра сглаживания сплайна из условия минимума ошибки вычисления первых производных при различных распределениях шумов измерений. Поэтому данная работа посвящена исследованиям, результаты которых позволят дать ответы на эти вопросы.
Ключевые слова: задача идентификации, уравнение Вольтера I рода, уравнение Вольтера II рода, сглаживающие кубические сплайны, выбор параметра сглаживания при различных распределениях шумов измерений.
ВВЕДЕНИЕ И ПОСТАНОВКА ЗАДАЧИ
Интегральные модели достаточно часто используются для описания и моделирования технических систем [1]. Если система является стационарной (т.е. ее параметры не меняются во времени), то для ее описания применяется
интегральное уравнение Вольтера I рода [1,2]:
t
Jk(t - т)ф(т)Л = f (t), t е[0,Г], (1)
k (t)
- импульсная переходная функция
где
(ИПФ) системы; ф(т)f (t) - ее входной и выходной сигналы. Техническую реализуемость системы определяет условие:
k(t) = 0 при t < 0 . (2)
Применительно к модели (1), задача непараметрической идентификации [2, 3] системы сводиться к вычислению оценки для k (t)
по измеренным значениям входного и выходного сигналов идентифицируемой системы из решения интегрального уравнения I
рода с разностным ядром вида:
t
Jф^-т)к(x)dx = f(t), t e[0,T] (3)
0
Сформулированная задача идентификации некорректна, поскольку могут не удовлетворятся одно или несколько условий корректности по Адамару [4, 5]. Наиболее часто
нарушается условие устойчивости решения к погрешностям (ошибкам) исходных данных, когда малые погрешности могут вызвать сколь угодно большие ошибки решения. Для построения устойчивых решений некорректных задач используются специальные методы регуляризации, как детерминированные [4], так и статистические [3, 5], которые непосредственно применяются к интегральному уравнению (1) или его конечномерной аппроксимации. В работе авторов [6] был построен эффективный алгоритм, вычисляющий оценку ИПФ из решения интегрального уравнения Вольтера II рода
к(/) + —[ фа - т)к(тШ = , Ф(0) [ ф(0)
г е[0, Т]. (4)
Решение этого интегрального уравнения является уже корректно поставленной задачей, но необходимо вычислить первые производные входного и выходного сигналов, что также является некорректно поставленной задачей и для устойчивого вычисления производных также требуется применение специальных алгоритмов. Заметим, что с устойчивым вычислением первой производной связан метод идентификации, основанный на
дифференцировании выходного сигнала ^ ' идентифицируемой системы при подаче на ее
вход в момент *0 функции Хэвисайда (сигнала с
о
прямоугольным единичным скачком амплитуды
в момент 0), т. е.
d
к(г) = лг'° )' {] • (5)
Устойчивое вычисление второй производной и производных более высоких порядков требуется в алгоритмах идентификации нелинейных ядер Вольтера, при подаче на вход системы прямоугольных импульсов в разные моменты времени [1].
Приведенные примеры (4), (5) показывают необходимость вычисления первых
производных зашумленных сигналов с достаточно высокой точностью. Аппроксимация производных конечными разностями вносит, как правило, значительную систематическую ошибку и требует согласования шага дискретизации с уровнем шума измерения (для выполнения условия сходимости полученных приближений к значению точной производной), что в реальном эксперименте практически невозможно. Применение сглаживающих сплайнов для дифференцирования зашумленных данных представляется весьма перспективным и наиболее подходящим для этого является сглаживающий кубический сплайн (СКС) дефекта 1, т. е. имеющий непрерывные первую и вторую производные на всем интервале определения сплайна. Такой сплайн использовался в работе [6] для устойчивого дифференцирования входного и выходного сигналов системы, входящих в уравнение (4). Это позволило в случае нормально распределенных шумов измерения получить решение с приемлемой точностью даже при высоких (10 % и выше) уровнях шумов измерений сигналов идентифицируемой системы [6].
Однако, при применении СКС на практике возникают следующие затруднения. Во-первых, используемые алгоритмы выбора параметра сглаживания ориентированы на получение наименьшей ошибки сглаживания, т.е. ошибки оценивания точных значений дифференцируемой функции. Для уравнений (4), (5) требуется минимизировать ошибку
дифференцирования и этому минимуму может соответствовать другое значение параметра сглаживания. Во-вторых, в научных
публикациях по сглаживающим сплайнам в качестве шума (погрешностей) измерений принимается модель шума, имеющего нормальное распределение. Опять возникает вопрос о выборе параметра сглаживания при других распределениях шума, встречающихся на практике (например, равномерное распределение, распределение Пуассона и т.д.). Поэтому данная работа посвящена исследованиям, которые дадут ответы на эти вопросы.
СГЛАЖИВАЮЩИЕ КУБИЧЕСКИЕ СПЛАЙНЫ И ВЫБОР ПАРАМЕТРА СГЛАЖИВАНИЯ
Предположим, что на некотором интервале
[0, Т] N
заданы узлов
0 = ^ < г2 <... < ^ = т (6)
ф. = ф(7. ), , = 1,..., N и определены значения т' ' '
функции ^ в этих узлах. Функция ф ( ) называется кубическим сплайном с узлами (6), если:
а) на каждом интервале [ !, г+1 ^ функция
^ (t)
является кубическим многочленом вида:
5ф (х) = а + ь (х - х)+с>(х - х )2 + (х - х )3;
(7)
^Ч 4, С)
б) функция Ф дважды непрерывно
дифференцируема на всём интервале [0, Т ] • Сплайн называется интерполяционным, если выполняются условия интерполяции
5ф &) = ф,-, г = 1,..., N , (8)
т.е. интерполяционный сплайн проходит через
заданные значения функции ),
удовлетворяет определенным краевым условиям
в точках ^, и этот сплайн имеет наименьшую кривизну, определяемую
1" (г )|2 ¿г
функционалом к (подробнее см. [5, 7]).
Предположим, что точные значения ф., г = 1,...,N
т'' ' ' , искажены шумом измерений
л , г = 1,...,N г , с определенными числовыми
характеристиками и регистрируются
(измеряются) значения:
% = ф, +л,, , = 1,..., N .
Очевидно, что в этом случае для фильтрации
л
шумов г сплайн не должен проходить через зашумленные значения, т.е. удовлетворять условиям (8). Поэтому перейдем к сглаживающему кубическому сплайну (СКС)
(г)
Ф> , удовлетворяющему естественным краевым условиям на вторую производную
^Ф'- (г) в виде ^ = ^) = 0 и
доставляющему минимум функционалу [5, 8]:
N п
^ (5) = а1"(/)|2 л + X Р- (Ф, - 5(г,-))2, (9)
1 ,=1
Р-1,, = 1,..., N „
где , - весовые множители. Как
видно из (9), сглаживающий сплайн не
Ф,
проходит через значения , и поэтому он используется для фильтрации (сглаживания) зашумлённых значений. Параметр сглаживания
а «управляет» гладкостью сплайна (а, следовательно, и ошибкой сглаживания): а = 0
- при " и сглаживающий сплайн становится интерполяционным (т.е. проходит
Ф,
и отсутствует
через все заданные точки фильтрация шума);
- при а сглаживающий сплайн
вырождается в прямую линию (эффект переглаживания зашумленных данных).
Алгоритм вычисления коэффициентов сглаживающего сплайна при заданном а подробно изложен в [6, стр. 45-49] и здесь не приводится. Заметим только, что первоначально решается СЛАУ с квадратной пятидиагональной
(N - 2)х(Ж - 2) матрицей размером - находится
вектор вторых производных, используя который
затем вычисляются коэффициенты полиномов
(7).
От величины параметра сглаживания а существенно зависит ошибка сглаживания. Очевидно желание выбрать параметр таким, чтобы его значение минимизировало величину ошибки сглаживания, которую определим, как среднеквадратическую ошибку (СКО) сглаживания
Д(а) = M M [■]
I (ф(t,) - ф (t,.))
(10)
где л - оператор математического ожидания по распределению шума измерения. Для
а„
вычисления оптимального значения
ор1
необходима априорная информация о точных ф(*,-)
значениях , которая на практике,
естественно, отсутствует.
В работе [5] рассмотрены несколько алгоритмов выбора параметра сглаживания. Показано, что эффективной (наилучшей)
оценкой для ор является величина , вычисленная на основе статистического критерия оптимальности фильтрации зашумлённых данных. Обоснование этого критерия и его свойства подробно изложены в [5]. Здесь приведём только конечные соотношения, необходимые для вычисления
а
оценки
KW
Введём в рассмотрение статистику
1 N 0 Р№ (а) = — £ />„
° '=1 (11)
О
где е' ,а ) - невязка ' -ого измерения.
Для случая нормального распределения шумов
измерения доказано, что если Рж (а)
некотором значении неравенству
при удовлетворяет
О < (а) <Э р (12)
—,Ы 1--, N
г 2
то такое значение можно принять в качестве а
оценки для ор (обозначим это значение как
а
,0 в
ж ). В неравенстве (12) величины 2 2
X2 N
квантили л -распределения с 1У степенями в ,1-в
свободы уровней 2' 2 соответственно.
Величина Р определяет вероятность ошибки первого рода при проверке гипотезы об
а
и, как правило, ап
оптимальности оценки
р 0 05. Заметим, что вычисление сводится к решению нелинейного уравнения
Рж(а)= N
а
итерационными алгоритмами. В качестве ж принимается очередное приближённо решение
п(")
а , которое удовлетворяет неравенству (12).
Эффективный алгоритм вычисления ,
построенный на основе метода Ньютона, приведён в [5]. Выполненные исследования
точности оценки аж [5] показали, что сплайн, а = аж
построенный при ж имеет:
а) ошибку сглаживания, незначительно (на 5-8 %) превышающую ошибку сглаживания а = а„
при параметре
ор1
(который можно найти
только в вычислительном эксперименте);
б) ошибку сглаживания, значительно (на 15-35%) меньше по сравнению с другими способами выбора параметра (подробнее см. [5, 8, 9]).
Заметим, что эти свойства оценки аж получены для нормального распределения шума измерения. Для ответа на вопрос «сохранятся ли эти свойства при других распределениях шума» посвящены следующие исследования.
ИССЛЕДОВАНИЯ АЛГОРИТМА ВЫБОР ПАРАМЕТРА СГЛАЖИВАНИЯ
В качестве тестового сигнала была принята функция, график которой показан сплошной кривой на рис. 1. Первая производная показана на этом рисунке штриховой линией. В качестве шумов измерений генерировались некоррелированные между собой псевдослучайные числа, подчиняющиеся следующим распределениям: • нормальному распределению с нулевым средним и дисперсией а2м;
а
равномерному распределению с
нулевым средним и дисперсией
^;
экспоненциальному распределению (распределению Пуассона) с плотностью распределения
10, при х < 0;
Р(х,Х, Н, „-л,
I Pe
, при х > 0,
с математическим ожиданием
M (х ) =
(среднее не равно нулю) и дисперсией Б (х) = -1.
Характеристики этих распределений определялись через требуемый относительный уровень шума, который определялся выражением:
8 ИФ%-Ф|1
где ф, (%( - векторы составленные из значений
Ф,.,ф%, , = 1,...,N, Ц - евклидова норма вектора.
На рис. 1 точечной кривой показаны значения сигнала, искаженные шумом с равномерным распределением и 8 = 0.10. Точность
сглаживания и дифференцирования
зашумленного сигнала определялись
сигнала относительными ошибками:
8» =
ЦФ
8Ф' (а) =
-Ф'
где фа - вектор составленный из значений 5ф,а(^),, = 1,...,N, (',Ф' - векторы составленные из значений производных
Ф'(г,), (г,.), г = 1,..., N, соответственно.
Рис. 1. Тестовый сигнал и его первая производная
Перейдем к ошибкам сглаживания и дифференцирования при параметре
сглаживания аг. В качестве иллюстрации минимизации ошибки сглаживания при этом параметре на рис. 2 приведены следующие кривые: сплошная кривая - относительная ошибка сглаживания точечной кривой, приведенной на рис. 1 (равномерный шум с 8 = 0.10), штриховая кривая - относительная
ошибка дифференцирования, штрих-точечная кривая - значения рг(а); точечные прямые -
квантили
, ^0
0.025,180' 0.975,180
Значения а, для которых рг(а) находится между точечными прямыми (т.е. удовлетворяется неравенство (12)), могут быть приняты в качестве оценки аг для оптимального параметра сглаживания и эти значения находятся в области минимальной ошибки сглаживания. Таким образом параметр сглаживания а может использоваться для построения СКС с минимальной ошибкой сглаживания функций, искаженных шумами с равномерным распределением. Однако точка минимума ошибки дифференцирования немного смещена вправо от точки минимума ошибки сглаживания.
Рис.2. Выбор параметра сглаживания аг при равномерном распределении шума
На Рис. 3 приведены значения тех же кривых, что и на Рис. 2, но полученных для шума измерений, подчиняющегося
распределению Пуассона. Анализ этих кривых позволяет сделать те же выводы, что и в случае равномерного шума.
Рис.3. Выбор параметра сглаживания аг при пуассоновском распределении шума
На Рис. 4 показаны точные значения первой производной тестовой функции (сплошная
кривая) и значения первой производной СКС, построенного при а = аг. Несмотря на ненулевое среднее шума измерения (экспоненциальное распределение при 8 = 0.10) удалось вычислить первую
производную с приемлемой точностью 8ф'(аг) = 0.206. Заметим, что в случае равномерного распределения шума эта величина была немного меньше: 8 , (аг) = 0.192 .
Рис. 4. Значения точной и вычисленной производных при экспоненциальном распределении шума
Заметим, что приведенные результаты были определены по одной реализации зашумленного сигнала и являются одной реализации соответствующих случайных величин, что явно недостаточно для объективных выводов о точности сглаживания и дифференцирования при разных распределениях шумов. Поэтому был проведен вычислительный эксперимент в котором генерировалось 50 (это объем выборки Ы!ат) векторов шума измерения
), I = 1,..., , формировались зашумленные векторы измерений ф%7) = ф + л(1),I = 1,...,, и по каждому вектору измерений вычислялись соответствующие характеристики, которые затем усреднялись по всей выборке, т.е. вычислялось выборочное среднее каждой характеристики.
Таблица 1
8n 8 (a ф V op ) aU 8Ф'(а'opt )
0.02 2.7 -10-4 9.6-10" 3 9.3 -10-4 0.087
0.05 1.1 10-3 0.022 2.1-10-3 0.152
0.10 2.6 -10-3 0.040 4.5 -10-3 0.232
0.15 4.6 -10-3 0.055 7.2-10-3 0.298
оптимального параметра сглаживания для вычислении первой производной а'ор1,
относительной ошибки сглаживания 8ф(аор(),
относительной ошибки дифференцирования
8ф' (а^() при различных распределениях шумов
измерений (нормальное распределение -Табл. 1, равномерное распределение - Табл. 2, экспоненциальное распределение - Табл. 3).
Таблица 2
8n aopt 8Ф(аopt ) a'opt Ma'opt )
0.02 2.8 -10-4 9.7 -10-3 9.5 -10-4 0.088
0.05 1.14-10-3 0.023 2.2-10-3 0.153
0.10 2.4 -10-3 0.041 4.5-10-3 0.232
0.15 4.6 -10-3 0.056 7.5 -10-3 0.303
Таблица 3
8n aopt Ma opt apt 8 ,(a' ф' V op
0.02 1.1-10-4 0.014 6.3 -10-4 0.081
0.05 5.2 -10-4 0.034 1.9-10-3 0.137
0.10 1.6-10-3 0.066 3.9-10-3 0.209
0.15 2.8 -10-3 0.098 5.2 -10-3 0.271
В Таблицах 1, 2, 3 приведены средние значения следующих характеристик: оптимального параметра сглаживания а ,
ВЫВОДЫ
Анализ данных таблиц 1, 2, 3 и графиков Рис. 2, 3 позволяют сделать следующие выводы:
• оценка аш оптимального параметра
сглаживания а , полученная первоначально
для нормального распределения шума измерения может успешно использоваться в сглаживающих кубических сплайнах для сглаживания шумов с равномерным и экспоненциальным распределениями;
• оптимальный параметр сглаживания а'ф незначительно смещается от ^ в сторону
больших значений и поэтому оценка а может также использоваться в СКС при вычислении первой производной от зашумленного сигнала;
• в силу ненулевого значения математического ожидания шума с экспоненциальным распределением минимальная ошибка сглаживания выше (в 1.5 -2 раза) по сравнению с нормальным распрямлением, однако минимальная ошибка дифференцирования ниже по сравнению с тем же нормальным распрямлением;
• в большинстве случаев значение дисперсии шума измерений неизвестно и в этом случае целесообразно обратиться к оценке
)
дисперсии, предложенной в работе [10] и позволяющей с приемлемой точностью (3-8%) оценить дисперсию шума.
Обобщая вышеизложенное, можно сделать вывод о целесообразности использования сглаживающих кубических сплайнов с параметром сглаживания аг для
дифференцирования зашумленных входных и выходных сигналов в задачах
непараметрической идентификации
динамических систем.
ЛИТЕРАТУРА
[1] Сидоров Д. Н. Методы анализа интегральных динамических моделей: теория и приложения. Иркутск: Изд-во ИГУ, 2013. 293 с.
[2] Бойков И. В. Аналитические и численные методы идентификации динамических систем. Пенза: ПГУ, 2016. 396 с.
[3] Воскобойников Ю.Е. Устойчивые алгоритмы непараметрической идентификации динамических систем. Новосиб. гос. архитектур.-строит. ун-т (Сибстрин). - Новосибирск: НГАСУ (Сибстрин), 2019. - 160 с.
[4] Тихонов А. Н. Методы решения некорректных задач. М: Наука, 1986. 285 с.
[5] Воскобойников Ю. Е. Математическая обработка эксперимента в молекулярной газодинамике. Новосибирск: Наука. 1984. 238 с.
[6] Воскобойников Ю.Е., Боева В.А. Новый устойчивый алгоритм непараметрической идентификации технических систем. Современные наукоемкие технологии. № 5. 2019. с. 25-29.
[7] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука. 1980. 352 с.
[8] Y. Wang. Smoothing Splines Methods and Applications. Ser. Monographs on Statistics and Applied Probability v. 121. A Chapman & Hall book. 2011. 347 P.
[9] Wahba G. Smoothing noisy data with spline functions. Numer. Math. 1975. V. 24, № 2. P. 383393.
[10] Воскобойников Ю.Е., Крысов Д.А. Оценивание характеристик шума измерения в модели «сигнал+шум». Автоматика и программная инженерия. 2018. № 3(25). С. 54-61.
Юрий Евгеньевич
Воскобойников, выпускник кафедры автоматики НГТУ (НЭТИ), доктор физ.-мат. наук, профессор,
Заслуженный работник
Высшей школы РФ, Соросовский профессор, действительный член МАИ, РАЕ, МАН ВШ, заведующий кафедрой прикладной
математики Новосибирского государственного архитектурно-строительного университета (Сибстрин), профессор кафедры
автоматики НГТУ. E-mail: [email protected]
Василиса Андреевна Боева,
выпускница кафедры
автоматики НГТУ (НЭТИ), аспирантка кафедры
прикладной математики НГАСУ (Сибстрин). Автор девяти публикаций по идентификации динамических систем и обработке данных. Новосибирск, просп. К. Маркса, д.20, НГТУ
Статья поступила 14.11.2019.
Researches of Efficiency of Using Smoothing Cubic Splines in Nonparametric
Identification Problems
Yu.E. Voskoboinikov V.A. Boeva2
1 NSTU, 2 NGASU (Sibstrin), Novosibirsk, Russia
Abstract. Mathematical models of many real technical systems have the representation in form of Volterra integral equation of the first kind with a difference kernel. For such systems, the identification problem is to calculate the estimation of pulse transition characteristics of the system from the registered values of input and output signals. The problem is ill-posed due to the violation of one or several well-posedness conditions. To calculate the unique stable solution various regularization algorithms, both deterministic and statistical, are applied. In their recent paper, authors propose the stable identification algorithm based on transition from original Volterra integral equation of the first kind solution to the well-posed problem of second kind integral equation solution. Smoothing cubic splines calculation for first derivatives of the input and output signals in new integral equation allows to reduce the negative impact of measurement noises on total level of identification error. However, smoothing cubic splines algorithm underlines a number of questions related to smoothing spline parameter selection depending on minimal first derivatives error condition with different measurement noises distributions. Current studies are dedicated to investigations provided all of answers stated.
Key words: identification problem, integral Volterra equation of the first kind, integral Volterra equation of the second kind, smoothing cubic splines, smoothing spline parameter choice with different measurement noises distributions.
REFERENCES dinamicheskikh modeley: teoriya i prilozheniya.
Irkutsk: Izd-vo IGU, 2013. 293 s. [1] Sidorov D. N. Metody analiza integral'nykh [2] Boykov I. V. Analiticheskiye i chislennyye metody
identifikatsii dinamicheskikh sistem. Penza: PGU, 2016. 396 s.
[3] Voskoboynikov Yu.Ye. Ustoychivyye algoritmy neparametricheskoy identifikatsii dinamicheskikh sistem. Novosib. gos. arkhitektur.-stroit. un-t (Sibstrin). - Novosibirsk: NGASU (Sibstrin), 2019. - 160 s.
[4] Tikhonov A. N. Metody resheniya nekorrektnykh zadach. M: Nauka, 1986. 285 s.
[5] Voskoboynikov Yu. Ye. Matematicheskaya obrabotka eksperimenta v molekulyarnoy gazodinamike. Novosibirsk: Nauka. 1984. 238 s.
[6] Voskoboynikov Yu.Ye., Boyeva V.A. Novyy ustoychivyy algoritm neparametricheskoy identifikatsii tekhnicheskikh sistem. Sovremennyye naukoyemkiye tekhnologii. № 5. 2019. s. 25-29.
[7] Zav'yalov Yu.S., Kvasov B.I., Miroshnichenko V.L. Metody splayn-funktsiy. M.: Nauka. 1980. 352 s.
[8] Y. Wang. Smoothing Splines Methods and Applications. Ser. Monographs on Statistics and Applied Probability v. 121. A Chapman & Hall book. 2011. 347 P.
[9] Wahba G. Smoothing noisy data with spline functions. Numer. Math. 1975. V. 24, № 2. P. 383393.
[10] Voskoboynikov YU.Ye., Krysov D.A. Otsenivaniye kharakteristik shuma izmereniya v modeli «signal+shum». Avtomatika i programmnaya inzheneriya. 2018. № 3(25). S. 54-61.
Yuri Evgenievich
Voskoboinikov, graduate of the Department of Automation NSTU (NETI), Doctor Phys.-Math. Sci., Professor, Honored Worker of the Higher School of the Russian Federation, Soros Professor, Full Member of MAI, RAE, MAN VSH, Head of the Department of Applied Mathematics, Novosibirsk State University of Architecture and Civil Engineering (Sibstrin), Professor, Department of Automation, NSTU. E-mail: [email protected]
Vasilisa Andreevna Boeva, a
graduate of the Department of Automatics of NSTU (NETI), PhD student of the Department of Applied Mathematics of NGASU (Sibstrin). The author of nine publications on the identification of dynamic systems and data processing. Novosibirsk, Prosp. of K. Marksa, h.20, NSTU.
The paper has been received on 14/11/2019.