CC BY
118
УДК 62-506.1
О. В. Шестернева, Т. В. Мальцева
МЕТОД ПОСТРОЕНИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ЛИНЕЙНОГО ДИНАМИЧЕСКОГО ОБЪЕКТА
Рассматривается авторский метод построения математической модели линейного динамического объекта, описываемого обыкновенным дифференциальным уравнением. Предлагаемый подход основан на предварительном определении порядка дифференциального уравнения объекта методами непараметрического моделирования и позволяет строить адекватные модели динамических объектов с меньшими затратами средств вычислительной техники.
Ключевые слова: модель линейного динамического объекта, оценка Надарая-Ватсона, коэффициент размытости, адекватность модели.
В современном мире множество технологических процессов, а также экономических и социальных систем могут быть представлены как некий линейный динамический объект, поэтому проблема создания адекватных математических моделей линейных динамических объектов является весьма актуальной. Существуют различные методы, решающие эту задачу, однако сам процесс моделирования зачастую достаточно трудоемкий (в основном за счет сложности выбора структуры модели) и требует больших затрат. В настоящей статье предложен подход, позволяющий упростить процесс выбора структуры модели объекта путем сочетания непараметрических и параметрических методов математического моделирования. Идея заключается в предварительном определении порядка дифференциального уравнения, описывающего объект, и последующем использовании полученной информации в создании параметрической модели. Порядок уравнения предлагается определять путем построения регрессионной непараметрической модели между входными и выходными сигналами объекта, после чего задача моделирования сводится к определению значений параметров параметрической модели известными методами, например, методом наименьших квадратов.
Построение непараметрической оценки регрессии. Прежде всего, необходимо указать, какой априорной информацией об объекте располагает исследователь. В данной работе рассматриваются объекты (процессы), относящиеся к классу линейных динамических (эти сведения априори имеются). Какой-либо другой информации о структуре объекта нет. Предполагается также, что существует возможность измерения входного сигнала и, поступающего на объект (объект представляет собой лишь некоторую часть более сложного процесса, в который исследователь не вмешивается, таким образом, поступающие на вход объекта данные не зависят от воли экспериментатора и могут быть лишь измерены), а также сигнала, полученного на выходе х (реакция объекта на входное воздействие). Измерения производятся в моменты времени t (/) со случайными помехами. Данные измерений формируют обучающую выборку {и (/'), х(/')}, /' = 1, п некоторого объема п, которая в дальнейшем будет являться единственной априорной информацией об объекте.
Остановимся кратко на вопросе построения непараметрической оценки кривой регрессии. Принципиальное отличие таких методов оценивания от параметри-
ческих заключается в том, что последние требуют знание структуры исследуемого объекта (процесса) с точностью до набора параметров и направлены на определение неизвестных параметров (при этом используются выборочные данные). Непараметрический подход позволяет отказаться от выбора структуры объекта и требует только наличия адекватной информативной выборки. В данной работе используется понятие непараметрической оценки регрессии, аппроксимирующей неизвестные стохастические зависимости по наблюдениям. Ставится задача построения оценки неизвестной зависимости между входным и выходным сигналами объекта при любом входном значении сигнала (априори вид стохастической зависимости не задан, предполагается, что она однозначная).
Непараметрическая (ядерная) оценка регрессии, основанная на использовании широко известной оценки плотности распределения Розенблатта-Парзена [1; 2], носит имя Надарая-Ватсона и имеет вид [3; 4]
Ли) = Х х-ф
и -и(г) И
/ "
Тф
и -и(г) И
(1)
Идея, лежащая в основе выражения (1), состоит в придании относительно большего веса наблюдениям, ближайшим к оцениваемой точке в смысле расстояния, определяемого ядром Ф(-).
Функция Ф(-) - ядро (колоколообразная, дельтаобразная функция) - удовлетворяет некоторым условиям сходимости [5], влияние же вида ядра на точность оценивания незначительно. В данной работе использовалось параболическое ядро:
Ф( I) =
(0,75 • (1 -1)2
¡0, Ы > 1.
< 1,
Параметр И в формуле (1) - коэффициент размытости, настройка которого производится согласно условию минимума среднеквадратичного критерия методом скользящего экзамена (по всей выборке, за исключением одной точки, вычисляется оценка регрессии (1), а в этой точке осуществляется проверка качества оценки). Заметим, что с ростом И сглаживающие свойства оценки нарастают, по И для каждого конечного объема выборки
существует некоторый оптимум (при малых И оценка представляет собой набор непересекающихся или слабо пересекающихся дельтаобразных функций и теряет свой смысл, а при больших И оценка становится сильно сглаженной и не отражает индивидуальных особенностей оцениваемой зависимости).
Доказана состоятельность оценки (1), а в случае выполнения двух условий для параметра размытости [1; 6; 7]:
И(п) ® 0, п • И(п) ® ¥,
п®¥ п®¥
приведенная оценка является также асимптотически несмещенной и асимптотически нормально распределенной. Зависимость коэффициента размытости от объема выборки п имеет вид
И(п) = с • п- 9, (2)
где 0 < 9 < 1; 0 < с < ¥.
При использовании квадратичного критерия наилучшего соответствия [6] было получено, что 9 = 1/ (т + 4), где т - размерность вектора и (число входных воздействий). Таким образом, настройка коэффициента размытости сводится к настройке параметра с.
Для многомерного входа оценка (1) имеет вид [1; 6]
\{и ) = Х Xі (0ПФ
i=1 1=1
(3)
ше данных, полученных на предыдущих шагах, необходимо учитывать.
Первоначально по выборке строится оценка регрессии (3), где в качестве аргументов используются входной (входные) и выходной сигналы объекта, а также значение выхода объекта в предыдущий момент времени (двумерная оценка регрессии). Минимизация критерия (4) по двум параметрам дает оптимальные значения коэффициентов размытости, а значение критерия является минимальной среднеквадратичной ошибкой и может быть использовано как показатель адекватности построенной непараметрической модели объекта (если это значение устраивает исследователя, то модель принимается).
Далее строятся непараметрические оценки регрессии вида (3), учитывающие все большее и большее число 5 предыдущих выходных сигналов (которые выступают в (3) в качестве аргументов X). В общем виде такой подход может выглядеть следующим образом:
x(t, u, x(t - j )) =
XXi -Ф
i =1
хП Ф
j=1
и -u(i)
_h
x(t - j) - x(t(i) - j)
(t - j)
/хпф1
/ ,=1 M
Настройка значений коэффициентов размытости в выражении (3) осуществляется одним из методов оптимизации путем минимизации среднеквадратичного критерия, который с учетом формулы (2) принимает следующий вид:
1 n
W(с) = - • X (х(1) - х‘„ (U, с ))2 ® min. (4)
n ,=i с
В настоящей работе использовался метод случайного спуска, где в качестве алгоритма поиска локального минимума был выбран последовательный симплексный метод [7; 8].
Определение порядка дифференциального уравнения объекта. Первая часть предложенного в работе метода построения математической модели линейного динамического объекта основана на определении порядка дифференциального уравнения, описывающего объект исследования, используя методы непараметрической аппроксимации стохастической зависимости входного (в общем случае входных) и выходного (выходных) сигналов. По выборочным данным строится непараметрическая оценка регрессии (2), где в качестве аргументов используется как входное воздействие на текущем шаге, так и значения выходного сигнала на предыдущих шагах. Такой подход позволяет учитывать динамику объекта, так как значения выходного сигнала объекта на нескольких шагах, являясь аргументами оценки регрессии (2) на последующих шагах, влияют на оценку выхода. Число предыдущих шагов, включаемых в модель, является аналогом порядка дифференциального уравнения: чем выше порядок, тем длиннее период функционирования объекта, влияющий на последующее его поведение, и тем боль-
Хф
i =1
хП ф
j=i
и -u(i)
_h
(5)
x(t - j) - x(t(i) - j)
h
(t - j)
Кроме того, в данной работе была сделана попытка оценить не только порядок производной 5 в левой части дифференциального уравнения, описывающего объект, но также порядок производной по управлению г (правая часть дифференциального уравнения), с учетом этого формула (5) будет преобразована к виду
х(/,и, х(/ - ]),и(/ - к)) =
'и -и(})
X x(i) -Ф
к,
хПФ
j=1
r
хПФ
x(t - j) - x(t(i) - j)
hx (t - j )
u(t - і) - u(t(i) - і)
h
u (t- і)
Хф( і-р-у.
i =1 hu
хПФ
j=1
x(t - j) - x(t(i) - j)
h
(t - j)
ф| и(ґ - к) - и(/ (і) - к)
к=1 ^ К (/-к)
В конечном счете, вход и выход объекта зависит от времени, таким образом, можно переписать непараметрическую модель в виде
Ёх-ф
і=1
<П ф
І=1
Г
хПф
х(ї) = -
и(ї) - и(і)
х(ї - І) - х(ї(і) - І) к (ї-І ) и(ї - к) - и(ї(і) - к)
к
и (ї- к )
ЁФ
і=1
<П ф
І=1
<П ф
и(ї) - и(і) х(ї - І) - х(ї(і) - І)
кх (ї- І )
и(ї - к) - и(ї(і) - к)
(6)
к
и (ї- к )
и с учетом этого выводить график выхода модели: сама модель описывается множественной регрессией, однако, все ее входные воздействия зависят от времени, поэтому поведение модели будет рассматриваться в различные периоды времени, без указания ее входных воздействий.
Для проверки предлагаемого алгоритма была проведена имитация линейного динамического объекта, в результате чего получена выборка зашумленных значений входного и выходного сигналов некоторого объема. Помеха накладывалась следующим образом: измерялся интервал изменения сигнальной части Д, задавался уровень помех р (от 0 до 1). С помощью генератора случайных чисел формировался вектор (размерность вектора совпадала с объемом выборки) значений равномерно распределенной на интервале [-Др; Д р] случайной величины, который впоследствии складывался с вектором значений сигнальной части.
Имитируемый объект описывался дифференциальным уравнением третьего порядка
йъ х(ї) й2 х(ї)
йї йх(ї) йї
3 + а2
йї2
+ а0 • х(ї) = Ь0 • и(ї),
Согласно данным табл. 2 в случае достаточно часто снимаемых измерений (выборка в этом случае является более информативной и позволяет лучше прослеживать динамику объекта моделирования) наименьшее значение среднеквадратичного критерия достигается при включении в модель трех предыдущих измерений выхода объекта. Это соответствует третьему порядку дифференциального уравнения (напомним, что сымитированный объект описывался дифференциальным уравнением третьего порядка). При увеличении шага дискретизации измерений динамика прослеживается хуже, что влияет на точность определения порядка дифференциального уравнения. Тем не менее, определяемый порядок близок к истинному, и полученная модель даже в таких случаях является адекватной.
Сравнение выхода объекта и выхода модели (с разным числом предыдущих шагов)
(7)
а = [1; 1; 2; 1], Ь0 = 2,
что априори предполагалось неизвестным. Информация об уровне помех отсутствовала. По выборочным данным было проведено непараметрическое исследование порядка дифференциального уравнения путем построения моделей (6), графические результаты которого приведены на рис. 1, а численные значения критерия (4) сведены в табл. 1.
Наилучшее (минимальное) значение среднеквадратичного критерия достигнуто при учете двух предыдущих шагов, что соответствует второму порядку дифференциального уравнения, однако, следует заметить, что это значение практически такое же, как и в случае трех предыдущих шагов, включаемых в модель (аналог дифференциального уравнения третьего порядка).
Были проведены исследования работы алгоритма для разного уровня помехи, действующей в каналах измерений, а также для различного объема и информативности выборки (частота дискретных измерений входного и выходного сигналов). Эти данные представлены в табл. 2, 3.
Рис. 1. Результаты непараметрического моделирования
Работоспособность предлагаемого метода проверена также для случая разного уровня помех, действующих в каналах измерений. Был рассмотрен случай отсутствия помех, а также случаи незначительной (10 % от полезного сигнала) и большой (80 %) помехи (табл. 3).
В тех случаях, когда помеха небольшая либо вовсе от -сутствует, определяемый порядок дифференциального уравнения, описывающего поведение исследуемого объекта, совпадает с истинным. При большой помехе качество выборочных данных снижается, что приводит к снижению точности определения порядка (отметим, что большие помехи создают дополнительные трудности построения адекватной модели объекта и в случае других часто применяемых методов моделирования).
Таким образом, основным фактором, влияющим на работоспособность предлагаемого метода определения порядка, является качество (информативность, точность) выборочных данных, что естественно, так как выборка -это единственная априорная информация, которой обладает исследователь. Однако даже в случае несовпадения истинного и получаемого порядка дифференциального уравнения, построенная модель может оказаться адекватной, например, за счет проверки значимости коэффициентов. Следует также отметить, что на практике истинный порядок дифференциального уравнения неизвестен, что зачастую приводит к ошибкам выбора структуры модели, тем не менее, если результаты моделирования устраивают заказчиков, модель принимается.
Построение параметрической модели объекта. Проверка адекватности. Последним этапом предлагаемого
метода является построение параметрической модели объекта. В частности, так как структура модели определена ранее методами непараметрического моделирования, задача сводится к нахождению оценок неизвестных параметров модели. Представим структуру модели (дифференциальное уравнение, известное с точностью до параметров) в виде разностного уравнения, так как этот тип модели является наиболее простым, но имеет довольно общий характер [9]. Применяя метод наименьших квадратов (в работе рассматривался наиболее простой случай некоррелированных равноточных измерений), получим уравнение для вектора оптимальных оценок параметров модели [7; 9].
В частном случае, когда объект описывался дифференциальным уравнением (7), непараметрические модели приведены на рис. 1, а параметрическая структура модели имела вид (наименьшее значение среднеквадратичный критерий (4) принимал при включении в модель трех предыдущих выходов объекта, см. табл. 2 для шага 0,2)
й3 х(ї) й2 х(ї)
йї3
- + а
йї2
йх(ї)
+а1---+ а0 • х(ї) = в0 • и(ї),
йї
Наконец, ставится вопрос об адекватности полученной модели. Данная проверка может быть осуществлена по одному из многочисленных критериев адекватности модели регрессии (см. разделы регрессионного и дисперсионного анализа в работе [10]). В данной работе в качестве характеристики точности подбора параметрической модели регрессии являлось так называемое значение Я2, которое в общем случае определяется формулой [10]
Я 2=Ё(хШо<ш -х)2 / Ё (х- х)2,
I=1 / /=1
где числитель представляет собой сумму квадратов, обусловленную регрессией, знаменатель - общую сумму квадратов, скорректированную на среднее х. Понятие сумм квадратов вводится для учета разброса (дисперсии) данных [10].
Выход построенной параметрической модели и данные выборки
(8)
где а = [а0, а1; а2, а3]; р0- оценки истинных парамет-
ров объекта, процедура МНК дала следующие значения параметров модели:
а = [1, 1,01, 1,9, 1,002], р0 = 1,98. (9)
Выход полученной модели хшО(Ш при выбранной структуре (9) и найденных параметрах (10) приведен на рис. 2, в сравнении с выборочными данными.
Время
Рис. 2. Выход параметрической модели и данные выборки
Таблица 1
Зависимость величины среднеквадратичной шибки моделирования от числа учитываемых предыдущих измерений выхода объекта
Число предыдущих измерений выхода объекта, включаемых в модель 1 2 3 4 5
Значение среднеквадратичного критерия 1,319 3 0,337 3 0,339 1 0,417 2 0,451 1
Таблица 2
Зависимость среднеквадратичной ошибки моделирования от числа учитываемых в модели шагов
Число измерений выхода Шаг в модели дискретизации измерений 1 2 3 4 5
0,1 2,584 8 0,180 2 0,104 5 0,235 3 0,273 0
0,2 1,319 3 0,339 3 0,312 3 0,417 2 0,451 1
0,4 1,190 9 0,609 8 0,617 8 0,818 9 0,893 6
0,8 1,043 5 1,144 7 1,425 2 1,510 6 1,526 3
1 1,901 3 1,656 6 1,983 6 1,865 2 1,769 5
2 1,912 0 2,134 4 2,172 0 1,659 4 1,652 1
Таблица 3
Поведение среднеквадратичной ошибки моделирования в зависимости от уровня помехи
———_____Число измерений выхода — модели Уровень помехи, % —— 1 2 3 4 5
0 2,025 5 0,150 1 0,145 7 0,209 3 0,250 2
10 2,584 8 0,180 2 0,164 5 0,235 3 0,273 0
80 2,775 3 1,421 6 1,360 9 1,304 0 1,303 0
Значение Я2 представляет собой меру «вклада кривой модели регрессии в общее отклонение от среднего» [10] и часто выражается в процентах (чем ближе значение Я2 к 1, тем лучше подобранная модель описывает объект).
Во всех рассмотренных случаях модель оказывалась адекватной (в случае, изображенном на рис. 2, Я2 = 0,903 3). Однако, при больших помехах, когда порядок был определен как пятый, коэффициенты при старших производных были близки к нулю (по сравнению с остальными), и исключение из модели соответствующих членов не повлияло на ее адекватность. Осуществление проверки значимости коэффициентов в некоторых случаях позволяет упростить получаемую модель (существует множество статистических пакетов обработки данных, позволяющих осуществлять проверку значимости коэффициентов и адекватности регрессии различными методами).
Проверка также может быть сделана на основании ^-критерия или /-критерия значимости регрессии [10], где выбор уровня значимости предоставляется исследователю.
Таким образом, в работе исследован метод построения параметрической модели линейного динамического объекта, базирующийся на непараметрическом подходе предварительного определения структуры модели. По сравнению с наиболее распространенным способом выбора структуры (перебор всех возможных вариантов) предлагаемый метод требует значительно меньше машинного времени, а лежащие в его основе расчеты гораздо проще (построение непараметрической модели во многих случаях является задачей менее сложной, нежели выбор структуры модели и настройка ее параметров, особенно в случае высоких порядков дифференциальных уравнений). Кроме того, не требуется постановка эксперимента (создания оптимальных планов покачивания), которая на многих реальных объектах является дорогостоящей либо вовсе невозможна. Измерение значений входных и выходных сигналов (являющееся единственной необходимой информацией для построения модели) в настоящее время не представляет особой проблемы, даже в случае малоинформативной обучающей выборки (небольшой объем) построенная модель оказывалась адек-
ватной. Проверена работоспособность алгоритма и в случае сильного зашумления снимаемых данных. Вопрос о выбросах не рассматривался, тем не менее, предполагается возможность использования робастных непараметрических оценок для более точного определения порядка дифференциального уравнения, описывающего объект [7]. Неточности в определении порядка зачастую могут быть устранены проверкой адекватности модели и значимости коэффициентов, что является во всех отношениях более простой задачей, чем построение многочисленных моделей для выбора наилучшей из них.
Библиографический список
1. Parzen, E. On estimation of probability Density Function / E. Parzen // Ann. Math. Stat. 1962. Vol. 33. P. 1065-1076.
2. Rosenblatt, M. Remarks on some non-parametric estimates of a density function / M. Rosenblatt // Ann. Math. Stat. 1956. Vol. 27. P. 832-837.
3. Надарая, Э. А. О непараметрических оценках плотности вероятности и регрессии / Э. А. Надарая // Теория вероятностей и ее применение. 1965. Т. 10(1). С. 199-203.
4. Watson, G. Smooth regression analysis / G. Watson // Sankhya. Ser. A. 1965. Vol. 26. Part 4. P. 359-372.
5. Лагутин, М. Б. Наглядная математическая статистика : учеб. пособие / М. Б. Лагутин. М. : БИНОМ. Лаборатория знаний, 2007.
6. Рубан, А. И. Идентификация стохастических объектов на основе непараметрического подхода / А. И. Рубан // Автоматика и телемеханика. 1979. N° 11. С. 106-118.
7. Рубан, А. И. Методы анализа данных : учеб. пособие / А. И. Рубан. 2-е изд., исправл. и доп. Красноярск : ИПЦ КГТУ 2004.
8. Химмельблау, Д. Прикладное нелинейное программирование / Д. М. Химмельблау. М. : Мир, 1975.
9. Пащенко, Ф. Ф. Введение в состоятельные методы моделирования систем : учеб. пособие : в 2 ч. Ч. 2. Идентификация нелинейных систем / Ф. Ф. Пащенко. М. : Финансы и статистика, 2007.
10. Дрейпер, Н. Прикладной регрессионный анализ : пер с англ. / Н. Дрейпер, Г. Смит. 3-е изд. М. : Изд. дом «Вильямс», 2007.
O. V. Shesterneva, T. V. Maltseva ON MATHEMATICAL LINEAR DYNAMIC OBJECT MODELLING
The construction method of the linear dynamic object model described with an ordinary differential equation is considered. The suggested approach is based on a prior degree of differential equation estimate with nonparametric modelling methods and allows to get adequate models using smaller cost of techniques than.
Keywords: linear dynamic object model, Nadaraya-Watson estimation, bandwidth, model adequacy.
