Адаптивное параметрическое оценивание квадратно-корневыми информационными алгоритмами

АДАПТИВНОЕ ПАРАМЕТРИЧЕСКОЕ ОЦЕНИВАНИЕ КВАДРАТНО-КОРНЕВЫМИ ИНФОРМАЦИОННЫМИ АЛГОРИТМАМИ

Содержание

Введение

1. Общие выкладки из теории

1.1. Общая постановка проблемы идентификации

1.2. Оценки максимального правдоподобия

1.3. Методы минимизации функций многих переменных

1.4. Метод квадратно-корневого информационного фильтра(ККИФ)

2. Оценивание параметров по методу максимального правдоподобия с использова-нием квадратно-корневых информационных фильтров

2.1. Постановка задачи

2.2. Функция правдоподобия и ее представление в терминах ККИФ

2.3. Градиент функции максимального правдоподобия

2.4. Значения производных переменных ККИФ

2.5. Описание алгоритма

3. Эксперименты и выводы

Заключение

Список используемой литературы

 

Задача идентификации формул и ру е тся следующим образом: по результатам наблюде н ий н ад вход н ыми и выходными переменным и системы должна быть построена оптимальная в некотором смысле модель, т. е. формал и зован н ое пр е дставление этой с и стемы

В зав и с и мост и от апр и орной информац и и об объекте управления различают задач и идентификац ии в узком и ш и роком смысле. Задача идент и фикаци и в у з ком смысле состоит в оцен и ван ии параметров и состоя н ия системы по результатам н аблюден и й над вход н ыми и выходными п е ременными, получен н ыми в условиях функционирования объ е кта. Пр и этом и звестна структура системы и задан класс мод е лей, к которому данный объ е кт относ и тся. Априорная информация об объект е достаточно велика

Апр и орная ин формация об объекте п р и идент и фикации в ш и роком смысл е отсутствует или оче н ь бедная, п оэтому пр и ходится предварительно р е шать большое ч и сло до п олн и тельных задач, такие как выбор структуры системы и задание класса модел е й, оце н ивание линейности объекта и действующих п ерем ен ных, оценивани е степени и формы влия н ия вход н ых переменных на выход н ые и др

Целью данной дипломной работы является исследование нового метода параметрической идентификации основанного на синтезе метода максимального правдоподобия и метода квадратно-корневого информационного фильтра, а также сравнение методов минимизации, использованных для минимизации выбранного функционала, с точки зрения сходимости, вычислительной точности, сложности, а также реализация данного метода на ЭВМ

 

Задача оценивания может быть сформулирована как задача нахождения наибольшего (наименьшего) значения некоторого функционала. Но т.к. значения параметров непосредственному наблюдению не доступны, то критерием выбора оптимума должен быть функционал от выходных значений . Примером такого функционала может служить либо функция правдоподобия, либо ее логарифм. Т.е. если - это выход объекта, -соответствующий выход модели и, также когда, невязки ошибок предсказания являются независимыми и имеют гауссовское совместное распределение с нулевым средним и матрицами ковариаций , тогда выражение для обратного логарифма функции максимального правдоподобия имеет следующий вид:

(1)

Тогда критерием выбора оптимума выберем выражение (2), которая является функцией многих переменных и для ее минимизации будем использовать наиболее известные и часто применяемые методы минимизации функций многих переменных: градиентный метод, метод Ньютона, метод сопряженных направлений.

Оценкой максимального правдоподобия является такое значение оцениваемых параметров , которое максимизирует вероятность события, при котором наблюдения, сгенерированные с подстановкой оцениваемых параметров, совпадают с действительными значениями наблюдений . Эта процедура эквивалентна минимизации обратного логарифма функции плотности условной вероятности невязок, представленный в формуле (2)

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

Для эффективного вычисления оценки максимального правдоподобия при использовании ККИФ, величины, входящие в выражение для (2) и его градиента, непосредственно выражаются либо через выходные значения КИИФ, либо легко находятся путем решения треугольных систем

 

Факт сходимости алгоритма максимального правдоподобия к оптимальным значениям параметров теоретически является недоказанным, поэтому в качестве основного метода исследования будем считать вычислительные эксперименты

Стоит заметить, что метод является достаточно сложным в вычислительном отношении, поскольку метод максимального правдоподобия с использованием ККИФ требует больших объемов вычислений: для перемножения, обращения, ортогональных преобразований матриц и поэтому для проведения экспериментов данный метод был реализован на ЭВМ

Модель, используемая в экспериментах, представленных на графиках, имеет следующий вид:……

В данной дипломной работе проведены эксперименты на сходимость метода максимального правдоподобия, используя различные алгоритмы минимизации. При этом варьировалось количество и расположение оцениваемых параметров в матрице перехода из состояния в состояние , которая, в данной случае, является устойчивой, а модель – наблюдаемой (на рисунках 1-3 представлены изменения оцениваемых параметров, используя при минимизации функционала градиентный метод; на рисунках 4-6 – изменение компонент градиента обратного логарифма функции правдоподобия; на рисунке 7 – нормализованная ошибка оценки параметров; на рисунках 6-21 – соответствующие графики для других методов минимизации). Также проведены эксперименты на выявление зависимости количества времени для одной итерации от количества измерений (рисунок 26), количества времени для одной итерации от размерности вектора оцениваемых параметров (рисунок 27), точности оценивания от количества наблюдений (рисунок 28), на выявления влияния соотношения сигнал/шум на точность оценивания (рисунок 29). На рисунках 22-25 представлена зависимость обратного логарифма функции максимального правдоподобия от параметров

 

После проведения серии вычислительных экспериментов были получены следующие результаты:

  • Как видно из графиков для обратного логарифма функции максимального правдоподобия по параметрам, минимум функции является не единственным, и как следствие этого возникают ситуации, когда методы минимизации сходятся не к истинному значению оцениваемых параметров. Так же стоит заметить, что график функционала, при больших отклонениях от истинных значениях параметров, идет практически параллельно горизонтальной оси координат. Из выше сказанного можно сделать вывод, что выбор начального приближения для параметров может оказать существенное влияние как на сходимость алгоритмов, так и на истинность полученных оценок.
  • На оценки параметров особенно сильное влияние оказывает наблюдаемость динамической системы объекта (наблюдаемость динамической системы является необходимым условием сходимости методов параметрической идентификации), а также соотношение сигнал/шум, причем с ростом соотношения точность оценок увеличивается.
  • Из исследованных алгоритмов наилучшей сходимостью обладает метод сопряженных направлений, а более точным является метод Ньютона, при этом он тоже обладает достаточно хорошей сходимостью. Поэтому предпочтительней использовать метод Ньютона, т.к. при использовании ККИФ матрица вторых производных функционала (в нашем случае это информационная матрица Фишера) вычисляется естественным путем из выходных данных.
  • Установлено, что в общем случае скорость сходимости с ростом размерности вектора параметров и количества наблюдений сильно падает, однако с увеличением количества входных данных растет точность оценок параметров. Но существует некоторый предел, при котором рост точности приостанавливается (при количестве наблюдений более 2500). Поэтому следует искать компромисс между скоростью и точностью.

 

 

Введение

С давних пор человеч е ство затрач и вает огромные усилия на установлен ие зако н омер н остей происходящих в пр и роде явл ен ий. П е рв и ч н ым в проц е сс е познания вс е гда являются результаты наблюдений. Они п ред с тавляют собо й отправ н ой пункт к модели, к абстрактному мышлению, а уже от модели осуществляется переход к практической деятельност и . Очевид н о, что эта схема позна н ия пр и ме н има н е зав и симо от того, и дет ли р е чь о е стестве нн ом или искусстве нн ом объекте. Созда н ие абстракт н ой модели обыч н о связано со “сжат и ем ” информации, содер ж ащейся в результатах н аблюде ний . Это объясняется тем, что каждый отдельный результат наблюден и й является случай н ым, поэтому построение адекватной модели реального объекта может быть осущ е ствле н о только на основе многократных наблюде н ий. Случайность каждого результата наблюде ни й объясняется, с од н ой сторо н ы, при н ципиальной невозможностью учесть все многообразие факторов, действующих на данный ко н крет н ый объект, каким бы простым он н и казался на первый в з гляд, и сложным и взаимосвязями этих факторов, а с другой стороны, несоверше н ством естестве нн ых ил и искусственных средств наблюден и я

Построе н ие модел и по результатам наблюдений представляет собой формализацию, необходимую для опр е делен и я основных приз н аков, связ е й, законом е р н остей, присущ и х объекту-ор и г и налу , и отсеива ни я второстепенных признаков. Во многих случаях модель, принятая при проектировании, существенно отличается от реального объекта, что значительно уменьшает или сводит на нет эффективность разработанной системы управления. В связи с этим возникло одно из новых и важных направлений в теории управления, связанное с построением модели на основании наблюдений, полученных в условиях функционирования объекта по его входным и выходным переменным. Это направление известно в настоящее время как идентификация систем

Задача идентификации формулируется следующим образом: по результатам наблюдений над входными и выходными переменными системы должна быть построена оптимальная в некотором смысле модель, т.е. формализованное представление этой системы

В зависимости от априорной информации об объекте управления различают задачи идентификации в узком и широком смысле. Задача идентификации в узком смысле состоит в оценивании параметров и состояния системы по результатам наблюдений над входными и выходными переменными, полученными в условиях функционирования объекта. При этом известна структура системы и задан класс моделей, к которому данный объект относится. Априорная информация об объекте достаточно велика

Априорная информация об объекте при идентификации в широком смысле отсутствует или очень бедная, поэтому приходится предварительно решать большое число дополнительных задач. К эт и м задачам относятся: выбор структуры системы и задание класса моделей, оценивание степени стационарност и и линейности объекта и действующих переменных, оценивание степени и формы влияния входных переменных на выходные, выбор и н форматив н ых перем енн ых и др. К н астоящ е му време н и накоплен большой о п ыт реш ен ия задач идентификации в узком смысле. М е тоды же реше ни я задач ид е нт и фикац ии в ш и роком смысле н ачали разрабатываться только в последн и е годы, и здесь результаты з н ач и тельно скромнее, что в первую оч е редь можно объясн и ть чр е звычайной трудностью задачи

Целью данной дипломной работы является исследование нового метода параметрической идентификации основанного на синтезе метода максимального правдоподобия и метода квадратно-корневого информационного фильтра (ККИФ), сравнение его с другими существующими алгоритмами с точки зрения вычислительной точности, быстродействия и сложности, а также реализация данного метода на ЭВМ

В первой главе данного дипломного проекта дана общая постановка задачи параметрической идентификации, а также общие сведения об оценке максимального правдоподобия и методах минимизации функций многих переменных. Также в этой главе рассмотрены ключевые моменты метода квадратно-корневого информационного фильтра, показаны его преимущества и недостатки перед стандартным фильтром Калмана

Во второй главе показано, как вычислить оценку максимального правдоподобия итеративным способом при помощи характеристического уравнения, которое включает в себя градиент обратного логарифма функции правдоподобия и информационной матрицы Фишера. В разделе 2.2 второй главы выведено выражение для функции правдоподобия, используя выходные значения естественным образом генерирующиеся ККИФ. В разделах 2.3 и 2.4 показан способ получения градиента обратного логарифма функции правдоподобия и формулы для информационной матрицы Фишера, а полный алгоритм для их подсчета представлен в разделе 2.5

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

В заключении данного дипломного проекта будут подведены итоги проделанной работы

1. Общие выкладки из теории

 

1.1. Общая постановка проблемы идентификации

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

которая позволяла бы приписать неизвестному параметру рассматриваемого объекта некоторое числовое значение (оценку) , причем эта оценка зависит от последовательности наблюдений , где - есть управление. Рассматриваемая ситуация изображена на рис.1

Часто подразумевается, что идентификация начинается из ”ничего” без всякой априорной информации об объекте. Но в большинстве технических задач такое предположение не реалистично; из структуры объекта и, по крайней мере, частичного понимания его функционирования можно извлечь определенную априорную информацию и, в частности, вид структуры модели. В этом случае остается только получить информацию о числовых значениях ряда параметров (коэффициентов дифференциальных уравнений, описывающих динамику объекта, и т. д.). В результате задача идентификации сводится к задаче оценивания параметров. Под оцениванием параметров понимается экспериментальное определение значений параметров, характеризующих динамику поведения объекта, в предположении, что структура модели объекта известна

При оценивании используют различные виды оценок, которые различаются объемом исходной информации об объекте. Например, при нахождении оценки по методу наименьших квадратов предполагается, что динамика объекта может быть аппроксимирована выбранной моделью. При получении ”марковских” оценок считается также известной ковариационная матрица шума. Для вычисления оценок максимального правдоподобия необходимо знание плотности вероятности измеряемого случайного процесса. Байесовские оценки, или оценки с минимальным риском, требуют знания априорных плотностей вероятности неизвестных параметров и величины штрафа за ошибки

В данной работе мы будем рассматривать частный случай показанной выше ситуации, т.е. рассмотрим параметрическое оценивание параметров динамической системы без управления, а неизвестные параметры будем оценивать методом максимального правдоподобия

1.2. Оценка максимального правдоподобия

Известно, что при оценивании существенный интерес представляет информация о рассматриваемых параметрах, в частности плотность вероятности или . Эта функция зависит от продолжительности интервала наблюдений или, что то же самое, от размера выборки и представляет собой наиболее полную информацию, которую можно получить, применяя статистические методы

Пусть известна плотность вероятности шума . По ней можно рассчитать плотность вероятности измерений , которая зависит от параметров объекта и обозначается через . Также положим, что априорное распределение вероятностей параметра в области возможных значений параметра. Тогда разумно выбирать , совпадающее со значением , максимизирующим . В соответствии с формулой Байеса

для произвольного имеем

Априорная плотность вероятности имеет вид

Апостериори после измерений выборочных значений , или , совместная плотность вероятности обозначается как

и называется функцией правдоподобия. По функции правдоподобия находиться оценка , причем логично выбирать такое значение , которое максимизирует . Тогда необходимое условие максимума имеет вид

или в силу монотонности логарифма

Эту формулу называют уравнением правдоподобия и, отыскивая решение этой системы уравнений, которое обеспечивает наибольшее значение или , находим оценку максимального правдоподобия , а, находя оценку максимального правдоподобия, мы тем самым находим неизвестные параметры системы

Таким образом, задача оценивания может быть сформулирована как задача нахождения наибольшего (наименьшего) значения некоторого функционала. Но т.к. значения параметров непосредственному наблюдению не доступны, то критерием выбора оптимума должен быть функционал от выходных сигналов или от математического ожидания ошибок оценок параметров. Примером такого функционала может служить либо функция правдоподобия, либо ее логарифм. Т.е. если - это выход объекта, -соответствующий выход модели и, также когда, невязки ошибок предсказания являются независимыми и имеют гауссовское совместное распределение с нулевым средним и матрицами ковариаций , тогда выражение для обратного логарифма функции максимального правдоподобия имеет следующий вид:

Тогда нахождение оценки максимального правдоподобия эквивалентно минимизации следующего функционала:

(1.2.1)

Тогда

1.3. Методы минимизации функций многих переменных

Критерием выбора оптимума, в нашем случае этим критерием есть выражение (1.2.1), является функция (функционал) многих переменных и для ее минимизации будем использовать наиболее известные и часто применяемые методы минимизации функций многих переменных

В общем случае будем рассматривать задачу

(1.3.1)

предполагая, что функция непрерывно дифференцируема на . Согласно определению дифференцируемости функции

(1.3.2)

где . Если , то при достаточно малых главная часть приращения (1.3.2) будет определяться дифференциалом функции . Справедливо неравенство Коши-Буняковского

,

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

Это замечательное свойство градиента лежит в основе ряда итерационных методов минимизации функций. Одним из таких методов является градиентный метод, к описанию которого мы переходим

  1. Градиентный метод. Будем считать, что некоторая точка уже выбрана. Тогда метод заключается в построении последовательности по правилу:

Существуют различные способы выбора в данном методе, но на практике нередко довольствуются нахождением какого-либо , обеспечивающего условие монотонности: . С этой целью задаются какой-либо постоянной и на каждой итерации метода берут . При этом для каждого проверяют условие монотонности, и в случае его нарушения дробят до тех пор, пока не восстановится монотонность метода

  1. Метод Ньютона. Градиентный метод является методом первого порядка, поскольку использует лишь первые производные минимизируемой функции. Однако если минимизируемая функция дважды непрерывно дифференцируема и производные вычисляются достаточно просто, то возможно применение метода минимизации второго порядка, которые используют квадратичную часть разложения этой функции в ряд Тейлора. Широко известный метод Ньютона представляет собой итерационный процесс:

  1. Метод сопряженных направлений. Метод сопряженных направлений является методом использующим лишь градиент функционала и описывается следующим образом:

,

где величина находится из условия

,

а вычисляется из следующих формул

Точное определение величины возможно лишь в редких случаях, поэтому реализация каждой итерации метода будет сопровождаться неизбежными погрешностями. Эти погрешности, накапливаясь, могут привести к тому, что векторы перестают указывать направление убывания функции, и сходимость метода может нарушиться. Чтобы бороться с этим явлением, метод сопряженных направлений время от времени обновляют, полагая

1.4. Метод квадратно-корневого информационного фильтра(ККИФ)

Пусть динамическая система с дискретным временем дана в следующем виде:

(1.4.1)

где вектор состояния , матрица переходов из состояния в состояние , и - случайный процесс, представляющий собой белый гауссовый шум с нулевым средним и ковариацией , т.е.

Пусть множество наблюдений задается уравнением:

(1.4.2)

где вектор наблюдения , матрица наблюдений , и - шум:

Предположим, что мы имеем априорный информационный массив , где - невырожденная априорная ковариация и оценки относятся к и следующим образом:

; (1.4.3)

где есть экстраполяционная оценка вектора состояния, а - матрица ковариации ошибки экстраполяционной оценки

Замечание: Для отфильтрованной оценки и ее матрицы ковариаций связь с информационным массивом аналогична

Предполагая, что мы можем записать следующее равенство:

(1.4.4)

Тогда модель наблюдений и предсказания выглядит следующим образом:

(1.4.5)

(1.4.6)

где , , - ортогональные матрицы такие, что матрицы и являются верхнетреугольными, а - остаточная ошибка, соответствующая решению наименьших квадратов

Заметим, что уравнения (1.4.5) и (1.4.6) эквивалентны уравнениям фильтра Калмана (доказательство данного факта см. в (2)). Но в отличие от традиционного фильтра Калмана, ККИФ позволяет избежать численной неустойчивости, являющейся результатом вычислительных погрешностей, поскольку вместо матриц ковариаций ошибки оценок на этапах экстраполяции и обработки измерений, по своей природе положительно определенных, ККИФ оперирует с их квадратными корнями. А это значит, что вычисления квадратного корня равносильно счету с двойной точностью для ковариации ошибок и кроме того устраняется опасность утраты матрицей ковариаций свойства положительно определенности. Недостатком данного метода является присутствие операций извлечения квадратного корня

 

2. Оценивание параметров по методу максимального правдоподобия с использованием квадратно-корневых информационных фильтров

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

2.1. Постановка задачи

Пусть динамическая система с дискретным временем дана в следующем виде:

(2.1.1)

где вектор состояния , матрица переходов из состояния в состояние , и - случайный процесс, представляющий собой белый гауссовый шум с нулевым средним и ковариацией , т.е. . Пусть множество наблюдений задается уравнением:

(2.1.2)

где вектор наблюдения , матрица наблюдений , и - шум:

Пусть также, матрицы и являются функциями неизвестного векторного параметра

Оценкой максимального правдоподобия является такое значение оцениваемых параметров , которое максимизирует вероятность события, при котором наблюдения, сгенерированные с подстановкой оцениваемых параметров, совпадают с действительными значениями наблюдений . Эта процедура эквивалентна минимизации обратного логарифма функции плотности условной вероятности, т.е. обратный логарифм функции правдоподобия представляется как:

(2.1.3)

где и - вычисляются согласно схеме фильтра Калмана следующим образом:

(2.1.4)

где есть ковариационная матрица ошибки экстраполяционной оценки

Запишем другие характеристики фильтра Калмана, которые нам понадобятся в дальнейшем:

Матрица Калмана:

(2.1.5)

Матрица ковариаций измененной по последним данным ошибки:

(2.1.6)

Невязка:

(2.1.7)

Измененная оценка:

(2.1.8)

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

(2.1.9)

где - оцениваемый векторный параметр; - индекс, определяющий номер итерации; - информационная матрица Фишера; - градиент обратного логарифма функции максимального правдоподобия

Стоит заметить, что итеративные алгоритмы, подобные (2.1.9), в среднем сходятся за меньшее число шагов, чем те алгоритмы, которые включают в себя только вычисления . С другой стороны, алгоритмы, содержащие и , требуют больше вычислений на каждом шаге

Модель наблюдений, в случае ККИФ, выглядит следующим образом:

(2.1.10)

где - ортогональная матрица такая, что - верхнетреугольная. Также, согласно (2.1.2) имеют вид:

(2.1.11)

где

(2.1.12)

тогда шум наблюдения имеет единичную ковариацию, что удовлетворяет ККИФ

Шаг предсказывания ККИФ, описывается следующим образом:

(2.1.13)

где матрица представляется уравнением (1.4.4), - ортогональная матрица такая, что матрица является верхнетреугольной

Информационным массивом ККИФ является массив данных . Он соотносится с оценкой состояния фильтра Калмана и матрицей ковариации ошибки оценивания следующими соотношениями:

(2.1.14)

(2.1.15)

2.2. Функция правдоподобия и ее представление терминах ККИФ

Для эффективного вычисления функции максимального правдоподобия при использовании ККИФ в фильтрации данных, необходимо выразить величины, входящие в выражение для , непосредственно через значения, которые вычисляются ККИФ-ом. Таким образом, две части (2.1.3):

(часть, зависящая от данных) (2.2.1)

(часть, зависящая от модели) (2.2.2)

выраженные через переменные, входящие в формулы ККИФ (2.1.10) и (2.1.13), приобретают следующий вид:

(2.2.3)

(2.2.4)

Доказательство (2.2.3) основано на следующем уравнении:

,

где - ортогональное преобразование такое, что матрица - верхнетреугольная. Находя нормы от обеих частей равенства, получим:

Уравнение (2.1.14) влечет за собой, следующее выражение:

Следовательно:

Далее, используя (2.1.8), чтобы переписать данное выражение в следующем виде:

После раскрытия скобок, опять используя уравнение (2.1.14), имеем:

Далее следует:

(2.2.5)

Наконец, используя (2.1.4), (2.1.5) и (2.1.15), получаем, что матрица, находящаяся внутри квадратных скобках выражения (2.2.5), является просто матрицей , что и доказывает (2.2.3)

Доказательство (2.2.4) основано на существовании ортогональных преобразований между определенными переменными ККИФ и остаточной ковариационной матрицей. Для упрощения записи положим, что

и

Используя обычную матричную алгебру и уравнения (2.1.4), (2.1.5), (2.1.6) и (2.1.15) можно показать, что выполняется следующее равенство:

Следовательно, если положить

(2.2.6)

получаем, что

- ортогональная матрица. Переписав (2.2.6) в виде

имеем

Из выражений для и получаем значения для детерминантов:

тогда имеем

Взяв логарифм от обеих частей выше записанного выражения и используя тот факт, что

получаем верное тождество для (2.2.4)

Простое преобразование формул, полученных выше, где шум наблюдения или измерения также является зависимым параметром, дает следующую формулу обратного логарифма функции правдоподобия в терминах ККИФ:

2.3. Градиент функции максимального правдоподобия

Для вычисления градиента , прежде всего, заметим, что градиент части, зависящей от модели (см. (2.2.2)), записывается следующим образом:

Так как матрицы и треугольные, а точнее верхнетреугольные, то только их диагональные элементы должны быть вычислены. Диагональные элементы первых трех матриц могут быть вычислены недорогим частичным обращением соответствующих величин ККИФ. Диагональные элементы последних двух матриц могут быть вычислены с использованием метода, описанного в разделе 2.4

Для градиента части, зависящей от данных, функции максимального правдоподобия, мы используем соотношение для изменения уравнений измерения ККИФ:

,

где - ортогональная матрица. Находя нормы от обеих частей равенства, получаем что:

Из последнего равенства имеем, что:

где значения может быть получено дифференцированием ККИФ, как показано в разделе 2.4. Значения матрицы получаются путем дифференцирования соотношения:

Таким образом:

(2.3.1)

Между тем, матрица - верхнетреугольная и должна равняться , где - верхнетреугольная часть матрицы на левой стороне (2.3.1) и - диагональная часть. И тогда, находится с помощью метода обратной подстановки решения треугольных систем

Подводя итог выше сказанному, имеем, что градиент обратного логарифма функции максимального правдоподобия приобретает вид:

,

где все входящие величины являются либо входными значениями КИИФ, либо легко находятся путем решения треугольных систем

Для выражения информационной матрицы Фишера в терминах ККИФ, вспомним, что - ый элемент матрицы Фишера записывается как:

Т.к. - случайный процесс с нулевым средним, то

(2.3.2)

где - -ая величина во временной последовательности, представляющей . Переписывая (2.3.2), используя ККИФ-форму представления , имеем, что - ый элемент матрицы Фишера приобретает вид:

Эта формула может быть использована и при замене ожидаемых значений переменных и вычисленными

2.4. Значения производных переменных ККИФ

Теперь дадим численно эффективный и точный метод для вычисления значений, , которые необходимы в формулах раздела 2.3

Для упрощения понимания положим, что преобразования ККИФ (2.1.10) и (2.1.13) имеют вид:

(2.4.1)

где - прямоугольная матрица, - ортогональная, которая при умножении с дает верхнюю трапециевидную матрицу . Элементы матрицы дифференцируемые функции параметра . Тогда, при заданных значениях производных , мы хотим определить матрицу . Явно прослеживается обобщение на случай ККИФ - преобразований и параметр заменяется вектором . Вдобавок, мы потребуем, чтобы и, следовательно, были квадратными и невырожденными

Для более полного представления ситуации, вначале остановимся на достаточно очевидном решении этой проблемы, однако, которое может вызвать определенные вычислительные трудности и, поэтому, не рекомендуется. Данный подход основан на использовании (2.4.1) и решении следующего уравнения:

,

далее прямое дифференцирование дает:

(2.4.2)

Используя тот факт, что и, следовательно, - верхнетреугольные, (2.4.2) может быть использовано для вычисления элементов . Для этого применяется алгоритм прямой подстановки, который является стандартным алгоритмом решения линейных систем. Вычислительные сложности могут возникнуть, когда матрица и, следовательно, плохо обусловлены (т.е. почти вырождены). Алгоритм прямой подстановки в этом случае может давать, мягко говоря, неточные результаты

Выбранный метод базируется на том наблюдении, что если матрица - ортогональная, то выполняется следующее равенство: . Дифференцируя это уравнение, находим, что

или (2.4.3)

Матрица , которая удовлетворяет соотношению - кососимметричная. Очевидно, что такая матрица имеет следующий вид: , где - строго нижнетреугольная. Поэтому (2.4.3) можно переписать в следующем виде:

(2.4.4)

для некоторой нижнетреугольной матрицы

Лемма: Нижнетреугольная матрица в соотношении (2.26), где удовлетворяет

(2.4.1), причем - невырождена, является нижнетреугольной частью матрицы , причем

(2.4.5),

где соответственно нижнетреугольная, диагональная и верхнетреугольная матрицы и удовлетворяет (2.4.4)

Доказательство:

Продифференцировав равенство (2.4.1), получаем:

(2.4.6)

из которого имеем:

Используя тот факт, что

получаем, что

(2.4.7)

Далее, так как и верхнетреугольные матрицы, то и и также верхнетреугольные. Поэтому, нижнетреугольная часть должна в точности соответствовать нижнетреугольной части с обратным знаком. Следовательно, если , то

Уравнение (2.4.7) дает метод для вычисления . Подстановкой (2.4.5) и (2.4.4) в (2.4.7) получаем, что

(2.4.8)

и, таким образом, имеем

Тем самым получаем путь нахождения :

    • вычислить ;
    • привести к виду ;
    • вычислить

и результатом будет

Уравнение (2.4.8) показывает всю опасность применения (2.4.6) непосредственно. Из соотношения (2.4.8) находим:

Первым слагаемым является , а вторым - . Записанное в таком виде соотношение (2.4.8), включает и в первое, и во второе слагаемое матрицу , но с разными знаками. Более того, матрица полная, то есть исключение происходит по всей матрице . Если значение матрицы небольшие, в сравнении с элементами матриц и , тогда сумма вычисляется достаточно точно. Однако если некоторые из элементов велики, тогда соответствующие элементы суммы будут неточны вследствие большого исключения элементов и результат в целом будет неточен

2.5. Описание алгоритма

Идеи, высказанные в разделе 2.4, могут быть использованы для создания очень небольшого и сжатого алгоритма вычисления обратного логарифма функции правдоподобия, его градиента и величин, входящих в выражение для информационной матрицы Фишера

Пусть определяет собой вектор параметров, по которым функция правдоподобия должна быть продифференцирована. Для вычисления , для каждого i, мы преобразуем уравнение изменения параметров наблюдающей модели следующим образом:

М1: Заменим (2.1.10) на

,

где первые два столбца повторяются для каждого

М2: Вычислить для каждого :

М3: Вычислить для каждого :

Заметим, что на шаге М2, матрица, которая должна быть обращена верхнетреугольная. Следовательно, результат обращения может быть получен с помощью метода обратной подстановки

Метод раздела 2.4, для нахождения величин , мы не можем применить напрямую к уравнению предсказания состояния (2.1.13), т.к. матрица, которая должна быть приведена к треугольному виду неквадратная, следовательно, необратимая. Однако этот алгоритм может быть применен при работе с подматрицами

Т1: Заменим (2.1.13) следующим соотношением:

где первые два столбца повторяются для каждого и , как в (2.1.13)

Т2: Вычислить для каждого :

здесь * - обозначение для первых q столбцов, не представляющих для нас интереса

Т3: Вычислить для каждого :

Соотношение для определяется дифференцированием уравнения:

и заменой производной , используя (2.4.4), значением вычисленным на шаге Т1. Заметим также, что , необходимое для этого последнего уравнения, легко вычисляется, т.к

где - верхнетреугольная

Но самое интересное, что значения , необходимые для вычисления матрицы Фишера, вычисляются автоматически на шаге М3

Отзыв

научного руководителя

на дипломную работу М.Ю. Кудрявцева

Тема работы ”Адаптивное параметрическое оценивание квадратно-корневым информационным алгоритмом” продиктована необходимостью проведения широкого спектра исследований по сравнительным оценкам различных подходов к проблеме идентификации моделей систем

Новизна и привлекательность рассматриваемого подхода обусловлены соединением в нем известного критерия максимума правдоподобия с алгоритмом квадратно-коневого информационного фильтра. Последний отличается высокой численной устойчивостью к погрешностям вычислений и к случаям плохой обусловленности схемы наблюдений. Однако факт совпадения оценок максимального правдоподобия с параметрами оптимального фильтра в общем случае не доказан. Этим объясняется актуальность экспериментальных исследований, устанавливающих условия, в которых данный критерий и соответствующий алгоритм идентификации могут считаться практически пригодными. М.Кудрявцев проанализировал заданный алгоритм, отличающийся сложным и разнообразным математическим аппаратом, сочетающим методы математической статистики, оптимального оценивания и матричных вычислений. Он обосновал вычислительные компактные схемы, включая вычисление обратного логарифма функции правдоподобия, его градиента по параметру неопределенности и информационной матрицы Фишера, и разработал необходимые программы

Основательность и настойчивость позволили М.Кудрявцеву выполнить эту работу без лишней торопливости, с глубоким пониманием существа вопроса и доведением всего исследования до наглядных экспериментальных результатов. Работа демонстрирует высокую подготовку ее автора по специальности ”Прикладная математика”, его способности к проведению самостоятельных исследований и заслуживает отличной оценки

 

  Д.т.н., профессор И.В. Семушин

Введение

Спектральный анализ - это один из методов обработки сигналов, который позволяет охарактеризовать частотный состав измеряемого сигнала. Преобразование Фурье является математической основой, которая связывает временной или пространственный сигнал (или же некоторую модель этого сигнала) с его представлением в частотной области

К обработке сигналов в реальном масштабе времени относятся задачи анализа аудио, речевых, мультимедийных сигналов, в которых помимо трудностей, связанных непосредственно с анализом спектрального содержания и дальнейшей классификацией последовательности отсчетов (как в задаче распознавания речи) или изменения формы спектра - фильтрации в частотной области (в основном относится к мультимедийным сигналам), возникает проблема управления потоком данных в современных вычислительных системах. Реальность накладывает отпечаток как на сами вычислительные алгоритмы, так и на результаты экспериментов, поднимая вопросы, с которыми не сталкиваешься при обработке всей доступной информации

При обработке сигналов обычно приходится решать задачи двух типов - задачу обнаружения и задачу оценивания. При обнаружении нужно дать ответ на вопрос, наблюдаем ли в данное время некоторый сигнал с априорно известными параметрами. Оценивание - это задача измерения значений параметров, описывающих сигнал (1)

Сигнал часто зашумлен, на него могут накладываться мешающие сигналы. Поэтому для упрощения указанных задач сигнал обычно разлагают по базисным составляющим пространства сигналов. Для многих приложений наибольший интерес представляют периодические сигналы. Вполне естественно, что используются Sin и Cos . Такое разложение можно выполнить с помощью классического преобразования Фурье

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

Постановка проблем, формулировка задач

На настоящее время существует большое количество алгоритмов и групп алгоритмов, которые так или иначе решают основную задачу спектрального анализа: оценивание спектральной плотности мощности, с тем чтобы по полученному результату судить о характере обрабатываемого сигнала .Основной вклад сделан такими исследователями как : Голд Б. (Gold B.), Рабинер Л. (Rabiner L.R.) , Бартлетт M. (Bartlett M.S.) Однако каждый из алгоритмов имеет свою область приложения. Например, градиентные адаптивные авторегрессионные методы не могут быть применены к обработке данных с быстро меняющимся во времени спектром. Классические методы имеют широкую область применения, но проигрывают авторегрессионным и методах, основанных на собственных значениях, по качеству оценивания. Но в реальном масштабе времени использование последних затруднено из-за вычислительной сложности

Более того, применение каждого из методов обычно требует выбора значений параметров (выбор окна данных и корреляционного окна в классических методах, порядка модели в авторегрессионном алгоритме и алгоритме линейного предсказания, предполагаемого числа собственных векторов в пространстве шума в метода Писаренко) и правильный выбор требует экспериментальных результатов с каждым классом алгоритмов

Таким образом, имеется следующая задача :

На основе существующих алгоритмов проанализировать возможность применения как к последовательной обработке сигналов в реальном времени, так и к блочной обработке и оценить качество получаемых результатов

Из вышесказанного можно сформулировать следующие подзадачи:

I. теоретическое и практическое исследование алгоритмов блочной обработки

II . анализ классических алгоритмов блочной обработки всей последовательности в части применения окон данных и корреляционных окон

  • анализ алгоритмов обработки сигналов в реальном масштабе времени
  • Существует несколько проблем, специфичных для обработки сигналов в реальном времени:

    · Необходимость в “одновременном” выполнении следующих основных этапов обработки данных:

  • Непосредственное получение последовательности входных данных (цифровые отсчеты аудио-сигнала, речевого сигнала).
  • Обработка получаемых отсчетов сигнала.
  • Представление обработанной информации
  • Возможность контролировать процесс обработки информации
  • · Ограничение длительности интервала выборки поступающих данных вычислительными ресурсами

    · Ограничение длительности интервала выборки характером сигнала

    Если первая проблема очевидна в рамках обработки данных в реальном времени, то вторая и третья проблемы требуют осмысления причин этих ограничений

    К сформулированным выше задачам добавляются :

  • задача построения схемы управления обработкой данных в реальном времени, основанной, в силу первой проблемы, на параллельных вычислениях и протоколах взаимодействия и синхронизации;
  • экспериментальный анализ по второй проблеме , то есть исследование влияния вычислительных ресурсов и методов оцифровки данных на максимально допустимую длину интервала выборки;
  • анализ длительности, исходя из характера сигнала .
  • Из постановки основной задачи вытекает необходимость в проведении большого количества экспериментов. Экспериментальные входные данные формируются следующим образом

    · для задачи анализа алгоритмов блочной обработки всей последовательности отсчетов формируются дискретизированные отсчеты данных тест-сигнала из суммы комплексных синусоид и аддитивных окрашенных шумовых процессов, сформированные посредством пропускания белого шума через фильтр с частотной характеристикой типа приподнятого косинуса или окна Хэмминга. Таким образом, в этом случае эксперимент определяется набором , где - последовательность комплексных синусоид с амплитудами дБ и частотами Гц , а - последовательность шумовых процессов с параметрами : центральная частота Гц. , динамический диапазон перекрываемых частот Гц. , мощность шума дБ

    · для анализа классических алгоритмов блочной обработки всей последовательности в части применения окон данных и корреляционных окон эксперимент и подсчет основных характеристик окон производится над дискретизированными отсчетами соответствующих функций

    · для анализа алгоритмов обработки сигналов в реальном масштабе времени данными являются аудио и речевой сигналы

    Выходными данными экспериментов являются :

    · для задачи анализа алгоритмов блочной обработки всей последовательности отсчетов :

    1.) оценка спектральной плотности мощности , полученная с помощью того или иного метода спектрального анализа, по которой можно судить о качестве применяемого метода, сравнивая истинную спектральную плотность мощности сформированного сигнала с полученной оценкой

    2.) вычислительные и временные затраты метода

    · для анализа окон данных и корреляционных окон - расчетные основные характеристики такие как : максимальный уровень боковых лепестков, эквивалентная ширина полосы, ширина полосы по уровню половинной мощности, степень корреляции и т.д

    · для анализа сигналов в реальном масштабе времени : спектральная плотность мощности (функция, зависящая в этом эксперименте также и от времени). Для оценки составляющих в спектре сигнала в данный момент времени

    Из ... Теоретический анализ существующих алгоритмов спектрального анализа

    Спектральная оценка, получаемая по конечной записи данных, характеризует некоторое предположение относительно той истинной спектральной функции, которая была бы получена, если бы в нашем распоряжении имелась запись данных бесконечной длины. Именно поэтому поведение и характеристики спектральных оценок должны описываться с помощью статистических терминов. Общепринятыми статистическими критериями качества оценки являются ее смещение и дисперсия

    Из формального определения спектра, следует, что спектр является некоторой функцией одних лишь статистик второго порядка, относительно которых в свою очередь предполагается, что они остаются неизменными, или стационарными во времени. Следовательно, такой спектр не передает полной статистической информации об анализируемом случайном процессе, а значит, дополнительная информация может содержаться в статистиках третьего и более высокого порядка. Кроме того, многие обычные сигналы, которые приходится анализировать на практике, не являются стационарными. Однако короткие сегменты данных, получаемые из более длинной записи данных, можно считать локально стационарными. Анализируя изменения спектральных оценок от одного такого сегмента к другому, можно затем составить представление и об изменяющихся во времени статистиках сигналов, то есть нестационарных

    Определение: Непрерывно-временным преобразованием Фурье называется функция

    Определение: Обратное преобразование Фурье определяется выражением

    Дискретное преобразование Фурье

    , где

    Спектральная Плотность Мощности

    Хотя выборочный спектр не является состоятельной оценкой истинной спектральной плотности мощности, эта оценка может быть использована если выполнять некоторого рода усреднение или сглаживания. На использовании этой оценки основан классический периодограммый метод определения спектральной плотности мощности

    Дальше идут методы ..... 1,2,3,4,.... вместе с экспериментальным анализом алгоритмов спектрального анализа

    Особенности реализации

    Для реализации алгоритмов был реализован язык проектирования алгоритмов, включающий в себя средства межзадачного обмена данными, то есть построение распределенных по процессам вычислительных алгоритмов, определенные части которого исполняются параллельно несколькими процессам. Дальнейшим развитием этого подхода является построение сетевых распределенных схем алгоритмов. Существует большое количество приложений этого подхода

    Заключение

    В данной работе :

    1. Tеоретически проанализированы методы спектрального анализа, а также возможность применения этих методов в современных вычислительных системах для обработки данных в реальном масштабе времени.
    2. Получены результаты поставленных экспериментов и на их основе выбран наиболее подходящий метод оценивания спектральной плотности мощности в аддитивной смеси комплексных синусоид и окрашенного стационарного шумового процесса.
    3. Дано описание и выполнена реализация схемы управления процессом обработки данных в реальном времени, использующая преимущества параллельной архитектуры вычислительных систем.
    4. Cформулирован ряд требований по вычислительным ресурсам при реальной обработке сделан анализ длины выборки данных при различном представлении входного сигнала.
    5. Получены результаты по эксперименту вычисления характеристик окон и на их основе выбрано оптимальное решение в каждом эксперименте по оцениванию спектральной плотности мощности тест-сигнала.

    Выводы

     

    Предзащита

    Введение

    Проблема идентификации линейной динамической системы заключается в создании модели процесса по его наблюдаемым входным и выходным сигналам в детерминистской или стохастической обстановке. Процесс идентификации включает в себя две независимые процедуры, а именно, структурную идентификацию и идентификацию параметров

    Когда неизвестны структура объекта и соответствующие физические законы, которым подчиняется его поведение, проводятся эксперименты, направленные на выявление структуры объекта и законов его поведения методами структурной идентификации. В случае, когда известна структура объекта (т.е. существует модель характеризующая его свойства), а неизвестными являются некоторые его характеристики, описываемые конечномерным вектором, последние определяются методами параметрической идентификации

    Постановка задачи

    Целью данной дипломной работы является исследование нового метода параметрической идентификации основанного на синтезе метода максимального правдоподобия и метода квадратно-корневого информационного фильтра (ККИФ), сравнение его с другими существующими алгоритмами с точки зрения вычислительной точности, быстродействия и сложности, а также реализация данного метода на ЭВМ

    Метод

    Как известно, оценкой максимального правдоподобия является значение оцениваемых параметров, которое максимизирует вероятность события, при котором наблюдения, сгенерированные с подстановкой оцениваемых параметров, совпадают с действительными значениями наблюдений. Вычисление оценки максимального правдоподобия может быть итеративно выполнено при помощи характеристического уравнения, которое включает в себя градиент обратного логарифма функции правдоподобия и информационную матрицу Фишера. Вычисления функции правдоподобия и информационной матрицы Фишера требуют применения фильтра Калмана (а также его производных для каждого параметра оценивания), который, как известно, не обладает достаточной устойчивостью. Бирман, занимавшийся построением численно устойчивых алгоритмов фильтрации, предложил для вычисления оценки максимального правдоподобия итеративным образом использовать квадратно-корневой информационный фильтр. В отличие от традиционного фильтра Калмана, ККИФ позволяет избежать численной неустойчивости, являющейся результатом вычислительных погрешностей, поскольку вместо ковариации ошибки оценок на этапах экстраполяции и обработки измерений, по своей природе положительно определенных, ККИФ оперирует с их квадратными корнями. Это значит, что вычисление квадратного корня равносильно счету с двойной точностью ковариации ошибок, кроме того устраняется опасность утраты матрицей ковариаций свойства положительно определенности. Недостатком данного метода является присутствие операций извлечения квадратного корня

    Таким образом, вычисление оценки максимального правдоподобия может быть осуществлено итеративно по следующей формуле:

    (1)

    где - конечномерный вектор оцениваемых параметров; - индекс, определяющий номер итерации; - информационная матрица Фишера; - градиент функции максимального правдоподобия

    Стоит заметить, что итеративные алгоритмы, подобные (1), в среднем сходятся за меньшее число шагов, чем те алгоритмы, которые включают в себя только вычисления . С другой стороны, алгоритмы, содержащие и , требуют больше вычислений на каждом шаге

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

    Заметим, что нахождение оценки максимального правдоподобия эквивалентно минимизации обратного логарифма функции правдоподобия, тогда критерием для метода является выражение:

    (2)

    где - невязка, - остаточная ковариация (т.е. ковариация невязок), подразумевается, что значения невязок в каждый момент времени независимы. Независимость же невязок обеспечивается при оптимальном фильтре, т.е. при точно известных значениях параметра . Из этого предположения следует, что начальные значения для параметра должны быть достаточно близкими к истинным его значениям

    Выводы

    Факт сходимости алгоритма максимального правдоподобия к оптимальным значениям параметров теоретически является недоказанным, поэтому в качестве основного метода исследования будем считать вычислительные эксперименты

    В рамках данного дипломного проекта были проведены следующие эксперименты:

    • Выявление зависимости точности оценивания от количества измерений.
    • Выявление зависимость точности оценивания от начальных условий для оцениваемых параметров.
    • Выявление зависимости времени оценивания от размерности задачи.
    • Проверка на сходимость метода с полностью наблюдаемой и ненаблюдаемой моделью системы.
    • Сравнение точности оценивания данного метода с другими существующими методами.
    • Сравнение времени оценивания данного метода с другими существующими методами.

    После проведения серии вычислительных экспериментов были получены следующие результаты:

    • Вышеописанный метод требует значительного количества времени для одной итерации по сравнению с другими методами параметрической идентификации, поскольку требуется вычисление градиента обратного логарифма функции правдоподобия и информационной матрицы Фишера. Данный факт показывает, что метод является достаточно сложным в вычислительном отношении.
    • Сходимость метода в значительной степени зависит от устойчивости матрицы перехода из состояния в состояние, от наблюдаемости динамической системы объекта, а также от количества оцениваемых параметров (наблюдаемость динамической системы является необходимым условием сходимости методов параметрической идентификации).
    • Метод критичен к начальным оценкам параметров.

    Заключение

    В данном дипломном проекте была проведена следующая работа:

    • Теоретически проанализирован алгоритм параметрической идентификации основанный на методе максимального правдоподобия с использованием квадратно-корневых информационных фильтров.
    • Данный метод программно реализован на ЭВМ.
    • Для проведения сравнительных экспериментов программно реализованы другие известные методы параметрической идентификации.
    • Поставлены эксперименты на выявление основных преимуществ и недостатков выше описанного метода по сравнению с другими реализованными методами.
    • Получены результаты поставленных экспериментов и на их основе сделаны выводы.

    Эксперименты

    Сходимость метода:

    ,

    , ,

    , ,

    1. Неизвестен один параметр

    начальные условия:

     

     

     

     

    

    Подобные работы:

    Актуально: