Научный журнал

Современные наукоемкие технологии

ISSN 1812-7320

"Перечень" ВАК

ИФ РИНЦ = 0,940

• Авторы
• Резюме
• Файлы
• Ключевые слова
• Литература

Воскобойников Ю.Е. 1, 2 Боева В.А. 1

---

1 ФГБОУ ВО «Новосибирский государственный архитектурно-строительный университет (Сибстрин)»

2 ФГБОУ ВО «Новосибирский государственный технический университет»

Математические модели многих технических систем имеют вид интегрального уравнения Вольтерра I рода с разностным ядром. Для таких систем задача идентификации заключается в построении оценки для импульсной переходной функции системы по зарегистрированным значениям входного и выходного сигнала. Такая задача является некорректно поставленной, в ней могут нарушаться одно или несколько условий корректности. Для нахождения единственного и устойчивого решения используют различные (как детерминированные, так и статистические) методы регуляризации. Для применения этих методов к рассматриваемой задаче исходное уравнение Вольтерра аппроксимируют системой линейных алгебраических уравнений с плохо обусловленной матрицей, для решения которой уже используют регуляризирующие алгоритмы. Такой подход имеет ряд недостатков. Во-первых, погрешности регистрации входного сигнала не учитываются при построении регуляризирующего алгоритма (только косвенно в некоторых алгоритмах выбора параметра регуляризации – в так называемом обобщенном методе невязки и его вариантах). Во-вторых, при невозможности задать требуемый шаг дискретизации велика ошибка аппроксимации исходного интегрального уравнения. Эти два момента (и не только они) приводят к возрастанию общей ошибки идентификации. В данной работе предлагается другой подход к построению устойчивого алгоритма идентификации, основанный на преобразовании исходного уравнения I рода к интегральному уравнению Вольтерра II рода, решение которого является уже корректно поставленной задачей. Использование сглаживающих кубических сплайнов для вычисления первых производных входного и выходного сигналов, входящих в новое интегральное уравнение, позволяет не только учесть, но и существенно уменьшить влияние погрешностей измерений сигналов на ошибку идентификации. Выполненные исследования показали преимущества предлагаемого подхода.

Статья в формате PDF

0 KB

задача идентификации

уравнение Вольтерра I рода

уравнение Вольтерра II рода

алгоритм решения уравнения Вольтерра II рода

сглаживающие кубические сплайны

выбор параметра сглаживания

1. Сидоров Д.Н. Методы анализа интегральных динамических моделей: теория и приложения. Иркутск: Изд-во ИГУ, 2013. 293 с.

2. Бойков И.В. Аналитические и численные методы идентификации динамических систем. Пенза: ПГУ, 2016. 396 с.

3. Первушин В.Ф. О непараметрической идентификации линейных динамических объектов // Вестник Томского государственного университета. Сер. Управление, вычислительная техника и информатика. 2013. № 4 (25). С. 95–104.

4. Тихонов А.Н. Методы решения некорректных задач. М.: Наука, 1986. 285 с.

5. Воскобойников Ю.Е. Устойчивые алгоритмы решения обратных измерительных задач / Научная монография. Новосибирск: Изд-во НГАСУ. 2007. 184 с.

6. Воскобойников Ю.Е. Математическая обработка эксперимента в молекулярной газодинамике. Новосибирск: Наука. 1984. 238 с.

7. Wang Y. Smoothing Splines Methods and Applications. Ser. Monographs on Statistics and Applied Probability v. 121. A Chapman & Hall book. 2011. 347 р.

8. Kubo L.H., J. de Oliveira. Smoothing by cubic spline modified applied to solve inverse thermal problem. Computational and Applied Mathematics. 2018. v. 37. no. 2 Р. 1162–1174.

9. Крянев А.В., Лукин Г.В. Метрический анализ и обработка данных. М.: ФИЗМАЛИТ, 2010. 279 с.

Интегральные динамические модели широко используются для описания и моделирования технических систем [1, с. 12–18]. В случае стационарной системы применяется уравнение Вольтерра I рода вида [2, с. 25–29]

(1)

где k(t) – импульсная переходная функция (ИПФ) системы (ядро интегрального уравнения); φ(τ), f(t) – входной и выходной сигналы системы. При этом выполняется условие, которое определяется технической реализуемостью системы

при . (2)

Задача непараметрической идентификации [2, 3] заключается в построении оценки для k(t) по зарегистрированным значениям входного и выходного сигналов идентифицируемой системы. Заметим, что найденная оценка ИПФ в дальнейшем используется как для вычисления количественных характеристик системы (время переходного процесса, величина перерегулирования и т.д.), так и для моделирования системы в целом с использованием модели «вход-выход».

Задача идентификации является некорректно поставленной (НПЗ), в которой может нарушаться одно или несколько требований корректности задачи по Адамару [4, с. 15–18]. Для вычисления устойчивого единственного решения НПЗ используются различные (детерминированные или статистические) регуляризирующие алгоритмы (РА). Для применения РА исходное уравнение или аппроксимируют системой линейных алгебраических уравнений с плохо обусловленной матрицей, или заменяют дискретной сверткой и строят РА на основе дискретного преобразования Фурье (подробнее см. [5, с. 102–106]). И в том, и в другом случаях возникают принципиальные трудности, приводящие к увеличению ошибки идентификации (отсутствует возможность в РА учесть случайный характер погрешностей измерений входного сигнала, сложность формирования периодической дискретной свертки [5, с. 103–106], значительная ошибка аппроксимации исходного интегрального уравнения при невозможности задать требуемый шаг дискретизации и другие).

В данной работе предлагается новый устойчивый алгоритм идентификации, основанный на сведении уравнения (1) к интегральному уравнению II рода (решение которого является корректно поставленной задачей), в которое входят первые производные от входного и выходного сигналов (вычисление которых является НПЗ). Для устойчивого вычисления производных предлагается использовать кубические сглаживающие сплайны, а ошибки дифференцирования минимизируются за счёт выбора параметра сглаживания на основе излагаемого критерия оптимальности. Все это позволяет существенно уменьшить влияние погрешностей регистрации входа и выхода системы на точность получаемых решений задачи идентификации ИПФ стационарной системы.

Метод исследования

Для построения устойчивого алгоритма решения задачи идентификации приведём уравнение (1) к виду уравнения задачи идентификации. При выполнении условия (2) получаем

. (3)

Дифференцируя это уравнение по переменной t и выполнив несложные преобразования, приходим к интегральному уравнению Вольтерра II рода

,

(4)

решение которого уже является корректно поставленной задачей. Однако построение алгоритма вычисления k(t) из этого уравнения сталкивается со следующими основными трудностями: дифференцирование зашумлённых входного и выходного сигналов идентифицируемой системы; вычисление интеграла свёртки с наименьшими ошибками интегрирования для уменьшения общей систематической ошибки алгоритма идентификации. Для преодоления этих трудностей привлечём математический аппарат сглаживающих кубических сплайнов (СКС). Приведём только основные понятия СКС, необходимые для дальнейшего изложения алгоритма идентификации.

Предположим, что на некотором интервале заданы N узлов

(5)

и определены значения функции φ(t) в этих узлах. Функция Sφ(x) называется кубическим сплайном с узлами (5), если:

а) на каждом интервале функция Sφ(t) является кубическим многочленом вида

(6)

б) функция Sφ(t) дважды непрерывно дифференцируема на всём интервале .

Сплайн называется интерполяционным, если выполняются условия интерполяции

. (7)

Естественным сглаживающим кубическим сплайном называется функция S _{φ, α}(t), удовлетворяющая естественным краевым условиям на вторую производную : и доставляющая минимум функционалу [6, с. 42–44; 7, с. 124–132]:

, (8)

где – весовые множители. Как видно из (8), сглаживающий сплайн не проходит через значения φi и поэтому он используется для фильтрации (сглаживания) зашумлённых значений. Параметр сглаживания α «управляет» гладкостью сплайна (а следовательно, и ошибкой сглаживания):

- при α = 0 сглаживающий сплайн становится интерполяционным (т.е. проходит через все заданные точки φ_i);

- при α→∞ сглаживающий сплайн вырождается в прямую линию (эффект переглаживания зашумленных данных).

Алгоритм вычисления коэффициентов сглаживающего сплайна при заданном α подробно изложен в [6, с. 45–49] и здесь не приводится. Для построения алгоритма решения уравнения (4) временно будем считать, что по зашумлённым значениям входного и выходного сигналов

, (9)

где ξi, ηi – случайные погрешности регистрации сигналов с дисперсиями соответственно построены сглаживающие кубические сплайны S _{φ, α}(t),  S _{f, α}(t) с производными соответственно.

Обратимся к интегралу свёртки, входящему в уравнение (4). Для вычисления оценки для k(t) предположим, что на каждом интервале k(t) постоянно и равно . Тогда, используя производную сглаживающего сплайна, интеграл свёртки будем аппроксимировать суммой

. (10)

Можно показать, что

(11)

где k = j – i. Подставляя (11) в (10), получаем

(12)

где Фj,i – элементы матрицы Ф, определяемые через коэффициенты сплайна S _{φ, α}(t):

(13)

где k = j – i. Таким образом, матрица Ф имеет размер (N – 1)×(N – 1), и если i > j, то Фi,j = 0. Заметим, что полученная квадратурная формула (13) позволяет достаточно точно и эффективно вычислять значения интегралов, тем самым уменьшая методическую ошибку предложенного алгоритма. Сформируем векторы

Тогда уравнение (4) можно аппроксимировать следующей системой линейных алгебраических уравнений (СЛАУ) вида

(14)

где I – единичная матрица размером . Решая эту систему, получаем вектор , проекции которого являются оценками для значений k(ti) идентифицируемой ИПФ системы.

Заметим, что СЛАУ (14) имеет хорошо обусловленную матрицу (число обусловленности не больше 3–5), и точность идентификации определяется только ошибками дифференцирования входного (элементы матрицы Ф) и выходного (проекции вектора ) сигналов.

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

(15)

где – оператор математического ожидания. Для вычисления оптимального значения α_opt необходима априорная информация о точных значениях f(ti), которая на практике, естественно, отсутствует.

В работе [6, с. 60–73] рассмотрены несколько алгоритмов выбора параметра сглаживания. Показано, что эффективной (наилучшей) оценкой для α_opt является величина α_W, вычисленная на основе статистического критерия оптимальности фильтрации зашумлённых данных. Обоснование этого критерия и его свойства подробно изложены в [6, с. 62–67]. Здесь приведём только конечные соотношения, необходимые для вычисления оценки α_W.

Введём в рассмотрение статистику

(16)

где – невязка i-го измерения. Доказано, что если ρ_W(α) при некотором значении α удовлетворяет неравенству

(17)

то такое значение можно принять в качестве оценки для α_opt (обозначим это значение как α_W). В неравенстве (17) величины – квантили χ2-распределения с N степенями свободы уровней соответственно. Величина β определяет вероятность ошибки первого рода при проверке гипотезы об оптимальности оценки α_W и, как правило, β = 0,05.

Заметим, что вычисление αW сводится к решению нелинейного уравнения

(18)

итерационными алгоритмами. В качестве α_W принимается очередное приближённое решение α⁽ⁿ⁾ , которое удовлетворяет неравенству (17). Эффективный алгоритм вычисления α_W, построенный на основе метода Ньютона, приведён в [6, с. 65–66]. Выполненные исследования точности оценки α_W [6, с. 72–76] показали, что сплайн, построенный при α = α_W , имеет:

а) ошибку сглаживания, незначительно (на 5–8 %) превышающую ошибку сглаживания при параметре α = α_opt (который можно найти только в вычислительном эксперименте);

б) ошибку сглаживания, значительно (на 15–35 %) меньше по сравнению с другими способами выбора параметра (подробнее см. [6, с. 73–76; 8, 9]).

Всё это позволяет сделать вывод о целесообразности использования сглаживающего сплайна с α = α_W для устойчивого вычисления производных входного и выходного сигналов при формировании матрицы Ф и вектора системы (14).

Результаты численных исследований

Для проверки работоспособности и эффективности предложенного алгоритма идентификации был проведен многочисленный вычислительный эксперимент. Остановимся только на результатах одной серии исследований. В качестве k(t) была взята ИПФ колебательного звена второго порядка (на рис. 1 – сплошная кривая). Значения входного и выходного сигналов искажались случайными нормально распределенными погрешностями с относительными уровнями: где φ, f – векторы, составленные из точных значений сигналов, – векторы зашумленных значений, – евклидова норма вектора. Ошибка идентификации определялась относительной величиной . На рис. 2 точечной кривой показана оценка для ИПФ, построенная по производным интерполяционных сплайнов (т.е. без сглаживания) при уровнях шумов δφ = 0,05, δ_f = 0,05, N = 190, сплошной кривой – значения точной ИПФ. Видно, что полученное решение существенно отличается от точной ИПФ, что иллюстрирует высокую неустойчивость такого решения к погрешностям регистрации входного и выходного сигналов. На рис. 1 точечной кривой показано решение, построенное по сглаживающим кубическим сплайнам (с параметрами сглаживания α_W) при тех же уровнях шумов. Видно, что решение достаточно хорошо совпадает с точной ИПФ системы (относительная ошибка идентификации δ_k = 0,077).

Рис. 1. Импульсная переходная функция и ее оценка

Также было выполнено сравнение предложенного алгоритма идентификации с регуляризирующим алгоритмом, построенным на основе ДПФ (подробно этот РА изложен в [5, с. 104–105]. Относительная ошибка идентификации предложенным алгоритмом была на 25–50 % меньше по сравнению со вторым РА. Это отчасти объясняется ошибкой аппроксимации исходного уравнения системой линейных алгебраических уравнений во втором РА. Заметим, что выигрыш по точности идентификации увеличивался по мере возрастания отношения , например при .

Рис. 2. Оценка ИПФ, построенная по интерполяционным сплайнам

Заключение

Высокая точность предлагаемого алгоритма идентификации объясняется следующими его преимуществами:

1. Матрица СЛАУ (14) хорошо обусловлена, и поэтому ошибки вычисления вектора правой части (ошибки дифференцирования выходного сигнала) гораздо слабее влияют на построенное решение, чем погрешности правой части в РА, построенном на основе ДПФ.

2. Наличие в предлагаемом алгоритме интеграла от производной входного сигнала (при вычислении элементов матрицы Ф (12), (13) в значительной степени «сглаживает» ошибку вычисления первой производной, что обуславливает уменьшение общей ошибки идентификации.

Библиографическая ссылка