Метод наименьших квадратов в задачах вариационного оценивания состояния нелинейных динамических систем Текст научной статьи по специальности «Математика»

обработка информации и управление X

УДК 629.191

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

В. И. Миронов,

доктор техн. наук, профессор Ю. В. Миронов,

доктор техн. наук, старший научный сотрудник Р. М. Юсупов,

член-корреспондент РАН, доктор техн. наук, профессор Санкт-Петербургский институт информатики и автоматизации РАН

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

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

Введение

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

Для решения данного круга задач часто применяется известный метод наименьших квадратов (МНК). Этот метод также находит широкое применение при обработке количественных результатов естественно-научных опытов, технических данных, астрономических и геодезических наблюдений и измерений.

Распространенность МНК во многом обусловлена тем, что при решении задач оценивания данным методом не требуется знания статистических характеристик ошибок измерений, которые

во многих случаях неизвестны или известны с невысокой точностью.

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

Вместе с тем МНК может быть реализован на основе использования условий оптимальности оценок вариационного типа. Некоторые вопросы обоснования и разработки соответствующей вариационной технологии рассматривались в работе авторов [10] для оценивания параметров орбитального движения КА.

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

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

Рассмотрим задачу оценивания параметров движения динамического объекта, которая заключается в наилучшем в некотором смысле определении ге-мерного вектора его исходного состояния x0 на заданный начальный момент времени г = t0 по результатам измерений, проводимых в N точках г, заданных на интервале измерений т = Т - г0.

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

X = ф(х, г), х(^) = Хо, г е [to, Т]. (1)

Измерениям подвергается т-мерный вектор

= у[хф].

Измеренное значение вектора у в момент ^ обозначим как у(^) = у; и представим модель измерений в виде

у(^) = \\1[х(Ъ)] + 8;,

I = 1(1)^; Ч е [to, Т]. (2)

Здесь 5; — т-мерный вектор случайных ошибок измерений.

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

1=]С ^ [*(4)]}, (3)

;=1

где

Рг = {У (Ч) - ^[х(Ч )]}Г ^ {у(#£) - у[х(Ч)]}, 1 =1(1)^;

Wi — симметрические матрицы весовых коэффициентов.

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

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

Вариационные условия оптимальности оценок

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

р{у(г), ^[х(#)]}={у(#)-^[х(#)]}Г W(^){y(^)-Y[x(^)]},

где у(£) и W(^) — произвольные непрерывные дифференцируемые функции, принимающие в моменты ^ соответственно значения у; и Wi (например, полиномы Лагранжа).

Тогда функционал (3) принимает вид

т N

I = / р{у(<), ¥ [х(<)]}^5(< - < )а t

<0 г=1

где о(г - ^i) — импульсная дельта-функция.

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

п N

Н(х, к, t) = ^Xj^j(х, t) + p[y(t), 1|/(х, t)]^5(# - ^)

j=l

1=1

и составим расширенный функционал

т

I*=/

#п

-V хj (#)дх-+н а#.

и 1 д# .

Рассмотрим далее первую вариацию этого функционала

87* = /-и X, №№■ 1+и

дХ:

,дН2

и -----ОХ;

• , у дt ) дХ; 1

1=1 1=1 1

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

81* = и Х j (<0 (<о) +£ * j (Т^ (Т) +

j=l j=l

Т п \дХ j дН'

-/-и

<0 j 1

дt

дxj

8xjdt.

(4)

Необходимые условия оптимальности оценки вектора x0 определяются из условия равенства нулю первой вариации функционала

51 = 0.

(5)

Поэтому из (4) и (5) получаем следующие необходимые условия оптимальности траектории x(x0, г) в виде известных уравнений Эйлера

а X j дН

- = -—, / = 1(1)га

а г дх}

и граничных условий

X j (#0) = 0, X j (Т) = 0, / = 1(1)п.

Эти условия вместе с уравнениями движения (1) образуют систему уравнений, решение которой относительно неизвестного значения x0

т

t

и определяет как оптимальную оценку x0, так и порождаемую ей траекторию x(x0, г).

Сформулируем данный результат в виде следующей теоремы.

Теорема 1. Оптимальная оценка вектора x0 и соответствующая ей оптимальная траектория доставляют решение краевой задаче для следующей системы дифференциальных уравнений:

X = ф(х, #); X = “7Х1 “^Х[у’ §(*_)

1=1

при граничных условиях

А(*о) = Х(Т) = 0.

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

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

и$)=нч)+тч),

где

АХ^) = 7Г )]}’1 =1(1)^-

ох

С учетом скачков сопряженных переменных теорему 1 можно переформулировать в следующем эквивалентном виде.

Теорема 2. Оптимальная оценка вектора x0 и соответствующая ей оптимальная траектория доставляют решение краевой задаче для следующей системы дифференциальных уравнений:

х = ф(х, t); X = -^— х (6)

дх

при граничных условиях

Х^о) = Х(Т) = 0; т

Х{4) = ЪЩ) + ^) Щ {у; - у [х(*;)]}>

ОХ ;

I = Ч№. (7)

Согласно этим условиям, для получения оптимальной оценки вектора x0 необходимо решить краевое уравнение

1^, Т) = 0,

заданное неявно на процедурах интегрирования системы (6). Для этого можно применить извест-

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

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

Так, например, если помимо дискретных измерений (2) проводятся и непрерывные измерения согласно модели

У^) = У!^)] + 51(г),

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

х = ф(х, Г); ^ = -^—X+^-W(^){yl(^) -^[х^)]};

дх дх

1(^) = 1(Т) = 0;

т

х)=х(ч )+-^т_^2щ {у. _¥)]}) 1=1(1)АГ.

с>х1

Регуляризация оптимальных статистических оценок

Как известно, многие задачи статистического оценивания могут быть отнесены к некорректным обратным задачам. Мощным средством решения таких задач является метод регуляризации, созданный А. Н. Тихоновым и развитый во многих работах.

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

1а = 1(х0) + а F(xo) (8)

принимает экстремальное значение.

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

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

Действительно, применение функционала (8) эквивалентно добавлению в функционале (3) слагаемого

Ро [х(#0)] = а F(xo).

Тогда по теореме 2 в системе условий (7) появляется дополнительное граничное условие

Х((0+) = ВД) -Шї°і.

0 дх0

Поскольку, согласно (7), А(^) = 0, то при регуляризации оптимальных оценок сопряженная система должна интегрироваться при начальном условии

>>( > = _Зро(Хр> = _а. дк0 9Хр

Соответствующий результат сформулируем в виде следующей теоремы.

Теорема 3. Регуляризированная оценка вектора Xo в исходной задаче и порождаемая ей регуляризированная траектория x(x0, Ї) доставляют решение двухточечной краевой задаче для канонической системы

х = ф(х, г); X = -^- X

дх

при граничных условиях

А(^) = -а ; А(Т) = 0; А(^+) =

дхо

т

= Жг) + ^) Щ {у; - )]}, І = 1(1)^. (9)

ОТ ;

При решении задач навигационного оценивания обычно для стабилизирующего функционала принимается выражение вида [8]

аґ(хо) = (хь -хо)тС(хь -хо),

где xb — заданный опорный вектор, близкий к истинному значению x0; С — некоторая симметрическая положительно определенная матрица.

В этом случае для определения начального значения сопряженного вектора А^0) в условиях (9) теоремы 3 получаем расчетное соотношение

А(^) =

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

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

Рассмотрим особенности применения вариационного МНК на примере решения задачи статистического оценивания параметров движения

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

где г = [х, у, z] и V = [их, иу, ] — векторы

координат и скорости движения КА соответственно; п0 = 398600,44 км3/с2 — постоянная притяжения Земли.

Проводятся прямые полные измерения элементов векторов г и v, так что модель измерений принимает вид

у! = г + 5г; у2 = v + 8^

где 8г, 8v — ошибки измерений.

Составим соответствующую сопряженную систему дифференциальных уравнений:

І — Я0

г -~Ь Г

К —у (ггГ Ж

Г

А = \АГ, А„]Т =\Хх, Ху, Хг, Хих, Хиу, Хиг]т.

В соответствии с теоремой 2 для оптимального оценивания вектора начального состояния КА % = [г, у]Т необходимо решить двухточечную краевую задачу для приведенных уравнений движения и сопряженных уравнений с учетом следующих граничных условий:

А(г0) = 0; А(Т) = 0; Хг ) = Хг (Ц) + [ух - г(^)];

XV & ) = Х„ ) + [у 2 - )], I = 1(1)

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

А^, Т) = 0.

Применение метода Ньютона приводит к следующему итерационному алгоритму:

Чк+1 = qk - Т) A(qk, Т), k = 0, 1, 2, .

дч Jk

Некоторые результаты расчетов приведены в табл. 1 и 2.

Расчеты проводились для спутника, находящегося на орбите с высотой h = 1000 км и эксцентриситетом е = 0,003. С помощью датчика случайных величин по нормальному закону распределения на мерном интервале Т = 100 с модели-

■ Таблица 1. Результаты оптимального оценивания

Оцени- ваемые параметры Точные значения Начальное приближе- ние Опти- мальные оценки Ошибки оценивания

х0, км 0 50 0,003 0,003

^0, км -7349,636 -7299,636 -7349,647 -0,011

^0, км 0 50 -0,011 -0,011

, км/с хо’ 0,89879 0,94879 0,89864 -0,00015

иУо’ км/с 0,00571 0,05571 0,00590 0,00019

v2o, км/с 7,32007 7,37007 7,32024 0,00017

■ Таблица 2. Сходимость вычислительного процесса

Номер итера- ции Ошибки оценивания

8x0, м 8^0, м 8г0, м 8 V м/с Чс’ м/с 5^ м/с

0 50000 50000 50000 50 50 50

1 6 -9 -8 0,32 0,13 0

2 3 -11 -11 0,15 0,19 0,17

3 3 -11 -11 0,15 0,19 0,17

ровалась с шагом Дt = 1 с статистическая выборка прямых измерений вектора текущего состояния КА. При этом предельные ошибки измерений задавались значениями 100 м по элементам вектора координат и 1 м/с — по элементам вектора скорости.

Литература

1. Аким Э. Л., Энеев Т. М. Определение параметров движения космических аппаратов по данным тра-екторных измерений // Космические исследования. 1963. Т. 1. № 1. С. 5-50.

2. Брандин Н. К., Разоренов Г. Н. Определение траекторий КА. — М.: Машиностроение, 1978. — 216 с.

3. Космические траекторные измерения / Под ред. П. А. Агаджанова, В. Е. Дулевича, А. А. Коростелева. — М.: Сов. радио, 1969. — 504 с.

4. Линник Ю. В. Метод наименьших квадратов и основы теории обработки наблюдений. — М.: Физмат-гиз, 1958. — 350 с.

5. Мудров В. И., Кушко В. П. Методы обработки измерений. — М.: Сов. радио, 1976. — 190 с.

6. Репин В. Г., Тартаковский Г. П. Статистический синтез при априорной неопределенности и адаптация информационных систем. — М.: Сов. радио, 1977. — 432 с.

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

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

Заключение

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

Работа выполнена при поддержке РФФИ (проект № 09-08-00259).

7. Статистические методы обработки результатов наблюдений / Под ред. Р. М. Юсупова / МО СССР, 1984. — 563 с.

Степанов М. Г. Введение в теорию смещенного оценивания параметров движения космических аппаратов по ограниченным данным / ВИККА им. А. Ф. Можайского. — СПб., 1993. — 135 с.

9. Эльясберг П. Е. Определение движения по результатам измерений. — М.: Наука, 1976. — 416 с.

10. Миронов В. И., Миронов Ю. В. Вариационный подход к статистическому оцениванию параметров орбитального движения космических аппаратов. — СПб.: ВИКУ им. А. Ф. Можайского, 2002. — 166 с.

11. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. — М.: Наука, 1979. — 288 с.