Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,899

NEW STABLE ALGORITHM OF NONPARAMETRIC IDENTIFICATION OF TECHNICAL SYSTEMS

Voskoboynikov Yu.E. 1, 2 Boeva V.A. 1
1 Novosibirsk State University of Architecture and Civil Engineering (Sibstrin)
2 Novosibirsk State Technical University
Mathematical models of many technical systems have the form of Volterra integral equation of the first kind with a difference kernel. For such systems, the task of identification is to build an estimate for the impulse transition function of the system from the registered values ??of the input and output signal. Such a task is incorrectly posed, in which one or several correctness conditions may be violated. To find a unique and sustainable solution, various (both deterministic and statistical) regularization methods are used. To apply these methods to the problem under consideration, the original Volterra equation is approximated by a system of linear algebraic equations with a poorly conditioned matrix, for the solution of which regularization algorithms are already used. This approach has several disadvantages. First, the errors of registration of the input signal are not taken into account when constructing a regularizing algorithm (only indirectly in some algorithms for choosing a regularization parameter – in the so-called generalized residual method and its variants). Secondly, if it is impossible to specify the required discretization step, the approximation error of the original integral equation is large. These two points lead to an increase in the total identification error. In this paper, we propose a different approach to constructing a stable identification algorithm, based on the transformation of the original equation of the first kind to the Volterra integral equation of the second kind, whose solution is already a correctly posed problem. The use of smoothing cubic splines for calculating the first derivatives of the input and output signals included in the new integral equation allows not only to take into account, but also to significantly reduce the influence of measurement errors of signals on the identification error. Completed studies have shown the advantages of the proposed approach.
identification problem
Volterra equation of the first kind
Volterra equation of the second kind
algorithm of solution of the Volterra equation of the second kind
smoothing cubic splines
selection of the smoothing parameter

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

vosk03.wmf (1)

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

vosk04.wmf при vosk05.wmf. (2)

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

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

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

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

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

vosk06.wmf. (3)

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

vosk07.wmf,

vosk08.wmf (4)

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

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

vosk10.wmf (5)

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

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

vosk13.wmf (6)

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

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

vosk15.wmf. (7)

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

vosk18.wmf, (8)

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

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

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

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

vosk20.wmf, (9)

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

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

vosk26.wmf. (10)

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

vosk27.wmf (11)

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

vosk28.wmf (12)

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

vosk29.wmf (13)

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

vosk30.wmf

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

vosk31.wmf (14)

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

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

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

vosk35.wmf (15)

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

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

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

vosk37.wmf (16)

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

vosk39.wmf (17)

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

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

vosk42.wmf (18)

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

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

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

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

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

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

vosk1.tif

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

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

vosk2.tif

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

Заключение

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

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

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