«АКТИВНАЯ ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ СТОХАСТИЧЕСКИХ ДИНАМИЧЕСКИХ СИСТЕМ НА ОСНОВЕ ПЛАНИРОВАНИЯ ЭКСПЕРИМЕНТА ...»

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ

УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ»

На правах рукописи

Чубич Владимир Михайлович

АКТИВНАЯ ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ

СТОХАСТИЧЕСКИХ ДИНАМИЧЕСКИХ СИСТЕМ НА ОСНОВЕ

ПЛАНИРОВАНИЯ ЭКСПЕРИМЕНТА

05.13.17 – Теоретические основы информатики Диссертация на соискание ученой степени доктора технических наук

Научный консультант Денисов В.И.

Новосибирск –

СОДЕРЖАНИЕ

Стр.

ОСНОВНЫЕ ОБОЗНАЧЕНИЯ И СОКРАЩЕНИЯ…...……….……….... ВВЕДЕНИЕ……………………………………………………………...……… 1 Проблема активной параметрической идентификации стохастических динамических систем и задачи диссертационного исследования……………………………………………………………………………. Теоретические и методологические основы активной параметрической идентификации………………...……………………………... 1.1.1 Процедура активной идентификации………………………... 1.1.2 Оценивание неизвестных параметров………………………... 1.1.3 Исходные понятия и результаты теории оптимального эксперимента………………………………………………………. 1.1.4 Прямая градиентная процедура синтеза непрерывных оптимальных планов…………………………………………….. 1.1.5 Двойственная градиентная процедура синтеза непрерывных оптимальных планов…………...…………………………........ 1.1.6 Построение дискретных оптимальных планов……………… 1.1.7 Схема процедуры активной параметрической идентификации систем с предварительно выбранной модельной структурой………………………………...……………………….…. Анализ современного состояния проблемы активной параметрической идентификации стохастических динамических систем…… Структурно-вероятностное описание моделей……………………...

1.3 1.3.1 Модели дискретных систем…………………………………... 1.3.2 Модели непрерывно-дискретных систем……………………. Цель и задачи исследования………………………………………… 1.4 Выводы………………………………………………………………… 1.5 2 Оценивание параметров моделей стохастических динамических систем………………………………….…………………………………….. Оценивание параметров моделей дискретных систем………..…… 2.1 Критерий максимального правдоподобия и алгоритм 2.1. вычисления его значения для линейных нестационарных моделей……………………………………………………….. Критерий максимального правдоподобия и алгоритм вычисления его значения для линеаризованных моделей……. Алгоритм вычисления градиента критерия максимального 2.1. правдоподобия для линейных нестационарных моделей…. Алгоритм вычисления градиента критерия максимального 2.1. правдоподобия для линеаризованных моделей……………. Оценивание параметров моделей непрерывно-дискретных систем 2.2 Особенности алгоритмов вычисления значений критериев 2.2. максимального правдоподобия для линейных нестационарных и линеаризованных моделей

Особенности алгоритмов вычисления градиентов критериев максимального правдоподобия для линейных нестационарных и линеаризованных моделей …………………... Выводы………………………………………………………………… 2.3 3 Планирование эксперимента для моделей стохастических дискретных систем…………………………………………………………………... Вычисление информационной матрицы Фишера…………………...

3.1 3.1.1 Вывод информационной матрицы Фишера для линейных нестационарных моделей……………………………………... 3.1.2 Алгоритм вычисления информационной матрицы Фишера для линейных нестационарных моделей…………………….. 3.1.3 Вычисление информационной матрицы Фишера для линейных нестационарных моделей, полученных в результате Планирование входных сигналов…………………………………….

Фишера по компонентам входного сигнала для линейных Алгоритм вычисления производных информационной матрицы Фишера по компонентам входного сигнала для Вычисление производных информационной матрицы Фишера по компонентам входного сигнала для моделей, полученных в результате временной линеаризации………….. Фишера по компонентам входного сигнала для моделей, полученных в результате статистической линеаризации…. Алгоритм вычисления производных информационной матрицы Фишера по компонентам входного сигнала для моделей, полученных в результате статистической линеаризации………………………………………………………... Планирование эксперимента как задача дискретного оптимального управления……………………………………… Планирование эксперимента в установившемся режиме Планирование начальных условий…………………………………... 3. Фишера по компонентам вектора начальных условий для Алгоритм вычисления производных информационной матрицы Фишера по компонентам вектора начальных условий для линейных нестационарных моделей………….... Планирование начальных условий на примере модели процесса изменения температуры в двухкомнатной квартире………………………………………………………........ Выводы………………………………………………...……………… 4 Планирование входных сигналов для моделей стохастических непрерывно-дискретных систем…………...……………………………….. Вычисление информационной матрицы Фишера…………………..

Вывод информационной матрицы Фишера для линейных нестационарных моделей……………………………………. Алгоритм вычисления информационной матрицы Фишера для линейных нестационарных моделей…………………… Вычисление информационной матрицы Фишера для линейных нестационарных моделей, полученных в результате линеаризации……………………………………………… Вычисление производных информационной матрицы Фишера по Дифференцирование информационной матрицы Фишера по компонентам входного сигнала для линейных нестационарных моделей………………...……………………….. Алгоритм вычисления производных информационной матрицы Фишера по компонентам входного сигнала для Вычисление производных информационной матрицы Фишера по компонентам входного сигнала для моделей, полученных в результате временной линеаризации……… Выводы………………………………………………………………...

5 Описание программных комплексов активной параметрической идентификации стохастических динамических систем на основе планирования входных сигналов………………………………………... Назначение и общие сведения о программных комплексах……….

Характеристика возможностей и организация программных комплексов………………………………………………………………... Описание интерфейса программных комплексов………………….. Выводы………………………………………………………………...

6 Примеры активной параметрической идентификации стохастических динамических систем……………………………………………….. Активная параметрическая идентификация дискретных систем….

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

ЗАКЛЮЧЕНИЕ………………………………………………………………... ПРИЛОЖЕНИЕ А. Справки о внедрении результатов диссертационной работы…………………….………………………………………………..

ОСНОВНЫЕ ОБОЗНАЧЕНИЯ И СОКРАЩЕНИЯ

- вектор неизвестных параметров размерности s;

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

- вектор истинных значений параметров;

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

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

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

- непрерывный нормированный план эксперимента - оптимальный по некоторому критерию непрерывный нормированный план эксперимента;

- информационная матрица плана;

- критерий оптимальности;

- информационная матрица Фишера одноточечного плана - r-мерный вектор управления (входа) в соответствующий - оценка одношагового прогнозирования состояния x(t k 1), x tk 1 | tk xij (t k 1 | t k 1) - оценка фильтрации состояния x(t k 1), соответствующая - p-мерный вектор шума системы в соответствующий момент - m-мерный вектор измерения (выхода) в момент времени y tk 1 | tk j-я реализация выходного сигнала, соответствующая входYij - относительная ошибка оценивания в пространстве откликов - относительная ошибка оценивания в пространстве откликов L( ;Y1 ) - m-мерный вектор шума измерения в момент времени t k 1 ;

v tk - m-мерный вектор обновления в момент времени t k 1 ;

- дискретный нормированный план эксперимента - область допустимых входных сигналов;

- область допустимых начальных условий;

- градиент скалярной функции g по аргументу - оператор математического ожидания;

- вещественное линейное пространство, состоящее из nRn - множество вещественных матриц, содержащих m строк и n - матрица, транспонированная к матрице А;

- матрица, обратная к невырожденной матрице А;

SpA - определитель матрицы А;

det A - дельта-функция Дирака;

ДУ - дифференциальное уравнение;

ИМФ - информационная матрица Фишера;

МАВ - максимум апостериорной вероятности;

ММП - метод максимального правдоподобия;

МНК - метод наименьших квадратов;

ОМП - оценка максимального правдоподобия;

ПК-I - программный комплекс активной параметрической идентификации стохастических нелинейных дискретных систем;

ПК-II - программный комплекс активной параметрической идентификации стохастических нелинейных непрерывнодискретных систем;

ПК-III - программный комплекс активной параметрической идентификации стохастических нестационарных линейных дискретных систем;

ВВЕДЕНИЕ

Актуальность исследования.

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

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

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

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

Степень разработанности проблемы.

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

Анализ современного состояния проблемы активной параметрической идентификации стохастических динамических систем показал, что наиболее значительный прогресс в ее решении достигнут применительно к линейным стационарным моделям и к моделям (в общем случае нелинейным) с детерминированными уравнениями состояний. Этому способствовали, в частности, труды таких признанных специалистов, как А.Ж. Абденов, Ю.П. Адлер, В.Г. Горский, В.И. Денисов, Э.К. Лецкий, В.Н. Овчаренко, А.А. Попов, А.М. Талалай в нашей стране и Г. Гудвин, М. Зейроп, Л. Льюнг, Р. Мехра, Р. Пейн, Л. Пронзато, Э. Уолтер за рубежом. Указанная проблема не рассматривалась для стохастических линейных нестационарных и нелинейных моделей с вхождением неизвестных параметров в уравнения состояния и измерения, в начальные условия и в ковариационные матрицы шумов системы и измерений.

В данной диссертационной работе решается проблема активной параметрической идентификации преимущественно для таких моделей.

Предмет исследования.

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

Цель и задачи исследования.

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

Для достижения поставленной цели в диссертационной работе ставятся и решаются следующие основные задачи:

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

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

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

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

Теоретическая и методологическая база исследования.

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

Научная новизна. Получены следующие новые результаты, которые выносятся на защиту:

1. Впервые выражения ИМФ для гауссовских линейных нестационарных и линеаризованных дискретных и непрерывно-дискретных моделей с параметрами в уравнениях состояния и измерения, в начальных условиях и в ковариационных матрицах шумов системы и измерений.

2. Разработаны алгоритмы вычисления ИМФ для гауссовских линейных нестационарных и линеаризованных дискретных и непрерывно-дискретных моделей.

3. Разработаны алгоритмы вычисления производных ИМФ по компонентам входного сигнала для линейных нестационарных дискретных и непрерывно-дискретных моделей, дискретных моделей, полученных в результате временной или статистической линеаризации и непрерывно - дискретных моделей, 4. Разработаны прямые и двойственные градиентные процедуры синтеза A- и D- оптимальных входных сигналов для гауссовских линейных нестационарных и линеаризованных дискретных и непрерывно-дискретных моделей с параметрами в уравнениях состояния и измерения, в начальных условиях и в ковариационных матрицах шумов системы и измерений.

5. Разработан алгоритм вычисления производных ИМФ по компонентам вектора начальных условий для линейных нестационарных дискретных моделей.

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

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

8. Разработаны и программно реализованы прямая и двойственная процедуры синтеза A- и D- оптимальных входных сигналов для установившегося режима гауссовских линейных стационарных дискретных моделей с неизвестными параметрами в уравнениях состояния и измерения, в начальных условиях и 9. Разработаны снабженные пользовательским интерфейсом программные комплексы ПК-I и ПК-II, предназначенные для активной параметрической идентификации стохастических нелинейных дискретных и непрерывнодискретных систем соответственно.

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

Программный комплекс ПК-I создан совместно с доцентом кафедры прикладной математики Новосибирского государственного технического университета О.С. Черниковой. При этом автором разработаны программы вычисления ИМФ и их производных по компонентам входного сигнала, программы нахождения значений критериев максимального правдоподобия и программы построения А- и D-оптимальных входных сигналов. О.С. Черникова разработала программы вычисления градиентов в процедурах оценивания параметров и синтеза А- и D- оптимальных входных сигналов.

Программный комплекс ПК-II создан совместно с Е.В. Филипповой. В нем автором разработаны программа вычисления ИМФ, программа нахождения значения критерия максимального правдоподобия и программы построения Аи D-оптимальных входных сигналов. Е.В. Филиппова разработала программу вычисления производных ИМФ по компонентам входного сигнала и программы вычисления градиентов в процедурах планирования А- и D-оптимальных входных сигналов.

Проектирование и реализация интерфейса к программным комплексам ПК-I и ПК-II осуществлялись совместно с О.С. Черниковой и Е.В. Филипповой.

Практическая ценность и реализация результатов исследования.

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

Разработанные процедуры и алгоритмы реализованы в программных комплексах активной параметрической идентификации стохастических нелинейных дискретных (ПК-I) и непрерывно-дискретных (ПК-II) систем (Свидетельство о государственной регистрации программы для ЭВМ № 2011612716. – М.: Федеральная служба по интеллектуальной собственности (Роспатент). Свидетельство о государственной регистрации программы для ЭВМ № 2011612718. – М.: Федеральная служба по интеллектуальной собственности (Роспатент). - 2011), в программном комплексе (ПК-III) активной параметрической идентификации стохастических нестационарных линейных дискретных систем (Свидетельство о государственной регистрации программы для ЭВМ № 2012612281. – М.: Федеральная служба по интеллектуальной собственности (Роспатент). - 2012).

Диссертационная работа выполнялась в рамках тематических планов НИР НГТУ по заданию Министерства образования и науки Российской Федерации на 2002-2005 гг. «Математическое моделирование многофакторных объектов на основе наблюдений» (№ 1.1.02), на 2006-2008 гг. «Моделирование статических и динамических многофакторных объектов стохастической природы и исследование вероятностных закономерностей» (№ 1.1.06), на 2009-2010 гг. «Методы и технологии моделирования и планирования экспериментов для исследования сложных многофакторных объектов» (№ 1.1.09), а также являлась частью исследований по ведомственной целевой программе «Развитие научного потенциала высшей школы» в 2006-2008 гг. (код проекта РНП.2.1.2.43) и федеральной целевой программой «Интеграция науки и высшего образования на 2002гг.» (код проекта Б0097/1376).

Проведение диссертационных исследований было поддержано грантами Федерального агентства по образованию (государственный контракт от «18»

ноября 2009 г. № П2365, научный руководитель Чубич В.М.) и Министерства образования и науки Российской Федерации (государственный контракт от «05» октября 2010 г. № 14.740.11.0587, научный руководитель Чубич В.М.) в рамках федеральной целевой программы «Научные и научно-педагогические кадры инновационной России на 2009-2013 гг.»

Соответствие диссертации паспорту научной специальности.

Содержание диссертации соответствует п.5 области исследований «Разработка и исследование моделей и алгоритмов анализа данных, обнаружения закономерностей в данных и их извлечениях, разработка и исследование методов и алгоритмов анализа текста, устной речи и изображений» паспорта специальности научных работников 05.13.17 – «Теоретические основы информатики» по техническим наукам.

Апробация работы.

Основные положения и результаты диссертационной работы докладывались и обсуждались на следующих конференциях и симпозиумах: Российская научно-техническая конференция «Информатика и проблемы телекоммуникаций» (г. Новосибирск, 1994 г.); Международные конференции «Актуальные проблемы электронного приборостроения» (г. Новосибирск, 1994 г., 2010 г.);

Международные научно-технические конференции «Информатика и проблемы телекоммуникаций» (г. Новосибирск, 1995 г., 1997г.); Российско-Корейские Международные Симпозиумы «Наука и технологии» KORUS’2003, KORUS’2004, KORUS’2005 (г. Ульсан, Корея, 2003 г., г.Томск, 2004 г., г. Новосибирск, 2005 г.); Всероссийская научно-практическая конференция «Перспективы развития информационных технологий» (г. Новосибирск, 2010 г.); Международная научно-техническая конференция «Системный анализ и информационные технологии» SAIT’2010 (г. Киев, Украина, 2010г.); Международная конференция IASTED по автоматике, управлению и информационным технологиям «Управление, диагностика и автоматика» ACIT-CDA’2010 (г. Новосибирск, 2010 г.), а также на научных сессиях факультета прикладной математики и информатики ФГБОУ ВПО «Новосибирский государственный технический университет».

Публикации.

Всего по результатам выполненных исследований опубликованы 42 работы, в том числе монография, 20 статей в журналах из Перечня ВАК ведущих рецензируемых научных изданий для опубликования основных результатов диссертаций на соискание учёной степени доктора и кандидата наук, 4 статьи в других журналах и сборниках научных трудов, 12 публикаций в материалах и сборниках трудов Международных и Российских конференций, 3 свидетельства о государственной регистрации программы для ЭВМ.

Структура и объем работы.

Диссертационная работа состоит из списка основных обозначений и сокращений, введения, шести разделов, заключения, списка использованных источников из 212 наименований и приложения. Общий объем работы составляет 247 страниц, включая 34 рисунка и 7 таблиц.

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

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

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

Третий и четвертый разделы посвящены теоретическим и прикладным аспектам планирования эксперимента для моделей стохастических дискретных и непрерывно-дискретных систем соответственно.

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

В шестом разделе рассмотрены примеры активной параметрической идентификации стохастических динамических систем.

Автор считает своим приятным долгом выразить глубокую благодарность д.т.н., профессору, заслуженному деятелю науки Российской Федерации В.И.Денисову за постоянную поддержку и внимание к работе.

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

Различают задачи идентификации в узком и широком смыслах [1]. В соответствии с [2], можно выделить следующие этапы решения задач идентификации в широком смысле:

- определение класса и характеристик входных воздействий;

- выбор критерия близости модели к системе;

- выбор и определение структуры модели на основании имеющейся априорной информации об исследуемом процессе и определенных эвристических соображений;

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

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

По способу проведения эксперимента, существующие методы идентификации можно разделить на пассивные и активные. При пассивной идентификации для построения математической модели используются реально действующие в системе сигналы и нормальный режим эксплуатации не нарушается. Методы пассивной идентификации достаточно полно описаны, например, в [4-15].

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

Эффект достигается в результате сочетания традиционных приемов параметрического оценивания с концепцией планирования эксперимента (см. например, [19-28]).

Процедура активной идентификации (см. рисунок 1.1) предполагает выполнение следующих основных этапов:

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

2. Определение структуры математической модели Это – весьма важная и ответственная часть процедуры идентификации. Определение общей структуры модели и класса уравнений, которыми предполагается описывать наблюдаемый процесс, является задачей структурной идентификации. Подчеркнем, что выбранная модельная структура должна быть идентифицируемой, т.е. позволять по имеющимся экспериментальным данным однозначно находить оценки неизвестных параметров. Исследованию идентифицируемости детерминированных линейных стационарных моделей в пространстве состояний посвящена монография [29].

Математические модели динамических систем можно классифицировать по следующим признакам [2]: одномерность – многомерность, линейность – нелинейность, стационарность – нестационарность, непрерывность – дискретность, детерминированность – стохастичность, характер возмущений, форма описания и т.д.

Описание динамической системы в пространстве состояний [1,30-32] позволяет учесть имеющиеся физические представления о механизмах работы системы. В отличие от моделей передаточных функций, которые используются при описании моделей линейных стационарных систем, методы пространства состояний позволяют создать компактную форму представления любых систем:

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

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

4. Планирование эксперимента При построении моделей динамических систем могут использоваться различные способы управления экспериментом. В простейшем случае управление экспериментом сводится к выбору оптимальных моментов измерений, в более сложном – к планированию оптимальных входных сигналов и начальных условий [24,33,34]. Возможны также варианты смешанных схем.

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

5. Проверка адекватности модели В результате выполнения предыдущих этапов процедуры активной идентификации получается определенная модель из некоторого класса, которая в соответствии с выбранным критерием качества наилучшим образом воспроизводит экспериментальные данные. Далее необходимо проверить соответствие модели данным наблюдений, априорной информации и поставленной прикладной цели. Проверка адекватности построенной модели осуществляется подачей на ее вход и вход самой системы тестирующих сигналов, которые не использовались для идентификации, и сравнением спрогнозированных по модели и реальных выходных данных (например, по сумме квадратов норм невязок). В случае необходимости осуществления оптимального выбора из нескольких конкурирующих модельных структур следует воспользоваться такими информационными критериями [35,36], как критерий Акаике, критерий Такеучи, критерий Шварца (байесовский критерий).

Предметом исследования диссертационной работы является проблема активной параметрической идентификации стохастических динамических систем с предварительно выбранной модельной структурой (структурно-вероятностное описание рассматриваемых моделей выполнено в п. 1.3). В этом случае процедура активной параметрической идентификации предполагает выполнение следующих этапов [37-45]:

а) вычисление оценок неизвестных параметров по измерительным данным, соответствующим некоторому плану эксперимента;

б) синтез на основе полученных оценок оптимального плана эксперимента;

в) пересчет оценок параметров по измерительным данным, соответствующим синтезированному плану.

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

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

Критерий идентификации формируется в соответствии с выбранным методом статистического оценивания. Здесь имеет смысл, прежде всего, выделить метод наименьших квадратов [6,7,10-15,46-62], не требующий знания закона распределения измерительных данных; метод максимального правдоподобия [4,6-11,13,46-57,60-69], использующий знание закона распределения выборочных данных; метод максимума апостериорной вероятности или байесовский подход к оцениванию [4,6,9-11,13,47,49-51,53,56,64-70], предполагающий случайность оцениваемых параметров и знание законов распределения оцениваемых параметров и измерительных данных.

Структурно-вероятностное описание рассматриваемых в диссертационной работе моделей и детерминированная природа подлежащих оцениванию неизвестных параметров обусловили выбор в качестве метода статистического оценивания метод максимального правдоподобия. Известно, что при выполнении некоторых общих условий (условий регулярности), накладываемых на функцию правдоподобия, оценки максимального правдоподобия (ОМП) обладают такими важными для практики асимптотическими свойствами как асимптотическая несмещенность, состоятельность, асимптотическая эффективность и асимптотическая нормальность [46-48,57,62-65,67-69].

полученная по выборке Y1, называется асимптотически несмещенной, если оценки (ОМП не всегда являются несмещенными оценками).

Здесь и далее обозначает векторную или матричную норму в зависимости от контекста.

Если N - оценка вектора параметров со смещением b N, то справедливо следующее неравенство информации (неравенство Рао – Крамера Фреше) для нижней границы ковариационной матрицы N [47,66,67-69]:

Здесь M 1 - дисперсионная матрица, являющаяся обратной к информационной матрице Фишера. Информационная матрица Фишера (ИМФ) определяется равенствами [6,11,47,48,67-69] где L( ;Y1 ) - функция правдоподобия (плотность совместного распределения Таким образом, в асимптотике ОМП обладают наименьшими ковариационными матрицами, совпадающими с M 1.

Оценка N вектора параметров называется асимптотически нормальN N нормальному распределению с нулевым математическим ожиданием и дисперM При заданном критерии идентификации задача нахождения оценок неизвестных параметров заключается в решении задачи нелинейного программирования с ограничениями:

где Методы нелинейного программирования [71-83] включают в себя большую группу численных методов, многие из которых ориентированы на решение оптимизационных задач соответствующего класса. Выбор того или иного метода обусловлен сложностью вычисления критерия идентификации, необходимой точностью решения, мощностью компьютера и т.д.

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

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

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

Эффективными методами нулевого порядка считаются метод вращающихся координат или метод Розенброка [71,72,75,80], метод деформируемого многогранника или метод Нелдера – Мида [71,75,80,83], метод сопряженных направлений или метод Пауэлла [71,77,80], метод конфигураций или метод Хука – Дживса [71,72,80,83] и разнообразные методы случайного поиска [71,73, 78,80].

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

Методы первого порядка используют кроме значений целевой функций значения ее производных. Наиболее эффективными среди методов первого порядка являются метод сопряженных градиентов или метод Флетчера-Ривса [71и методы переменной метрики или квазиньютоновские методы (в последних используются аппроксимации матрицы вторых производных или обратной к ней). Хронологически первым квазиньютоновским методом является метод Дэвидона-Флетчера-Пауэлла [71-83], но в настоящее время, по мнению многих исследователей, лучшим квазиньютоновским методом является метод Бройдена-Флетчера-Голдфарба-Шенно [73,74,76,77,79,81-83]. Методы первого порядка целесообразно применять в тех случаях, когда есть возможность вычислять производные минимизируемой функции.

Методы второго порядка дополнительно требуют вычисления вторых производных и быстро сходятся для выпуклых целевых функций. К этой группе относятся метод Ньютона и его многочисленные модификации [71-83], среди которых заслуживают внимание методы доверительной области [74,79,82-84], оказавшиеся весьма эффективными для задач высокой размерности, насчитывающих многие сотни и тысячи переменных.

Другой общепринятый подход к решению оптимизационной задачи (1.2) предполагает либо непосредственный (в этом случае мы приходим к методам возможных направлений, из которых наиболее популярным является метод проекции градиента [71-75,77,78,80,82]), либо опосредованный (через функцию Лагранжа) учет ограничений (1.3). Метод проекции градиента представляет практический интерес лишь в тех случаях, когда проектирование выполняется легко, что определяется простотой устройства допустимого множества (например, представляет собой шар, координатный параллелепипед, задается системой ограничений в виде линейных равенств и неравенств и т.д.). Методы второй группы, самым передовым и эффективным из которых, по мнению многих современных исследователей, является метод последовательного квадратичного программирования (общепринятая аббревиатура – SQP, от английского Sequential quadratic programming), в последнее время развиваются наиболее интенсивно. Этот метод [79,81-83] реализован в рамках пакета Optimization Toolbox [84,85] программной системы MATLAB и используется в диссертационной работе для численного нахождения оценок неизвестных параметров, а также для синтеза непрерывных оптимальных планов см. п. 1.1.4 и п. 1.1.5).

Метод последовательного квадратичного программирования позволяет организовать решение оптимизационной задачи (1.2) с ограничениями (1.3) в соответствии со следующей схемой:

Шаг 2. Вычислить где в противном случае вычислить отвечающие найденному значению g k множители Лагранжа k и k.

Шаг 3. Положить Шаг 4. Увеличить k на единицу и перейти на шаг 2.

На шаге 2 решается задача квадратичного программирования (1.4), (1.5), пользование в целевой функции матрицы H k (матрицы вторых производных для функции Лагранжа L задачи (1.2), (1.3)) позволяет учесть важную информацию о членах второго порядка в аппроксимации ограничений и приводит к высокой скорости сходимости рассматриваемой процедуры в целом. В пакете Optimization Toolbox для решения задачи (1.4), (1.5) используется метод проекции градиента с аппроксимацией матрицы вторых производных H k по формуле Бройдена-Флетчера-Голдфарба-Шенно.

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

Весьма полезны при численном нахождении оценок неизвестных параметров методы глобальной оптимизации [86,87].

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

Изложим некоторые основополагающие понятия и результаты теории планирования эксперимента.

и соответствующих им долей повторных запусков (весов):

Множество планирования определяется ограничениями на условия проведения эксперимента.

Под непрерывным нормированным планом эксперимента условимся понимать совокупность величин В отличие от дискретного нормированного плана в непрерывном нормированном плане снимается требование рациональности весов pi.

В общем случае непрерывный нормированный план соответствует вероятностной мере (d ), заданной на области и удовлетворяющей условиям неотрицательности и нормировки:

При этом нормированная информационная матрица ( ) плана определяется соотношением Для плана (1.7) интеграл в (1.8) переходит в сумму, т.е.

где M( i ) – информационные матрицы точек спектра плана, вычисляющиеся для каждого значения i в соответствии с равенствами (1.1).

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

Теорема 1.1 [34,88-94]. Нормированные информационные матрицы обладают следующими свойствами:

1. Для любого непрерывного нормированного плана информационная матрица M - вещественная, симметричная, неотрицательно - определенная матрица порядка s s (s – количество неизвестных параметров).

2. Множество матриц M, соответствующее всем возможным нормированным планам (1.7), является выпуклым. Если множество планирования замкнуто, то и множество информационных матриц замкнуто.

3. Для любого непрерывного нормированного плана всегда найдется дисs(s 1) точек и нормированная информационная матрица M( ) которого совпадает с информационной матрицей M плана.

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

Подчеркнем, что в [88,89,91,93,94] приведено универсальное доказательство теоремы 1.1, применимое как к статическим, так и к динамическим моделям.

Выбранный определенным образом план эксперимента позволяет повысить точность оценивания неизвестных параметров. Будем судить о качестве плана по значению некоторого функционала от информационной матрицы или соответствующей ей дисперсионной матрицы D. Перечислим наиM более важные критерии, отражающие точность оценивания неизвестных параметров [19-25].

План * называется А – оптимальным, если его дисперсионная матрица имеет наименьший след, т.е.

А – оптимальный план позволяет найти оценки неизвестных параметров с минимальной средней дисперсией. При этом эллипсоид рассеивания имеет минимальную сумму квадратов длин осей и наименьшую длину диагоналей параллелепипеда, описанного около этого эллипсоида.

План * называется D – оптимальным, если Для D – оптимального плана объем эллипсоида рассеивания оценок параметров наименьший.

План * называется Е – оптимальным, если он минимизирует (максимизирует) максимальное (минимальное) собственное значение дисперсионной (информационной) матрицы, т.е.

где i - собственное значение матрицы D( ) или M( ). Е – оптимальный план минимизирует длину максимальной оси эллипсоида рассеивания оценок параметров.

При активной параметрической идентификации динамических систем в качестве критериев оптимальности используются также (см., например, [95,96]) след информационной матрицы SpM или след взвешенной информационной матрицы Sp WM.

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

Теорема 1.2 [34,88-94,97,98]. Следующие утверждения:

1. План минимизирует X M( ) ;

эквивалентны между собой. Информационные матрицы планов, удовлетворяющих условиям 1 – 3, совпадают. Любая линейная комбинация планов, удовлетворяющих 1 –3, также удовлетворяет 1-3.

представлены в таблице 1.1.

Таблица 1.1 - Соответствие значений параметров обобщенной теоремы эквивалентности критериям оптимальности А – оптимальности Подчеркнем, что в работах [88,89,91,93,94,98] приведено универсальное доказательство теоремы 1.2, применимое как к статическим, так и к динамическим моделям.

Укажем два следствия, удобные для проверки планов на А - и D - оптимальность.

ет своего максимального значения SpM 1( *).

В точках D – оптимального плана * функция (, *) достигает своего максимального значения s.

Теорема 1.2 дает еще одну статистическую интерпретацию D – оптимальному планированию экспериментов, утверждая эквивалентность D – и G – оптимальных планов.

Критерий G – оптимальности относится к группе критериев, направленных на повышение точности прогнозируемых по модели выходных данных [19-25].

План * называется G – оптимальным (минимаксным), если он удовлетворяет условию G – оптимальный план обеспечивает не слишком высокую обобщенную дисперсию ошибки прогнозирования за счет минимизации максимальной по 1.1.4 Прямая градиентная процедура синтеза непрерывных Оптимизационную задачу можно решать непосредственно (напрямую) с помощью общих методов численного поиска экстремума [71-83], в том числе и с применением методов глобальной оптимизации [86,87]. Характерной особенностью этого подхода является большая размерность экстремальной задачи. Для критериев А- и D- оптимальности получается задача выпуклого программирования, возможные варианты решения которой представлены, например, в [21,91,99-103]. Остановимся на следующем варианте прямой градиентной процедуры построения непрерывных оптимальных планов из [39,40,44,94,98,104-106] :

Шаг 1. Зададим начальный невырожденный план в котором q 1, s – количество неизвестных параметров. Вычислим информационные матрицы M одноточечных планов для i 1,2,...,q и по формуле (1.9) информационную матрицу всего плана 0. Положим k 0.

Шаг 2. Считая веса p1,pk,...,pq фиксированными, для задачи Соответствие значений параметров X M( ), прямой процедуры критериям Аи D- оптимальности такое же, как в таблице 1.1.

выполним одну итерацию метода последовательного квадратичного программирования (см. п. 1.1.2). Составим план где i 1 - точки, найденные на шаге 2, и вычислим M Шаг 3. Зафиксируем точки спектра полученного плана и для задачи выполним одну итерацию метода последовательного квадратичного программирования. Составим план Шаг 4. Если для малого положительного числа выполняется неравенство перейдём на шаг 5, в противном случае для k k 1 повторим шаги 2 и 3.

Шаг 5. Проверим необходимое условие оптимальности плана Если требуемое условие оптимальности выполняется (см. следствия из теоремы 1.2), закончим процесс. В противном случае повторим всё сначала, скорректировав начальный план 0.

Приведённый алгоритм требует вычисления градиентов Начнем с критерия А – оптимальности (1.10). Для него с учетом формулы (1.9) имеем:

Перейдем теперь к критерию D- оптимальности (1.11). В этом случае Таким образом (см. соотношения (1.15)-(1.18)), применение прямой градиентной процедуры синтеза непрерывных А – или D - оптимальных планов предполагает разработку алгоритмов вычисления ИМФ и их производных по компонентам точек спектра плана эксперимента (1.7).

1.1.5 Двойственная градиентная процедура синтеза непрерывных Другой подход (его называют двойственным) к решению оптимизационной задачи (1.12) основан на обобщенной теореме эквивалентности. В этом случае рассматриваемая задача уже не является задачей выпуклого программирования, но размерность пространства варьируемых параметров может оказаться существенно меньше, чем при прямом подходе. Возможные варианты двойственных процедур представлены, например, в [20,25,91,97,100-103]. Остановимся на следующей двойственной градиентной процедуре построения непрерывных оптимальных планов из [39-41,43,44,94,98,104-106] :

Шаг 1. Зададим начальный невырожденный план 0 и по формуле (1.9) вычислим нормированную матрицу M 0 плана. Положим k 0.

Шаг 2. Найдём локальный максимум методом последовательного квадратичного программирования (см. п. 1.1.2).

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

Шаг 3. Вычислим k по формуле Здесь где произведём его «очистку» в соответствии с рекомендациями из [20], положим k k 1 и перейдём на шаг 2.

Соответствие значений параметров X M( ), двойственной процедуры критериям А- и D- оптимальности такое же, как в таблице 1.1.

Приведённый алгоритм требует вычисления градиента Для критерия А – оптимальности (1.10) получаем:

В случае критерия D – оптимальности (1.11) имеем:

Соотношения (1.20),(1.21) показывают, что, как и в случае прямой градиентной процедуры, для применения двойственной градиентной процедуры синтеза непрерывных А - или D - оптимальных планов необходимо разработать алгоритмы вычисления ИМФ и ее производной по компонентам точек спектра плана эксперимента.

1.1.6 Построение дискретных оптимальных планов Предположим, что мы при помощи прямой, двойственной или комбинированной прямой – двойственной процедуры синтезировали непрерывный план Его практическое использование затруднено тем обстоятельством, что веса p* представляют собой произвольные вещественные числа, заключенные в интервале от нуля до единицы, и в случае заданного числа возможных запусков системы величины k* p* могут не являться целыми числами. Проведение эксперимента требует округления величин k* до целых чисел. Очевидно, что полученный в результате такого округления план будет отличаться от оптимального непрерывного плана, причем приближение будет тем лучше, чем больше число возможных запусков.

Приведем следующий алгоритм «округления» непрерывного плана до точного из [25]:

Шаг 1. Вычислим числа i и i по формулам (здесь z - ближайшее к z целое число, большее z; z - целая часть числа z).

Шаг 2. Вычислим ' и '', воспользовавшись выражениями их значений. Положим j=1.

рядоченном наборе, то положим s j 1, в противном случае s j 0.

Шаг 5. Если j q, увеличим j на единицу и перейдем на шаг 4. В противном случае сформируем приближённый дискретный план 1.1.7 Схема процедуры активной параметрической идентификации систем с предварительно выбранной модельной структурой Изложенные в данном пункте вопросы оценивания неизвестных параметров и планирования эксперимента позволяют в случае предварительно выбранной структуры математической модели представить процедуру активной параметрической идентификации в виде схемы, помещенной на рисунке 1.2.

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

Книга В.Г. Горского, Ю.П. Адлера, А.М. Талалая [24] представляет ставшее классическим введение в указанную проблему и содержит универсально изложенные постановки задач планирования эксперимента и возможные подходы к их решению, применимые для различных способов задания моделей динамических систем.

Вопросы построения стохастических дискретных моделей в форме передаточных функций с использованием планирования входных сигналов во временной и частотной областях всесторонне рассматривались в классических монографиях Г. Гудвина и Р. Пейна [107], М. Зейропа [108], Л. Льюнга [11], а также, например, в работах Э. Уолтера и Л. Пронзато [13], Л. Пронзато [109], Рисунок 1.2 - Схема процедуры активной параметрической идентификации В.И. Денисова, И.Л. Еланцевой и В.М. Чубича [103], В.И. Денисова, И.Л. Еланцевой и И.А. Полетаевой [110,111].

Различные аспекты проблемы планирования входных сигналов во временной и частотной областях при построении дискретных линейных стационарных моделей со стохастическими уравнениями состояний и измерений обсуждались, например, в классических статьях Р. Мехры [88,89], монографиях Г.Гудвина, Р. Пейна [107] и В.И. Денисова, В.М. Чубича, О.С. Черниковой, Д.И. Бобылевой [94], а также в работах А.Ж. Абденова [92,112,113], А.Ж. Абденова и В.Ю. Пейсаховича [114], А.А. Попова [115], В.И. Денисова, В.М. Чубича и О.С. Черниковой [104,105,116-118]. Активная параметрическая идентификация динамических систем, описываемых стохастическими дискретными линейными нестационарными моделями в пространстве состояний, рассматривалась в статье В.И. Денисова, И.Л. Еланцевой и В.М. Чубича [103]. Подчеркнем, что в публикациях [94,104,105,116-118] результаты представлены для математических моделей с неизвестными параметрами в матрицах состояния, управления, возмущения, измерения, в начальных условиях и в ковариационных матрицах возмущений (шума системы) и ошибок измерений (шума измерений).

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

Вопросы планирования входных сигналов при построении непрерывнодискретных моделей с детерминированными (в том числе и нелинейными) уравнениями состояний рассматривались, например, в монографии В.И. Денисова и А.А. Попова [33], в статьях В.Н. Овчаренко [120-124], Э. Морелли [125], К. Жобертье, Ф. Бурнонвилля, П. Котон и Ф. Ренделя [126], К. Жобертье, Л. Денис-Видаль, П. Котон и Г. Джоли-Бланшар [127]. Различным аспектам проблемы планирования входных сигналов во временной области при построении непрерывно-дискретных линейных стационарных моделей со стохастическими уравнениями состояний и измерений посвящены, например, соответствующие разделы монографии В.И. Денисова, В.М. Чубича, О.С. Черниковой и Д.И. Бобылевой [94], статьи А.Ж. Абденова [128], А.Ж. Абденова, В.И. Денисова, В.М. Чубича [129], А.А. Попова [115], В.И. Денисова, В.М. Чубича, Д.И. Бобылевой [98, 130,131]. Подчеркнем, что в публикациях [94,98,130,131] результаты представлены для математических моделей с неизвестными параметрами в матрицах состояния, управления, возмущения, измерения, в начальных условиях и в ковариационных матрицах возмущений и ошибок измерений в различных комбинациях.

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

1.3 Структурно-вероятностное описание моделей 1.3.1 Модели дискретных систем Рассмотрим управляемую, наблюдаемую, идентифицируемую модель динамической системы вида где x t k – n - вектор состояния; u t k – детерминированный r - вектор управления (входа); w t k – p - вектор шума системы (возмущения); y t k 1 – m - вектор измерения (выхода); v t k 1 – m - вектор шума (ошибки) измерения.

Априорные предположения:

случайные векторы w t k и v t k 1 образуют стационарные белые гауссовские последовательности, для которых (здесь и далее ki - символ Кронекера);

начальное состояние x t 0 имеет нормальное распределение с параметрами и не коррелирует с w t k и v t k 1 при любых значениях k ;

включающий в себя элементы вектор - функций f x t k,u t k, t k, комбинациях.

Частным случаем модели (1.22), (1.23) являются модели линейной нестационарной системы и линейной стационарной системы При активной параметрической идентификации нелинейных систем (1.22), (1.23) с указанными априорными предположениями будем применять временную [132-135] и статистическую [3,136-142] линеаризации, в результате сводя исходную задачу к соответствующей задаче для модели вида (1.24), (1.25) со специальным образом определенными векторами a(t k ), A(t k 1) и матрицами F(t k ), H(t k 1).

Линеаризация во временной области возможна в предположении, что вектор – функции f x t k,u t k, t k и h x t k 1, t k 1 непрерывны и дифференцируемы по x t k, u t k и x t k 1 соответственно [39,44,143].

Зададим номинальную траекторию xн t k 1, k 0,1,..., N 1 в соответствии с равенством и, отбросив члены второго и более высоких порядков, запишем уравнения линеаризованной модели:

Положив приходим к модели линейной нестационарной системы (1.24), (1.25).

В отличие от временной линеаризации, статистическая линеаризация применима к неоднозначным функциям и к существенным нелинейностям, имеющим характеристики с угловыми точками и разрывами. Статистическая линеаризация заключается в замене нелинейной характеристики эквивалентной (в вероятностном смысле) линеаризованной функциональной зависимостью [41,106,144,145]. В соответствии с [141] получим:

Отметим, что вычисление интегралов в формулах (1.35), (1.36) можно существенно упростить с помощью представленных в [3,136,138,139,141] выражений для коэффициентов статистической линеаризации типовых одномерных нелинейностей, встречающихся в системах автоматического управления.

Подставив выражения (1.33), (1.34) в уравнения (1.22), (1.23) соответственно, с учетом обозначений приходим к модели линейной нестационарной системы (1.24), (1.25).

1.3.2 Модели непрерывно-дискретных систем Рассмотрим управляемую, наблюдаемую, идентифицируемую модель динамической системы вида:

где x t – n - вектор состояния; u t – детерминированный r - вектор управления (входа); w t – p - вектор шума системы (возмущения); y t k 1 – mвектор измерения (выхода); v t k 1 – m- вектор шума (ошибки) измерения.

Априорные предположения:

случайные векторы w t и v t k 1 являются стационарными белыми гауссовскими шумами, для которых (здесь и далее (t ) - дельта-функция Дирака);

начальное состояние x t 0 имеет нормальное распределение с параметрами и не коррелирует с w t и v t k 1 при любых значениях переменной k;

комбинациях.

Частным случаем модели (1.43), (1.44) являются модели линейной нестационарной системы и линейной стационарной системы Непрерывность и дифференцируемость вектор – функций f x t,u t, t и вить линеаризацию во временной области [40,43,146] нелинейной модели (1.43), (1.44).

венства Разложив вектор - функции f x t,u t, t и h x t k 1, t k 1 в ряды Тейлора в окрестностях точек x н t,u н t и x н t k 1 соответственно и отбросив члены второго и более высоких порядков, запишем уравнения линеаризованной модели Положив приходим к модели линейной нестационарной системы (1.45),(1.46).

Выполним статистическую линеаризацию нелинейной модели (1.43), (1.44). В соответствии с [141] получим:

где следующих дифференциальных уравнений:

Подчеркнем, что при вычислении интегралов в формулах (1.56), (1.57) может оказаться полезным использование выражений из [3,136,138,139,141] для коэффициентов статистической линеаризации типовых одномерных нелинейностей.

Подставив выражения (1.54), (1.55) в уравнения (1.43), (1.44) соответственно, с учетом обозначений приходим к модели линейной нестационарной системы (1.45), (1.46).

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

Для достижения поставленной цели в диссертационной работе ставятся и решаются следующие основные задачи:

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

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

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

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

1.5. Выводы 1. Определены предмет, теоретическая и методологическая база исследований.

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

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

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

2 Оценивание параметров моделей стохастических Оценивание параметров моделей стохастических динамических систем в пространстве состояний можно осуществлять либо путем обработки всех накопленных к моменту окончания идентификационного эксперимента измерительных данных, либо в реальном масштабе времени, когда оценки неизвестных параметров рекуррентно пересчитываются по очередному поступившему измерению. В первом случае используются методы наименьших квадратов (МНК), максимального правдоподобия (ММП) и максимума апостериорной вероятности (МАВ), упоминавшиеся в п. 1.1.2. Во втором случае (см., например, [1,13,70,129,147-150]) фильтр Калмана [1,141,151-156], что предполагает модификацию вектора состояния рассматриваемой модели за счет включения в него оцениваемых параметров.

В выполнялось оценивание параметров непрерывнодискретных моделей с нелинейными детерминированными уравнениями состояний при помощи МНК. В [94,104,158] и в [94] данный метод использовался для оценивания параметров моделей, соответственно, дискретных (1.26), (1.27) и непрерывно-дискретных (1.47), (1.48) линейных стационарных систем во временной области.

В [49,51,159] и [150] вычислялись ОМП моделей нелинейных непрерывно-дискретных систем с детерминированными и стохастическими уравнениями состояний соответственно, причем во втором случае выполнялась линеаризация модели во временной области и функция правдоподобия строилась на основе уравнений расширенного фильтра Калмана. Оценивание параметров линейных стационарных дискретных и непрерывно-дискретных моделей вида (1.26), (1.27) и (1.47), (1.48) с применением ММП во временной области рассматривалось в [70,94,104,160] и [62,94,130,161,162] соответственно. Как и в [160], в [163,164] при записи функции правдоподобия использовались фильтры квадратного корня [156,165,166], но применительно к линейным нестационарным дискретным моделям.

Оценивание параметров моделей нелинейных непрерывно-дискретных систем с детерминированными уравнениями состояний при помощи МАВ рассматривалось в [167].

К передовым методам параметрической идентификации, использующимся при построении стохастических моделей линейных стационарных дискретных систем и реализованным в рамках пакета System Identification Toolbox [168,169] программной системы MATLAB, следует отнести семейство так называемых подпространственных методов (в английском варианте Subspace methods), к которым относятся методы N4SID (Numerical algorithms for subspace state space system identification), MOESP (Multivariable output error state space) и CVA (Сanonical variate analysis), описанные в [14,149,170-173], [14,149,172-175] и [172,176,177] соответственно. Указанные методы не предполагают численного решения задачи нелинейного программирования с ограничениями (1.2), что устраняет присущую МНК, ММП и МАВ проблему нахождения глобального экстремума, и работают наиболее эффективно при использовании канонических форм представления многомерных модельных структур [11,173,178].

В данном разделе мы рассмотрим вопросы оценивания параметров моделей стохастических динамических систем, взяв за основу подход и алгоритмы из [94] и модифицировав их с точки зрения вычислительной рациональности и универсальности, обеспечивая возможность их применения как к линейным нестационарным, так и к нелинейным дискретным и непрерывно-дискретным моделям.

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

запусков системы, причем сигнал U1 он подает на вход системы k1 раз, сигнал U 2 - k 2 раз и т.д., наконец, сигнал U q - k q раз. В этом случае дискретный нормированный план эксперимента (1.6) представляет собой совокупность точек U1, U 2,…, U q и соответствующих им долей повторных запусков:

Множество планирования U определяется ограничениями на условия проведения эксперимента.

Обозначим через Yij j-ю реализацию выходного сигнала, соответствуюyij (t1) проведения по плану идентификационных экспериментов будут подготовлены данные наблюдений Ui,Yij, i 1,2,...,q, j 1,2,...,ki. Подчеркнем, Априорные предположения, изложенные в п. 1.3, позволяют воспользоваться ММП для оценивания неизвестных параметров (полезные свойства ОМП перечислены в п. 1.1.2). В соответствии с ММП необходимо найти такие значения параметров, для которых где ln L ; - логарифмическая функция правдоподобия, которую необходимо составить.

Поскольку запуски системы выполняются независимо, критерий максимального правдоподобия ( ; ) связан с логарифмическими функциями правдоподобия для одного запуска ln L( ;Ui,Yij) соотношением где ln L( ;Ui,Yij) выписывается для каждой модельной структуры по-своему.

Для численного решения задачи нелинейного программирования (2.2) воспользуемся методом последовательного квадратичного программирования (см. п. 1.1.2), что предполагает необходимость разработки алгоритмов вычисления значений критериев идентификации и их градиентов.

2.1 Оценивание параметров моделей дискретных систем 2.1.1 Критерий максимального правдоподобия и алгоритм вычисления его значения для линейных нестационарных моделей Следуя [4,94,161,162], для математической модели (1.24), (1.25) запишем что позволяет представить соотношение (2.3) в виде где ij (t k 1) и B(t k 1) вычисляются по ниже приведенным уравнениям дискретного фильтра Калмана [1,66,141,151-156,179-181]:

с начальными условиями: xij (t 0 | t 0 ) x(t 0 ), (t 0|t 0 ) xij (t k 1 | t k ) - оценка одношагового прогнозирования состояния x(t k 1), соответствующая паре Ui,Yij ;

(t k 1 | t k ) - ковариационная матрица ошибок одношагового прогнозирования;

xij (t k 1 | t k 1) - оценка фильтрации состояния x(t k 1), соответствующая паре Ui,Yij ;

(t k 1 | t k 1) - ковариационная матрица ошибок фильтрации.

Эквивалентная выражению (2.4) запись позволяет предложить следующий алгоритм вычисления значения критерия Шаг 1. Определить Q, R, x(t 0 ), (t 0 ).

Шаг 3. Определить F(t k ), (t k ), H(t k 1) и по формулам (2.6), (2.8), (2.9), (2.11) Шаг 5. Определить a ui t k,t k.

Шаг 6. Положить j=1.

Шаг 7. Выбрать xij (t k | t k ) x(t 0 ), если k=0.

Шаг 9. Положить Шаг 10. Увеличить j на единицу. Если j k i, перейти на шаг 7.

Шаг 11. Увеличить i на единицу. Если i q, перейти на шаг 5.

2.1.2 Критерий максимального правдоподобия и алгоритм вычисления его значения для линеаризованных моделей В случае применения временной или статистической линеаризации (подробно см. п. 1.3.1) к математической модели (1.22), (1.23) что при подстановке в соотношение (2.3) приводит к следующей записи критерия максимального правдоподобия [41,44,106] где ij (t k 1) и Bi (t k 1) вычисляются по уравнениям (см. формулы (2.5)-(2.11)):

с начальными условиями: xij (t 0 | t 0 ) x(t 0 ), i (t 0|t 0 ) Уточним, что при временной линеаризации причем При статистической линеаризации где Эквивалентная выражению (2.13) запись позволяет предложить следующий алгоритм вычисления значения критерия Шаг 1. Определить Q, R, x(t 0 ), (t 0 ).

Шаг 3. Определить (t k ). Положить 1 0, i 1.

Шаг 4. Выбрать i (t k | t k ) P(t 0 ), если k=0.

Шаг 5. Определить Fi (t k ), Hi (t k 1) в случае временной линеаризации при помощи выражений (2.22), (2.24); в случае статистической линеаризации необходимо воспользоваться равенствами (2.26), (2.28). По формулам Pi (t k 1 | t k 1). Определить a ui t k,t k, Ai (t k 1) в случае временной линеаризации при помощи соотношений (2.21), (2.23); в случае статистической линеаризации необходимо воспользоваться формулами (2.25), Шаг 7. Выбрать xij (t k | t k ) x(t 0 ), если k=0.

Шаг 9. Положить Шаг 10. Увеличить j на единицу. Если j k i, перейти на шаг 7.

Шаг 11. Положить Шаг 12. Увеличить i на единицу. Если i q, перейти на шаг 4.

2.1.3 Алгоритм вычисления градиента критерия максимального правдоподобия для линейных нестационарных моделей и симметричности матрицы B(t k 1), получим Шаг 3. Определить F(t k ), (t k ), H(t k 1), по формулам (2.6), (2.8), (2.9), (2.11) следующим из равенств (2.6), (2.8), (2.9), (2.11), Шаг 7. Положить j=1.

Шаг 8. Выбрать xij (t k | t k ) Шаг 10. По формулам, следующим из равенств (2.5), (2.7), (2.10), Шаг 11. Положить Шаг 12. Увеличить j на единицу. Если j k i, перейти на шаг 8.

Шаг 13. Увеличить i на единицу. Если i q, перейти на шаг 6.

Шаг 15. Увеличить k на единицу. Если k N 1, перейти на шаг 3, иначе закончить процесс.

2.1.4. Алгоритм вычисления градиента критерия максимального правдоподобия для линеаризованных моделей (2.31) и симметричности матрицы Bi (t k 1), получим случае статистической линеаризации необходимо воспользоваться равенствами (2.26),(2.28). По формулам (2.15), (2.17), (2.18), (2.20) найти Шаг 6. По формулам, следующим из равенств (2.15), (2.17), (2.18), (2.20), в случае временной линеаризации при помощи соотношений (2.21),(2.23);

в случае статистической линеаризации необходимо воспользоваться формулами (2.25), (2.27).

Шаг 8. Выбрать xij (t k | t k ) x(t 0 ), Шаг 10. По формулам, следующим из равенств (2.14), (2.16), (2.19), Шаг 11. Положить Шаг 12. Увеличить j на единицу. Если j k i, перейти на шаг 8.

Шаг 13. Положить Шаг 14. Увеличить i на единицу. Если i q, перейти на шаг 4.

2.2 Оценивание параметров моделей непрерывно-дискретных систем Для моделей непрерывно-дискретных систем точки спектра плана (2.1) 2.2.1 Особенности алгоритмов вычисления значений критериев максимального правдоподобия для линейных нестационарных и линеаризованных моделей Для линейных нестационарных моделей непрерывно-дискретных систем критерий максимального правдоподобия определяется, как и в п. 2.1.1, выражением (2.4), но ij (t k 1) и B(t k 1) в этом случае вычисляются по следующим уравнениям непрерывно-дискретного фильтра Калмана [151,154,161,162]:

с начальными условиями: xij (t 0 | t 0 ) x(t 0 ), (t 0|t 0 ) Данное обстоятельство позволяет говорить о том, что большинство шагов алгоритма вычисления значения критерия максимального правдоподобия для линейных нестационарных моделей при некотором фиксированном совпадут с соответствующими шагами алгоритма, представленного в п. 2.1.1.

Уточним, какие именно шаги требуют корректировки, не повторяя шаги, оставшиеся без изменений:

Шаг 3. Определить H(t k 1) и по формулам (2.35), (2.37), (2.38), (2.40) найти Шаг 5. Определить a ui t,t.

В случае применения временной или статистической линеаризации (подробно см. п. 1.3.2) к математической модели (1.43), (1.44) критерий максимального правдоподобия по-прежнему записывается в виде (2.13), с той разницей, что ij (t k 1) и Bi (t k 1) вычисляются теперь по следующим уравнениям (см. формулы (2.34)-(2.40)):

с начальными условиями: xij (t 0 | t 0 ) x(t 0 ), i (t 0|t 0 ) Уточним, что при временной линеаризации причем При статистической линеаризации где Логика предыдущих рассуждений позволяет сделать вывод о том, что большинство шагов алгоритма вычисления значения критерия максимального правдоподобия для линеаризованных моделей при некотором фиксированном совпадут с соответствующими шагами алгоритма из п. 2.1.2. Уточним, какие именно шаги требуют корректировки, не повторяя шаги, оставшиеся без изменений:

Шаг 3. Положить 1 0, i 1.

Шаг 5. Определить Fi (t), t t k, t k 1, Hi (t k 1) в случае временной линеаризации при помощи выражений (2.49), (2.51); в случае статистической линеаризации необходимо воспользоваться равенствами (2.53), (2.55).

случае временной линеаризации при помощи соотношений (2.48), (2.50);

в случае статистической линеаризации необходимо воспользоваться формулами (2.52), (2.54).

2.2.2 Особенности алгоритмов вычисления градиентов критериев максимального правдоподобия для линейных нестационарных и Для линейных нестационарных моделей градиент критерия максимального правдоподобия определяется, как и в п. 2.1.3, выражением (2.32). Это позволяет заметить, что в алгоритме вычисления градиента критерия максимального правдоподобия при некотором фиксированном большинство шагов совпадают с соответствующими шагами алгоритма из п. 2.1.3. Уточним, какие именно шаги требуют корректировки, не повторяя шаги, оставшиеся без изменений:

Шаг 3. Определить H(t k 1), по формулам (2.35), (2.37), (2.38), (2.40) найти Шаг 10. По формулам, следующим из равенств (2.34), (2.36), (2.39), Точно также большинство шагов алгоритма вычисления градиента критерия максимального правдоподобия для линеаризованных моделей при некотором фиксированном совпадает с соответствующими шагами алгоритма из п. 2.1.4. Уточним, какие именно шаги требуют корректировки, не повторяя шаги, оставшиеся без изменений:

при помощи выражений (2.49), (2.51); в случае статистической линеаризации необходимо воспользоваться равенствами (2.53), (2.55). По формулам (2.42),(2.44),(2.45),(2.47) найти Pi (t k 1 | t k ), i (t i Шаг 6. По формулам, следующим из равенств (2.42), (2.44), (2.45), (2.47), соотношений (2.48), (2.50); в случае статистической линеаризации необходимо воспользоваться формулами (2.52), (2.54).

Шаг 10. По формулам, следующим из равенств (2.41), (2.43), (2.46), В соответствии с целью исследования представлены алгоритмы вычисления значений критериев максимального правдоподобия и их градиентов для гауссовских линейных нестационарных и линеаризованных дискретных и непрерывно-дискретных моделей с неизвестными параметрами в уравнениях состояния и измерения, в начальных условиях и в ковариационных матрицах шумов системы и измерений. Данные алгоритмы обеспечивают возможность нахождения ОМП неизвестных параметров рассматриваемых математических моделей методом последовательного квадратичного программирования с вычислением соответствующих градиентов по рекуррентным аналитическим формулам и реализованы в рамках программных комплексов ПК-I и ПК-II активной параметрической идентификации стохастических динамических систем, описанных в разделе 5.

3 Планирование эксперимента для моделей стохастических дискретных систем 3.1 Вычисление информационной матрицы Фишера ИМФ занимает одно из центральных мест в теории планирования эксперимента, поскольку именно она используется для вычисления информационной матрицы всего плана и фигурирует в алгоритмах численного построения оптимальных планов.

3.1.1 Вывод информационной матрицы Фишера для линейных нестационарных моделей Первым аналитическое выражение ИМФ для гауссовских линейных стационарных дискретных систем, описываемых моделями в пространстве состояний, получил Р.К. Мехра в [88,89], для случая вхождения неизвестных параметров в матрицы состояния и управления. Это выражение в дальнейшем было уточнено в [34,90,91,182] и обобщено в [183] на случай вхождения неизвестных параметров также и в матрицу измерения. Отметим, что при выводе выражений ИМФ в [34,88-91,182,183] было использовано представление фильтра Калмана в форме нормированной обновляющей последовательности.

Другой подход к вычислению ИМФ был предложен А.А. Поповым в [184]. Как и в [34,88-91,182], он рассмотрел случай вхождения неизвестных параметров в матрицы состояния и управления, но применил уравнения фильтра Калмана в форме ненормированной обновляющей последовательности, что позволило оптимизировать процедуру вычисления информационной матрицы за счет исключения необходимости извлечения квадратного корня из матрицы в каждый момент времени. Эта идея впоследствии была использована в [94,185] при выводе ИМФ для гауссовских линейных стационарных систем в случае вхождения неизвестных параметров в матрицы состояния, управления, возмущения, измерения, в начальные условия и в ковариационные матрицы возмущений и ошибок измерений.

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

Теорема 3.1 [143] Для математической модели (1.24), (1.25) с априорными предположениями из п. 1.3.1 и неизвестными параметрами, входящими в x t 0 в различных комбинациях, элементы ИМФ определяются выражением кретного фильтра Калмана.

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

Введем в рассмотрение обновляющую последовательность где x t k 1 | t k определяется по следующим уравнениям дискретного фильтра Калмана [1,66,141,151-156,179-181]:

Запишем логарифмическую функцию правдоподобия:

Элементы ИМФ связаны с логарифмической функцией правдоподобия (3.9) соотношением [6,11,47,48,67-69]:

в котором усреднение берется по выборочному пространству.

Подставим последнее равенство в (3.10). Воспользовавшись тем, что получим:

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

Поскольку можем записать:

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

формулой которая выводится по индукции из уравнений (3.3), (3.7) фильтра Калмана и в которой Положив в выражении (3.17) получим Принимая во внимание, что будем иметь:

Для вычисления второго математического ожидания продифференцируем x t k 1 | t k из (3.18) по i :

В силу того, что с учетом соотношений (3.19), (3.20) получим Для вычисления третьего математического ожидания продифференцируx tk 1 | tk Поскольку с учетом соотношений (3.19), (3.20), (3.23), (3.26) получим Подставив соотношения (3.19), (3.21), (3.24), (3.27) в равенства (3.14) и (3.16), приходим к тому, что Соотношения (3.15), (3.28), (3.29) позволяют записать формулу (3.11) в виде (3.1). Таким образом, теорема доказана.

Следствие. Для математической модели стационарной системы (1.26), R, P t 0 и вектор x t 0 в различных комбинациях, формула (3.1) принимает вид что совпадает с результатом из [94,185].

Заметим, что если в модели стационарной системы (1.26), (1.27) неизвестные параметры входят только в матрицы F и, выражение (3.30) записывается следующим образом и совпадает с соответствующим аналогом из [184].

Непосредственно соотношение (3.1) непригодно для программной реализации в силу того, что оно содержит нераскрытые математические ожидания.

Перейдем к вычислению этих математических ожиданий.

Введем в рассмотрение расширенный вектор состояния размером n (s 1) и матрицы Ci [,...,,,,..., ], i 0,1,...,s размером n n(s 1). С учётом этих обозначений Соотношения (3.31) – (3.36) позволяют записать выражение (3.1) в виде Обозначим Тогда Получим рекуррентные соотношения для x A t k 1 и A t k 1.

Дифференцируя (3.42) по i, с учетом формулы (3.12) приходим к тому, что С учетом этих обозначений равенства (3.42) и (3.43) можно объединить следующей компактной формулой При k 0 из соотношения (3.3) имеем Поэтому и, следовательно, Вычисляя математические ожидания от (3.47) и (3.48), получаем:

Из равенства (3.49) сразу же вытекает, что При k 1,2,..., N 1 в силу соотношений (3.39), (3.47), (3.49) имеем:

Принимая во внимание выражения (3.19), (3.21) и (3.24), находим Таким образом, второе и третье слагаемые в формуле (3.51) равны нулю.

С учетом равенства (3.50) получаем следующее рекуррентное соотношение:

3.1.2 Алгоритм вычисления информационной матрицы Фишера для линейных нестационарных моделей Приведенные в предыдущем пункте аналитические выкладки позволяют разработать алгоритм вычисления ИМФ для математической модели (1.24), (1.25) с априорными предположениями из п. 1.3.1 при некотором фиксированном значении вектора оцениваемых параметров в уравнениях состояния, наблюдения, в начальных условиях и в ковариационных матрицах шумов системы и измерений. Представим возможный вариант такого алгоритма, следуя с незначительными изменениями [143]:

Шаг 1. Определить Шаг 4. Если k 0, вычислить x A t k 1 и A t k 1 по формулам (3.49) и (3.52) соответственно, после чего перейти на шаг 8.

Шаг 6. Сформировать матрицы FA t k, K A t k и вектор a A t k, воспользовавшись равенствами (3.44) - (3.46).

по формуле Шаг 10. Используя выражения (3.37), с учетом (3.38), (3.40), получить приращение, отвечающее текущему значению k.

3.1.3 Вычисление информационной матрицы Фишера для линейных нестационарных моделей, полученных в результате линеаризации В случае применения линеаризации во временной области (подробно см.

п. 1.3.1) к математической модели (1.22), (1.23) алгоритм вычисления ИМФ из п. 3.1.2 претерпевает незначительные изменения. Следуя [143], уточним, какие именно шаги требуют корректировки и представим их в новой редакции, не повторяя шаги, оставшиеся без изменений:

Шаг 3. Определить uн (t k ), u(t k ) и при помощи равенства (1.29) найти a u(t k ), t k В случае применения статистической линеаризации (подробно см.

п. 1.3.1) к математической модели (1.22), (1.23) алгоритм вычисления ИМФ из п. 3.1.2 претерпевает следующие изменения (шаги, оставшиеся без изменений, опущены):

Шаг 3. Определить u(t k ) и при помощи равенства (1.39) найти Используя выражение (1.40), получить F x(t k ), P(t k ), u(t k ), t k и мул (1.37) и (1.38) соответственно. При помощи равенств (1.41), (1.42) 3.2 Планирование входных сигналов Планирование входных сигналов является, по-видимому, наиболее эффективным способом управления экспериментом, использующимся при построении моделей стохастических динамических систем. Отметим, что свобода в выборе входных характеристик существенно зависит от приложений. В экономических и экологических системах у экспериментатора нет возможности воздействовать на систему с целью проведения идентификационных экспериментов, в то время как в лабораторных условиях и на стадиях разработки нового оборудования выбор входных величин имеет лишь амплитудные и мощностные ограничения.

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

Следуя [24], будем считать, что непрерывный нормированный план (1.7) в данном случае может быть задан в виде При этом каждая точка спектра плана представляет собой последовательность импульсов, «развернутую во времени». Применительно к моделям, структурновероятностное описание которых представлено в п. 1.3.1, это означает, что Замкнутое ограниченное множество U представляет собой область допустимых входных сигналов. На практике чаще всего используются два типа ограничений: на амплитуду и на мощность [24,88,89]. В первом случае область допустимых входных сигналов определяется соотношением и представляет собой координатный параллелепипед. Во втором случае и является шаром.

Нормированная информационная матрица плана (3.53) определяется соотношением в котором ИМФ одноточечных планов* вычисляются в соответствии с подразделом 3.1.

Построение оптимальных планов будем осуществлять при помощи прямой и двойственной градиентных процедур, что предполагает вычисление следующих градиентов (см. формулы (1.13), (1.14) и (1.19)):

Соотношения (1.15) - (1.18) и (1.20), (1.21) показывают, что для этого необходимо получить выражения для производных от ИМФ по компонентам входного сигнала и разработать соответствующие вычислительные алгоритмы.

Напомним, что ИМФ зависят от подлежащих оцениванию неизвестных параметров.

3.2.1 Нахождение производных информационной матрицы Фишера по компонентам входного сигнала для линейных нестационарных Будем считать, что в модели состояния управляемой линейной нестационарной системы (1.24), (1.25) Q, R, P t 0 и векторы b t k, A t k 1, x t 0 и выполнены априорные предположения из п. 1.3.1.

Для вычисления производных представим ИМФ в виде суммы двух слагаемых, одно из которых зависит от входного сигнала, а другое – нет:

Чтобы получить такое разложение, подставим равенства (3.38), (3.40) в формулу (3.37). В результате получим следующие выражения для элементов матриц Из соотношений (3.56), (3.57) следует, что этого, воспользовавшись формулами (3.45) и (3.55), представим выражение для a A (t k ) в виде суммы двух слагаемых, одно из которых не зависит, а другое – зависит от u(t k ) :

где Из соотношения (3.61) следует, что Подставляя равенство (3.55) в формулу (3.60) при k=0, с учетом выражения (3.63) получаем 3.2.2 Алгоритм вычисления производных информационной матрицы Фишера по компонентам входного сигнала для линейных нестационарных моделей Приведенные в предыдущем пункте аналитические выкладки позволяют разработать алгоритм вычисления производных ИМФ по компонентам входноU; ) ванном значении вектора неизвестных параметров в уравнениях состояния, наблюдения, в начальных условиях и ковариационных матрицах шумов системы и измерений. Представим возможный вариант такого алгоритма, следуя с незначительными изменениями [44]:

Шаг 4. Если k 0, вычислить (см. (3.49)) Шаг 5. Найти t k по формуле (см. (3.41)) Шаг 6. Сформировать матрицу FA t k и вектор-столбец a A t k в соответствии с равенствами (см. (3.44), (3.45)):

Шаг 7. Вычислить x A t k 1 по формуле (см. (3.49)) шения (см. (3.4), (3.5), (3.6), (3.8)):

ва (3.65) и перейти на шаг 14.

Шаг 14. Используя выражение (3.59), получить приращения Шаг 17. Увеличить k на единицу. Если k N 1, перейти на шаг 3. В противном случае закончить процесс.

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

п. 1.3.1) к математической модели (1.22), (1.23) алгоритм вычисления производных от ИМФ по компонентам входного сигнала из п. 3.2.2 претерпевает незначительные изменения. Следуя [44], уточним, какие именно шаги требуют корректировки и представим их в новой редакции, не повторяя шаги, оставшиеся без изменений:

Шаг 2. Положить Шаг 3. Определить uн (t k ), u(t k ). Воспользовавшись равенствами (см. (1.29), (3.55), (1.30)) 3.2.4 Нахождение производных информационной матрицы Фишера по компонентам входного сигнала для моделей, полученных в результате статистической линеаризации В случае применения статистической линеаризации (подробно см. п.

1.3.1) к уравнениям (1.22), (1.23) в получившейся линеаризованной модели вида (1.24), (1.25) a x t k,P t k,u t k,t k, F x t k,P t k,u t k,t k, A x t k 1,P t k 1,t k и H x t k 1,P t k 1,t k 1 зависят (см. соотношение (1.37)) от u t 0,u t1,....,u t k.

A t k 1 также зависят от указанных переменных, разложение (3.56) становится невозможным и существенно усложняется вычисление производных ИМФ по компонентам входного сигнала.

Остановимся на варианте с нелинейной моделью (1.22) для вектора состояния и линейной моделью (1.23) для вектора измерения. Это устранит зависимость A t k 1 и H t k 1 от входного сигнала и в некоторой степени упростит вычисление производных от ИМФ по его компонентам.

Подставив равенства (3.38), (3.40) в формулу (3.37), получим следующее выражение для элементов ИМФ:

Из соотношения (3.66) следует, что (3.49). В результате получим:

xA tk формулой (3.44):

где выражения (3.45) следует, что ния (3.52). Получим:

3.2.5 Алгоритм вычисления производных информационной матрицы Фишера по компонентам входного сигнала для моделей, полученных в результате статистической линеаризации Приведенные в предыдущем пункте аналитические выкладки позволяют разработать алгоритм вычисления производных ИМФ по компонентам входноM U;

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

с использованием равенств (1.39) и (1.40) соответственно.

Шаг 4. Если k 0, то определить x A t k 1 и A t k 1 по формулам (см.

(3.49), (3.52)) соответствии с равенствами (3.44), (3.46) и (3.45):

Шаг 7. Вычислить x A t k 1 и A t k 1 по формулам (см. (3.49), (3.52)) ношения (3.4), (3.5), (3.6), (3.8):

ношений (см. шаг 10):

формулам (1.37), (1.38) и перейти на шаг 3, иначе закончить процесс.

3.2.6 Планирование эксперимента как задача дискретного оптимального Ограничимся построением одноточечных планов и выберем в качестве критерия оптимальности след ИМФ, т.е. по аналогии с работами [95,96], в которых рассматривались детерминированные линейные непрерывные модели, будем решать оптимизационную задачу но применительно к стохастическим нелинейным дискретным моделям. Отметим, что в отличие, например, от критериев D-, A-, Е- оптимальности, критерий (3.73) не имеет четкой статистической интерпретации, но позволяет свести задачу планирования эксперимента к задаче дискретного оптимального управления и применить соответствующие методы для ее решения. Покажем, как это можно сделать, следуя материалам статьи [45].

Предположим, что выполнены все априорные предположения из п. 1.3.1 и математическая модель динамической системы (1.22), (1.23) позволяет осуществить временную линеаризацию. Тогда, в силу разложения (3.56) для ИМФ линеаризованной модели, задачу (3.73) можно заменить эквивалентной:

С учетом выражения (3.57) после некоторых преобразований получим, что где причем Равенство (3.61) позволяет представить формулу (3.49) в виде FA t k A t k и bA t k определяются выражениями (3.63), (3.62) соответственно.

Подставив равенство (3.75) в формулу (3.74), с учетом симметричности матрицы ZA1 t k 1, получим, что Здесь Поскольку d0 t k в соотношении (3.76) не зависит от U, при максимизации критерия J1 это слагаемое можно не учитывать. В результате приходим к критерию который вместе с системой (3.75) позволяет рассматривать задачу синтеза оптимальных входных сигналов как задачу дискретного оптимального управления с суммарным показателем качества [186]. Для ее решения воспользуемся методом последовательного квадратичного программирования с использованием конечно-разностной аппроксимации градиента.

Сформулированную задачу можно свести к следующей задаче оптимизации конечного состояния [186]:

Для решения этой задачи воспользуемся методом последовательного улучшения управлений Л.И. Шатровского [187,188], адаптировав его к дискретной задаче.