авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 |

«ISSN 2075-6836 Ф е д е ра л ь н о е го с уд а р с т в е н н о е б юд ж е т н о е у ч р е ж д е н и е н а у к и институт космических исследований российской академии наук (ики ран) ...»

-- [ Страница 2 ] --

a(1- e 2 ) p = a(1- e) » 147,1106 км rp = = 1+ e 1+ e — расстояние до перигелия Земли. Аналогично находим, что расстояние до афелия Земли равно ra = a(1 + e) » 152,1106 км.

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

4.3. Параметры орбиты спутника p ra =. (4.19) a rp Д о к а з а т е л ь с т в о. Отношение расстояния до апоцентра к расстоянию до перицентра равно p (1- e) (1 + e) ra = =.

p (1 + e) (1- e) rp С другой стороны, µ (1 + e) p (1 + e) p = =, (1- e) a µ (1- e) p что и доказывает теорему.

Пример 4.4. Космический аппарат на высоте 300 км имел скорость 10 км/с, причем вектор скорости был направлен па раллельно поверхности Земли. Рассчитать скорость аппарата в апогее, если известно, что расстояние от центра Земли до апо гея составляет около 370 тыс. км.

Р е ш е н и е. Тот факт, что на высоте 300 км вектор скоро сти спутника был направлен параллельно поверхности Земли, означает, что в этой точке орбита спутника имеет перигей, т. е.

расстояние до перигея равно 6671 км. Тогда по правилу рычага величина скорости спутника в апогее составляет r a = a p = 554,6 км с-1.

rp 4.3.2. круговые орбиты Круговые орбиты являются частным случаем эллиптических.

Выше было показано, что для круговых орбит выполняются ус ловия µ e = 0, h 0, = кр (r ) =.

r В отличие от эллиптической, круговая орбита полностью определяется одним параметром — чаще всего используется 52 4. невоЗмущенное движение сПутника в Плоскости орбиты радиус круговой орбиты или высота спутника над поверхностью притягивающего тела.

4.3.3. гиперболичеСкие орбиты При движении спутника вокруг притягивающего центра по ги перболической орбите выполняются соотношения 2µ e 1, h 0,.

r Рассмотрим такое движение подробнее (рис. 8).

Пусть, как и ранее, S — текущее положение спутника, дви жущегося вокруг притягивающего центра C;

O — центр гипер болы;

P — перицентр орбиты. Истинная аномалия J, так же как и для эллиптических орбит, представляет собой угол между текущим радиус-вектором спутника r =CS и направлением CP из притягивающего центра на перицентр. Однако для ги перболических орбит истинная аномалия изменяется в пределах -J lim J J lim, где J lim — предельное значение угла J, при котором знамена тель правой части уравнения орбиты спутника (4.3) обращается в нуль, т. е. расстояние от спутника до притягивающего центра неограниченно увеличивается:

1 1 + e cos J lim = 0 J lim = arccos - = - arccos.

e e Фокальный параметр p и эксцентриситет e гиперболической орбиты определяются так же, как и для эллиптических орбит:

OC p = CD, e =.

OP Для гиперболических орбит определяют действительную полуось a и мнимую полуось b:

a = OA = OP, b = OB.

Расстояние до перицентра равно p rp = = a(e -1).

1+ e 4.3. Параметры орбиты спутника Рис. 8. Параметры гиперболической орбиты Можно доказать справедливость следующих соотношений, аналогичных соотношениям (4.15)–(4.17):

p a=, (4.20) e - b2 = a2 (e 2 -1), (4.21) b p = a(e 2 -1) =. (4.22) a Эти соотношения позволяют полностью определить гипер болическую орбиту по любым двум из следующих параметров:

• фокальный параметр p;

• эксцентриситет e;

• действительная полуось a;

• мнимая полуось b;

• расстояние до перицентра rp.

4.3.4. параболичеСкие орбиты Параболические орбиты являются предельным случаем эллип тических и гиперболических орбит, соответствующим выполне нию равенств µ e =1, h = 0, =.

r 54 4. невоЗмущенное движение сПутника в Плоскости орбиты Уравнение параболической орбиты имеет вид p r=, 1 + cos J откуда непосредственно следует, что p rp =.

Из уравнения (4.15) или (4.20) следует, что для параболиче ских орбит можно считать a=.

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

4.4. время ПрОхОждения сПутника Через заданную тОЧку Орбиты 4.4.1. эллиптичеСкие орбиты Рассмотрим движение спутника S вокруг притягивающего цен тра C по эллиптической орбите с центром в точке O, большой полуосью a и эксцентриситетом e (рис. 9). Проведем окруж ность радиуса a с центром в точке O. Проведем через точку S перпендикулярно линии апсид отрезок прямой от линии апсид до пересечения с окружностью. Точку пересечения обозначим через S, а величину угла COS — через E.

Угловая величина E называется эксцентрической аномалией, она наряду с истинной аномалией J характеризует смещение спутника относительно перицентра. Отличие заключается в том, что для эксцентрической аномалии угловое смещение спутника отсчитывается из центра орбитального эллипса O, а для истинной аномалии — из притягивающего центра C.

Можно показать [Сихарулидзе, 2011], что истинная аномалия J связана с эксцентрической аномалией E соотношением 1+ e E tg = tg.

2 1- e 4.4. время прохождения спутника через заданную точку орбиты Рис. 9. Иллюстрация к выводу уравнения Кеплера для эллиптических орбит Рассмотрим в плоскости орбиты декартову систему коорди нат Oxo yo с центром, совпадающим с центром орбитального эл липса O, и осью абсцисс, направленной на перицентр P (см.

рис. 9). Пусть текущее положение спутника определяется векто ром ro = OS, имеющим в этой системе координаты r x ro = o.

r y o Тогда нетрудно убедиться, что rx = a cos E, o ry = b sin E.

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

Обозначим через tp время прохождения спутником перицен тра P. Тогда можно показать [Охоцимский, Сихарулидзе, 1990], что моменты времени t и tp связаны соотношением a t -t p = (E - e sin E ), (4.23) µ которое называется уравнением Кеплера для эллиптических ор бит. В этом уравнении — гравитационный параметр притяги вающего центра.

56 4. невоЗмущенное движение сПутника в Плоскости орбиты Рассмотрим также величину µ n=, a называемую средним движением, и величину M = n(t – tp), назы ваемую средней аномалией. Тогда уравнение (4.23) примет вид E - e sin E = M. (4.24) Теорема 4.2. При известных значениях средней аномалии M и эксцентриситета e уравнение Кеплера (4.24) всегда имеет един ственное действительное решение относительно эксцентрической аномалии E.

Д о к а з а т е л ь с т в о. Рассмотрим функцию F(E) = E – – sin E – M. Эта функция, очевидно, является непрерывной.

Кроме того, dF F (-) = -, F (+) = +, =1- e cos E 0, dE т. е. функция F(E) — монотонно возрастающая, и, следователь но, уравнение F(E) = 0 всегда имеет единственное решение.

Решение уравнения (4.24) может быть найдено приближен но — например, с использованием схемы метода последователь ных приближений [Эльясберг, 1965]:

E n = M + e sin E n-1, где En – 1 и En — предыдущее и последующее приближения иско мой величины E. В качестве первого приближения можно взять значение E1 = M.

Из уравнения (4.23) следует, что время перелета спутника между двумя точками орбиты с эксцентрическими аномалия ми E1 и E2 равно a ( E2 - E1 - e (sin E2 - sin E1 )).

t2 - t1 = µ Определение 4.5. Периодом обращения спутника при невоз мущенном движении вокруг притягивающего центра называет ся промежуток времени между двумя последовательными про хождениями спутником одной и той же точки своей орбиты.

Пусть tp — момент прохождения спутника через пери центр, t p — момент следующего прохождения через перицентр, 4.4. время прохождения спутника через заданную точку орбиты соответствующий значению эксцентрической аномалии E = 2.

Тогда период обращения спутника T, согласно уравнению (4.23), будет равен a T = t p - t p =. (4.25) µ Как следует из последнего соотношения, период обраще ния спутника зависит только от величины большой полуоси ор биты и от гравитационного параметра притягивающего центра.

Пример 4.5. Космический аппарат «Восток-2» имел на пер вых витках высоту перигея 244 км и высоту апогея 183 км. Най ти период обращения этого спутника.

Р е ш е н и е. Принимая радиус Земли равным RE = 6371 км, найдем расстояние до перигея и расстояние до апогея орбиты спутника:

rp = RE + hp = 6371 +183 = 6554 км, ra = RE + ha = 6371 + 244 = 6615 км.

Тогда величина большой полуоси орбиты данного спутника равна ra + rp a= = 6584,5 км.

Теперь, принимая значение гравитационного параметра Земли равным = 3,986·105 км3·с–2, найдем искомый период об ращения:

a T = 2 = 5317 с = 1 ч 28 мин 37 с.

µE Пример 4.6. Рассчитать радиус орбиты и скорость, которые должен иметь спутник, вращающийся по круговой орбите во круг Земли в плоскости экватора, если требуется, чтобы этот спутник все время находился над одной и той же точкой поверх ности Земли.

Р е ш е н и е. Для того чтобы спутник находился постоян но над одной и той же точкой поверхности Земли, необходимо, чтобы период T его обращения по орбите совпадал с периодом вращения Земли вокруг своей оси. Примем вначале для про стоты значение этого периода равным 24 ч или 86 400 с. Тогда 58 4. невоЗмущенное движение сПутника в Плоскости орбиты величина большой полуоси орбиты спутника, которая для кру говой орбиты совпадает, очевидно, с ее радиусом, должна со ставлять T a = 3 µ E = 42 241,1 км, т. е. высота такой орбиты будет равна 35 870 км. Скорость дви жения спутника по этой орбите, согласно формуле (4.13), составит µE » 3,072 км с-1.

= кр (a) = a Заметим, что если рассмотреть более точное значение пе риода обращения Земли, равное 23 ч 56 мин 4 с или 86 164 с, то радиус круговой орбиты спутника получится равным 42 164 км, т. е. высота орбиты будет равной 35 793 км. Скорость движе ния спутника по круговой орбите такого радиуса составит 3,075 км·с–1. Как видим, изменение требуемого значения пери ода обращения на 236 с привело к изменению высоты орбиты спутника на 77 км.

Замечание 4.3. Если рассмотреть два спутника, движущиеся вокруг общего притягивающего центра по эллиптическим ор битам со значениями большой полуоси a1 и a2 и имеющие пе риоды обращения T1 и T2 соответственно, то из формулы (4.25) следует соотношение T12 a =, T22 a выражающее третий закон Кеплера.

4.4.2. гиперболичеСкие орбиты Теперь рассмотрим случай, когда спутник S движется вокруг притягивающего центра C по гиперболической орбите с цен тром в точке O, действительной полуосью a и эксцентрисите том e (рис. 10). Рассмотрим новую переменную H такую, что E = iH, H = -iE, 4.4. время прохождения спутника через заданную точку орбиты Рис. 10. Иллюстрация к выводу уравнения Кеплера для гиперболических орбит где i — мнимая единица, определяемая как решение уравнения i 2 = –1. Величина H является аналогом эксцентрической анома лии E для случая гиперболического движения. Истинная анома лия J связана с H соотношением [Сихарулидзе, 2011] e +1 H J tg = th.

2 e -1 Рассмотрим, как и в предыдущем случае, в плоскости орби ты декартову систему координат Oxo yo с центром, совпадающим с центром орбитального эллипса O, и осью абсцисс, направлен ной из точки O на перицентр P (см. рис. 10). Обозначим через ro радиус-вектор, определяющий текущее положение спутника в этой системе координат. Тогда компоненты rx, ry этого ра диус-вектора будут определяться по формулам o o rx = -a ch H, ry = b sh H o o (эти формулы справедливы для случая, когда рассматривается движение по левой ветви гиперболы, как показано на рис. 10).

Если в некоторый момент времени t спутник находится в точке орбиты, эксцентрическая аномалия которой равна H, и tp — время прохождения спутника через перицентр P, то мож но показать [Охоцимский, Сихарулидзе, 1990], что моменты времени t и tp связаны соотношением 60 4. невоЗмущенное движение сПутника в Плоскости орбиты a t -t p = (e sh H - H ), (4.26) µ которое называется уравнением Кеплера для гиперболических орбит.

Последнее соотношение приводимо к виду e sh H - H = M, (4.27) где M = n(t – tp) — средняя аномалия;

n — среднее движение.

4.5. ФОрмула ламберта Рассмотрим движение спутника S по эллиптической орбите с большой полуосью a вокруг притягивающего центра C с из вестным гравитационным параметром (рис. 11).

Пусть S1 и S2 — две различные точки орбиты, s — расстоя ние между ними, r1 и r2 — расстояния от этих точек до притя гивающего центра. Тогда время перелета спутника от точки S до точки S2 может быть найдено по формуле Ламберта* [Балк, 1965;

Иванов, Лысенко, 2004] * Иоганн Генрих Ламберт (1728–1777) — выдающийся немецкий астроном, физик, математик, философ. Разработал принципы фото метрии — основного метода наблюдательной астрономии. С помощью своих методов Ламберт смог точно оценить относительную яркость Луны, ослабление света в земной атмосфере, также на основе фотоме трии им проводилась оценка межзвёздных расстояний. Занимался ис следованиями орбит комет. Ламберту принадлежит и космологическая концепция развития Вселенной, во многом предвосхитившая научные открытия будущего.

Ламберт внес также огромный вклад в развитие математики, осо бенно тригонометрии, теории конических сечений и гиперболических функций, стал одним из основателей неевклидовой геометрии. Он впервые доказал иррациональность чисел и e, составил таблицу про стых чисел до 102 000. Значительный вклад Ламберт внес и в развитие философии, он состоял в переписке с И. Кантом, высоко ценившим его ум и способности.

Ламберт был членом нескольких научных обществ — литератур ного и физического в Куре, членом-корреспондентом научного обще ства в Гёттингене, членом Баварского научного общества и Берлинской Академии наук. Скончался в возрасте 49 лет.

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

4.5. Формула ламберта Рис. 11. Иллюстрация к формуле Ламберта a ( - sin 1 ) - ( 2 - sin 2 ), t2 - t1 = µ где 1 и 2 — корни уравнений r1 + r2 + s cos 1 =1-, 2a r1 + r2 - s cos 2 =1-.

2a Очевидно, что эти уравнения имеют бесконечное множе ство решений. Приведем правила, позволяющие однознач но выбрать числа 1 и 2 для подстановки в формулу Ламберта [Балк, 1965].

1. Если спутник, двигаясь по эллиптической орбите, проходит дугу S1S2 после прохождения перицентра и до прохожде ния следующего апоцентра (т. е. так, как показано на рис. 11), то следует выбирать числа 1 и 2 из условия 0 1, 0 2.

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

62 4. невоЗмущенное движение сПутника в Плоскости орбиты Таблица 1. Выбор чисел 1 и 2 при различных случаях расположения точек S1 и S2 на орбитальном эллипсе Возможные случаи 1 1. Сегмент не содержит 1 ни C, ни C 2. Сегмент содержит C 2 - 1 - и C 3. Сегмент содержит C - и не содержит C 4. Сегмент не содержит 2 - 1 C, но содержит C Пример 4.7. Космический корабль совершает перелет между орбитами Земли и Марса. Можно принять, что орбита перелета является эллиптической, причем в одном из ее фокусов распо ложено Солнце, а расстояния перигелия и афелия данной ор биты составляют 120 и 240 млн км соответственно. Точка нача ла перелета расположена на расстоянии 150 млн км от Солнца, точка окончания перелета — на расстоянии 228 млн км от Солн ца. Сколько времени должен занять такой перелет?

Р е ш е н и е. Для нахождения времени перелета между дву мя точками орбиты удобно воспользоваться формулой Ламбер та. Для этого по известным значениям расстояний до афелия и перигелия сначала найдем значения a, e и p:

ra + rp ra + rp a= =180 млн км, e = =, 2 ra - rp p = a (1- e 2 ) =160 млн км.

4.5. Формула ламберта Теперь найдем истинные аномалии J1 и J 2 точек S1 и S2.

Из уравнения орбиты (4.3) находим:

p-r J = arccos, er откуда J1 = 1,369 рад, J 2 = 2,679 рад. Тогда расстояние s между точками S1 и S2 равно s = r12 + r22 - 2r1r2 cos(J 2 -J1 ) » 238,3 млн км.

Найдем теперь множители 1 и 2:

r1 + r2 + s r1 + r2 - s cos 1 =1- » -0,71, cos 2 =1- » 0,61, 2a 2a откуда 1 = 2,363 рад, sin 1 = 0,702;

2 = 0,912 рад, sin 2 = 0,791.

Подставляя найденные величины в формулу Ламберта, найдём (180 106 ) (2,363 - 0,702) - (0,912 - 0,791) = t2 - t1 = 132,5 =10,21406 106 с =118,22 сут.

5. невозМущенное движение сПутниКа в инерциальной систеМе Координат 5.1. Элементы Орбиты До сих пор нами рассматривалось движение спутника в плоско сти орбиты и ставилась задача определения положения спутни ка на орбите по тем или иным известным величинам. Перейдем теперь к рассмотрению задачи определения пространственного положения спутника в какой-либо заданный момент времени относительно некоторой инерциальной системы координат.

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

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

5.1. Элементы орбиты Рассмотрим движение спутника S по эллиптической орбите вокруг притягивающего центра C (рис. 12). Зададим в простран стве некоторую инерциальную систему отсчета Cxyz, начало ко торой примем совпадающим с притягивающим центром C. Ор бита спутника будет пересекать плоскость Cxy в двух точках.

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

Определение 5.2. Долготой восходящего узла орбиты на зывается угол между положительным направлением оси аб сцисс Cx и направлением из притягивающего центра на вос ходящий узел (см. рис. 12). Долгота восходящего узла меняется в пределах от 0 до 2.

Определение 5.3. Наклонением орбиты называется угол i между плоскостью орбиты и координатной плоскостью Cxy. На клонение орбиты меняется в пределах от 0 до. Если для задан ной орбиты i, то движение спутника называется прямым, если же i, то движение называется обратным. При i = 2 орбита спутника называется экваториальной, при i = — по лярной.

Рис. 12. Элементы орбиты спутника: наклонение i, долгота восходящего узла и аргумент перигея 66 5. невоЗмущенное движение сПутника в инерциальной системе координат Определение 5.4. Аргументом перицентра орбиты называет ся угол между направлениями из притягивающего центра на восходящий узел и на перицентр орбиты. Аргумент перицентра отсчитывается в направлении движения спутника и меняется в пределах от 0 до 2.

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

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

• наклонение орбиты i, долготу восходящего узла и аргу мент перицентра, определяющие пространственное поло жение плоскости орбиты;

• большую полуось a и эксцентриситет e, характеризующие форму и размер орбиты;

• время = tp прохождения спутника через перицентр.

Определение 5.5. Элементы орбиты a, e, i,,, образуют полную систему и называются кеплеровыми элементами невозму щенного движения.

Замечание 5.1. Полную систему элементов орбиты, очевид но, можно выбрать множеством способов. В частности, вместо большой полуоси орбиты a можно рассмотреть фокальный па раметр p, вместо аргумента перицентра или вместо времени прохождения перицентра — угол u = J +, называемый аргу ментом широты, вместо параметра также часто рассматривают среднюю аномалию M и т. д. Особенностью кеплеровой систе мы элементов является то, что с ее помощью движение спутника определяется шестью постоянными величинами, не зависящими от времени.

5.2. ОПределение ПОлОжения и скОрОсти сПутника ПО Элементам Орбиты Покажем, как по известным значениям кеплеровых элемен тов орбиты определить компоненты rx, ry и rz радиус-вектора r 5.3. определение элементов орбиты спутника по положению и скорости и компоненты x, xy, z вектора скорости спутника в заданный момент времени t [Эльясберг, 1965].

1. Находим фокальный параметр p = a(1 – e2 ) и среднюю ано µ малию M = (t - ).

a 2. Из уравнения Кеплера E – e sinE = M какой-либо численной процедурой находим эксцентрическую аномалию E.

1+ e E 3. По формуле = tg находим значение истинной 1- e аномалии J, а затем — аргумент широты u = J +.

p 4. Из уравнения орбиты r = находим расстояние 1 + e cos J до спутника r.

5. Координаты спутника определяем из следующих соотно шений:

rx = r (cos cos u - sin sin u cos i ), ry = r (sin cos u + cos sin u cos i ), rz = r sin u sin i.

6. Определяем величину радиальной и трансверсальной компонент вектора скорости:

µ µ r = e cos J, n = (1 + e cos J).

p p 7. Компоненты вектора скорости определяем по формулам x = r (cos cos u - sin sin u cos i ) - n (cos sin u + sin cos u cos i ), y = r (sin cos u + cos sin u cos i ) - n (sin sin u - cos cos u cos i ), z = r sin u sin i + u cos u sin i.

5.3. ОПределение ЭлементОв Орбиты сПутника ПО ПОлОжению и скОрОсти Теперь покажем, как вычислить кеплеровы элементы орбиты, если в некоторый заданный момент времени t известны поло 68 5. невоЗмущенное движение сПутника в инерциальной системе координат жение спутника, определяемое радиус-вектором r = (rx, ry, rz ), и вектор его скорости = ( x, y, z ) [Эльясберг, 1965].

1. Вычислим вспомогательные величины (компоненты век торной постоянной площадей c) cx = ry z - rz y, c y = rz x - rx z, cz = rx y - ry x, 2 2 c = c = cx + c y + cz и найдем долготу восходящего узла, а также наклонение орбиты i по формулам c c = - x, cos i = z, 0 i.

cy c 2. Введем вспомогательный угол (угол между вектором ско рости и нормалью к радиус-вектору спутника) и найдем его из условия r sin =, -, 2 r 2 2 2 2 2 где r = r = rx + ry + rz, = = x + y + z.

Тогда, очевидно, r 2 2 - (r ) cos =.

r r Вычислим также величину k =.

µ 3. Большую полуось a и эксцентриситет орбиты e найдем по формулам r a=, e = 1- k(2 - k )cos.

2-k 4. Истинную аномалию J определим из условий k cos 2 - k sin cos sin =, cos =, e e аргумент широты u — из условий rx cos + ry sin rz cos u =, sin u =, r sin i r 5.3. определение элементов орбиты спутника по положению и скорости и затем аргумент перицентра — из равенства = u -J.

5. Время прохождения перицентра найдем по формуле a =t- (E - e sin E ), µ где эксцентрическая аномалия E определяется по формуле 1- e J E tg = tg.

2 1+ e 6. возМущенное движение сПутниКа 6.1. ПОнятие О вОзмущеннОм движении.

метОд Оскулирующих ЭлементОв Рассмотренная выше ограниченная задача двух тел описы вает, вообще говоря, идеальный случай, когда спутник движется при наличии только гравитационного взаимодействия с притя гивающим телом. При описании реального движения спутника приходится учитывать целый ряд других факторов, оказываю щих влияние на характер его орбиты. Такие факторы носят на звание возмущающих факторов орбиты, возмущающих воздей ствий или, короче, возмущений. К основным типам возмущений относятся [Суханов, 2010]:

• гравитационные возмущения: влияние несферичности притягивающего тела, влияние притяжения других небес ных тел;

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

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

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

6.1. Понятие о возмущенном движении. метод оскулирующих элементов Уравнение возмущенного движения спутника вокруг при тягивающего тела в самом общем виде записывается следую щим образом:

µ = - r +, (6.1) r r где — суммарное ускорение, получаемое спутником в резуль тате действия всех возмущающих факторов и называемое воз мущающим ускорением. Если возмущающее ускорение отлич но от нуля, то движение, очевидно, уже не будет совершаться по одному из конических сечений. Однако в каждый момент времени t по известным значениям r(t) и (t) при помощи опи санной выше процедуры можно подобрать значения кепле ровых элементов орбиты, соответствующих невозмущенному движению в данной точке пространства с данной скоростью.

Иначе говоря, в каждый момент времени реальная орбита дви жения спутника будет касаться некоторой орбиты невозму щенного движения, соответствующей найденным значениям кеплеровых элементов орбиты. Таким образом мы можем полу чить зависимости кеплеровых элементов от времени, т. е. рас сматривать функции a(t), e(t), i(t), (t), (t), (t), описывающие возмущенное движение спутника. Если эти функции известны, то положение и скорость спутника в любой момент времени t могут быть определены аналогично случаю невозмущенного движения.

Определение 6.2. Рассмотренные функции a(t), e(t), i(t), (t), (t), (t), представляющие собой элементы непрерывно меняющейся со временем орбиты невозмущенного движения, называются оскулирующими элементами, а сама такая орбита — оскулирующей орбитой. Описание движения спутника при по мощи оскулирующей орбиты называется методом оскулирующих элементов*.

Введем обозначения:

r x =, f ( x ) = µ, g( x, t ) =, r r где x — вектор состояния, f ( x ) — правая часть векторной фор мы уравнения невозмущенного движения спутника (3.1).

* Идея и разработка метода оскулирующих элементов принадле жат Ж. Лагранжу.

72 6. воЗмущенное движение сПутника Тогда уравнение возмущенного движения в векторной фор ме запишется в виде x = f ( x ) + g( x, t ).

Обозначим через q вектор оскулирующих элементов орбиты q = q( x, t ) = (a(t ), e(t ), i(t ), (t ), (t ), (t )).

Тогда q q f ( x ) + g( x ).

q= x= x x Если g(x, t) = 0 т. е. движение невозмущенное, то q = const и, следовательно, q f (x) = 0.

x Тогда q q q= g( x, t ) =.

x v Введем орбитальную систему координат O (рис. 13).

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

Пусть векторы и в орбитальной системе имеют коор динаты S =, = T.

W Очевидно, что = r — радиальная скорость, = n — трансверсальная скорость, = 0. Также нетрудно заметить, что =W (в правой части уравнения возмущенного движения (6.1) первое слагаемое представляет собой правую часть уравне ния невозмущенного движения (3.1), т. е. вектор, лежащий в плоскости орбиты;

следовательно, третья компонента этого вектора в орбитальной системе координат равна нулю).

6.1. Понятие о возмущенном движении. метод оскулирующих элементов Рис. 13. Векторы r и в орбитальной системе координат Тогда q q q q= W.

S+ T+ vr vn v Подставляя в это соотношение выражения для кеплеровых элементов орбиты, в результате можно получить следующую си стему уравнений для оскулирующих элементов, которая назы вается уравнениями Ньютона – Лагранжа:

p 2a p a= e sin JS + T, r µ 1- e r p r sin JS + 1 + cos JT + eT, e= p p µ r i = cos uW, µp (6.2) r sin u = W, µp sin i 1 r p cos J r sin u S + 1 + sin JT = W, - p tg i e e p µ r p = eN sin JS - cos JS + NT, µe r 74 6. воЗмущенное движение сПутника J 2 p2 cos J d J где u = J + — аргумент широты, N = (1 + e cos J)2.

r2 Часто на практике вместо элемента орбиты рассматрива ют истинную аномалию J. Тогда последнее из уравнений Нью тона – Лагранжа* (6.2) заменяется более простым уравнением p µ cos J 1 r J= S - 1 + sin JT.

2+ e e p µ r Если возмущающее ускорение не зависит явно от време ни, то правые части уравнений Ньютона – Лагранжа также не зависят явно от t. Поэтому в качестве независимой перемен ной в этом случае как правило выбирают аргумент широты u = + J учитывая, что r3 sin u µp W.

u= 1 tg i r µp * Жозеф Луи Лагранж (1736–1813) — великий французский ма тематик, астроном и механик итальянского происхождения. Наряду с Эйлером признан лучшим математиком XVIII века. Особенно про славился исключительным мастерством в области обобщения и синтеза накопленного научного материала.

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

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

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

Лагранж состоял членом нескольких Академий наук, причем в Берлинскую был избран в возрасте 20 лет (!) по рекомендации Л. Эй лера и Д’Аламбера.

Имя Лагранжа внесено в список 72 величайших учёных Франции, помещённый на первом этаже Эйфелевой башни.

6.1. Понятие о возмущенном движении. метод оскулирующих элементов Обозначим =. (6.3) r sin u 1- W µp tg i µp Тогда u = 2, и система уравнений (6.2) преобразуется r к виду da a r 2 2a p == e sin JS + T, r µ 1- e du u de e r r r sin JS + 1 + cos JT + eT, == p du u p µ r di i == cos uW, du u µp (6.4) d r 3 sin u == W, µp sin i du u d r r r sin u - cos JS + e 1 + sin JT - e == W, p tgi du u µe p d r 4 p == eN sin JS - cos JS + NT, du u µe µp r где определяется формулой (6.3).

Начальные условия для уравнений Ньютона – Лагранжа определяются по начальному положению спутника согласно формулам, приведенным в п. 5.3.

7. задача оПределения движения 7.1. мОдель Оценивания Каждый человек в своей жизни и общество в целом ежеминут но решают задачи нахождения оценок (или, короче, оценивания) параметров тех или иных объектов (систем). При этом иссле дуемая система описывается с помощью совокупности некото рых параметров, характеризующих ее и изменяющихся по ка ким-либо законам. Бытовые оценки параметров носят обычно описательный характер: богатый – бедный, высокий – невысо кий, умный — неумный и т. п. Способ нахождения оценок па раметров называют методом оценивания. При использовании статистических методов оценивания рассматривается матема тическая модель, описывающая поведение системы с некоторой точностью, и определяются количественные оценки одного или нескольких параметров этой модели.

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

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

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

В соответствии с этим будем полагать, что модель измерений представляется в виде y = f ()+. (7.1) n m n Здесь y — вектор измерений, f : — вектор-функ ция, определяющая модель, — суммарный вектор ошибок мо дели и измерений.

7.2. несмещенный алгОритм Оценивания Если n = m, то естественно искать оценку вектора как реше ние системы уравнений y = f ().

В случае, когда число измерений равно числу оцениваемых параметров (n = m) и существует однозначная обратная функция f -1( y ), решение последних уравнений имеет вид = f -1 ( y ).

Если положить = 0, то определенный таким образом оце ниватель при n = m обладает свойством =, т. е. оценка пара метра совпадает с его истинным значением при отсутствии ошибок исходных данных. Такое свойство будем называть не смещённостью оценивателя.

Рассмотрим теперь понятие несмещенного оценивателя при n m.

Теорема 7.1 (о несмещенном оценивателе). Пусть — неко торая векторная норма. Оцениватель = arg min y - f ( x ) (7.2) x является несмещённым.

78 7. Задача оПределения движения Доказательство. Если в (7.1) положить = 0, то минимум в (7.2) достигается при x = — в этом случае норма равна нулю, а при x норма положительна. Поэтому для рассмотренного оценивателя при = 0 имеем =, что и доказывает его несме щённость.

Рассмотрим два вида нормы.

1. Пусть норма евклидова, т. е.

m xi2, x=x x x = i = m где xi — компоненты вектора x. Тогда оцениватель (7.2) имеет вид МНК = arg min y - f ( x ) (7.3) x и называется методом наименьших квадратов (МНК). Оценку МНК вектора, полученную с использованием оценивате ля (7.3), будем называть оценкой метода наименьших квадратов или, короче, оценкой наименьших квадратов. В более общем слу чае, когда норма определяется выражением x=x x Wx w для некоторой положительно определенной матрицы W, кото рую далее будем называть весовой, оцениватель имеет вид = arg min y - f ( x ) W y - f ( x ) ( ) ( ) МНК W x и называется обобщенным методом наименьших квадратов (ОМНК). Соответствующую оценку МНК W будем называть оценкой обобщенного метода наименьших квадратов или, короче, оценкой ОМНК.

2. Пусть норма имеет вид m xi.

x=x i = Тогда оцениватель (7.2) имеет вид n МНК = arg min yi - fi ( x ) x i = и называется методом наименьших модулей (МНМ), а оценка МНК — оценкой наименьших модулей.

7.4. одномерная линейная модель Качество оценки характеризуется величиной ее ошибки -, а точнее — статистическими характеристиками этой ошибки. Обычно эти характеристики удается найти лишь для линейной модели измерений.

7.3. линеаризация мОдели измерений Пусть 0 — некоторое номинальное значение вектора. Пред полагая разность – 0 достаточно малой, проведем линеариза цию модели измерений, оставляя только слагаемые первого по рядка. Тогда она запишется в виде y - f (0 ) = H ( - 0 ) +, где f () H = — матрица размерности nm. Окончательно линеаризованная модель записывается в виде y = H +, 0 где y = y - f ( ) + H — приведенный вектор измерений. Эта приближенная модель обычно используется для вычисления оценки =, относительно которой модель снова линеаризу ется. Процесс последовательных приближений продолжается до тех пор, пока оценки s и s+1 не окажутся достаточно близ s+ кими друг к другу. Тогда оценка = принимается за оценку в первоначальной нелинейной модели.

7.4. ОднОмерная линейная мОдель Одномерная линейная модель имеет вид yi = + i, i =1,, n, где — неизвестный скалярный параметр.

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

80 7. Задача оПределения движения n = y yi, n i = а оценка наименьших модулей — выборочную медиану вектора y = ( y1, y2, yn ).

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

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

l = b, b m (далее параметр l будем называть контролируемым), а также со вокупности подобных параметров l = B, s где l — вектор контролируемых параметров;

B — матрица размерности ms;

b и B известны. Если по результатам измере ний найдена оценка вектора, то оценки указанных параме тров можно найти по формулам l l = b, = B.

Для линейной модели оценка ОМНК, соответствующая ве совой матрице W, имеет вид = arg min ( y - H x ) W ( y - H x ) = Y y, x где Y (HWH )-1 HW.

Как видим, если выполняется условие det(HWH ) 0, то ОМНК представляет собой линейный относительно вектора y несмещенный алгоритм (доказательство несмещённости предо ставляется читателю). При этом = B = X y, X = B Y = B (HWH )-1 HW.

l 7.5. линейный несмещенный алгоритм и метод наименьших квадратов Отвлечемся от метода наименьших квадратов и рассмотрим сразу линейную оценку вектора l:

= X y, l где X — матрица размерности ns определяющая эту оценку.

Условие несмещённости линейного алгоритма имеет вид = l HX = B.

l Тогда ошибка оценки будет равна l - l = X.

l Таким образом, ОМНК дает линейную несмещенную оцен ку любого векторного параметра l = B (а значит, и любого скалярного параметра l = b ). Точнее, каждой весовой матри це W однозначно соответствует одна матрица X для линейного оценивания вектора l = B (или один вектор x для оценивания параметра l = b ).

Рассмотрим соотношение между оценкой ОМНК и линей ной несмещенной оценкой с другой точки зрения. Пусть из вестна матрица X размерности ns, удовлетворяющая условию несмещённости HX = B (в скалярном виде это условие имеет упомянутый выше вид Hx = b). Как найти весовую матрицу W такую, чтобы оценка ОМНК вектора l = B, соответствующая этой матрице, имела вид = X y ? Эта задача эквивалентна задаче нахождения неот l рицательно определенной матрицы W, являющейся решением матричного уравнения X = B (HWH )-1 HW.

Такая задача рассматривалась различными авторами (см., например, [Rao, 1975], [Бахшиян, 1983]). Было показано, что существует целое множество весовых матриц, для каждой из ко торых оценка наименьших квадратов совпадает с заданной не смещенной оценкой. В явном виде это множество в настоящем курсе мы приводить не будем.

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

82 7. Задача оПределения движения Замечание 7.1. Рассмотренный выше метод наименьших модулей не является линейным несмещенным оценивателем.

Поэтому он, вообще говоря, не эквивалентен никакой оценке ОМНК.

7.6. выЧисление тОЧнОсти ОценОк При известнОй кОвариациОннОй матрице ОшибОк измерений. теОрема гаусса – маркОва Пусть вектор ошибок представляет собой случайный вектор с математическим ожиданием E() равным нулевому вектору, и известной матрицей ковариаций K E ( ). Тогда дисперсия ошибки оценки l = x y скалярного параметра l = b равна D(l ) = x Kx, где l = l - l. Данная формула является универсальной, так как она справедлива для любого оценивателя x и используется для вычисления не только дисперсии оценок скалярных параме тров, но и для определения ковариационной матрицы оценок векторных параметров. В этом случае возникает вопрос о на хождении оптимального оценивателя, для которого дисперсия оценки будет минимальной. Можно показать, что минимум дисперсии по оценивателю x при условии несмещённо сти Hx = b достигается при x = b (HK-1H )-1 HK-1.

Это устанавливается с помощью метода множителей Лагранжа.

Можно заметить, что оцениватель получился такой же, как если бы использовался ОМНК с весовой матрицей W = K –1. Таким образом, с точки зрения минимизации дисперсии оценки лю бого скалярного параметра оптимальным является ОМНК с ве совой матрицей, обратной матрице ковариаций ошибок измере ний. Оптимальная дисперсия при этом равна D * = b K W b, где K W (HK-1H )- — ковариационная матрица оценки вектора.

7.6. вычисление точности оценок при известной ковариационной матрице ошибок… Этот известный результат называется теоремой Гаусса*– Маркова**.

Замечание 7.2. Можно показать, что матрица K - K W является неотрицательно определенной, если K — ковариаци онная матрица оценки, полученной любым линейным несме щенным алгоритмом.

** Карл Фридрих Гаусс (1777–1855) — немецкий математик, меха ник, физик и астроном. Считается одним из величайших математиков всех времён, «королём математиков». Число научных открытий Гаусса настолько велико, что ему приходилось вести для их записи специаль ный дневник.

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

Им же был изобретён электрический телеграф. В трактате Гаусса «Тео рия движения небесных тел» изложена каноническая теория учёта воз мущений орбит.

Гаусс знал в совершенстве множество языков, в 62 года стал из учать русский, чтобы прочитать труды Лобачевского. До старости при вык большую часть вычислений производить в уме.

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

** Андрей Андреевич Марков (1856–1922) — русский математик, академик, внёсший большой вклад в развитие теории вероятностей, математического анализа и теории чисел. С его именем связаны, в част ности, такие понятия теории случайных процессов как марковский про цесс, цепь Маркова, неравенство Маркова.

8. КлассичесКие задачи Планирования эКсПериМента 8.1. усреднение мОдели измерений m Пусть — вектор неизвестных параметров системы, и при каждом значении i = 1, …, n проводится ri измерений yij, m j = 1, …, ri некоторой скалярной функции h i, где h i — из вестный вектор. Обычно считают, что измерения yij при задан ном i проводятся в заданный момент времени ti (возможен и другой параметр привязки всех измерений yij при заданном i, например дальность полета при движении точки по траекто рии). Тогда можно измерения в этот момент считать сеансом измерения. Например, космические измерения проводятся се ансами, соответствующими небольшим интервалам времени, в течение которых есть радиовидимость космического объекта.

Эти интервалы часто можно при планировании космического эксперимента заменять отдельными моментами ti.

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

Согласно теореме Гаусса – Маркова, в классе линейных несме щенных алгоритмов оценивания минимальную дисперсию оценки любого параметра l = b дает МНК, соответствующий весовой матрице W = I. Учитывая это, осредним модель оценивания.

При использовании МНК исходная модель Лемма 8.1.

измерений yij = h i + ij, i =1,, n, j =1,, ri эквивалентна усредненной модели yi = h i + i, i =1,, n, (8.1) 8.1. усреднение модели измерений где r r 1i 1i yij, i = r ij yi = ri j =1 i j = есть средние арифметические указанных выше ri измеренных зна чений yij и ошибок ij. Эквивалентность данных моделей означает, что оценки параметра, получаемые при использовании этих мо делей, совпадают. При этом усредненные ошибки i являются не коррелированными между собой случайными величинами с нулевы ми математическими ожиданиями и дисперсиями Di =.

ri Для доказательства леммы, согласно теореме Гаусса – Мар кова, запишем дисперсию оценки МНК некоторого скалярного параметра l = b, которая соответствует весовой матрице W = I, как минимальную дисперсию в классе всех линейных несме щенных алгоритмов. Получаем ri n lмнк = xij yij, i =1 j = где коэффициенты xij оценивателя находятся из условия мини мума дисперсии оценки при условиях несмещённости:

n ri ri n xij hi = b.

DlМНК = min xij : (8.2) xi i =1 j =1 i =1 j = Обозначим ri xi = xij. (8.3) j = ri xij Минимум по xij суммы в (8.2) при условии (8.3) и фикси j =1 x рованных xi достигается при xij = i для всех j, т. е. все повторя ri ющиеся измерения в оптимальном случае оцениваются с одним оценивателем (это ясно и из соображений симметрии). Тогда (8.2) переписывается в виде n x2 n DlМНК = min i : xi h i = b. (8.4) xi i =1 ri i = 86 8. классические Задачи Планирования ЭксПеримента Согласно теореме Гаусса – Маркова равенство (8.4) доказывает лемму 8.1, так как линейный оцениватель с коэффициентами xi соответствует модели измерений (8.1).

Замечание 8.1. Так как в (8.4) суммирование записано по всем i, то, чтобы исключить суммирование по тем элемен там, где ri = 0, целевая функция должна быть доопределена сле дующим образом:

xi2 0, если xi = 0 и ri = 0, = ri +, если xi 0 и ri = 0.

8.2. ПОстанОвка задаЧи Пусть задано общее число измерений:

n ri = N.

i = Вместо чисел ri будем использовать вектор p с координатами r pi = i, i =1,, n.

N Далее будем, как это обычно делают, пренебрегать целочислен ностью ri и считать, что p — непрерывный план эксперимента, т. е.

любой вектор из симплекса n n = p n : pi 0, pi =1.

i = m Пусть b j, j = 1, …, s — заданные векторы, а l j = b j — контролируемые параметры (обычно s m ). При заданном пла не p рассмотрим линейные оценки каждого из этих параметров:

l j = ij yi, p {i : pi 0}, (8.5) i p где числа ij, i = 1, …, n определяют линейный оцениватель па раметра lj. Этот оцениватель дает несмещенные оценки параме тров b j тогда и только тогда, когда он удовлетворяет условиям несмещённости 8.2. Постановка задачи hi ij = b j, j =1,, s, (8.6) i p т. е. bj есть линейная комбинация тех векторов hi, для которых pi 0. Далее будем обозначать символом n множество тех p из n, для которых уравнения (8.6) разрешимы относительно ij.

Для каждого p n рассмотрим наилучшие линейные не смещенные оценки, для которых оцениватель определяется из условия минимальных дисперсий при линейном несмещенном оценивании:

ij = min : h i ij = b j, j =1,, s.

(8.7) NDl j ij i p pi i p Задача оптимального планирования эксперимента состоит в нахождении плана p n, доставляющего нижнюю грань за данному критерию оптимальности L(p) [Бахшиян, Соловьев, 1998;

Бахшиян, 2012]:

L = inf {L(p): p n }. (8.8) Будем рассматривать два критерия оптимальности, предпола гая далее, что оценки (8.5) являются наилучшими нелинейными несмещенными. Первый критерий 1 s L(p) = (8.9) N Dl j j = называется критерием L-оптимальности. Второй критерий { } L(p) = N max Dl1,, Dl s (8.10) будем называть MVs-критерием (при s = m и lj = j, j = 1, …, m, критерий (8.10) называют просто MV-критерием). В зависимо сти от критерия (8.9) или (8.10) задачу (8.8) будем соответствен но называть L-задачей или MVs-задачей.

При решении задачи (8.8) будем предполагать, что rang(h1,, h n ) = m (8.11) (т. е. ранг составной матрицы максимален), что равносильно до пущению о возможности линейного несмещенного оценивания всех компонент вектора в случае использования всех n групп измерений.


88 8. классические Задачи Планирования ЭксПеримента Замечание 8.2. Нетрудно показать, что суммирование по i p в соотношениях (8.5)–(8.7) можно заменить суммирова нием по i = 1, …, n, если доопределить слагаемые в (8.7) при pi = 0 следующим образом:

0, если = 0 и p = 0, ij = ij i pi +, если ij 0 и pi = 0.

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

Замечание 8.3. Согласно теореме Гаусса – Маркова, при вы полнении условия (8.6) наилучшая нелинейная несмещенная оценка совпадает с оценкой взвешенного метода наименьших квадратов l j = b j, где — любое решение нормальных урав нений n M(p) = pi h i yi, (8.12) i = а n M(p) = pi h i h i i = — нормированная на число N информационная матрица. Эти уравнения всегда совместны (доказательство этого факта предо ставляется читателю), но решение будет единственно только в случае невырожденной матрицы M(p). План p в этом случае называется невырожденным. В противном случае план p являет ся вырожденным, и решение уравнений (8.12) неоднозначно.

Однако условие (8.6) эквивалентно тому, что оценка b j каж дой функции b j, j = 1, …, s, одна и та же для всех решений уравнения (8.12). При этом параметры b j называются оценива емыми. Таким образом, условия несмещённости и оцениваемо сти эквивалентны. Далее, условия несмещённости (8.6) можно трактовать как принадлежность векторов bj линейному много образию, порожденному векторами p1 h1,, pn h n или, эк вивалентно, столбцами информационной матрицы. Это означа m ет, что найдутся векторы j такие, что выполняется условие оцениваемости параметров b j :

M(p) j = b j, j =1,..., s. (8.13) 8.3. скалярная задача планирования эксперимента и алгоритм ее решения 8.3. скалярная задаЧа ПланирОвания ЭксПеримента и алгОритм её решения При s = 1 L-задача (а также MVs-задача) сводится к минимиза ции дисперсии Dl, где l = l1 и называется C-задачей планирова ния эксперимента. Согласно (8.4), C-задача записывается в виде n x2 n ND0 = min i : xi h i = b, x n, p n.

* (8.14) pi, xi pi i = i =1 Проведем далее минимизацию дисперсии по pi, пренебре гая тем, что ri должно быть целым. Это оправдано при достаточ но большом N. Получим, используя метод множителей Лагран жа, что оптимальный план измерений есть xi pi =.

n xk k = При этом для заданного оценивателя с коэффициентами xi n 1 ND0 =. (8.15) xi N i = Замечание 8.4. Указанная оптимизация по ri здесь проведе на не совсем строго (так же как и в основополагающей работе [Elfving, 1952]). Действительно, мы проводили оптимизацию по pi только при условии p1 + … + pn = 1. Однако, условие несме щённости может выполняться не при всех таких pi, т. е. оценка параметра l может быть получена не всегда. В связи с этим ус ловие несмещённости также называют условием оцениваемости.

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

Оптимизация величины D0, полученной при оптималь ных ri, по оценивателю x приводит к задаче линейного програм мирования n n N * D1 = min xi : xi h i = b.

* xi i =1 i = 90 8. классические Задачи Планирования ЭксПеримента Введем новые переменные zi = xi, zi 0. Тогда xi hi = zi ai, где h, x 0, ai = i i -h i, xi 0.

При этом оптимальную задачу можно записать в виде n n 1 = min zi : zi a i = b, zi0, a i i, * (8.16) zi, a i i =1 i = где i = {h i, - h i } — множества, состоящие из двух векторов (рис. 14).

Сформулированную задачу можно интерпретировать как задачу линейного программирования, состоящую в минимиза ции суммы коэффициентов zi 0 разложения вектора b по век торам ±hi. Геометрический способ ее решения (в двумерном случае) обсуждался в пп. 1.3.4. У этой задачи существует реше ние, содержащее не более m отличных от нуля коэффициентов, т. е. вектор b лежит внутри многогранного угла, образованно го некоторыми m векторами ± hi (базисом), причем гиперпло скость, проведенная через их концы, отсекает от вектора b наименьшую (считая от конца вектора) часть по сравнению с другими комбинациями из m векторов. Легко заметить, что указанная гиперплоскость, проходящая через концы векторов оптимального базиса, принадлежит выпуклой оболочке всех 2n векторов ± hi, i = 1, …, n (см. рис. 2).

Рис. 14. Построение выпуклой оболочки conv( 1, 2,, n) множеств i при решении задачи (8.16) 8.4. оценивание параметров параболической траектории по измерениям дальности Таким образом, можно сделать следующий важный вывод.

При решении скалярной задачи планирования эксперимента следует выбрать оптимальные m измерений и найти оценива тель для усредненной модели из соотношения H 0 x 0 = b x 0 = (H 0 )-1 b, где H0 — матрица размерности mm, составленная из векто ров hi, входящих в оптимальный базис. Это соответствует тому, что вектор находится из m уравнений yi = h i, i =1,, m, где нумерация по i условна. После этого оптимальные значения величин ri находятся из условия их пропорциональности вели чинам xi.

8.4. Оценивание ПараметрОв ПарабОлиЧескОй траектОрии ПО измерениям дальнОсти Рассмотрим в качестве примера плоское движение материаль ной точки, брошенной под углом к горизонту с начальной скоростью 0. Введем следующие обозначения: r, — векторы положения и скорости точки в момент времени t;

r0, 0 — векто ры положения и скорости в начальный момент времени t0 = 0;

r(t ) x(t ) = — вектор состояния.

(t ) Движение точки описывается радиус-вектором r (t ) x r(t ) =, ry (t ) где rx (t ) = x t, ry (t ) = y t - gt 2, T y x = 0 cos, y = 0 sin, t [0, T ], = 0 — время подъ 2 g 0 ёма на максимальную высоту. Пусть на интервале [0, tk ] могут производиться измерения дальности l = r(t ) до движущейся 92 8. классические Задачи Планирования ЭксПеримента точки. Ошибки дальности будем считать случайными величина ми с нулевым математическим ожиданием и стандартным от клонением l = 1 м.

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

Требуется:

• решить задачу оптимального оценивания дальности на интервале [0, tk ], т. е. найти оптимальные значения моментов измерений t1, t2 и соответствующую среднеквадратичную ошиб ку оценки 1 = D1 ;

• доказать, что при любом tk справедливо t2 = tk;

• интерпретировать результат как C-задачу планирования эксперимента, найти соответствующие доли p1, p2 общего числа измерений в моменты времени t1 и t2.

Запишем линеаризованную модель измерений в виде l (t ) = h r r0 + h 0 + (t ) = h + (t ), где h (h r, h ) — составной вектор-строка.

Нахождение оптимальных моментов измерений, согласно изложенной выше теории, сводится в данном случае к решению двумерной задачи линейного программирования (используется только вектор h(t)) k k min xi : xi h (ti ) = b, (8.17) xi i =0 i = b = (T, 0). Будем искать решение этой задачи геометрически где в соответствии с подходом, изложенным в пп. 1.3.4. При этом для простоты и наглядности будем пренебрегать дискретностью моментов измерений, рассматривая вместо моментов ti непре рывную переменную t.

Компоненты вектора h(t) вычисляются из соотношений l (t ) l (t ) r(t ) l (t ) (t ) h (t ) = =.

+ r(t ) (t ) = 0 = 0 = 8.4. оценивание параметров параболической траектории по измерениям дальности Производные измеряемой функции вычисляются l(t) по формулам l (t ) r l (t ) = =,.

r(t ) (t ) r Окончательно векторы h(t) для каждого t будут иметь вид ry (t ) t rx (t ) t h (t ) =,, (8.18) rx2 (t ) + ry2 (t ) 2 rx (t ) + ry (t ) где rx (t ) = x t, ry (t ) = y t - gt.

0 Годограф функции h(t) имеет вид, показанный на рис. 15.

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

Построение линейной оболочки и оптимальный базис показа ны на рис. 16. Нетрудно заметить, что при любом значении tk [0, T ] вектор h(tk) входит в состав оптимального базиса, т. е.

выполняется условие t2 = tk.

Рис. 15. Годограф функции h(t) (см. формулу (8.18)) 94 8. классические Задачи Планирования ЭксПеримента Рис. 16. Построение оптимального базиса в задаче (8.17) Момент времени t1 может быть найден из условия касания годографа и прямой, проходящей через его крайнюю точку, со ответствующую моменту времени tk. После этого из (1.12) опре деляются x1 и x2, а по ним — p1 и p2. Значение среднеквадрати ческой ошибки, согласно (8.16), равно оптимальному значению целевой функции задачи (8.17).

Возьмем, например, tk = 80 c. Для этого случая найдем:

t2 = tk = 80 с, t1 = 38 с, x1 = –1,49, x2 = 1,96, p1 = 0,43, p2 = 0,57, 1 = 3,46.

9. задача оПтиМальной линейной иМПульсной КорреКции и алгоритМ её решения Определение 9.1. Коррекцией движения будем называть зада чу исправления траектории движения космического аппарата по результатам измерений и составления на их основе прогноза движения.

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

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


9.1. задаЧа ОПтимальнОй кОррекции При известнОм кОрректируемОм вектОре Задача оптимальной коррекции при известном векторе коррек тируемых параметров имеет большой самостоятельный интерес s и разнообразные практические приложения. Пусть l — не который вектор параметров системы, который может быть из менен путем импульсной коррекции ее траектории. Такая кор рекция заключается в мгновенном изменении компонент вектора скорости в моменты времени t1, …, tn. Изменение в каж дый момент ti будем характеризовать корректирующим векто ром (импульсом) ui размерности si, принадлежащим евклидово s му пространству i. Например, в случае коррекции траектории 96 9. Задача оПтимальной линейной имПульсной коррекции… космического аппарата, если импульс ui может производиться в произвольном направлении, то si = 3, i = 3 — трехмерное s евклидово пространство;

если же заранее задана плоскость или прямая, которым должен принадлежать импульс, то si равно s двум или единице, а i есть плоскость или прямая. Будем рас сматривать идеальную линейную коррекцию, т. е. предполагать, что импульсы производятся без ошибок и изменение вектора l в момент ti под действием импульса ui равно Ui ui, где Ui — ма трица размерности ssi, характеризующая влияние импульса и называемая матрицей влияния. Тогда для изменения вектора l s на векторную величину b нужно импульсы ui произвольно выбирать из условия n Ui ui = b.

i = Пусть затраты на коррекцию в момент ti характеризуются s величиной pi (ui), где pi(·) — некоторая норма в i. Например, в случае коррекции траектории космического аппарата путем ориентации в пространстве с использованием одного двигателя, жестко связанного с аппаратом, затраты на коррекцию малого импульса пропорциональны величине 1 3 pi (ui ) = ui ui ui = uij.

ui 2 j = Если же для коррекции используются двигатели, действующие вдоль трех осей координат, то = uij, pi (ui ) = ui j = где uij — составляющие вектора ui вдоль направления действия двигателей.

При сделанных допущениях поставим задачу оптимизации суммарных затрат на коррекцию. Эта задача имеет вид n n L = min pi (ui ): U i ui = b.

ui i =1 i = Решение данной задачи определяет и оптимальные моменты приложения импульса, так как, если ui = 0, то в момент ti им пульс не проводится. Представим каждый вектор ui в виде ui = xi i, xi pi (ui ) 0, 9.1. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения s где i i принадлежит множеству {i : pi (i ) = 1}. Тем самым, xi есть длина вектора ui в метрике, определяемой нормой pi (·), а i есть единичный вектор в этой метрике, направленный вдоль импульса.

Используя введенные новые переменные, запишем опти мальную задачу коррекции в виде n n L = min xi : xi a i = b, xi 0, a i i, i = 1,, n, xi, a i i =1 i = { }. Эта задача по виду si где i = a i : a i = U i i, pi ( i ) =1, i напоминает задачу линейного программирования, однако здесь переменными оптимизации являются не только числа xi, но и векторы ai в линейных условиях-равенствах, выбираемые из множеств i. Такая задача называется обобщенной задачей ли нейного программирования [Данциг, 1966, Бахшиян и др., 1980].

Ее можно также назвать многопараметрической задачей линейно го программирования. Приведение задачи коррекции к обобщен ной задаче линейного программирования и алгоритм ее реше ния были предложены М. Л. Лидовым* [Лидов, 1971].

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

Геометрический способ решения обобщенной задачи ли нейного программирования аналогичен рассмотренному выше геометрическому способу решения обычной задачи линейного программирования. Для случая, когда нормы pi являются евкли довыми, множества i представляют собой линейное преобра зование сфер i =1, т. е. являются эллипсоидами. Такой слу чай изображен на рис. 17. Если вектор b находится внутри многогранного угла, образованного некоторыми k векторами a i i, i = 1, …, k, то векторы ai, i = 1, …, k дают решение задачи в следующем смысле (см. рис. 17).

* Михаил Львович Лидов (1926–1993) — советский и российский учёный в области прикладной небесной механики, известен результата ми в теории эволюции орбит и в космической баллистике.

98 9. Задача оПтимальной линейной имПульсной коррекции… Рис. 17. Геометрический способ решения оптимальной задачи коррекции Нужно провести k s корректирующих импульсов, вели чина каждого импульса равна xi, направление определяется век тором i. При этом случай k s уже не является исключитель ным. Например, если вектор b совпадает по направлению с вектором ai, принадлежащим выпуклой оболочке i, то k = 1.

Симплексный алгоритм также сводится к вводу в базис векто ра as, который пересекает гиперплоскость, проходящую через концы векторов базиса. Проверка достаточного условия опти мальности, таким образом, сводится к проверке условия непе ресечения этой гиперплоскости любым из бесконечного множе ства векторов a i i, i = 1, …, n. Если это условие проверяется эффективно, то алгоритм симплекс-метода также эффективен.

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

Пусть B = (a1, …, as ) — матрица базиса на текущей итера ции, а вектор однозначно находится из системы уравнений B = (1,,1) j =1- a j = 0, j =1,..., s. (9.1) Таким образом, относительные разности j для векторов бази са равны нулю (см. шаг 2 алгоритма, приведённого в пп. 1.3.2).

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

Достаточное условие оптимальности базиса состоит в том, что все векторы из множеств i лежат ниже указанной гипер плоскости. Математически это записывается в виде равенства нулю минимума относительной разности на множестве всех векторов a i, i = 1, …, n, min min min(1- a) = 0. (9.2) i ai Из (9.2) получаем min =1- max max a =1- max i. (9.3) i i ai Здесь i = max a = max U i = U i, (9.4) = ai U i и максимум реализуется при =.

U i Если условие (9.2) не выполняется, т. е. существует номер k такой, что k = max i 1, то в базис вводится любой вектор i a k k, для которого не выполняется условие (9.2). Вывод век тора из базиса осуществляется по обычным правилам симплекс метода.

Замечание 9.1. Если в рассмотренной задаче коррекции все импульсы одномерные, т. е. si = 1 для всех i, то эта задача сводит ся к задаче линейного программирования, которая возникает при рассмотрении C-задачи планирования эксперимента.

Замечание 9.2. Если полученный указанным алгоритмом оптимальный базис содержит не более одного вектора с ненуле вым весом xi из каждого множества i, то этот базис дает реше ние задачи оптимальной коррекции. Данное условие может не выполняться только в случае, когда некоторое множество i содержит плоскость, проходящую через совокупность aij, j = 1, …, ri векторов из i, имеющих ненулевые веса xij (подроб нее об этом см. также в статье [Бахшиян, 1989], где исследуется 100 9. Задача оПтимальной линейной имПульсной коррекции… обобщенная задача линейного программирования). Но в этом случае один из векторов указанной совокупности векторов за меняется некоторой их выпуклой комбинацией — вектором ri a i = j a ij, j = где xj j =.

ri xil l = После этого множество i будет содержать лишь один вектор ai с ненулевым весом, т. е. новый базис дает решение задачи опти мальной коррекции.

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

Замечание 9.3. В отличие от обычной задачи линейного программирования, в обобщенной задаче линейного програм мирования изложенный алгоритм требует, вообще говоря, бес конечного числа итераций (см. рис. 17). Однако за конечное число итераций удается найти приближенное решение с задан ной точностью с использованием на каждой итерации оценки для оптимума n xi n L xi, i = i = где max {1,, n }. Данная оценка является следствием со отношения (1.22).

9.3. коррекция параболической траектории летательного аппарата 9.2. ПрОектная задаЧа идеальнОй кОррекции Пусть вектор требуемой величины коррекции не задан, а из вестно лишь, что он лежит в заданной области:

b.

Такой случай имеет место, например, когда при расчетах траек тории движения летательного аппарата заранее, еще до запуска, необходимо рассчитать возможный расход топлива. Естествен но поставить задачу нахождения максимального значения сум марного импульса на оптимальную коррекцию при условии, что корректируемый вектор может принимать любые значения из. Таким образом отыскивается n n L = max min pi (ui ): U i ui = b. (9.5) i =1 b ui i = Можно показать, что если множество представляет со бой выпуклый многогранник, то максимум в (9.5) достигается в одной из вершин этого многогранника.

9.3. кОррекция ПарабОлиЧескОй траектОрии летательнОгО аППарата Рассмотрим задачу коррекции траектории летательного аппа рата (ЛА), движущегося в гравитационном поле по параболиче ской орбите в некоторой плоскости с известной начальной ско ростью. Такое движение может рассматриваться как частный случай движения материальной точки, брошенной под углом к горизонту (см. п. 8.4). Рассмотрим для данного случая зада чу идеальной коррекции. Предположим, что коррекция может производиться либо в произвольном, либо в заданном направ лении с помощью одного двигателя, жестко связанного с ЛА.

Пусть известна начальная скорость ЛА 0 = x, y, 0 и проводятся измерения двух параметров движения — дально сти и высоты полета (т. е. вектор контролируемых параметров l двумерный). Допустим, что в момент времени T ЛА должен прийти в требуемую точку, которая является целью полета, од нако из-за ошибок исполнения начального импульса реальная траектория ЛА отличается от требуемой и заканчивается в какой-то другой точке. Рассмотрим вектор b, компоненты 102 9. Задача оПтимальной линейной имПульсной коррекции… которого равны разности контролируемых параметров ЛА для требуемой и реальной траекторий в конечный момент времени (т. е. представляют собой требуемую величину коррекции по каждому из параметров).

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

Пусть ti, i = 1, …, n — моменты времени, в которые может быть осуществлена коррекция, ui — корректирующий импульс в момент ti, p (ui ) = ui — затраты на проведение импульса в момент ti. Тогда задача оптимизации суммарных затрат на коррекцию имеет вид n n L = min ui : U i ui = b, ui i =1 i = где Ui — матрица влияния импульса, произведенного в мо мент ti, на вектор контролируемых параметров l.

Представим ui в виде ui = xi i, где xi ui — длина вектора импульса, определяемая нормой ;

i — единичный вектор, определяющий направление импульса.

В результате приходим к следующей задаче:

n n L = min xi : xi U i i = b, xi 0, i =1,..., n.

i =1 xi, i i = { } Введем векторы ai = Ui i, a i i = U i i : i, i =1.

Тогда задача оптимизации принимает вид n n L = min xi : xi a i = b, xi 0, a i i, i =1,..., n.

xi, a i i =1 i = Уравнения движения ЛА в произвольный момент времени t записываются в виде r (t ) = t, x x r (t ) = t - gt, y y (t ) =, x x y (t ) = y - gt.

9.3. коррекция параболической траектории летательного аппарата Вектор b определяется как r (T ) - r (T ) b = x x, r (T ) - r (T ) y y где T — время полета ЛА, rx (T ), ry (T ) — требуемые значения координат в момент времени T.

Запишем уравнения движения ЛА в момент времени T че рез значения координат и скорости в текущий момент времени t:

r (T ) = r (t ) + (t )(T - t ), x x x r (T ) = r (t ) + (t )(T - t ) - g (T - t ).

y y y (9.6) (T ) = (t ), x x y (T ) = y (t ) - g(T - t ).

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

1. Импульсы направлены произвольно. В данной задаче зави симость r(T) от (t) нелинейна, поэтому нахождение матрицы влияния сопровождается линеаризацией модели:

r (T ) r (T ) x x (t ) (t ) r(T ) x T - t y U(t ) = = = = (T - t ) E 2, T - t (t ) ry (T ) ry (T ) x (t ) y (t ) где E2 — единичная матрица размерности 22, r(t ) = (rx (t ), ry (t )), (t ) = ( x (t ), y (t )).

Тогда в дискретном случае U i = U(ti ) = (T - ti )E 2, множества i имеют вид { } i = (T - ti ) i : i 2, i =1, т. е. при каждом i концы векторов из множества i лежат на окружности радиуса (T – ti). Если коррекция может проводиться непрерывно в произвольные моменты времени из некоторого 104 9. Задача оПтимальной линейной имПульсной коррекции… интервала [t1, tn], то система множеств i переходит в систему множеств (t ), t1 t tn.

На рис. 18 показан геометрический способ нахождения ре шения. Несложно заметить, что оптимальным будет приложе ние одного корректирующего импульса в направлении векто ра b в момент времени t = t1.

Найдем величину этого импульса:

x1 a1 = x1 (T - t1 )1 = b, откуда b x1 = u1 =.

T - t 2. Импульсы направлены вдоль вектора скорости. В этом слу чае, как и в предыдущем, корректирующие импульсы ui будут иметь размерность ki = 2. Представим вектор скорости в виде произведения (t ) = (t ) (t ), (t ) где (t ) — длина вектора скорости;

(t ) = — единичный (t ) вектор, определяющий направление вектора скорости.

Из (9.6) следует, что 1 r(T ) = r(t ) + (t )(T - t ) - g (T - t )2.

2 Рис. 18. Геометрический способ решения задачи коррекции 9.3. коррекция параболической траектории летательного аппарата Рис. 19. Годограф множеств i (слева, коррекция проводится в задан ные дискретные моменты времени) и множеств (t ) (справа, коррек ция проводится в произвольные моменты времени) Тогда матрица влияния определяется следующим образом:

2 r(T ) 1 U(t ) = = r(t ) + (t )(T - t ) - g (T - t ) = (t ) (t ) 2 1 = r(t ) + (t ) (t )(T - t ) - g (T - t ) = (T - t ) (t ).

(t ) В дискретном случае отсюда следует, что U i = U(ti ) = (T - ti ) i, а множества i имеют вид { } i = (T - ti ) i i : i 1, i = ±1, или A i = {(T - ti ) i, - (T - ti ) i }, т. е. каждое из множеств i состоит из двух векторов. На рис. показаны годографы этих множеств в случае, если коррекция может производиться в заданные дискретные моменты време ни, и множеств (t ) при t1 t tn, если коррекция может про водиться в любой момент времени.

Возможны следующие случаи расположения вектора b.

1. Вектор b лежит в области годографа множеств i (рис. 20). В этом случае оптимальным решением будет приложе ние двух корректирующих импульсов (коррекция производится 106 9. Задача оПтимальной линейной имПульсной коррекции… в заданные дискретные моменты времени (см. рис. 20, слева), оптимальные моменты приложения импульсов — t = t1 и t = t2, импульсы производятся в положительном направлении, т. е.

при = +1). Если же коррекция может проводиться в любой мо мент времени (см. рис. 20, справа), то оптимальным будет при ложение одного корректирующего импульса, направление ко * торого определяется вектором a(t ).

2. Вектор b лежит вне области годографа множеств i (рис. 21). В случае, изображенном на левой части рис. 21, для коррекции ЛА необходимо приложить два корректирующих им пульса, один — в направлении вектора a1, т. е. в направлении вектора скорости при t = t1, а второй — в направлении векто ра –a4, т. е. в момент времени t = t4 в направлении, противопо ложном направлению вектора скорости.

Рис. 20. Построение решения в задаче оптимальной коррекции для дискретного (слева) и непрерывного (справа) случаев: вектор b лежит в области годографа Рис. 21. Построение оптимального решения в задаче коррекции для дискретного (слева) и непрерывного (справа) случаев: вектор b лежит вне области годографа 9.3. коррекция параболической траектории летательного аппарата В случае, изображенном на правой части рис. 21, для кор рекции ЛА также необходимо приложить два корректирующих импульса. Моменты времени и направление импульсов опреде ляются аналогично.

10. гарантированные хараКтеристиКи точности 10.1. критика классиЧескОгО ПОдхОда к Оцениванию тОЧнОсти При довольно естественных условиях, налагаемых на матрицу плана измерений и ковариационную матрицу ошибок, диспер сия оценки Гаусса – Маркова стремится к нулю с увеличением числа измерений n [Эльясберг, 1976]. Однако на практике та кой картины не наблюдается. Простой пример приведен в книге [Эльясберг, 1983]. Пусть измеряются показания различных ча сов одной точности и находится среднее арифметическое всех таких показаний. Такая оценка есть оценка МНК, которая со впадает с оценкой Гаусса – Маркова при отсутствии корреляции между наблюдениями и одинаковой дисперсии. При увеличе нии числа показаний должно было бы получиться сколь угод но точное время, но на практике этого не произойдет, так как существует корреляция между ошибками в точности хода, об условленная преемственностью способов изготовления часов, и эта корреляция может меняться в зависимости от времени их изготовления. Поэтому классический подход к вычислению точности оценивания может давать слишком оптимистические результаты.

10.2. гарантирующий ПОдхОд к выЧислению тОЧнОсти Оценивания Этот подход может дать, наоборот, более пессимистические прогнозы точности, но для практики иногда важно «перестрахо ваться». Поясним сущность этого подхода на рассмотренном выше примере вычисления дисперсии оценки. Пусть известно 10.2. гарантирующий подход к вычислению точности оценивания лишь множество, которому принадлежит ковариационная матрица ошибок измерений. Элементы K этого множества под чиняются условию неотрицательной определенности, которое будем записывать в виде K 0. Тогда в наихудшем с точки зре ния величины дисперсии случае получим гарантированную дис персию Dmax = max D(l).

K Отметим, что гарантированная дисперсия вычисляется при заданном оценивателе x, т. е., например, при заданной весовой матрице, которая уже не может выбираться в соответствии с те оремой Гаусса – Маркова, так как не известна ковариационная матрица ошибок. Вопрос об оптимизации оценивателя рассмо трим ниже. Сейчас выпишем Dmax для простых случаев задания множества. Допустим, что множество определяется следу ющим образом:



Pages:     | 1 || 3 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.