В работе рассматривается авторский метод построения математической модели линейного динамического объекта, описываемого обыкновенным дифференциальным уравнением. Предлагаемый подход основан на предварительном определении порядка дифференциального уравнения объекта методами непараметрического моделирования и позволяет строить адекватные модели динамических объектов с меньшими затратами средств вычислительной техники.
ВВЕДЕНИЕ
В современном мире множество технологических, процессов, а также экономических и социальных систем могут быть представлены как некий линейный динамический объект, поэтому проблема создания адекватных математических моделей линейных динамических объектов является весьма актуальной. Существуют различные методы, решающие эту задачу, однако сам процесс моделирования зачастую достаточно трудоемкий (в основном за счет сложности выбора структуры модели) и требует больших затрат. В настоящей статье предложен подход, позволяющий упростить процесс выбора структуры модели объекта путем сочетания непараметрических и параметрических методов математического моделирования. Идея заключается в предварительном определении порядка дифференциального уравнения, описывающего объект и последующем использовании полученной информации в создании параметрической модели. Порядок уравнения предлагается определять путем построения регрессионной непараметрической модели между входными и выходными сигналами объекта, после чего задача моделирования сводится к определению значений параметров параметрической модели известными методами, например, методом наименьших квадратов.
- ПОСТРОЕНИЕ НЕПАРАМЕТРИЧЕСКОЙ ОЦЕНКИ РЕГРЕССИИ
Прежде всего, необходимо указать, какой априорной информацией об объекте располагает исследователь. В данной работе рассматриваются объекты (процессы), относящиеся к классу линейных динамических (эти сведения априори имеются). Какой-либо другой информации о структуре объекта нет. Предполагается также, что существует возможность измерения входного сигнала 
, поступающего на объект (объект представляет собой лишь некоторую часть более сложного процесса, в который исследователь не вмешивается, таким образом, поступающие на вход объекта данные не зависят от воли экспериментатора и могут быть лишь измерены), а также сигнала, полученного на выходе x (реакция объекта на входное воздействие). Измерения производятся в моменты времени 
со случайными помехами. Данные измерений формируют обучающую выборку 
некоторого объема n, которая в дальнейшем будет являться единственной априорной информацией об объекте.
Остановимся кратко на вопросе построения непараметрической оценки кривой регрессии. Принципиальное отличие таких методов оценивания от параметрических заключается в том, что последние требуют знание структуры исследуемого объекта (процесса) с точностью до набора параметров и направлены на определение неизвестных параметров (при этом используются выборочные данные). Непараметрический подход позволяет отказаться от выбора структуры объекта и требует только наличия адекватной информативной выборки. В данной работе используется понятие непараметрической оценки регрессии, аппроксимирующей неизвестные стохастические зависимости по наблюдениям. Ставится задача построения оценки неизвестной зависимости между входным и выходным сигналами объекта при любом входном значении сигнала (априори вид стохастической зависимости не задан, предполагается, что она однозначная).
Непараметрическая (ядерная) оценка регрессии, основанная на использовании широко известной оценки плотности распределения Розенблатта-Парзена [4, 8], носит имя Надарая-Ватсона и имеет вид [3, 10]:
. (1)
Идея, лежащая в основе (1), состоит в придании относительно большего веса наблюдениям, ближайшим к оцениваемой точке в смысле расстояния, определяемого ядром Ф(.).
Функция Ф(.) – ядро (колоколообразная, дельтаобразная функция) –удовлетворяет некоторым условиям сходимости [2], влияние же вида ядра на точность оценивания незначительно. В данной работе использовалось параболическое ядро:
Параметр h в формуле (1) – коэффициент размытости, настройка которого производится согласно условию минимума среднеквадратичного критерия методом скользящего экзамена (по всей выборке, за исключением одной точки, вычисляется оценка регрессии (1), а в этой точке осуществляется проверка качества оценки). Заметим, что с ростом h сглаживающие свойства оценки нарастают, по h для каждого конечного объема выборки существует некоторый оптимум (при малых h оценка представляет собой набор непересекающихся или слабо пересекающихся дельтаобразных функций и теряет свой смысл, а при больших h оценка становится сильно сглаженной и не отражает индивидуальных особенностей оцениваемой зависимости).
Доказана состоятельность оценки (1), а в случае выполнения двух условий для параметра размытости [4, 6, 7]:
приведенная оценка является также асимптотически несмещенной и асимптотически нормально распределенной. Зависимость коэффициента размытости от объема выборки n имеет вид:
(2)
где 0 < q < 1, 0 < c < ∞.
При использовании квадратичного критерия наилучшего соответствия [6], было получено, что q = 1/ (m+4), где m – размерность вектора u (число входных воздействий). Таким образом, настройка коэффициента размытости сводится к настройке параметра c.
Для многомерного входа оценка (1) имеет вид [4, 6]:
. (3)
Настройка значений коэффициентов размытости в (3) осуществляется одним из методов оптимизации путем минимизации среднеквадратичного критерия, который с учетом (2) принимает вид:
. (4)
В настоящей работе использовался метод случайного спуска, где в качестве алгоритма поиска локального минимума был выбран последовательный симплексный метод [7, 9].
Отметим, что в случае больших помех, действующих в каналах измерений, целесообразно для настройки параметра c коэффициента размытости для каждого входного воздействия использовать другой критерий качества. Этот критерий получается из квадратичного интегрального критерия качества раскрытием квадрата разности и последующей заменой x на оценку xn в удвоенном произведении:
. (5)
В критерии (5) в отличие от критерия (4) не участвует исходная выборка, а значит, помехи не будут напрямую влиять на его значение, что позволит получить аппроксимирующую кривую, которая не «тянется» за сильно зашумленными выборочными данными. Однако следует отметить, что критерий (4) является более наглядным (так как непосредственно производится сравнение выхода объекта и модели), поэтому в дальнейшем будут приводиться значения именно этого критерия (см. таблицы 1 – 3), хотя оптимизацию параметров размытости в случае больших помех (более 50 %) целесообразно проводить по критерию качества (5).
- ОПРЕДЕЛЕНИЕ ПОРЯДКА ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ОБЪЕКТА
Первая часть предложенного в работе метода построения математической модели линейного динамического объекта основана на определении порядка дифференциального уравнения, описывающего объект исследования, используя методы непараметрической аппроксимации стохастической зависимости входного (в общем случае входных) и выходного (выходных) сигналов. По выборочным данным строится непараметрическая оценка регрессия (2), где в качестве аргументов используется как входное воздействие на текущем шаге, так и значения выходного сигнала на предыдущих шагах. Такой подход позволяет учитывать динамику объекта, так как значения выходного сигнала объекта на нескольких шагах, являясь аргументами оценки регрессии (2) на последующих шагах, влияют на оценку выхода. Число предыдущих шагов, включаемых в модель, является аналогом порядка дифференциального уравнения: чем выше порядок, тем длиннее период функционирования объекта, влияющий на последующее его поведение, и тем больше данных, полученных на предыдущих шагах, мы должны учитывать.
Первоначально по выборке строится оценка регрессии (3), где в качестве аргументов используются входной (входные) и выходной сигналы объекта, а также значение выхода объекта в предыдущий момент времени (двумерная оценка регрессии). Минимизация критерия (4) по двум параметрам дает оптимальные значения коэффициентов размытости, а значение критерия является минимальной среднеквадратичной ошибкой и может быть использовано как показатель адекватности построенной непараметрической модели объекта (если это значение устраивает исследователя, то модель принимается).
Далее строятся непараметрические оценки регрессии вида (3), учитывающие все большее и большее число s предыдущих выходных сигналов (которые выступают в (3) в качестве аргументов 
). Так в случае первого порядка дифференциального уравнения, описывающего объект, непараметрическая модель будет иметь вид:
.
Если объект описывается дифференциальным уравнением второго порядка, то соответствующая ему непараметрическая модель:
.
Наконец, модель, описывающая поведение объекта, имеющего третий порядок дифференциального уравнения:
В общем виде такой подход может быть представлен в виде:
. (6)
Кроме того, в данной работе была сделана попытка оценить не только порядок производной s в левой части дифференциального уравнения, описывающего объект, но также порядок производной по управлению r (правая часть дифференциального уравнения), с учетом этого формула (6) будет преобразована к виду:
.
В конечном счете, вход и выход объекта зависит от времени, таким образом, можно переписать непараметрическую модель в виде:
(7)
и с учетом этого выводить график выхода модели: сама модель описывается множественной регрессией, однако, все ее входные воздействия зависят от времени, поэтому поведение модели будет рассматриваться в различные периоды времени, без указания ее входных воздействий.
Для проверки предлагаемого алгоритма была проведена имитация линейного динамического объекта, в результате чего получена выборка зашумленных значений входного и выходного сигналов некоторого объема. Помеха накладывалась следующим образом: измерялся интервал изменения сигнальной части Δ, задавался уровень помех ρ (от 0 до 1). С помощью генератора случайных чисел формировался вектор (размерность вектора совпадала с объемом выборки) значений равномерно распределенной на интервале [-Δ.ρ; Δ.ρ] случайной величины, который впоследствии складывался с вектором значений сигнальной части.
Имитируемый объект описывался дифференциальным уравнением третьего порядка
(8)
что априори предполагалось неизвестным. Информация об уровне помех отсутствовала. По выборочным данным было проведено непараметрическое исследование порядка дифференциального уравнения путем построения моделей (7), графические результаты которого приведены на рисунке 1, а численные значения критерия (4) сведены в таблицу 1.
Рисунок 1. Результаты непараметрического моделирования
Таблица 1. Зависимость величины среднеквадратичной шибки моделирования от числа учитываемых предыдущих измерений выхода объекта
|
Число предыдущих измерений выхода объекта, включаемых в модель |
1 |
2 |
3 |
4 |
5 |
|
Значение среднеквадратичного критерия |
1.3193 |
0.3373 |
0.3391 |
0.4172 |
0.4511 |
Наилучшее (минимальное) значение среднеквадратичного критерия достигнуто при учете двух предыдущих шагов, что соответствует второму порядку дифференциального уравнения, однако, следует заметить, что это значение практически такое же, как и в случае трех предыдущих шагов, включаемых в модель (аналог дифференциального уравнения третьего порядка).
Были проведены исследования работы алгоритма для разного уровня помехи, действующей в каналах измерений, а также для различного объема и информативности выборки (частота дискретных измерений входного и выходного сигналов). Эти данные представлены в таблицах 2, 3 (наилучшее значение критерия для каждого случая обозначено курсивом).
Таблица 2. Зависимость среднеквадратичной ошибки моделирования от числа учитываемых в модели шагов
|
Число измерений выхода Шагв модели дискретизации измерений |
1 |
2 |
3 |
4 |
5 |
|
0.1 |
2.5848 |
0.1802 |
0.1045 |
0.2353 |
0.2730 |
|
0.2 |
1.3193 |
0.3393 |
0.3123 |
0.4172 |
0.4511 |
|
0.4 |
1.1909 |
0.6098 |
0.6178 |
0.8189 |
0.8936 |
|
0.8 |
1.0435 |
1.1447 |
1.4252 |
1.5106 |
1.5263 |
|
1 |
1.9013 |
1.6566 |
1.9836 |
1.8652 |
1.7695 |
|
2 |
1.9120 |
2.1344 |
2.1720 |
1.6594 |
1.6521 |
Из данных таблицы 2 можно увидеть, что в случае достаточно часто снимаемых измерений (выборка в этом случае является более информативной и позволяет лучше прослеживать динамику объекта моделирования) наименьшее значение среднеквадратичного критерия достигается при включении в модель трех предыдущих измерений выхода объекта. Это соответствует третьему порядку дифференциального уравнения (напомним, что сымитированный объект описывался дифференциальным уравнением третьего порядка). При увеличении шага дискретизации измерений динамика прослеживается хуже, что влияет на точность определения порядка дифференциального уравнения. Тем не менее, определяемый порядок близок к истинному, и полученная модель даже в таких случаях является адекватной.
Работоспособность предлагаемого метода проверена также для случая разного уровня помех, действующих в каналах измерений. Был рассмотрен случай отсутствия помех, а также случаи незначительной (10% от полезного сигнала) и большой (80%) помехи (данные представлены в таблице 3).
Таблица 3. Поведение среднеквадратичной ошибки моделирования в зависимости от уровня помехи
|
Число измерений выхода в модели Уровень помехи, % |
1 |
2 |
3 |
4 |
5 |
|
0 |
2.0255 |
0.1501 |
0.1457 |
0.2093 |
0.2502 |
|
10 |
2.5848 |
0.1802 |
0.1645 |
0.2353 |
0.2730 |
|
80 |
2.7753 |
1.4216 |
1.3609 |
1.3040 |
1.3030 |
В тех случаях, когда помеха небольшая, либо вовсе отсутствует, определяемый порядок дифференциального уравнения, описывающего поведение исследуемого объекта, совпадает с истинным. При большой помехе качество выборочных данных снижается, что приводит к снижению точности определения порядка (отметим, что большие помехи создают дополнительные трудности построения адекватной модели объекта и в случае других часто применяемых методов моделирования).
Таким образом, основным фактором, влияющим на работоспособность предлагаемого метода определения порядка, является качество (информативность, точность) выборочных данных, что естественно, так как выборка – это единственная априорная информация, которой обладает исследователь. Однако даже в случае несовпадения истинного и получаемого порядка дифференциального уравнения, построенная модель может оказаться адекватной, например, за счет проверки значимости коэффициентов. Следует также отметить, что на практике истинный порядок дифференциального уравнения неизвестен, что зачастую приводит к ошибкам выбора структуры модели, тем не менее, если результаты моделирования устраивают заказчиков, модель принимается.
- ПОСТРОЕНИЕ ПАРАМЕТРИЧЕСКОЙ МОДЕЛИ ОБЪЕКТА. ПРОВЕРКА АДЕКВАТНОСТИ ПОЛУЧЕННОЙ МОДЕЛИ
Последним этапом предлагаемого метода является построение параметрической модели объекта. В частности, так как структура модели определена ранее методами непараметрического моделирования, задача сводится к нахождению оценок неизвестных параметров модели. Представим структуру модели (дифференциальное уравнение, известное с точностью до параметров) в виде разностного уравнения, так как этот тип модели является наиболее простым, но имеет довольно общий характер [5]. Применяя метод наименьших квадратов (в работе рассматривался наиболее простой случай некоррелированных равноточных измерений), получим уравнение для вектора оптимальных оценок параметров модели [5, 7].
В частном случае, когда объект описывался дифференциальным уравнением (8), непараметрические модели приведены на рисунке 1, а параметрическая структура модели имела вид (наименьшее значение среднеквадратичный критерий (4) принимал при включении в модель трех предыдущих выходов объекта, см. таблицу 2 для шага 0.2):
(9)
где 
– оценки истинных параметров объекта.
Процедура МНК дала следующие значения параметров модели:
(10)
Выход полученной модели 
при выбранной структуре (9) и найденных параметрах (10) приведен на рисунке 2 в сравнении с выборочными данными.
Рисунок 2. Выход параметрической модели и данные выборки
Наконец, ставится вопрос об адекватности полученной модели. Данная проверка может быть осуществлена по одному из многочисленных критериев адекватности модели регрессии (см. разделы регрессионного и дисперсионного анализа) [1]. В данной работе в качестве характеристики точности подбора параметрической модели регрессии являлось так называемое значение R2, которое в общем случае определяется формулой [1]:
где числитель представляет собой сумму квадратов, обусловленную регрессией, знаменатель – общую сумму квадратов, скорректированную на среднее 
. Понятие сумм квадратов вводится для учета разброса (дисперсии) данных [1].
Значение R2 представляет собой меру «вклада кривой модели регрессии в общее отклонение от среднего» и часто выражается в процентах (чем ближе значение R2 к 1, тем лучше подобранная модель описывает объект).
Во всех рассмотренных случаях модель оказывалась адекватной (в случае, изображенном на рисунке 2, R2 = 0.9033). Однако, при больших помехах, когда порядок был определен как пятый, коэффициенты при старших производных были близки к нулю (по сравнению с остальными), и исключение из модели соответствующих членов не повлияло на ее адекватность. Осуществление проверки значимости коэффициентов в некоторых случаях позволяет упростить получаемую модель (существует множество статистических пакетов обработки данных, позволяющие осуществлять проверку значимости коэффициентов и адекватности регрессии различными методами).
Проверка также может быть сделана на основании F-критерия или t-критерия значимости регрессии [1], где выбор уровня значимости предоставляется исследователю.
ЗАКЛЮЧЕНИЕ
В работе исследован метод построения параметрической модели линейного динамического объекта, базирующийся на непараметрическом подходе предварительного определения структуры модели. По сравнению с наиболее распространенным способом выбора структуры (перебор всех возможных вариантов) предлагаемый метод требует значительно меньше машинного времени, а лежащие в его основе расчеты гораздо проще (построение непараметрической модели во многих случаях является задачей менее сложной, нежели выбор структуры модели и настройка ее параметров, особенно в случае высоких порядков дифференциальных уравнений). Кроме того, не требуется постановка эксперимента (создания оптимальных планов покачивания), которая на многих реальных объектах является дорогостоящей, либо вовсе невозможна. Измерение значений входных и выходных сигналов (являющееся единственной необходимой информацией для построения модели) в настоящее время не составляет особой проблемы, даже в случае малоинформативной обучающей выборки (небольшой объем) построенная модель оказывалась адекватной. Проверена работоспособность алгоритма и в случае сильного зашумления снимаемых данных. Вопрос о выбросах не рассматривался, тем не менее, предполагается возможность использования робастных непараметрических оценок для более точного определения порядка дифференциального уравнения, описывающего объект [7]. Неточности в определении порядка зачастую могут быть устранены проверкой адекватности модели и значимости коэффициентов, что во всех отношениях задача более простая, чем построение многочисленных моделей для выбора наилучшей из них.
ЛИТЕРАТУРА
1. Дрейпер Н., Смит Г. Прикладной регрессионный анализ, 3-е изд.: Пер с англ. – М.: Издательский дом «Вильямс», 2007. – 912 с.
2. Лагутин М. Б. Наглядная математическая статистика: Учеб. пособие / М. Б. Лагутин. – М.: БИНОМ. Лаборатория знаний, 2007. – 427 с.
3. Надарая Э.А. О непараметрических оценках плотности вероятности и регрессии// Теория вероятностей и ее применение. 1965. Т.10(1). С.199-203.
4. Parzen E. On estimation of probability Density Function. // Ann. Math. Stat. -- 1962. -- Vol.33. -- P.1065--1076.
5. Пащенко Ф. Ф. Введение в состоятельные методы моделирования систем: Учеб. пособие: В 2-х ч. Ч. 2. Идентификация нелинейных систем. – М.: Финансы и статистика, 2007. – 288 с.
6. Рубан А.И. Идентификация стохастических объектов на основе непараметрического подхода// Автоматика и телемеханика. -- 1979. -- N 11. -- С.106-118.),
7. Рубан А. И. Методы анализа данных: Учеб. пособие. 2-е изд., исправл. и доп. / А. И. Рубан. – Красноярск: ИПЦ КГТУ, 2004. – 319 с.
8. Rosenblat M. Remarks on some non-parametric estimates of a density function // Ann. Math. Stat. -- 1956. -- Vol.27. -- P.~832--837.
9. Химмельблау Д. Прикладное нелинейное программирование. – М.: Мир, 1975. – 534 с.
10. Watson G. Smooth regression analysis //Sankhya, ser.A. 1965. Vol.26, part 4. P.~359--372
1
