«Давыдов Алексей Алексеевич. МАТЕМАТИЧЕСКИЕ МОДЕЛИ ДЛЯ АНАЛИЗА ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ МАЛЫХ КОСМИЧЕСКИХ АППАРАТОВ Специальность 01.02.01 – Теоретическая механика. ДИССЕРТАЦИЯ на соискание учёной степени кандидата ...»

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

Давыдов Алексей Алексеевич.

МАТЕМАТИЧЕСКИЕ МОДЕЛИ ДЛЯ АНАЛИЗА

ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ МАЛЫХ

КОСМИЧЕСКИХ АППАРАТОВ

Специальность 01.02.01 – Теоретическая механика.

ДИССЕРТАЦИЯ

на соискание учёной степени кандидата физико-математических наук

Научный руководитель:

доктор физико-математических наук, профессор В.В. Сазонов Москва – 2012 2 Содержание Введение

Глава 1. Исследование режима гашения угловой скорости космического аппарата в нештатной ситуации

1.1. Уравнения вращательного движения КА и их стационарные решения

1.2 Устойчивость стационарных решений

1.3. Области притяжения стационарных решений

1.4 Реконструкция фактических реализаций режима гашения угловой скорости

1.5 Исследование стационарных вращений КА с учётом инерционных характеристик маховиков

Глава 2. Определение параметров вращательного движения КА ДЗЗ по телеметрическим данным о токе солнечных батарей

2.1. Математическая модель вращательного движения КА

2.2. Метод определения вращательного движения КА

2.3. Реализация метода Гаусса-Ньютона

2.4. Поиск начального приближения

2.5. Результаты определения вращательного движения КА

Глава 3. Определение параметров вращательного движения малого спутника связи по данным измерений тока солнечных батарей.. 3.1. Математическая модель вращательного движения КА

3.2. Метод определения вращательного движения КА

3.3. Реализация метода Левенберга-Марквардта

3.4. Поиск начального приближения

3.5. Результаты определения вращательного движения КА

Глава 4. Разработка модели вращательного движения

4.1 Вывод уравнений движения

4.2 Расчёт собственных частот модели

5. Заключение

6. Приложение 1

6.1 Программа для мониторинга состояния КА

7. Приложение 2

7.1. Программная реализация модели КА

7.2 Методика использования

7.3 Реализация тестовой задачи

7.4 Пример использования программного интерфейса библиотеки........ 8. Список литературы

9. Таблицы и рисунки

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

Актуальность темы. Одним из направлений деятельности ФГУП ГКНПЦ им. М.В. Хруничева в последние годы является разработка и эксплуатация МКА. Результатом работы стало создание МКА дистанционного зондирования Земли «Монитор-Э», геостационарных спутников связи «КазСат-1», «Экспресс-МД1», «КазСат-2». В настоящее время в ГКНПЦ ведётся разработка ряда перспективных МКА. Основой для создания всех перечисленных МКА стала разработанная в ГКНПЦ космическая платформа «Яхта» [1]. Платформа спроектирована в негерметичном исполнении и является универсальной, то есть на ее основе может быть создан широкий спектр МКА различного назначения, функционирующих как на низких, так и на высокоэллиптических и геостационарных орбитах. Тематика производства космических аппаратов является новой для предприятия, поэтому в ходе разработки перспективных МКА, платформа существенно модернизируется: применяются новые конструктивные решения, корректируется аппаратный состав платформы, совершенствуется алгоритмическое обеспечение бортовой системы управления. Это обстоятельство определяет повышенное внимание ко всем аспектам проектирования и эксплуатации МКА, одним из которых является исследование динамических характеристик МКА и его вращательного движения. Задача анализа динамических характеристик МКА тесно связана с разработкой математической модели вращательного движения МКА. Другая, не менее важная задача, возникающая при проведении лётных испытаний МКА – исследование его фактического вращательного движения. Наряду с практическим подтверждением корректности разработанной математической модели МКА, особенный интерес данные исследования представляют в случае, когда штатные измерительные средства МКА по каким либо причинам недоступны. В этом случае исследование вращательного движения МКА является источником дополнительной информации, позволяющей, например, спрогнозировать дальнейшее движение МКА, дать оценку энергобаланса и температурного состояния МКА и т.д.

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

Содержание работы.

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

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

Измерения угловой скорости используются при управлении вращательным движением КА с помощью двигателей-маховиков. Возникает задача исследования функционирования штатных алгоритмов управления при отсутствии части необходимых измерений. В диссертации эта задача решена для алгоритма, обеспечивающего гашение угловой скорости КА. Исследование функционирования алгоритма сводится к исследованию устойчивости стационарных решений системы нелинейных дифференциальных уравнений. Подобный класс задач хорошо изучен и широко представлен в литературе [2-10]. В частности, вопросы исследования систем с вращающимися массами изложены в [9]. В настоящей работе, с использованием разработанной В.В. Румянцевым и А.С. Озиранером теории устойчивости по части переменных [11] и теоремы Барбашина-Красовского [4, 5] показано, что эффективная работа исследуемого алгоритма возможна не при всех начальных условиях движения. В общем случае реализуется один из двух возможных финальных режимов, описываемых устойчивыми стационарными решениями уравнений движения. В одном из них компонента угловой скорости КА относительно оси, для которой отсутствуют измерения, отлична от нуля. С помощью численных расчетов получены оценки областей притяжения этих стационарных решений. Предложен простой способ, позволяющий вывести начальные условия режима гашения угловой скорости из области притяжения нежелательного решения. Для проверки адекватности исследуемой модели была проведена реконструкция нескольких имевших место реализаций этого режима. Реконструкция выполнена посредством аппроксимации решениями уравнений вращательного движения КА телеметрических значений компонент угловой скорости и суммарного кинетического момента двигателей-маховиков. Аппроксимация выполнена методом наименьших квадратов [12-14] с помощью методики, разработанной на основе подхода, предложенного в [15, 16]. Практическое применение и работоспособность данного подхода показаны в большом числе работ, например [17-21] и др.

Во второй главе реконструируется неуправляемый полёт КА в нештатной ситуации. Этот полёт проходил при отсутствии штатной телеметрической информации о параметрах вращательного движения. Для определения вращательного движения КА была использована доступная косвенная информация – телеметрические значения электрического тока, снимаемого с солнечных батарей (СБ). Идея использовать данные о токе, снимаемом с СБ, как источник информации для определения вращательного движения КА не нова [22-25]. В приведённых публикациях, данные о токе СБ используются для непосредственного определения текущей ориентации КА, часто – с целью использования полученных данных в алгоритмах системы управления. При этом предъявляются определённые требования к пространственной конфигурации панелей СБ, а сами данные об упомянутом токе используются, как правило, в сочетании с какой-либо дополнительной информацией, например – измерений магнитометра. В настоящей работе представлена методика иного рода, предназначенная для использования при пост-обработке телеметрической информации на Земле. К достоинствам разработанной методики можно отнести отсутствие какихлибо требований к пространственной конфигурации СБ, отсутствие необходимости в дополнительных данных о движении КА. В качестве недостатков можно указать невозможность использовать данную методику в контуре управления КА вследствие ее вычислительной сложности и наличия режимов вращения, при которых определение этого вращения по данной методике становится невозможным.

Разработанная интегральная статистическая методика обработки телеметрических данных также является развитием подхода, предложенного в [15, 16]. Для ее реализации разработана математическая модель движения КА с учетом действия на последний гравитационного и восстанавливающего аэродинамического моментов. Данные, собранные на отрезке времени длиной несколько десятков минут обрабатываются совместно методом наименьших квадратов с помощью интегрирования уравнений вращательного движения КА. Методика позволила определить фактическое вращательное движение КА в нештатной ситуации, уточнить значения моментов инерции КА и углов, задающих положение СБ в связанной с КА системе координат.

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

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

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

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

Вопросы динамики подобных систем представлены в литературе [26-34].

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

Для построения процедуры расчёта матрицы кинетической энергии предложен специальный векторный аппарат. Применение данного аппарата для описания динамики роботов-манипуляторов приведено в работе [35]. Такой подход позволил обеспечить достаточную простоту программной реализации разработанной модели, и одновременно её «программную автономность», то есть возможность реализации на вычислительной машине без привлечения сторонних математических библиотек и пакетов программ. Указанное обстоятельство существенно, так как позволяет применять разработанную модель как для наземной отработки системы управления, так и для использования её в качестве «эталонной модели движения», реализуемой на бортовой вычислительной машине КА.

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

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

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

стационарные решения КА представляет собой гиростат. Он состоит из твердого главного тела, на котором установлены три двигателя-маховика [34]. Каждый маховик имеет относительно главного тела одну степень свободы – может вращаться вокруг собственной оси материальной симметрии.

Систему координат, образованную главными центральными осями инерции КА, обозначим x1 x2 x3. В этой системе тензор инерции КА задается матрицей I diag( I1, I 2, I 3 ), H (h1, h2, h3 ) – собственный кинетический момент маховиков (гиростатический момент КА), (1, 2, 3 ) – абсолютная угловая скорость главного тела. По физическому смыслу I i (i 1, 2, 3). Оси вращения маховиков параллельны осям xi, так что каждая компонента гиростатического момента hi создается собственным маховиком.

В режиме гашения угловых скоростей управление кинетическим моментом маховиков можно с приемлемой точностью описать уравнениями hi ki i (i 1, 2, 3). Здесь и ниже точка над символом означает дифференцирование по времени t, ki – положительные параметры. Для реализации такого управления измерялись компоненты угловой скорости i. Как уже говорилось, в нештатной ситуации измерения величины 3 отсутствовали, и маховик, создававший компоненту гиростатического момента h3, не управлялся. Значение этой компоненты оставалось неизменным:

h3 h30 const.

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

Уравнения (1.1) допускают первый интеграл выражающий постоянство модуля кинетического момента КА в его движении относительно центра масс. Ниже используется функция Её производная по времени в силу уравнений (1.1) имеет вид Вследствие последнего соотношения, произвольное решение уравнений (1.1), оставаясь на поверхности интеграла (1.2) с неизменным значением G, с течением времени стремится к решению, в котором 1 2 0.

Подставив последние соотношения в уравнения (1.1), получим Выписанным соотношениям удовлетворяют два стационарных решения уравнений (1.1). Точнее, два семейства таких решений. Одно из них Здесь h10, h20 и 30 – произвольные постоянные, которые с постоянh10 h20 h30 G 2, ( I330 h30 )2 G 2. В случае решения (1.6) примем G I330 h30.

1.2 Устойчивость стационарных решений Интерес представляет асимптотическая устойчивость решений (1.5) и (1.6), однако вследствие существования у системы (1.1) первого интеграла (1.2) в данном случае она невозможна. Речь здесь может идти только об условной асимптотической устойчивости или асимптотической устойчивости по части переменных. Начнем с решения (1.6) при G 0.

Будем говорить, что решение (1.6) условно асимптотически устойчиво, если любое решение системы (1.1), начальные условия которого лежат в достаточно малой окрестности точки (1.6) и на той же самой поверхности интеграла (1.2), стремится к (1.6) при t. Чтобы исследовать такую устойчивость, можно с помощью интеграла (1.2) при G I330 h30 исключить 3 из системы (1.1) и исследовать обычную асимптотическую устойчивость стационарного решения 1 2 0, h1 h2 0 получившейся системы. В данном случае нет необходимости выполнять это понижение порядка в явном виде. Достаточно исследовать поведение функции T на поверхности интеграла (1.2) в окрестности точки (1.6) и воспользоваться результатами Е.А. Барбашина и Н.Н. Красовского [4, 5].

Положим 3 30 и разрешим соотношение (1.2) относительно. Получим Здесь и ниже многоточие означает члены третьей и более высокой степени относительно 1, 2, h1 и h2. Подставим 3 30 с учетом выписанного выражения для в выражение для T. Будем иметь Возьмем функцию Ляпунова в виде Ее производная по времени в силу системы, полученной из (1.1) исключением 3, имеет вид (ср. (1.4)) V 2k11 2k22. Условие положительной определенности квадратичных слагаемых V выражается неравенством G30 0. Множество V 0 при 30 0 не содержит целых траекторий новой системы кроме ее тривиального стационарного решения. По теореме Барбашина – Красовского при G30 0, k1 0, k2 0 это тривиальное решение асимптотически устойчиво. Если же G30 0, k1 0, k2 0, то согласно теореме Красовского рассматриваемое тривиальное решение неустойчиво.

Исследуем теперь устойчивость положения равновесия (1.5). Воспользуемся обобщением теоремы Барбашина – Красовского для задачи устойчивости по части переменных [11]. Фазовые переменные системы (1.1) разобьем на две группы, называемые переменными y и z :

y (1, 2, 3 ), z (h1, h2 ). Функция (1.3) при некотором положительном коэффициенте I удовлетворяет условию T I || y ||2. Здесь и ниже || x || – евклидова норма вектора x. Все решения системы (1.1) ограничены. В самом деле, в силу последнего неравенства и неравенства T 0 норма || y || ограничена. Отсюда в силу интеграла (1.2) следует ограниченность || z ||.

Множество T 0 в случае, k1 0, k2 0, | h10 | | h20 | 0 не содержит целых траекторий уравнений движения кроме стационарного решения (1.5). Следовательно, это решение асимптотически y-устойчиво (теоремы 19.1 и 19. в [11]). Иными словами, в любом решении системы (1.1) с начальными условиями из достаточно малой окрестности точки (1.5) y(t ) 0 при Отыскание стационарных решений системы (1.1) и исследование их устойчивости допускают обычную для задач такого рода интерпретацию в терминах задачи на условный экстремум [8, 10]. Согласно соотношению (1.4) устойчивым стационарным решениям системы (1.1) отвечают точки минимума функции (1.3) на поверхности интеграла (1.2), т. е. точки условного минимума этой функции. Точки невырожденного условного экстремума этой функции должны отвечать стационарным решениям этой системы.

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

Есть две возможности их удовлетворить.

Первая возможность: 0. В этом случае из первых трех уравнений (1.7) получаем 1 2 3 0. Переменные h1 и h2 могут принимать любые значения. Это стационарное решение (1.5). Вторая возможность:

I11 h1 0, I 22 h2 0. Здесь из первых двух уравнений (1.7) получаем 3 (G h30 ) / I3 30. Далее имеем h1 h2 0. Пришли к стационарному решению (1.6). В обоих случаях надо полагать G 0. Если G 0, то оба найденных экстремума являются вырожденными.

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

1.3. Области притяжения стационарных решений Построение областей притяжения стационарных решений выполнялось посредством численного интегрирования уравнений (1.1) при различных начальных условиях. Некоторые результаты такого интегрирования приведены на рис. 1.1 – 1.4. Каждый рисунок естественным образом разбивается на три части: левую, среднюю и правую. Левые и средние части рисунков иллюстрируют проекции вычисленных траекторий системы (1.1) на плоскости (1, 2 ) и (1, 3 ) соответственно. В правых частях рисунков изображены проекции этих траекторий на плоскость (h1, h2 ). Начальные точки траекторий помечены маркером в виде кружка. На рисунках и далее в тексте угловые скорости выражены в град./с, компоненты гиростатического момента – в Нмс. Траектории построены при I1 643 кгм 2, I 2 720 кгм 2, I3 253 кгм 2, k1 k2 3 кгм 2 /с. Начальные условия траекторий указаны в подписях к рисункам. При t траектории на рис. 1.1, 1.2 стремятся к решению (1.5), траектория на рис. 1.3, 1.4 – к решению (1.6). Как видно из рисунков, траектории системы (1.1) весьма разнообразны.

Для исследования областей притяжения перепишем систему (1.1) в виде: y Y(y,z), z Z(y,z). Здесь использованы векторные обозначения предыдущего раздела. Ниже ограничимся рассмотрением решений системы (1.1), для которых в (1.2) G 0. Для таких решений из области притяжения решения (1.5) при t имеем lim || y(t ) || 0, lim || z(t ) || 0 ; для решений из области притяжения решения (1.6), при t имеем lim || y(t ) || 0, lim || z(t ) || 0. Положим lim || y(t ) || || z(t ) ||. Значения 0 отвечают решениям из области притяжения решения (1.5), а значения 0 – решениям из области притяжения решения (1.6). Решения, для которых 0, лежат вблизи границы этих областей. Что представляет собой эта граница, неизвестно. Однако, как показывают расчеты, она располагается в малой окрестности гладкой поверхности. Величина зависит от начальных условий решения и параметров системы (1.1). Она рассчитывалась следующим образом. Система (1.1) интегрировалась из начальной точки t 0 до тех пор, пока выполнялось хотя бы одно из неравенств || Y[y(t ), z(t )]|| 108 с 2, || Z[y(t ),z(t )]|| 10 8 Нм и t 3600 с. В точке окончания интегрирования t t принималось || y(t ) || || z(t ) ||.

Исследуя зависимость величины от начальных условий решения системы (1.1) и ее параметров, можно получить оценки областей притяжения стационарных решений (1.5), (1.6). В общем случае такие области следует рассматривать в 11-мерном пространстве, координатами в котором служат 5 начальных условий и 6 параметров. Однако, учитывая механический смысл задачи, размерность исследуемого пространства целесообразно уменьшить. Во-первых, поскольку рассматривается конкретный КА, параметры I1, I 2, I 3, k1, k2 следует фиксировать. В расчетах использовались их значения, указанные выше. Во-вторых, процесс гашения угловых скоростей в нештатной ситуации начинался после «выбега» маховиков, то есть после остановки последних под действием трения. По этой причине следует принять h1 (0) h2 (0) 0. В результате указанного сокращения числа параметров области притяжения будем строить в 4-мерном пространстве величин 1 (0), 2 (0), 3 (0) и h30.

«Граница» областей притяжения строилась так. Фиксировались какие-либо три параметра, оставшийся параметр варьировался. Находились два значения этого параметра, которым отвечали значения разных знаков. Методом деления отрезка пополам находился корень уравнения 0.

Соответствующая точка в пространстве R4[1 (0), 2 (0), 3 (0), h30 ] считалась принадлежащей искомой «границе».

Результаты построения «границы» представлены на рис. 1.5 – 1.8. На рис. 1.5 представлена часть «границы», лежащая в гиперплоскости 2 (0) 0. По существу это поверхность в пространстве R3[1 (0), 3 (0), h30 ].

Она симметрична относительно плоскости 1 (0) 0. Сечение этой поверхности плоскостью 1 (0) 0, т.е. пересечение исходной «границы» в четырехмерном пространстве с подпространством 1 (0) 2 (0) 0, представляет собой прямые 3 (0) 0 и 3 (0) h30 / I 3. Это следует из найденного выше условия G30 0 условной асимптотической устойчивости решения (1.6). В данном случае кавычки в слове «граница» можно опустить: полученный результат – точный.

На рис. 1.6 представлена часть «границы», лежащая в гиперплоскости 1 (0) 0. Это – поверхность в пространстве R3[2 (0), 3 (0), h30 ]. Рассматриваемая поверхность аналогична поверхности на рис. 1.5. В частности, ее сечение плоскостью 2 (0) 0 представляет собой те же прямые 3 (0) 0 и 3 (0) h30 / I 3. Рисунок 1.7 демонстрируют зависимость проекции «границы» на пространство R3[2 (0), 3 (0), h30 ] при 1 (0) 0, а рисунок 1.8 – ту же зависимость при h1 (0) 0 и h2 (0) 0.

Анализируя представленные поверхности, можно отметить достаточно большую область притяжения стационарного решения (1.6), в которой сохраняется ненулевая компонента угловой скорости 3. Это обстоятельство говорит о наличии достаточно большой вероятности того, что полное гашение угловых скоростей в сложившейся ситуации окажется невозможным. Однако (см. выше), если процесс гашения начать после выбега маховиков, т.е. при h1 (0) h2 (0) h30 0, то движение системы будет стремиться к стационарному решению (1.5).

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

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

одноименные величины, деленные на I1. При построении аппроксимации телеметрических данных параметры системы (1.8),, k1, k2 и h30 выступают в роли параметров согласования.

Телеметрическая информация об угловой скорости КА и гиростатическом моменте КА (суммарном кинетическом моменте маховиков) в обработанном виде представляет собой последовательности чисел Здесь i( n ) и h(j n ) – измерения величин i, h j, …, выполненные в момент времени tn, t1 t2... t N. Разности tn1 tn обычно не превышают нескольких секунд. Продолжительность t N t1 обрабатываемых отрезков данных составляет от 3 до 6 мин. На самом деле телеметрические значения компонент угловой скорости и гиростатического момента относятся к так называемой строительной системе координат КА, которая связана с элементами конструкции его корпуса. Однако оси этой системы были практически параллельны осям системы x1 x2 x3, поэтому ниже данные телеметрии будем относить к системе x1 x2 x3. В принятой модели h3 h30 const, поэтому последовательность данных h3( n ) будем аппроксимировать константой h30.

Аппроксимацию телеметрических данных (1.9) будем строить методом максимального правдоподобия. Примем, что ошибки в данных i( n ), h(j n ) независимы и имеют нормальные распределения с нулевыми средними значениями. Стандартные отклонения ошибок в данных i( n ) одинаковы и равны, стандартные отклонения ошибок в данных h(j n ) также одинаковы и равны h. Значения этих стандартных отклонений неизвестны. Аппроксимацией данных (1.9) будем считать соответствующие компоненты фазового вектора в решении системы (1.8), заданном на отрезке t1 t t N и доставляющем минимум функционалу [12, 36] Последний член в (1.10) служит для учёта априорной информации о некоторых уточняемых параметрах и регуляризации задачи поиска минимума. В расчётах принимались: 0 2.54 и 0 0.73, k10 k20 0.12 (проектные значения этих параметров), w1 105, w2 104. Минимизация проводится по начальным условиям решения системы (1.8) в точке t1 и параметрам,, k1, k2, h30. Для простоты письма все уточняемые величины объединим в вектор q R10 и будем считать функционал (1.10) функцией этого вектора: (q). Аппроксимирующее решение отвечает значению q arg min (q).

Минимизация функционала (1.4) выполнялась в несколько этапов.

Сначала методом случайного поиска находилась грубая оценка q*, которая затем уточнялась методом Гаусса-Ньютона [37]. Точность аппроксимации данных (1.9) и разброс в определении компонент q будем характеризовать соответствующими стандартными отклонениями. Оценки стандартных отклонений ошибок в значениях i( n ) и h(j n ) находятся по формулам [12, 13] стандартные отклонения компонент вектора q равны квадратным корням из соответствующих диагональных элементов матрицы C 1 (q ), где C (q) и Ch (q) – матрицы систем нормальных уравнений, возникающих при минимизации выражений (q) и h (q) методом ГауссаНьютона. Стандартные отклонения величин i (t1 ), h j (t1 ), и т. п. будем обозначать i, hj,,… Аппроксимация данных (1.9) была построена на трех интервалах времени. На всех этих интервалах представлена заключительная часть процесса гашения угловой скорости. Основные характеристики интервалов приведены в табл. 1.1. Здесь для каждого интервала указаны: дата и декретное московское время момента t1, длина интервала t N t1, число N включённых в обработку моментов времени с измерениями. Результаты аппроксимации представлены на рис. 1.9. Здесь сплошные кривые – графики решения системы (1.8) (на нижних графиках – прямые h3 h30 ), марt, и t,h метрические данные. Как видно из рисунка, на всех трех интервалов движение КА стремится к стационарному решению (1.5).

Стандартные отклонения ошибок в телеметрических данных и найденных оценок параметров приведены в табл. 1.2. Здесь, h, i, hi и kj выражены в с 1. Как видно из таблицы, наиболее точная аппроксимация построена для интервала 1. Оценки параметров системы (1.8) для этого интервала k2 0.131с 1. Для интервалов 2 и 3 оценки параметров не приводим, так как их точность невысока. Эти интервалы отвечают окончанию процесса гашения угловой скорости, и относительные ошибки в данных 1n ), (2n ) на них выше, чем на интервале 1. Основываясь на результатах, полученных для интервала 1, можно сделать вывод об адекватности уравнений (1.8) и (1.1).

1.5 Исследование стационарных вращений КА с учётом инерционных характеристик маховиков Уравнения вращательного движения КА с явным учетом моментов инерции маховиков можно записать в виде уравнений Эйлера-Лагранжа [27]. Кинетическую энергию движения КА относительно центра масс представим в виде Здесь и ниже без дополнительных пояснений используются введенные выше обозначения; новые обозначения: J – осевой момент инерции каждого маховика (маховики полагаем одинаковыми), i – угол поворота маховика, ось которого параллельна оси xi, относительно главного тела.

Уравнения Эйлера-Лагранжа имеют вид Здесь Qi – механические моменты, приложенные к маховикам относительно их осей. Следуя сделанным ранее предположениям, примем Q1 k11, Q2 k22, Q3 0. В силу последнего соотношения 3 – циклическая координата, и ее можно исключить из уравнений движения методом Рауса.

Соответствующий циклический первый интеграл имеет вид Введем функцию Рауса R Tr h303, исключив из нее переменную 3 с помощью циклического интеграла. Уравнения вращательного движения КА можно записать в виде выписанных выше уравнений ЭйлераЛагранжа с функцией R вместо Tr и без уравнения для 3. Явный вид этих уравнений Сделаем в последних уравнениях замену переменных h1 J 1, h2 J 2 и преобразуем их к виду Для рассматриваемого КА J min( I1, I 2, I3 ). При J 0 уравнения (1.11) переходят в уравнения (1.1).

Уравнения (1.11) допускают первый интеграл (интеграл кинетического момента) явный вид которого При J 0 этот интеграл переходит в интеграл (1.2). Введем функцию Ее производная по времени в силу уравнений (1.11) имеет вид В силу последнего соотношения произвольное решение уравнений (1.11), оставаясь на поверхности интеграла (1.12) с неизменным значением G, с течением времени стремится к решению, в котором 1 2 0. Подставив последние соотношения в уравнения (1.11), получим 3h1 3h2 0, 3 0, h1 h2 0. Выписанным соотношениям удовлетворяют два семейства стационарных решения уравнений (1.11). Они имеют вид (1.5) и (1.6).

Произвольные постоянные в этих решениях связаны с постоянной G в (1.12) соотношениями: h10 h20 h30 G 2, ( I3 J )30 h30 G. Исследование устойчивости решений (1.5), (1.6) системы (1.11) выполняется по изложенной выше схеме. Приведем результат, полагая k1 0, k2 0, Решение (1.6) условно асимптотически устойчиво при G30 0 и неустойчиво при G30 0. Решение (1.5) асимптотически устойчиво по переменным 1, 2, 3 при | h10 | | h20 | 0. Таким образом, переход к системе (1.11) не вносит ничего существенно нового в результаты, полученные для системы (1.1).

Заметим, что штатное гашение угловой скорости описывается приведенными выше уравнениями Эйлера-Лагранжа при Q3 k33, k3 0. В силу этих уравнений T k11 k22 k33, и они допускают устойчивое в целом по переменным 1, 2, 3 семейство стационарных решений в котором 1 2 3 0.

Глава 2. Определение параметров вращательного движения КА ДЗЗ по телеметрическим данным о токе солнечных батарей Космический аппарат (КА) «Монитор-Э» вследствие нештатной ситуации совершал неуправляемый полёт, при этом отсутствовала телеметрическая информация о параметрах его вращательного движения. Возникла задача определения вращательного движения КА по доступной косвенной информации – электрическому току, снимаемому с СБ. В данной работе описана интегральная статистическая методика, позволившая решить эту задачу. Значения тока, полученные на отрезке времени длиной несколько десятков минут, обрабатывались совместно методом наименьших квадратов с помощью интегрирования уравнений вращательного движения КА. В результате обработки оценивались начальные условия движения, уточнялись моменты инерции КА и углы, задающие положение СБ в связанной с КА системе координат. Приведены результаты обработки нескольких отрезков данных, позволившие реконструировать фактическое вращательное движение КА.

2.1. Математическая модель вращательного движения КА КА считаем твёрдым телом, геоцентрическое движение центра масс которого – кеплерово эллиптическое. Элементы этого движения находятся по данным радиоконтроля орбиты. Для записи уравнений движения КА относительно центра масс введём две правые декартовы системы координат – орбитальную OX1 X 2 X 3 и образованную главными центральными осями инерции КА Ox1x2 x3. Точка O – центр масс КА, оси OX 3 и OX направлены соответственно вдоль геоцентрического радиус-вектора точки O и по трансверсали к орбите в этой точке.

Положение системы Ox1x2 x3 относительно системы OX1 X 2 X 3 зададим углами, и, которые введём следующим образом. Система OX1 X 2 X 3 может быть переведена в систему Ox1x2 x3 тремя последовательными поворотами. Первый поворот выполняется вокруг оси OX 2 на угол. Второй поворот выполняется вокруг новой оси OX 3 на угол. Третий поворот производится на угол вокруг оси OX 1, полученной в результате первых двух поворотов и совпадающей с осью Ox1. Матрицу перехода от системы Ox1x2 x3 к системе OX1 X 2 X 3 обозначим A || aij ||i3, j 1, где aij – косинус угла между осями OX i и Ox j. Элементы этой матрицы выражаются через введённые углы с помощью простых тригонометрических формул.

Система уравнений вращательного движения КА образована динамическими уравнениями Эйлера для компонент угловой скорости и кинематическими соотношениями Пуассона для элементов первой и третьей строк матрицы A. В уравнениях Эйлера учитываются гравитационный и восстанавливающий аэродинамический моменты. Уравнения вращательного движения имеют вид Здесь i и i – компоненты векторов абсолютной угловой скорости КА и скорости его центра масс относительно поверхности Земли в системе Ox1x2 x3 ; 0 – среднее орбитальное движение КА; J i – моменты инерции КА относительно осей Oxi, – плотность атмосферы в точке O ; pi – параметры аэродинамического момента. Плотность атмосферы рассчитывается согласно модели [39, 40].

Переменные a1i и a3i связаны условиями ортогональности матрицы A, поэтому начальные условия для них выражаются через углы, и.

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

2.2. Метод определения вращательного движения КА Решения уравнений (2.1), аппроксимирующие фактическое вращательное движение КА, выбираются из условия наилучшего сглаживания с их помощью телеметрических данных о токе, снимаемом с СБ. При построении аппроксимации наряду с начальными условиями движения КА, уточняются и параметры уравнений (2.1), и pi (i 1,2,3).

Ток, вырабатываемый СБ, примерно пропорционален косинусу угла падения солнечных лучей на их светочувствительную поверхность. Последняя расположена в плоскости, неподвижной относительно системы Ox1x2 x3. Нормаль n к светочувствительной стороне СБ зададим углом S между вектором n и его проекцией на плоскость Ox1 x2 ( | S | 10 ) и углом S между осью Ox1 и проекцией n на плоскость Ox1 x2. Компоненты n в системе Ox1x2 x3 имеют вид: n1 cos S cos S, n2 sin S cos S, n3 sin S.

Орт S направления «Земля-Солнце» задается в орбитальной системе координат компонентами Si (t ) (i 1,2,3), которые рассчитываются по приближенным формулам [40] и элементам кеплеровой орбиты КА. Упомянутый косинус и снимаемый с СБ ток задаются формулами Здесь I 0 – ток, вырабатываемый СБ на орбите Земли при перпендикулярном падении солнечных лучей на их плоскость, I 0 45 А.

На самом деле расчёт тока более сложен и требует знания труднодоступной дополнительной информации, но и упрощённые формулы (2.2) позволяют получить приемлемые результаты. Эти формулы можно ещё более упростить, если учесть, что в те моменты времени, когда телеметрические значения тока превышали некоторый положительный предел I min, заведомо выполнялось условие 0. Для таких моментов расчётные значения тока можно находить по формуле I I 0. При обработке данных, полученных на КА Монитор-Э принималось I min 10 А.

Телеметрическая информация о токе СБ представляет собой последовательность чисел Здесь I n – приближённое значение тока в момент времени tn, t1 t2... t N. Разности tn1 tn, как правило, не превышают нескольких секунд. В обработку включаются отрезки данных, длина которых t N t составляет от 20 до 40 минут.

Обработка данных (2.3) выполняется методом наименьших квадратов. Предполагается, что ошибки в значениях I n независимы и имеют одинаковое нормальное распределение с нулевым средним значением и стандартным отклонением. Значение неизвестно. На решениях уравнений (2.1), заданных на отрезке t1 t t N, определим функционал Аппроксимацией фактического движения КА на этом отрезке будем считать решение, доставляющее такому функционалу минимум. Минимизация проводится по начальным условиям решения в точке t1 : 0 (t1 ), Для простоты письма объединим все уточняемые величины в один вектор и будем рассматривать (2.4) как функцию (q). Тогда q* arg min (q) – искомая оценка вектора q.

Минимизация функционала (2.4) выполнялась в несколько этапов.

Сначала методом случайного поиска [41] находилась грубая оценка q*, которая затем уточнялась методом Левенберга-Марквардта и – на заключительном этапе – методом Гаусса-Ньютона [37].

2.3. Реализация метода Гаусса-Ньютона Применение этого метода в задачах определения вращательного движения спутников по данным измерений бортовых датчиков описано в [15-21]. На каждой итерации этого метода поправка q, уточняющая имеющуюся оценку q q*, определяется системой Здесь C (q) – матрица нормальных уравнений, E13 – единичная матрица порядка 13, e jnm – символ Леви-Чивиты, m qk – псевдопроизводные [15], служащие для представления истинных производных / qk при k 11. Значения m (tn ) qk (k 1,...,11) определяются в процессе интегрирования уравнений совместно с уравнениями (3.1) и уравнениями в вариациях относительно m qk. Ненулевые начальные условия m (t1 ) qk и m (t1 ) qk имеют вид Ненулевые производные n j qk в уравнениях (2.5) определяются формулами:

Точность аппроксимации данных (2.3) и разброс в определении компонент q* arg min будем характеризовать, следуя методу наименьших квадратов, соответствующими стандартными отклонениями. Стандартное отклонение ошибок в данных (2.3) находится по формуле стандартные отклонения компонент вектора q* равны квадратным корням из соответствующих диагональных элементов матрицы 2C 1 (q* ).

Ниже стандартные отклонения величин 0, 10, p1 и т. п. будем обозначать Чтобы в результате минимизации функционала (2.4) значения параметров pi,, лежали в разумных с физической точки зрения пределах, в этот функционал вводились дополнительные слагаемые Здесь 1 и 2 – положительные числа, 0 2.54 и 0 0.73 – проектные значения параметров и. В расчётах принимались 1 2 10.

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

2.4. Поиск начального приближения Для надежной работы алгоритмов Левенберга-Марквардта и ГауссаНьютона необходимо иметь достаточно точное начальное приближение точки минимума. Поиск начального приближения осуществляется методом случайного поиска с обучением [41]. При этом сначала функционал (2.4) рассматривается на укороченном начальном отрезке данных (2.3) и считается функцией меньшего числа переменных, чем входит в вектор q. Вопервых, некоторые компоненты фиксируются: принимается p1 p2 p3 0, параметры,, S и S приравниваются проектным значениям. Во-вторых, матрица A(t1 ) параметризуется углами и, которые вводятся так. Угол 0 между направлением на Солнце и нормалью к рабочей поверхности СБ в момент t1 фиксируется и задается формулой 0 arccos[ I (t1 ) I 0 ]. Этот угол не меняется при произвольных поворотах системы Ox1x2 x3 вокруг ортов S и n. Угол поворота системы Ox1x2 x3 вокруг S обозначается, угол поворота вокруг n –. Приведем формулу, определяющую A(t1 ).

Орт S в орбитальной системе будем задавать углами S и S. В этой системе S (sin S cos S,sin S,cos S cos S ). Матрицу A(t1 ) определим формулой Матрица K определяется положением Солнца в орбитальной системе координат в момент t1, матрица M зависит от положения орта n в системе Ox1x2 x3, Q – матрица последовательных поворотов КА: сначала вокруг орта n на угол, затем вокруг орта S на угол.

Для удобства задания начального значения угловой скорости КА и контроля его последующего уточнения, величины i 0 ( i 1,2,3 ) заменяютarctan(20 10 ), arcsin( 30 0 ). В итоге функционал (2.4) на укороченном начальном отрезке данных (2.3) минимизируется по 5 параметрам:,, *,,.

После выполнения заданного числа попыток случайного поиска осуществляется переход к параметрам q и полному отрезку данных. Здесь опять используется случайный поиск, а после его завершения применяется метод Левенберга-Марквардта.

2.5. Результаты определения вращательного движения КА Определение движения КА по данным о токе СБ было выполнено на 12 интервалах времени. Основные характеристики этих интервалов приведены в табл. 2.1. Здесь для каждого интервала указаны: дата и декретное московское время первого измерения t1, длина интервала t N t1, число N включённых в обработку измерений. Полученные результаты представлены в табл. 2.22.4 и на рис. 2.22.6. В табл. 2.2 и 2.3 приведены результаты минимизации функционала (2.4) на интервалах из табл. 2.1. Здесь указаны значения некоторых компонент вектора q*, стандартные отклонения этих компонент и стандартное отклонение ошибок в данных (2.3). Параметры 0, 0, 0, pi были определены со значительной погрешностью и поэтому не приводятся. Погрешность определения начального углового положения КА на интервалах 1 6 характеризуется значениями,, 1.5, на интервалах 712 – значениями,, 0.5, причем доминируют значения, близкие верхним границам. Стандартные отклонения параметров pi также очень велики. Большая погрешность определения величин 0, 0, 0 вызвана малой информативностью данных (2.3), а параметров pi – еще и слабым влиянием аэродинамического момента на вращательное движение КА.

Рис. 2.2, 2.3 содержат примеры графиков зависимости от времени компонент угловой скорости КА i и расчетного тока СБ I. Графики построены на соответствующих отрезках t1 t t N ; началом отсчета времени служит точка t1. Маркерами рядом с графиками I (t ) указаны аппроксимируемые данные измерений (2.3). Хотя графики на рис. 2.2 построены для интервала 2, они типичны для интервалов 1 6 из табл. 2.1; графики на рис.

2.3 относятся к интервалу 8 и типичны для интервалов 7 12 (см. [42, 43], где аналогичные графики приведены для всех интервалов из табл. 2.1). Судя по графикам и значениям в табл. 2.2 построенные аппроксимации тока СБ достаточно точны.

Анализ табл. 2.2 и графиков функций i (t ) показывает, что движение КА на интервалах 16 было близко к стационарному вращению вокруг оси Ox2 – оси максимального главного центрального момента инерции. Модуль угловой скорости составлял примерно 10 град/с. На интервалах движение КА было близко движению Эйлера-Пуансо, полодии которого охватывали ось Ox2. Модуль угловой скорости составлял около 1 град/с.

Более медленные движения КА были и более сложными. Изменение тока СБ в случае сложного движения оказалось более информативным, поэтому оценки q* на интервалах 712 получились заметно точнее, чем на интервалах 16 (ср. графики функций I (t ) на рис. 2.2 и 2.3, а также стандартные отклонения в табл. 2.2, 2.3 для интервалов 16 и 712).

Для визуализации движения КА можно использовать графики зависимости от времени углов, и. Однако в данном случае такое описание не наглядно. Если КА имеет большую угловую скорость, то наибольший интерес представляет движение той его главной центральной оси инерции, проекция угловой скорости на которую сохраняет знак. В данном случае такой осью является ось Ox2. С учетом этого обстоятельства вращательное движение КА представим в виде суперпозиции медленной эволюции вектора L кинетического момента КА в его движении относительно центра масс и быстрого движения КА вокруг этого вектора [27]. Движение L будем рассматривать относительно инерциальной системы координат OZ1Z 2 Z3, которая совпадает с орбитальной системой, фиксированной на момент прохождения КА через восходящий узел орбиты. Положение вектора L в системе OZ1Z 2 Z3 зададим углом L между L и плоскостью орбиты, и углом L между проекцией L на плоскость орбиты осью OZ3. Угол L считаем положительным, если проекция L на ось OZ 2 положительна. Направление отсчёта L согласовано с направлением оси OZ 2. Аналогичные углы S и S задают в системе OZ1Z 2 Z3 положение вектора S. Угол между векторами L и S обозначим S.

Как показали результаты обработки полученных данных, на интервале 1 имели место неравенства L 0.012 град/с и L 0.012 град/с, на интервалах 26 – неравенства L 0.003 град/с и L 0.003 град/с, на интервалах 712 – неравенства L 0.12 град/с и L 0.12 град/с. Указанные величины существенно меньше угловой скорости движения КА на соответствующих интервалах.

Для представления движения КА вокруг вектора L, будем рассматривать проекции годографа орта оси Ox2 на координатные плоскости системы Ol1l2l3, связанной с этим вектором (рис. 2.4). Система Ol1l2l3 вводится следующим образом. Ось Ol2 направлена вдоль L ; ось Ol1 лежит в плоскости Ol2 Z1, перпендикулярна L и образует тупой угол с осью OZ1 ;

ось Ol3 дополняет систему до правой. Введем также угол нутации между осями Ox2 и Ol2.

Предложенное в [27] представление движения в виде суперпозиции двух движений учитывает тот факт, что движение КА относительно вектора L является достаточно простым. В частности, угловые скорости КА на интервалах 1 6 описываются приближенными формулами где и W – постоянные величины, Формулы (2.6) получены посредством решения линеаризованных уравнений Эйлера без учета внешних моментов. Для каждого интервала t1 t t N постоянные и W определяются соотношениями Среднеквадратичные отклонения характеризуют точность формул (2.6). При отсутствии действующих на КА внешних моментов должно выполняться соотношение W 0. Если к тому W для интервалов 1 6 приведены в табл. 2.4. Как видим, постоянство угловой скорости 2 выполнено с высокой точностью, постоянство величины 3 r 21 выполнено менее точно.

Для интервалов 712 формулы (2.6) менее точны, чем для интервалов 1 6, но и здесь величины,, W и W оказываются полезными характеристиками. Отношение W / W для интервалов 712 оказалось примерно таким же как, для интервалов 1 6. Это неудивительно, так как функция 3 r 22 является первым интегралом нелинейных уравнений Эйлера, записанных без учета внешних моментов. Отношение / | | изза увеличения роли нелинейных эффектов для интервалов 712 оказалось существенно больше, чем для интервалов 1 6.

В табл. 2.4 приведены также параметры, характеризующие положение векторов L и S в системе OZ1Z 2 Z3. Поскольку эти векторы менялись медленно, их положение на каждом интервале t1 t t N удобно описывать усреднёнными углами L, L, S, S и s. Усреднение углов выполнялось по той же схеме, что и усреднение величин 2, 3 r 21. В табл. 2.4 указаны значения этих усредненных углов, а также еще максимальные значения max угла на обработанных интервалах.

На рис. 2.4, 2.6 представлена реконструкция движения оси Ox2 на интервалах 2 и 8. Здесь приведены проекции годографа орта оси Ox2 на плоскости (l1, l2 ), (l1, l3 ) и графики зависимости от времени углов L, L.

Рис. 2.4 и 2.6 типичны соответственно для групп интервалов 1 6 и (для всех интервалов табл. 2.1 аналогичные рисунки приведены в [42, 43]).

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

Глава 3. Определение параметров вращательного движения малого спутника связи по данным измерений тока солнечных батарей Рассматривается малый КА – спутник связи «КазСат-1», находящийся на геостационарной орбите. Для управления вращательным движением КА используются двигатели-маховики. КА после нештатной ситуации перешёл в состояние неуправляемого вращения, при этом непосредственная телеметрическая информация о параметрах его вращательного движения отсутствовала. Вследствие этого возникла задача определения вращательного движения спутника по имеющейся косвенной информации – току, снимаемому с солнечных батарей. Телеметрические измерения тока солнечных батарей, полученные на интервале времени длиной несколько часов, обрабатывались совместно методом наименьших квадратов с помощью интегрирования уравнений вращательного движения спутника. Приведены результаты обработки 10 отрезков данных измерений, позволившие определить фактическое вращательное движение КА и оценить суммарный кинетический момент двигателей-маховиков.

3.1. Математическая модель вращательного движения КА КА считаем гиростатом, совершающим свободное вращательное движение. Для записи уравнений этого движения введём две правые декартовы системы координат – неподвижную в абсолютном пространстве X1 X 2 X 3 и образованную главными центральными осями инерции КА x1 x2 x3. Ось X 1 направлена на Солнце, положение которого в системе X1 X 2 X 3 полагаем неизменным на рассматриваемом отрезке времени, ось X 3 направлена в сторону Северного полюса мира, лежащего в плоскости X1 X 3.

Положение системы x1 x2 x3 относительно системы X1 X 2 X 3 зададим углами, и, которые введём с помощью следующего условия. Система X1 X 2 X 3 может быть переведена в систему x1 x2 x3 тремя последовательными поворотами: 1) на угол вокруг оси X 1, 2) на угол вокруг новой оси X 2, 3) на угол вокруг новой оси X 3, совпадающей с осью x3.

Матрицу перехода от системы x1 x2 x3 к системе X1 X 2 X 3 обозначим A || aij ||i3, j 1, где aij – косинус угла между осями X i и x j. Элементы этой матрицы выражаются через введённые углы с помощью формул Система уравнений вращательного движения образована уравнениями, выражающими теорему об изменении кинетического момента спутника в его движении относительно центра масс и кинематическими соотношениями Пуассона для элементов первой и третьей строк матрицы A. Эти уравнения имеют вид Здесь i – компоненты вектора абсолютной угловой скорости КА;

J i – моменты инерции КА относительно осей xi ; J1hi – компоненты собственного кинетического момента двигателей-маховиков. В уравнениях (3.1) и далее, если не оговорено другое, компоненты векторов относятся к системе x1 x2 x3.

Переменные a1i и a3i зависимы – связаны условиями ортогональности матрицы A, поэтому начальные условия для них выражаются через углы, и. В процессе интегрирования уравнений (3.1), элементы второй строки матрицы A вычисляются как векторное произведение её первой и третьей строк. Размерность величин i и hi – с 1, единица измерения времени – секунда. Параметры hi считаются постоянными на участке интегрирования. Значения параметров и : 2.765, 0.474.

3.2. Метод определения вращательного движения КА Решения уравнений (3.1), аппроксимирующие фактическое движение КА относительно центра масс, выбираются из условия наилучшего сглаживания с их помощью телеметрических данных о токе, снимаемом с солнечных батарей (СБ). При этом определяются начальные условия движения КА и параметры hi (i 1,2,3).

Ток, вырабатываемый СБ, примерно пропорционален косинусу угла падения солнечных лучей на их светочувствительную поверхность. Полагаем, что эта поверхность представляет собой плоскость, неподвижную в системе x1 x2 x3. Компоненты орта нормали n (n1, n2, n3 ) к светочувствиn1 0.9997, n2 0.0191, тельной поверхности имеют значения:

n 3 0.0165. Орт S направления «Земля-Солнце» в системе координат X1 X 2 X 3 имеет компоненты S (1, 0, 0). Расчётное значение снимаемого с СБ тока и косинус угла между ортами n и S задаются формулами Здесь I 0 – максимальный ток, вырабатываемый батареями на орбите Земли при перпендикулярном падении солнечных лучей на их плоскость, I 0 102 А; На практике точный расчёт тока солнечных батарей достаточно сложен, однако упрощённые формулы (3.2) позволяют получить приемлемые для рассматриваемой задачи результаты. Эти формулы можно ещё упростить, учитывая, что в обработку следует включать только те моменты времени, для которых телеметрические значения тока превышают некоторый положительный предел I min. Для таких моментов времени заведомо выполнялось условие 0, и расчётные значения тока можно находить по I I 0. При обработке полученных данных принималось формуле Поскольку в формулах (3.2) величина зависит только от параметров a1k, график (t ) инвариантен по отношению к углу. С учётом этого обстоятельства, уравнения для a3i можно исключить из системы (3.1). Однако эти уравнения оставлены для того, чтобы наряду с угловыми скоростями КА, найти угловое положение КА на исследуемых отрезках времени. Далее всюду начальное значение угла принято равным нулю.

Телеметрическая информация о токе, снимаемом с батарей, представляет собой последовательность чисел Здесь I n – приближённое значение тока в момент времени tn, t1 t2... t N. Разности tn1 tn, как правило, не превышают нескольких секунд. В обработку включаются отрезки данных, длина которых t N t составляет от 3 до 11 часов.

Сглаживание данных (3.3) выполняется методом наименьших квадратов. Пусть ошибки в значениях I n независимы и имеют одинаковое нормальное распределение с нулевым средним значением и стандартным отклонением. Значение неизвестно. На решениях уравнений движения, заданных на отрезке t1 t t N, определим функционал Аппроксимацией фактического движения КА на этом отрезке будем считать решение, доставляющее такому функционалу минимум. Минимизация проводится по начальным условиям движения КА в точке t1 :

0 (t1 ), 0 (t1 ), i 0 i (t1 ), и параметрам hi. Как было указано выше (t1 ) 0. Для простоты письма все уточняемые параметры объединим в q* arg min (q) – искомая оценка вектора q.

Система (3.1) и выражения (3.2) инвариантны относительно преобразований векторов q и n, приведенных в табл. 3.1. В первой строке таблицы указаны значения компонент этих векторов, которые назовём номинальными. Эти значения задают некоторую функцию (t ). Во второй и последующих строках приведены преобразованные компоненты векторов q и n, задающие ту же самую функцию (t ). Поскольку компоненты ni (i 1,2,3) неизменны и удовлетворяют условиям | n2 || n1 |, | n3 || n1 |, n1 0, при минимизации функционала (3.4) для каждого исследуемого отрезка времени следует рассматривать четыре значения q, приближённо соответствующие первым четырём строкам таблицы 3.1. Для выбора значения q, отвечающего фактическому движению КА, проводился анализ полученных оценок компонент суммарного кинетического момента h (h1, h2, h3 ) двигателей-маховиков. Так как на всём рассматриваемом интервале движения КА один из двигателей-маховиков был отключён, то область вариации указанного вектора была достаточно сильно сжата в направлении оси вращения неработающего маховика. В связи с этим, практически во всех рассмотренных случаях из четырёх найденных значений h только одно удовлетворяло условию нахождения внутри области вариации суммарного кинетического момента. Это значение и принималось в качестве искомого решения.

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

3.3. Реализация метода Левенберга-Марквардта Применение этого метода в задачах определения вращательного движения спутников по данным измерений бортовых датчиков описано в [15-21]. Метод является одним из вариантов метода Гаусса-Ньютона. На каждой итерации этого метода поправка q, уточняющая имеющееся приближённое значение q, определяется системой Здесь C (q) – матрица системы нормальных уравнений, – положительный параметр, E8 – единичная матрица порядка 8, e jnm – символ ЛевиЧивиты, k qi – псевдопроизводные [15], служащие для представления истинных производных qk.

Значения m (tn ) qk (k 1,...,8) определяются в процессе интегрирования уравнений совместно с уравнениями (3.1) и уравнениями в вариациях относительно m qk. Ненулевые начальные условия m (t1 ) qk и m (t1 ) qk имеют вид Точность аппроксимации данных (3.3) и разброс в определении компонент q* arg min будем характеризовать, следуя методу наименьших квадратов, соответствующими стандартными отклонениями. Стандартное отклонение ошибок в значениях I n находится по формуле стандартные отклонения компонент вектора q* равны квадратным корням из соответствующих диагональных элементов матрицы * C 1 (q* ).

Ниже стандартные отклонения величин,, i, hi будем обозначать, 3.4. Поиск начального приближения Чтобы алгоритм минимизации Левенберга-Марквардта был надёжен, необходимо иметь достаточно точное начальное приближение точки минимума. Последнее находилось методом случайного поиска с обучением [41]. При этом функционал (3.4) рассматривался на укороченном начальном отрезке данных (4.3) и считался функцией четырёх переменных. Для сокращения числа переменных использовались следующие соображения [44, 45]. Начальный угол 0 (рис. 3.1) между направлением на Солнце и нормалью к рабочей поверхности СБ в момент t1 фиксировался и вычислялся по формуле 0 arccos[ I (t0 ) I 0 ]. Этот угол инвариантен по отношению к произвольным поворотам системы x1 x2 x3 вокруг ортов S и n. Как указывалось выше, поворот на произвольный угол вокруг S не влияет на характер зависимости (t ). Напомним, что этот поворот соответствует углу, (t1 ) 0. Искомое начальное угловое положение КА, описываемое матрицей A(t1 ), можно параметризовать единственным углом поворота вокруг n. Параметризацию зададим в виде Здесь Q – матрица поворота КА вокруг орта n на угол. Элементы матрицы Q зависят также от начального угла 0 между S и n. Соотношение (3.5) позволяет выразить углы 0 и 0 в функции угла :

Для удобства выбора начального значения и последующего контроля величин угловой скорости КА и суммарного кинетического момента двигателей-маховиков, величины i 0 и hi ( i 1, 2, 3 ) заменялись величинами Было известно, что КА до возникновения нештатной ситуации находился в орбитальной ориентации, а кинетический момент двигателеймаховиков был близок к нулю. Вращение КА после аварии было обусловлено в основном раскруткой двигателей-маховиков, поэтому угол между векторами кинетического момента корпуса КА и вектором h был близок к 180°. С учётом этого обстоятельства, при подборе начальных значений принималось: h k 0, h, h, k 0.02. Таким образом, в процедуре случайного поиска для функционала (3.4) принималось упрощённое выражение ( z ), z (, 0,, ) R4. После окончания случайного поиска найденное значение z пересчитывалось в значение q, которое служило начальным приближением для вычислений по методу Левенберга-Марквардта.

3.5. Результаты определения вращательного движения КА Определение фактического движения КА относительно центра масс по данным (3.3) было выполнено на 10 интервалах времени. Основные характеристики этих интервалов приведены в табл. 3.2. Здесь для каждого интервала указаны: дата и декретное московское время первого измерения t1, длина интервала t N t1, число N включённых в обработку измерений.

Полученные результаты представлены в табл. 3.33.5 и на рис. 3.43.13. В табл. 3.33.5 приведены результаты минимизации функционала (3.4) на интервалах из таблицы 3.2. Здесь указаны значения компонент вектора q*, стандартные отклонения этих компонент и стандартное отклонение ошибок данных (3.3). Стандартные отклонения не превышали 3.6, стандартные отклонения не превышали 2.4. Максимальные значения i составили для 1 – 0.023 град/с, для 2 – 0.032 град/с, для 3 – 0. град/с. Максимальные значения hi составили для h1 – 0.00041 с-1, для h 2 – 0.00046 с-1, для h 3 – 0.00008 с-1.

Интервал 1 данных получен приблизительно через сутки после развития нештатной ситуации, интервалы 28 охватывают промежуток времени с 17 августа по 9 сентября – до момента выключения двигателеймаховиков. Интервалы 9, 10 получены во время полёта КА с остановленными двигателями-маховиками. Анализ найденных угловых скоростей КА и суммарного кинетического момента двигателей-маховиков показал следующее. На интервалах движения 18 суммарный кинетический момент системы корпус КА – двигатели-маховики не превышал 2Hмс, а кинетический момент корпуса КА находился в диапазоне 20 22Нмс. На интервалах 9 и 10 суммарный кинетический момент составлял ~ 1Hмс, а кинетический момент корпуса КА не превышал величины 2.2Hмс. При этом на всех интервалах движения КА вектор суммарного кинетического момента двигателей-маховиков занимал произвольное положение в системе x1 x2 x3.

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

Для наглядного представления углового движения КА можно было бы использовать графики углов,,. Однако для найденного движения КА такое представление не обладает наглядностью. В связи с этим, поступим следующим образом. Введём обозначения (рис. 3.2): L KA – вектор кинетического момента корпуса КА, L КУДМ – вектор кинетического момента двигателей-маховиков, L – вектор суммарного кинетического момента КА с маховиками. Для указанных векторов выполняется равенство L L КА L КУДМ. Вектор L КУДМ неподвижен в системе x1 x2 x3, вектор L неподвижен в системе X1 X 2 X 3. Иллюстрировать характер движения КА будем с помощью проекций годографа L KA на координатные плоскости системы l1l2l3 (рис. 3.3), связанной с вектором суммарного кинетического момента [27]. Эта система определяется следующим образом: ось l2 совпадает по направлению с L ; ось l1 лежит в плоскости l2 X 1, перпендикулярна L и образует острый угол с осью X 1 ; ось l3 дополняет систему координат до правой.

На рисунках 3.43.13 представлена реконструкция вращательного движения КА указанным способом. Все рисунки скомпонованы одинаково, они содержат проекции годографа L KA на плоскости системы координат l1l2l3, графики величин 1 (t ), 2 (t ), 3 (t ), I (t ). На графиках угловых скоростей КА применены условные обозначения: круглые маркеры чёрного цвета – 1, круглые маркеры белого цвета – 2, квадратные маркеры – 3. На графиках расчётного тока СБ I (t ) маркерами указаны аппроксимируемые данные измерений (3.3).

На всех рассмотренных интервалах, за исключением двух последних, вращательное движение КА носит одинаковый характер. Отличие состоит лишь в изменении ориентации вектора суммарного кинетического момента двигателей-маховиков в связанных осях КА и как следствие – изменение относительной ориентации векторов L KA и L. На последних двух интервалах движение КА более медленное, так как скорости вращения двигателей-маховиков были значительно снижены. В таблицах 3.6, 3.7 представлены некоторые величины, рассчитанные на основе полученных при минимизации параметров движения КА. В таблице 3.6 приведены угловые скорости КА и двигателей-маховиков. Можно отметить, что модуль угловой скорости КА находился в диапазоне 0.60.75 град/с, наибольшую угловую скорость КА имел по оси Ox3. В таблице отсутствуют данные по угловой скорости вращения третьего двигателя-маховика, так как этот маховик был остановлен. В таблице 3.7 приведены модули векторов L КУДМ и L, и компоненты соответствующих ортов. В таблице также приведены углы между векторами L КУДМ и L КА ( L КУДМ ^ L КА ), и между векторами L и S ( L ^ S ). Последний угол на всём рассматриваемом интервале времени (около трёх месяцев) медленно менялся в диапазоне от ~130 до ~145. Это изменение вызвано воздействием на вращающийся КА внешних факторов:

световое давление, гравитационные моменты и др.

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

Глава 4. Разработка модели вращательного движения Для разрабатываемого КА дистанционного зондирования Земли была выбрана конструкция привода панелей солнечной батареи, позволяющая в ходе полёта КА изменять конфигурацию его панелей солнечной батареи. Помимо этого, на рассматриваемом КА установлена целевая аппаратура (ЦА), масса подвижной части которой составляет существенную долю массы КА. Возникла задача разработать математическую модель вращательного движения КА с учётом его указанных конструктивных особенностей. При этом необходимо было обеспечить простоту программной реализации созданной модели и её независимость от каких-либо математических программных пакетов. В ходе разработки математической модели, уравнения вращательного движения КА получены в форме уравнений Эйлера-Лагранжа и приведены к виду, удобному для их численного интегрирования на ЭВМ. Исследованы малые колебания тел, составляющих КА, в окрестности положения равновесия, которое существует при отсутствии внешних моментов; найдены собственные частоты колебаний. Разработанная математическая модель, реализована на языке программирования Borland С++. Программная реализация представляет собой динамическую библиотеку, вместе с которой поставляется тестовая задача и пример использования указанной библиотеки в программе с графическим интерфейсом. Разработанная библиотека позволяет включить её как составную часть математической модели КА в состав стенда для отработки системы ориентации и стабилизации КА. Графический интерфейс реализует визуализацию результатов интегрирования уравнений движения КА и позволяет осуществлять настройку параметров модели, адаптируя её к изменениям параметров конструкции КА. Краткое описание разработанного программного обеспечения приведено в приложении 2. В качестве теста приведены результаты интегрирования уравнений движения в случае свободного движения КА с вращающимся зеркалом.

4.1 Вывод уравнений движения Уравнения движения КА [26, 28] записываются в виде уравнений Лагранжа 2 рода. Расчётная схема КА представлена на рис. 4.1. КА рассматривается как механическая система, состоящая из 7 элементов, включая собственно корпус КА, корневое звено СБ, четыре панели СБ и вращающееся зеркало ЦА. Соединение элементов системы, за исключением ЦА, осуществляется через двухстепенные шарниры, допускающие поворот одного элемента относительно другого вокруг двух взаимно перпендикулярных осей. Дополнительные степени свободы введены для моделирования упругих колебаний конструкции КА.

С каждым элементом конструкции связана своя система координат bi xi yi zi, ( i 1,...,7 ), начало которой лежит в точке bi, совпадающей с центром вращения соответствующего шарнира (рис. 4.1). При полностью раскрытых панелях солнечных батарей, соответствующие оси всех систем координат совпадают по направлению. В уравнениях будем учитывать три вращательные степени свободы КА, вращательную степень свободы зеркала ЦА, а также вращения остальных элементов системы в шарнирах.

Введём обозначения: 1, 2, 3 – угловые скорости КА вокруг осей ССК;

2,..., 7 – углы поворота в шарнирах корневого звена, панелей СБ и зеркала ЦА вокруг осей x2,..., x7 ; 2,..., 6 – углы поворота в шарнирах корневого звена и панелей СБ вокруг осей z2,..., z6. Нулевые значения всех угловых величин соответствуют положению КА с полностью раскрытыми недеформированными панелями СБ, а положительные направления угловых отклонений согласованы с соответствующими осями координат. В шарнирах, за исключением зеркала ЦА, присутствуют упругие элементы и вязкое трение. Далее будем использовать обозначение mi – масса i -го элемента, J i – тензор инерции i -го элемента в системе координат, связанной с центром масс i -го элемента, оси которой совпадают по направлению с осями системы bi xi yi zi. Массово-инерционные характеристики всех элементов системы известны. Обозначим bi (bi1, bi 2, bi 3 )T – радиус-вектор центра вращения i -го шарнира в системе координат bk xk yk zk элемента, являющегося базовым для i -го шарнира ( i 2,...,7 ). Здесь для i 2,...,6 : k i 1, для i 7 : k 1. Введём также обозначение сi (сi1, сi 2, сi 3 )T – радиус-вектор центра масс i -го элемента в системе координат bi xi yi zi. Численные значения компонент векторов ci и bi известны.

Введём обозначение []n, означающее, что некоторый вектор задан в системе координат bn xn yn zn n -го элемента. Для перехода от одной системы координат к другой введём матрицы перехода Ank таким образом, что [a]n Ank [a]k.

Выражение для кинетической энергии движения рассматриваемой системы относительно общего центра масс записывается в виде где [k ]k (1k, 2k, 3k )T – вектор угловой скорости k -го элемента в собственной системе координат, [vk ]1 [k ]1 [C ]1 – скорость центра масс k-го элемента относительно общего центра масс системы, Выражения для угловых скоростей элементов в собственных системах координат:

1. для k 2,...,6 : [k ]k [k 1 ]k k [ek1 ]k k [ek 3 ]k, где [ek1 ]k (1,0,0)T, [ek 3 ]k (0,0,1)T, 2. [7 ]7 [1 ]7 7 [e71 ]7, где, [e71 ]7 (1,0,0)T.

Здесь ek1, ek 3 – орты осей вращения шарниров, соединяющих k элемент с базовым, k, k – угловые скорости в шарнирах (рис. 4.1).

Выражения для линейных скоростей центров масс элементов в системе координат, связанной с корпусом КА:

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

Первый поворот осуществляется вокруг оси xk на угол k, второй поворот (для двухстепенных шарниров панелей СБ) осуществляется вокруг новой оси zk на угол k.

Согласно приведённым выражениям для i и i, связи, наложенные на систему, стационарны. Используя приведённые выражения, выразим матрица кинетической энергии, x (11, 21, 31, q)T – вектор обобщённых обобщённых координат системы.

Матрица M (q) строится на основе выражения для кинетической энергии движения системы относительно общего центра масс. Чтобы, располагая выражениями для i и i, найти элементы матрицы M (q), необходимо иметь возможность выполнять операции над трёхмерными векторами и линейными формами обобщённых скоростей. Для этого введём два типа векторов и скаляров [35].

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

Например, к первому типу относятся радиус-векторы отдельных точек системы и скалярные произведения таких векторов. Примером векторов первого типа являются вектора ei, bi, ci. Вектора (скаляры) первого типа представляются одномерными массивами, длины 3. Элементы этих массивов (числа), в общем случае, рассчитываются для заданных значений обобщённых координат и скоростей.

Ко второму типу относятся векторы (скаляры), зависящие как от обобщённых координат, так и от обобщённых скоростей, причём зависимость от последних должна быть линейной и её требуется указывать в явном виде. Примерами векторов второго типа могут служить скорости отдельных точек системы, угловые скорости элементов системы. Векторы второго типа представляются двумерными массивами – матрицами, размера 3 n, где n – число степеней свободы системы. Чтобы из такой матрицы получить трёхмерный вектор, её необходимо умножить справа на вектор обобщённых скоростей. Скаляр второго типа представляется одномерным массивом длинны n. Элементы указанных массивов рассчитываются для заданных значений обобщённых координат, значения обобщённых скоростей здесь не рассматриваются. При заданных значениях последних вектор второго типа можно преобразовать в вектор первого типа. Векторы второго типа используются только для получения матрицы M (q), построить которую можно, выполнив над векторами второго типа соответствующие операции.

Представим скаляр второго типа s, выражающий линейную скалярную функцию величины s от n мерного вектора производных обобщённых координат в виде строки s (s1,..., sn ), где координата si отражает зависимость s от i компоненты вектора производных обобщённых координат. Трёхмерный вектор второго типа v, соответствующий трёхмерной величине v представляется в виде строки обычных трёхмерных векторов (в виде матрицы 3 n ) – v (v1,..., vn ), где координата vi отражает зависимость v от i компоненты вектора производных обобщённых координат.

Далее для всех векторов и скаляров второго типа будем использовать введённое здесь обозначение..

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

сложение, вычитание, умножение на число, результат – скаляр второго типа, умножение на вектор первого типа, результат – вектор второго типа, умножение на вектор из пространства x, результат – число. С векторами второго типа разрешены следующие операции: сложение, вычитание, умножение на число, векторное умножение на вектор первого типа, результат – вектор второго типа, скалярное умножение на вектор в E 3, результат – скаляр второго типа, перевод в другую систему координат путём умножения на соответствующую матрицу 3 3.

Используя введённые типы векторов, найдём выражения для элементов матрицы M (q). Введём вектора второго типа k1, k 3 и 1 такие, что k ek1 k1x, k ek 3 k 3 x, 1 1x. Как и в случае с обычными векторами, введём обозначение [ ]n, означающее, что некоторый вектор второго типа задан в системе координат bn xn yn zn n -го элемента. Для перехода от одной системы координат к другой используются введённые ранее матрицы перехода Ank, что [ ]n Ank [ ]k. Введённые вектора второго типа представляют собой блочные матрицы вида:

[ k 3 ]k 0,..., 0, ek 3, 0,..., 0, для k 2,...,6. Здесь 0 – нулевой вектор длинны 3.

С учётом записанных ранее выражений для линейных и угловых скоростей, запишем выражения для векторов второго типа k, k и v k, таких, что k k x, k k x, vk vk x. В качестве примера выпишем ещё раз выражение для 2 : [2 ]2 [1 ]2 2[e21 ]2 2[e23 ]2. Используя введённые ранее вектора второго типа, перепишем последнее выражение в виде [2 ]2 [1 ]2 x [ 21 ]2 x [ 23 ]2 x. Учитывая что 2 2 x, запишем выражение для 2 : [2 ]2 [1 ]2 [ 21 ]2 [ 23 ]2.

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

Для k 2,...,6 :

Выражения для линейных скоростей относительно общего центра масс системы: [vk ]1 [k ]1 [C ]1, k 2,...,7.

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

Последнее выражение можно переписать в виде:

матрицы кинетической энергии системы в функции вектора q обобщённых координат системы.

Для управления вращательным движением КА используются 3 двигателя-маховика, жёстко установленные на КА таким образом, что положительные направления угловых скоростей роторов маховиков совпадают с положительными направлениями осей ССК. Чтобы учесть динамику собственного вращения роторов двигателей-маховиков и влияние этого вращения на движение КА, уравнение для кинетической энергии системы необходимо дополнить выражением вой момент инерции роторов двигателей-маховиков, – угол поворота ротора маховика вокруг оси собственного вращения. Коэффициентами d ki, ввиду их малости, можно пренебречь, кроме того, для полого учёта динамики двигателей-маховиков, их массово-инерционные характеристики должны также быть учтены в инерционных характеристиках КА.

Уравнения движения КА относительно центра масс запишем в виде:

Здесь M1,2,3 – компоненты момента внешних сил, действующих на КА (момент вычисляется относительно центра масс КА), M k – суммарные моменты, приложенные к осям вращения маховиков, сi и i – диагональные элементы матриц квадратичных форм потенциальной энергии и диссипативной функции взаимодействия элементов системы.

Для интегрирования выписанных уравнений вращательного движения на ЭВМ более удобна запись указанных уравнений относительно обобщённых импульсов:

h – вектор кинетического момента комплекса двигателей-маховиков, c (с1, с2,..., c10,0) – вектор коэффициентов упругости в шарнирных соединениях, q0 – вектор угловых отклонений в шарнирах в невозмущённом положении, (1, 2,..., 11,0)T – вектор коэффициентов вязкого трения в шарнирах, mмаг, mаэр, mграв, mшар – момент магнитных катушек, аэродинамический, гравитационный, шарнирный моменты.

Надо отметить, что в программной реализации модели величины k не включались в вектор q. Это сделано для обеспечения изолированности программных кодов моделей движения КА и моделей двигателеймаховиков для того, чтобы иметь возможность их раздельной компиляции.

В процессе интегрирования уравнений движения системы, вектор x угловых скоростей системы находится из соотношения x M 1P, угловое положение корпуса КА в ИСК находится по стандартным кинематическим соотношениям [46], вектор q угловых отклонений в шарнирах находится из соотношения q x.

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

4.2 Расчёт собственных частот модели Как уже отмечено выше, разработанная модель представляет КА как систему из 7 твердых тел. Она включает корпус КА, корневое звено, 4 панели солнечных батарей и зеркало. Зеркало представляет собой осесимметричный ротор, который может вращаться относительно своей оси симметрии, жестко закрепленной на корпусе. Корпус КА, корневое звено и панели батарей образуют линейную цепочку. Звенья цепочки соединены двухстепенными шарнирами. В шарнирах действуют упругие восстанавливающие моменты так, что в равновесной конфигурации корневое звено, панели батарей и ось зеркала лежат в одной плоскости. Малые собственные колебания такой системы относительно центра масс описываются dim 14, i – компоненты вектора бесконечно малого поворота корпуса в системе координат c1x1 y1z1, M – матрица кинетической энергии системы, вычисленная для равновесной конфигурации, C – матрица жесткостей системы, у которой 10 диагональных элементов, отвечающих обобщенным координатам 2, 2,, 6, 6, равны коэффициентам упругости в шарнирах, а остальные элементы равны нулю.

Общее решение выписанного уравнения имеет вид Здесь B0 и B1 – векторы размерности 14, у которых компоненты, отвечающие координатам 2, 2,, 6, 6 равны нулю, а остальные компоненты – произвольные постоянные, ck – произвольные постоянные, k и Ak – решения задачи на собственные значения 2 MA CA, 0. Величины k – резонансные частоты упругих колебаний системы.

Решения указанной задачи на собственные значения находились с помощью стандартной процедуры линейной алгебры [47]. Ее решения приведены в табл. 4.1, 4.2. Эти таблицы получены при одной и той же матрице C и разных матрицах M. Табл. 4.1 получена при бесконечно большой массе корпуса КА. Бесконечно большая масса корпуса означает заделку шарнира, связывающего корневое звено с корпусом, в неподвижную стенку. Такая заделка отвечает расчету, выполненному методом конечных элементов, для рассматриваемой конфигурации батарей. Проведение расчетов при указанном выборе матрицы M по нашей модели позволило выбрать правильные значения коэффициентов жесткости в шарнирах. Коэффициенты жесткости, отвечающие углам 2,, 6, были приняты равными 412.525 Нм, коэффициенты, отвечающие углам 2,, 6, были приняты равными 97.92 Нм. Табл. 4.1, 4.2 получены при этих значениях.

В соответствии с исходными данными, низшие тона колебаний панелей батарей составляют 0.227 Гц (изгибные колебания) и 0.729 Гц (крутильные колебания). Данные, приведённые в табл. 4.1, показывают соответствие низших частот предложенной модели этим значениям. Табл. 4. содержит собственные частоты в случае номинального значения массы корпуса КА. В этом случае изгибные колебания батарей происходят с частотой 0.493 Гц, крутильные колебания – с частотой 0.736 Гц. Существенное изменение частоты изгибных колебаний по сравнению с крутильными колебаниями объясняется меньшим отношением момента инерции КА относительно оси c1 x1 к моментам инерции панелей СБ вокруг осей bi xi по сравнению с аналогичным соотношением для моментов инерции в крутильном направлении (вдоль осей bi zi ).

5. Заключение Основные результаты, полученные в диссертации:

1. Разработана математическая модель управляемого вращательного движения МКА «Монитор-Э» в режиме гашения угловой скорости при условии отсутствия измерений компоненты угловой скорости относительно одной из связанных осей аппарата. Найдены стационарные решения модельных уравнений движения и исследована их устойчивость. Построена оценка областей притяжения стационарных решений модельных уравнений. Дана рекомендация, обеспечивающая успешное гашение угловой скорости МКА.

2. Разработаны и реализованы в виде программы для персонального компьютера две интегральные статистические методики реконструкции вращательного движения МКА по телеметрической информации определенного вида. Одна из этих методик позволила по данным измерений двух компонент угловой скорости МКА и суммарного кинетического момента двигателей-маховиков реконструировать вращательное движение МКА в нескольких реализациях указанного в п. 1 режима гашения угловой скорости. Реконструкция подтвердила адекватность разработанной математической модели и эффективность выданной рекомендации. С помощью второй методики выполнена реконструкция вращательного движения МКА «Монитор-Э» и «КазСат-1» по телеметрическим значениям тока, снимаемого с солнечных батарей. Созданное программное обеспечение использовано в инженерном сопровождении летных испытаний указанных МКА при парировании нештатных ситуаций.

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

6. Приложение 6.1 Программа для мониторинга состояния КА В соответствии с разработанной методикой на языке Pascal в среде Turbo Delphi была разработана программа для персонального компьютера.

Программа имеет графический интерфейс пользователя и состоит из одного исполняемого файла. Главное окно программы представлено на рисунках 6.1 и 6.2. В верхней части окна находятся закладки для перехода к страницам «Расчёт» и «Результат». На рисунке 6.1 показана страница «Расчёт», на рисунке 6.2 – страница «Результат». Страницу «Расчёт»

условно можно разделить на левую и правую части. Левая часть предназначена для работы с исходными данными, в правой части находятся инструменты для реализации начального поиска параметров движения и последующего уточнения этих параметров с помощью минимизации. В левой верхней части страницы «Расчёт» имеется область «Характеристики КА и параметры расчёта». Здесь пользователь, нажав на кнопку «Загрузить из файла», может из специально подготовленного текстового файла загрузить в программу параметры тензора инерции КА, параметры солнечной батареи, установленной на КА, признаки работоспособности двигателеймаховиков, а также некоторые параметры процедуры численного интегрирования уравнений движения КА.

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

1-я строка: шаг интегрирования модели; характерная точность; регуляризирующий множитель для кинетических моментов двигателеймаховиков (ДМ); минимальный учитываемый ток СБ; максимальный ток СБ.

2-я, 3-я и 4-я строки: тензор инерции КА.

5-я строка: признаки работы ДМ (4 признака, 1 – работает, 0 – не работает).

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

Ниже, в левой части окна, находится область для работы с исходными телеметрическими данными о токе СБ. Телеметрические данные представляют собой файл с расширением «mat» (в формате системы Matlab), версия формата 5.0. Данные в файле располагаются по столбцам. Первый столбец – время в формате TDateTime среды разработки Delphi. В этом формате дата и время кодируются вещественным числом, причём целая часть соответствует числу дней, разделяющих кодируемую дату и 31 декабря 1899 года, а дробная часть – кодирует долю суток, прошедшую с полуночи. Второй столбец – ток солнечной батареи в амперах. После загрузки данных, они начинают отображаться на графике под кнопкой «Загрузить» и кнопками выбора отображаемого интервала времени. По оси абсцисс отложено время в формате календарной даты и времени суток чч:мм:сс. По оси ординат – измеренные значения тока солнечных батарей в амперах. Вертикальные красная и синяя линии на графике показывают, соответственно, начало и конец отрезка данных, которые будут использованы для идентификации вращательного движения КА. После выбора интересующего отрезка данных, они принимаются к расчёту. Помимо выбора начальной и конечной точек, пользователь может выбрать также интервал прореживания, т.е. кратность выборки точек из исходных данных. После того, как данные выбраны и загружены, в нижней части окна программы появляется точечный график черного цвета с загруженным отрезком данных.

В правой части страницы «Расчёт» располагаются элементы управления для реализации методики поиска начального приближения (область «поиск начального приближения») и последующего уточнения решения с помощью метода Левенберга-Марквардта (область «минимизация»). В части поиска начального приближения пользователю предоставлено две таблицы, одна из которых отображает начальный вектор искомых величин, а также параметры, используемые при начальном подборе параметров методом случайного поиска с обучением. Вторая таблица отражает текущее значение искомого вектора в процессе работы алгоритма начального поиска, а также текущее значение интервала поиска. Желтым цветом выделены параметры, относящиеся к ДМ и к конструктивным параметрам КА. Справа от таблиц расположен элемент управления для задания начального угла между нормалью к СБ и направлением на Солнце, а также две вертикальные шкалы, отображающие степень завершения процесса начального подбора вектора состояния. Запуск алгоритма начального подбора осуществляется нажатием на кнопку «Start», расположенную под правой таблицей.