CC BY
72
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013 Управление, вычислительная техника и информатика № 2(23)
УДК 519.24
Н.А. Сергеева, М.В. Цепкова
О НЕПАРАМЕТРИЧЕСКОМ МОДЕЛИРОВАНИИ ДИНАМИЧЕСКИХ ПРОЦЕССОВ
Рассматривается задача идентификации линейной динамической системы (ЛДС) на основе интеграла свертки с использованием оценки весовой функции. Весовая функция при этом получена в результате реакции ЛДС на входное возмущающее воздействие в виде кусочно-постоянной аппроксимации 5-функции. Приводятся результаты численного моделирования непараметрических алгоритмов.
Ключевые слова: дискретно-непрерывные динамические процессы, весовая функция, непараметрическая идентификация, интеграл Дюамеля.
Настоящая работа посвящена малоизученной проблеме идентификации динамических процессов в «широком» смысле в отличие, от широко распространенной параметрической идентификации или идентификации в «узком» смысле. В частности, рассматривается построение модели линейной динамической системы, когда параметрическая структура объекта неизвестна. Непараметрическая модель ЛДС основывается на соответствующей оценке интеграла Дюамеля, построенного по результатам экспериментов на исследуемом процессе.
1. Идентификация в «узком» и «широком» смыслах
На рис. 1 представлена схема исследуемого процесса с запаздыванием. Особенность постановки задачи состоит в том, что процесс построения модели линейной динамической системы осуществляется в условиях, когда параметрическая структура объекта неизвестна. Ниже приведем общепринятую схему идентификации динамической системы с той лишь разницей, что параметрическая модель, описывающая процесс, неизвестна.
u(t)
ИУ
|4(t)
Процесс
x(t)
uh
ИУ
vh
xt
Рис. 1. Схема динамической системы
h
h
Введены следующие обозначения: и(/) - входная переменная и(/) е Я , х(/) -
выходная переменная процесса, которую без нарушения общности, можем считать скалярной, £(/) - случайное воздействие с ограниченной дисперсией:
D {|(t)} < Const и нулевым математическим ожиданием: м{¡#)} = о, hu,hx-
случайные помехи, действующие в каналах измерения входной и выходной пере-
h h ~ менных процесса, ut, xt - измерения входной и выходной переменных процесса
соответственно, ИУ - измерительное устройство. Для простоты записи будем обозначать измерения входных и выходных переменных Uh, Xh через {(u , xt ), i = 1, s} = {us , xs }, где xs = (Xj, X2...Xs ), us = (Uj, u2...us) - временные векторы [1].
Таким образом, переменная x(t) может быть представлена в виде объективно существующей зависимости
x(t) = A(u(t-0), §(t), t), (1)
где А - неизвестный оператор процесса, 0 - величина запаздывания.
Измерение переменных x(t) и u(t) осуществляется со случайными ошибками, имеющими нулевое математическое ожидание и ограниченную дисперсию, плотность вероятности их неизвестна. Обозначим эти наблюдения xt, ut, t = 1,2..., здесь t дискретное время.
Исследователь при моделировании подобных процессов преследует цель построения математической модели
x(t) = B(u(t -0), t), (2)
где B - класс операторов, который определяется на основании имеющейся априорной информации, x(t) - выход модели. Ясно, что в этом случае стремятся к близости x(t) к x(t) в смысле принятого критерия оптимальности. Очевидно, что модель (2), описывающая процесс (1), может быть удовлетворительной, если класс операторов B выбран удачно. Описанная ситуация является типичной не только для разнообразных технологических процессов, но и для многих других. Проблема моделирования подобных процессов усугубляется недостатком априорной информации об операторе A, что приводит к необходимости рассматривать задачи идентификации в «широком» смысле.
При моделировании разнообразных технологических процессов в настоящее время доминирует теория идентификации в «узком» смысле [1]. Ее содержание состоит в том, что на первом этапе каким-то образом на основе априорной информации определяется параметрический класс операторов B, например
xp (t) = Be (u (t -0), P), (3)
а на втором этапе осуществляется оценка параметров в на основе имеющейся выборки {xt, ut, t = 1, s}, s - объем выборки. Успех решения задачи идентификации в
этом случае существенно зависит от того, насколько «удачно» определен оператор (3).
Идентификация в «широком» смысле предполагает отсутствие этапа выбора параметрического класса оператора (3), если, конечно, для этого нет достаточных априорных сведений. Часто оказывается значительно проще определить класс операторов (2) на основе сведений качественного характера, например линейности процесса или типа нелинейности, однозначности либо неоднозначности и др.
В этом случае задача идентификации состоит в оценивании этого оператора на основе выборки {xt, ut, t = 1, s} в форме
*s (t) = Bs (u (t-0), xs, us). (4)
Для некоторых классов операторов теория идентификации в «широком» смысле была развита на основе методов непараметрической статистики [2]. В этом случае формулы типа (4) являются непараметрическими моделями процесса (1).
При моделировании дискретно-непрерывных процессов следует отличать запаздывание, присущее самому объекту, от задержки, необходимой для измерения переменных объекта. В дальнейшем предполагается, что величина запаздывания 0 известна.
2. Непараметрические модели
Известно, что ЛДС может быть описана в виде интеграла Дюамеля, который при ненулевых начальных условиях имеет следующий вид:
t
x(t) = k (0)u(t) + | h(t -t)u (r)d т, (5)
о
где k(t) и h(t) - соответственно переходная и весовая функции системы. Ниже приведен дискретный аналог интеграла свертки:
t / Дт
х(t) = k(0)u(t) + ^ h(t - Ti )u(Ti )Дт, (6)
i=1
где Дт - шаг дискретизации по времени, Ti = /Дт- значения времени дискретизации. Вычисление выхода x (t) при этом возможно, если известна переходная и импульсная переходная функции системы.
При непараметрическом подходе требуется построить непараметрические оценки переходной и весовой функций, связанных соотношением h(t) = dk / dt [3]. Основой для восстановления оценки переходной функции k(t) служат наблюдения реакции объекта на единичное возмущающее воздействие.
Пусть на вход ЛДС подано единичное возмущающее воздействие 1( t ),0 < t < T, где T - время окончания переходного процесса, а 1(t) - функция Хэ-
висайда. Обозначим {xi = kt, ut = 1, i = 1, s} - реализацию наблюдений «входа-
выхода» объекта, причем наблюдения выходной переменной x(t) осуществляются в дискретном времени через интервал Дt со случайной статистически независимой помехой с нулевым математическим ожиданием и ограниченной дисперсией.
Непараметрической оценкой x(t) (непараметрической моделью ЛДС) будет статистика [2]
t
xs (t) = ks (0)u (t) + | hs (t - r)u (r)d т,
о
где ks (t) и hs (t) - соответственно оценки переходной и весовой функций системы. Ранее [3] предлагалось восстанавливать весовую функцию объекта как производную переходной функции, полученной в результате измерения выхода объекта при подаче на вход единичной функции.
Ниже предлагается подход, при котором в качестве входного воздействия используется аппроксимация 5-функции. В этом случае выходная переменная объекта представляет собой некоторое приближение к весовой функции системы. В этой ситуации отсутствует необходимость оценивания производной переходной функции и тем самым несколько упрощается задача статистического оценивания весовой функции. В дальнейшем рассматриваются различные аппроксимации 5-функции, которые обозначим как 5д.
Для построения оценок весовой функции используется непараметрическая оценка регрессии, построенная по выборке выходных значений динамического процесса, полученных при входном воздействии и(/)=5д. Восстановление весовой функции ht по наблюдениям с помехами осуществляется следующим образом:
h(t)=C £hH (О (7)
здесь h - весовая функция системы, s - объем выборки, Cs - параметр размытости, B(z) - ядерная функция, обладающая известными свойствами сходимости [4], Т -время установления переходного процесса.
В качестве непараметрической модели x(t) предлагаются непараметрические статистики класса Н-аппроксимаций [4], вытекающих из оценок Надарая - Ватсона [5] при равномерной дискретизации аргумента t:
s 1 t/ Дт i t -т. -t Л
xs(t) = ks(0)M(t) + £— £ hHI C—L Iм(тг)Дт.
j=1 SCs i=1 V Cs )
В качестве ядерной функции выбрано усеченное параболическое ядро. Вопросы оптимизации модели по параметру размытости рассмотрены в [3].
3. Вычислительные эксперименты
В ходе численного исследования рассматривается подход к идентификации импульсной переходной (весовой) функции линейного динамического объекта при входном возмущающем воздействии в виде 5Д. Исследуются модели вида (6) при нулевых начальных условиях.
В вычислительном эксперименте исследуемый объект описывался уравнением второго порядка
ax” (t) + bx' (t) + cx(t) = u(t) (8)
или в разностном виде:
2a + bДt b Дt2
ut.
г _ 2 г-1 2 г-2 ~ 2 г
ь + адг+сдг ь + адг+сдг ь + адг + сдг
Далее на вход объекта подавался сигнал и(г) = 5Д, а на выходе наблюдался хг - суть - весовая функция Иг объекта.
На практике невозможно подать 5-функцию на вход объекта, поэтому в вычислительном эксперименте в качестве 5Д была принята система вида
(г) _ (N, г е [0, дт+дг],
и(г) =од = <
д (о, г г[0,дт+дг],
где N - высота ступени, ДТ - ширина интервала, Дг - шаг дискретизации.
На рис. 2, а изображен вид 5Д, на рис. 2, б - оценка весовой функции Н() и ее аналитический аналог (решение дифференциального уравнения (8), когда и(г) = 5(0).
и(г)
N
¿(г), ¿/г)
Рис. 2. Пример 5Д и весовой функции: а - входное воздействие 5Д; б - весовая функция
Представленная на рис. 2 оценка весовой функции в дальнейшем используются при моделировании объекта, так как при данной 5Д она наиболее приближена к аналитическому решению, что повышает точность моделирования.
В рамках данной работы рассматривалось два варианта приближения к 5-функции. Способ задания входного воздействия влияет на известные начальные параметры, на основе которых определяется 5Д.
Входное воздействие объекта можно задавать:
1. Через ширину интервала ДТ;
2. Через высоту ступени N.
Первый способ заключается в задании ширины интервала ДТ и в определении через нее высоты ступени N. В этом случае 5Д рассматривается в виде трапеции, например как на рис. 2. Первоначально определяется шаг сетки, для этого вводится количество интервалов под ступенькой п, оно равно числу отрезков, на которое делится ширина интервала:
Дг =—, (9)
где п - количество интервалов под ступенькой.
Высота ступени рассчитывается из площади трапеции:
N =---------2--------, (10)
(дт + (дт - 2дг))
где Дг - шаг сетки, определяемый по формуле (9).
Второй способ заключается в вычислении ширины интервала ДТ через задаваемую высоту ступени N. При данном способе 5-образная функция принята в виде прямоугольника, ширина интервала соответственно рассчитывается из площади прямоугольника:
ДТ = —.
N
После того как определена ширина интервала или высота ступени, строиться входное воздействие 5Д на всем временном интервале [0, Т], где Т - время установления переходного процесса.
п
После получения весовой функции она оценивается по формуле (7). На рис.3 показана способность непараметрической оценки сглаживать влияние помех в выборке. Показан случай, когда на объект действует помеха, распределенная по нормальному закону с нулевым математическим ожиданием и среднеквадратическим отклонением с = 0,03 (рис. 3, а) и с = 0,1 (рис. 3, б). Принято обозначение: х(г) - выход объекта, при входном воздействии 5Д, И,(г) - оценка весовой функции.
По рис. 3 видно, что оценка весовой функции дает хорошее приближение выхода модели и объекта.
Для проверки адекватности полученной оценки считается ее относительная ошибка на временном отрезке [0, Т]. Формула для расчета ошибки имеет вид
Е=1 £ \И(г) - И,(г )| 100%.
=1 \ Ишах - Итт\
Рассмотрим влияние способа задания входного воздействия на оценку весовой функции. Наблюдается величина относительной ошибки при различных значениях параметров 5Д.
Результаты исследований отражены в табл. 1 и 2, где значениям ширины интервала ДТ и высоты ступени N соответствует относительная ошибка оценки весовой функции е на временном отрезке [0;5]. Табл. 1 соответствует результатам исследования при задании 5Д через ширину интервала ДТ, в табл. 2 приведены результаты при другом способе задания входного воздействия. Объем выборки составляет , = 150 точек, на объект действует помеха, распределенная по нормальному закону с нулевым математическим ожиданием и среднеквадратическим отклонением с = 0,03.
Таблица 1 Таблица 2
N 110 13 10 3 1
е,% 0,62 1,86 2,15 3,84 5,6
ДТ 0,001 0,076 0,1 0,3 1
ДТ 0,01 0,076 0,1 0,3 0,9
е,% 0,52 0,85 1,11 2,57 5,45
N 111 14 13 3 1
Сравнивая полученные результаты, можно заметить, что при близких параметрах 5д величина относительной ошибки е в табл. 1 получается меньше, чем в табл. 2. Поэтому при восстановлении весовой функции входное возмущающее воздействие 5Д рекомендуется задавать через ширину интервала АТ.
На рис. 4 представлена графическая интерпретация полученных результатов. Введено обозначение И(Г) - истинное значение весовой функции, И^) - оценка весовой функции. Истинным значением весовой функции называем аналитическое решение дифференциального уравнения (8) при и(/)=5(/).
Рис. 4. Сравнение весовых функций: а - 5А, заданная через ширину интервала; б - 5а, заданная через высоту ступени
На рис. 4, а аппроксимация задавалась через ширину интервала, на рис. 4, б через высоту ступени. По графикам видно, что в первом случае (рис. 4, а) оценка весовой функции более приближена к истинному значению, эти же результаты подтверждает величина относительной ошибки, которая равна е = 0,52 и 0,82 для первого и второго случаев соответственно. При построении модели используется оценка весовой функции, приведенная на рис. 4, а.
После того как определена оценка весовой функции, строится модель объекта, основанная на дискретной форме интеграла Дюамеля при нулевых начальных условиях, которая имеет вид
/ / Дт
^ (0 = X (/ - т )и(т,. )Дт , (11)
/=1
где И^) - это оценка весовой функции, и(т) - входное воздействие, Ат - шаг дискретизации по времени, т = /Ат - значения времени дискретизации.
В вычислительных экспериментах был принят линейный динамический объект второго порядка вида
х( = 1,67 х(-1 - 0,76 х/-2 + 0,638^.
На первом этапе входное возмущающее воздействие объекта рассматривается в виде щ = 5д. Для восстановления весовой функции используется 5д (вида рис. 2, а), заданная через ширину интервала вышеописанным способом (10). После этого
строится модель объекта (11) на основании наблюдений весовой функции И. При построении модели используется только та оценка весовой функции, которая максимально приближена к ее истинному значению. Когда модель построена на вход объекта и модели подается входное воздействие в виде и (ґ) = -4^(3,6ґ), объем выборки 5 = 150 точек, на объект действует помеха, распределенная по нормальному закону со среднеквадратическим отклонением с = 0,1, объект моделируется на временном отрезке [0;5]. На рис. 5 приведена графическая интерпретация проведенного эксперимента. Введено обозначение х(ґ) - выход объекта, при произвольном входном воздействии и(ґ), х5(ґ) - модель объекта.
Рис. 5. Результат моделирования при и (Г) = -41(^(3,6Г): а - модель объекта; б - оценка весовой функции
По результатам моделирования (рис. 5) видно, что модель практически повторяет поведение объекта, что также подтверждает величина относительной ошибки 0,07.
Изменим условия моделирования. Увеличим объем выборки 5 и используем в качестве входного воздействия функцию и (Г) = 7 Бт(2,5Г) + 7 соб(Г). Объем выборки составит 5 = 200 точек, помеха, действующая на объект, имеет нулевое математическое ожидание и среднеквадратическое отклонение с = 0,07, объект моделируется на временном отрезке [0; 5],
Качество работы модели подтверждается величиной относительной ошибки приближения выхода модели и выхода объекта, которая для этого случая составляет
0,074. По результатам моделирования можно сделать вывод, что непараметрическая модель ЛДС (11) адекватна объекту при произвольных входных воздействиях.
В условиях реального функционирования объекта невозможно подать 5-функцию, так как она представляет собой бесконечный сигнал. При этом если использовать 5д в качестве входного воздействия объекта, то возможно снять на объекте аналог весовой функции. Из проведенных экспериментов можно сделать вывод, что модель, построенная с помощью аппроксимации весовой функции, снятой непосредственно на объекте, может быть адекватна ему.
Рис. 6. Результат моделирования при и(Г) = 78т(2,5г) + 7cos(г) : а - модель объекта; б - оценка весовой функции
Заключение
В заключение следует отметить, что настоящая работа посвящена малоизученному классу задач идентификации в «широком» смысле, т.е. в случае, когда априори параметрическая структура объекта не известна. Ранее непараметрическая идентификация ЛДС осуществлялась по оценкам переходных функций [4].
В настоящей работе сделана попытка исследования вопроса о возможности приближения 5-функции в виде 5-образной функции, подаваемой на объект. При этом непараметрическая модель линейного динамического объекта строится по наблюдениям с ошибками весовой функции объекта в отличие от [4, 5].
Важно отметить, что возможно построение непараметрической модели динамического объекта с запаздыванием, которое содержится в результатах наблюдения выходной переменной объекта х(Г). Модель динамического объекта строится на основании измерения весовой функции системы. Этот путь оказывается эффективным, так как нет необходимости в оценивании производной переходной функции, что с точки зрения статистического оценивания является более трудной задачей.
В данной работе приводились исследования получаемой весовой функции в зависимости от вида дельтаобразной функции 5Д. Исследования показали, что при задании 5Д через ширину интервала весовая функция наиболее приближена к истинному значению, что повышает эффективность непараметрического моделирования дискретно-непрерывных процессов. Многочисленные вычислительные эксперименты показали эффективность предложенных непараметрических моделей.
Приведенные непараметрические модели ЛДС могут служить основанием для настройки параметров ныне существующих и действующих стандартных регуляторов П, ПИ, ПИД [7], которыми оснащены многие технологические процессы и объекты в различных отраслях промышленности.
ЛИТЕРАТУРА
1. Эйкхофф П. Основы идентификации систем управления. Оценивание параметров и состояния. М.: Мир, 1875. 681 с.
2. Medvedev A.V. Identification and control for linear dynamic systems of unknown order // Techniques IFIP, Prague: Verlag, 1975. P. 48-55.
3. Sergeeva N.A. To nonparametric models of nonlinear dynamic systems // CDAM. V. 3. Minsk: BSU, 2001. P. 92-97.
4. МедведевА.В. Непараметрические системы адаптации. Новосибирск: Наука, 1983. 174 с.
5. Надарая Э.А. Непараметрическое оценивание плотности вероятностей и кривой регрессии. Тбилиси: Изд-во Тбилисского университета, 1983. 193 с.
6. Хардле В. Прикладная непараметрическая регрессия. М. Мир, 1993. 349 с.
7. Ротач В.Я. Теория автоматического управления. 5-е изд. М.: Изд. дом МЭИ, 2008. 396 с.
Сергеева Наталья Александровна Цепкова Мария Викторовна Сибирский федеральный университет
E-mail: sergena@list.ru; kidahime@mail.ru Поступила в редакцию 1 февраля 2013 г.
Tsepkova Maria V., Sergeeva Natalia A. (Siberian Federal University. Krasnoyarsk). About nonparametric modeling of dynamic processes.
Keywords: weight function, nonparametric identification, the Dirac delta function, the integral of Duhamel.
The problem of non-parametric modeling of linear dynamic processes is discussed. Discrete form of the Duhamel integral for zero initial conditions is used as the model:
t / At
xs (t) = Z hs (t-Ti )u (Ti)AT i=1
where hs(t) is the assessment of the weight function, «(т) - input action, Дт - time-step discretization, Ti = iДт - the value of time discretization.
As an estimate of the weight function is used nonparametric regression estimator. By definition: the weighting function is a reaction to the system input action in the S-function form. In this paper, as an approximation of S-function system is used:
« ) = Ш,t e [0,AT +At]; w (0,t g [0, AT + At].
Input effect can be given through the width of the interval ог through height of level. These methods have different input 5Д parameters, in the first case, the width of the interval is given and height of level is defined, in the second case, the height of level is fixed and the width of the interval is defined. The accuracy of the weighting function depends on these parameters, so when the input effect is defined with the width of the interval, weight function is closer to the true value than the other way.
If the estimate of the weight function is used in the simulation, which is as close to the true value, then the model is effective in determining the system response to arbitrary input action. This method is applicable for modeling linear dynamic systems, without reference to the order of equation of the describing system.
