Космические исследования, 2020, T. 58, № 4, стр. 321-330

“НЕПОЛНАЯ” КРАЕВАЯ ЗАДАЧА

В данной работе, в отличие от вышеупомянутых, краевая задача решается в исходной постановке, без осреднения уравнений движения, как и в [8]. В работе [8] задача рассматривалась без учета влияния тени (и других возмущающих факторов) в двух постановках: (А) с возможностью многократного включения и выключения тяги, и (Б) с постоянно работающей, постоянной по величине малой тягой. А именно, рассматривалась задача оптимального многовиткового перелета космического аппарата в центральном ньютоновском поле Земли с электроракетными двигателями малой тяги, с питанием от солнечных батарей. Управление КА осуществлялось с помощью вектора малой тяги. На направление тяги не накладывалось никаких ограничений. С использованием принципа максимума Понтрягина задача сводилась к двухточечной краевой задаче. Для этого вводились сопряженные переменные, формировался гамильтониан задачи, с его помощью находились выражения для оптимального направления вектора тяги, как функции фазовых и сопряженных переменных, и дифференциальные уравнения для фазовых и сопряженных переменных

где t – время, P – величина тяги, x = 〈φ, m, h, y, z, v, w〉 – вектор фазовых переменных, λ = 〈λ_φ, λ_m, λ_h, λ_y, λ_z, λ_v, λ_w〉 – вектор сопряженных переменных. Здесь m – масса КА, а φ, h, y, z, v, w – равноденственные переменные, которые определяются через гравитационную постоянную Земли μ_Е и оскулирующие переменные p, e, θ, ω, Ω, i формулами:

Задача сводилась к поиску начальных значений сопряженных переменных, которые обеспечили бы выполнение краевых условий в конечный момент времени. Таким образом получалась краевая задача в случае (Б). В случае (А), вдобавок к выражениям для оптимального направления тяги, была получена, также с помощью принципа максимума, т.н. функция переключения Sw, которая дает критерий для оптимальных моментов включения и выключения тяги: при Sw < 0 тяга нулевая, при Sw > 0 тяга максимальная (особые случаи не рассматривались).

В данной работе используется постановка (Б), но добавляется условие обнуления тяги в области тени: есть функция тени Sh, при Sh < 0 тяга нулевая, при Sh > 0 тяга максимальная (случаи движения по границе тени рассматривать не будем). Для модели “цилиндрической” тени функция Sh конкретизирована в [9]:

где r – радиус-вектор КА, R_S – радиус-вектор Солнца, r = |r|, R_S = |R_S|, R_Е – радиус Земли, h_min – заданная константа (например, высота атмосферы).

Получающаяся задача очень похожа на задачу (А), но вместо функции переключения Sw величина тяги “задается” функцией Sh. Может возникнуть впечатление, что простой заменой одной функции на другую мы получим новую краевую задачу, решая которую сможем получать искомые траектории. Такая “краевая” действительно дает возможность находить траектории перелета на конечную орбиту с обнулением тяги в области тени, но эти траектории не оптимальны, ибо не учитываются условия оптимального пересечения границ тени. Это т.н. условия скачка в сопряженных переменных [10]:

где G = Sh, λ^– и λ⁺ – значения сопряженных переменных слева и справа от границы тени, а μ – параметр, для определения которого необходимо численно решать уравнение, в явном виде относительно μ неразрешимое.

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

Основные причины, по которым было решено использовать для поиска траекторий “неполную” краевую задачу, следующие. Во-первых, такая задача существенно проще “полноценной” краевой задачи. На многовитковых траекториях могут быть десятки и сотни участков с тенью [9], и соответственно раза в два больше условий скачка. Во-вторых, предполагалось, что, как и в случае учета влияния других возмущений (гравитационного притяжения Луны, Солнца и вариаций геопотенциала [11–13]), ее решения будут не сильно отличаться от номинальных. А если получаемые решения по функционалу (расходу рабочего вещества) окажутся близки к номинальным траекториям – их можно будет рассматривать как приближенно-оптимальные и, возможно, использовать вместо оптимальных решений.

НАЧАЛЬНЫЕ ОРБИТЫ

Уточним постановку задачи. В качестве конечной орбиты рассматривается геостационарная орбита (ГСО), а в качестве начальной – орбита с наклонением i₀ = 13°, высотой перигея H_π = 9.2 и апогея H_α = 76.8 тыс. км. Эти три параметра задают конфигурацию орбиты. А ее положение в пространстве определяют параметры ω₀, Ω₀. Как указано в [14], в случае оптимальных по быстродействию перелетов между некомпланарными начальной эллиптической и конечной круговой орбитами со свободной линией апсид начальной орбиты, из условий трансверсальности следует, что линия апсид должна совпадать с линией пересечения плоскостей начальной и конечной орбит. В [4] также отмечено, что задача перелета между такими орбитами типична при выведении КА на высокие круговые орбиты. Поэтому зафиксируем угловое расстояние перицентра от узла ω₀ = 0°. Другие характеристики орбиты приведены в табл. 1. Дополнительно полагается, что КА стартует из перигея начальной орбиты, т.е. истинная аномалия θ₀ = 0° (следовательно, φ₀ = Ω₀), а конечный момент времени определяется заданной угловой дальностью перелета, а именно, целым количеством витков N в переменной φ (φ_к = φ₀ + 2πN). Как правило, значение N подбирается оптимальным, для чего приходится решать несколько краевых задач. В результате у нас остается два свободных параметра: долгота восходящего узла Ω₀ и время (дата) старта t₀ (измеряется в сутках от 00.00 1.I.2018).

НОМИНАЛЬНЫЕ ТРАЕКТОРИИ

Анализ влияния тени будем проводить сравнением многовитковых (с наилучшим количеством витков N) траекторий перелета на геостационарную орбиту с тенью и траекторий номинальных, т.е. оптимальных траекторий перелета без выключения тяги в тени (со своим наилучшим N). Для последних, в силу центральности ньютоновского поля и отсутствия возмущений, имеется симметрия по отношению к параметру Ω₀: изменение Ω₀ приводит к повороту начальной орбиты и траектории перелета как целого вокруг оси вращения Земли. Параметр t₀ также не влияет на траекторию перелета, поскольку от него не зависят ни положение начальной орбиты в пространстве, ни положение КА на орбите в начальный момент времени (в силу принятых выше допущений). Поэтому характеристики номинальных траекторий одинаковы для всех значений параметров Ω₀, t₀. Они приведены в табл. 2.

ПОЛОЖЕНИЕ УЧАСТКОВ С ТЕНЬЮ НА ТРАЕКТОРИЯХ

В отличие от номинальных, перелеты с нулевой тягой в области тени зависят от положения теневых участков на траектории, а значит, зависят от положения Солнца относительно Земли в небесной системе координат (и, следовательно, от t₀), и от положения траектории относительно Солнца (а следовательно, и от Ω₀).

Некоторое общее представление о том, как с изменением параметров меняется количество участков с тенью и на каких витках они расположены, можно получить из т.н. “портретов” номинальных траекторий с тенью (рис. 1, 2). По оси ординат отсчитываются номера витков, по оси абсцисс – либо параметр t₀ (при заданном значении Ω₀), либо параметр Ω₀ (при заданном значении t₀), и черными квадратиками отмечены витки с тенью.

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

РЕЗУЛЬТАТЫ РАСЧЕТОВ

В табл. 3 приведены результаты решения “неполной” краевой задачи для отдельных значений параметра t₀. Дата старта изменяется в интервале года с шагом 1.5 мес., угол Ω₀ – в диапазоне [0°–350°] с шагом 10°. Приведены величины “дополнительных” затрат рабочего вещества ∆М_РВ на перелет по сравнению с номинальными траекториями в процентах от номинальных затрат М_РВ (485.3 кг). Если считать, что дополнительные затраты, не превышающие 2% – это “хороший” результат, то, судя по данным табл. 3, высказанные выше ожидания оправдались на 2/3. В 45 вариантах (≈1/6 часть) затраты на перелет с тенью оказались меньше, чем на номинальной траектории. Но треть результатов превышает 2% (они выделены жирным), причем есть ∆М_РВ > 30%.

СЛАБОЕ ВЛИЯНИЕ ТЕНИ НА ВИТКАХ В КОНЦЕ ТРАЕКТОРИИ

В табл. 3 особо выделяются два столбца, с t₀ = = 91 и 274 суток, в которых все варианты “хорошие”. На первый взгляд, то, что ΔМ_РВ слабо меняется с изменением Ω₀ во всем диапазоне значений, выглядит странно, особенно в сравнении с другими столбцами, в которых изменения большие. Но это частично, а именно в отношении того, почему группа витков в конце траектории с изменением параметра Ω₀ почти не влияет на функционал, объясняется просто.

Даты t₀ = 91 и 274 суток (от начала года) дней на 10 позже весеннего и осеннего равноденствия. Продолжительность номинального перелета 180 суток – около полугода, поэтому завершение перелетов – лишь на неделю позже осеннего и весеннего равноденствия, когда Солнце находится вблизи экватора. Поэтому витки в конце траектории, имеющие малое наклонение, пересекают цилиндр тени. И это хорошо видно для t₀ = 91 на рис. 2: группа витков с тенью в конце траектории состоит из 23–32 витков, в зависимости от Ω₀ (для t₀ = 274 картина похожая).

Но эти витки имеют и малый эксцентриситет. Изменение параметра Ω₀ означает поворот номинальной траектории как целого вокруг оси вращения Земли. На это можно посмотреть по-другому: что Солнце вращается в плоскости экватора вокруг траектории, так что участки тени смещаются по кругу вдоль почти круговых витков. Угловые размеры теневых участков ≈16°. Насколько ситуация близка к симметричной, видно из рис. 3, где показаны “цилиндр” тени и каждый седьмой виток из группы последних витков с тенью со 129-го по 157-й на номинальной траектории с t₀ = 274 суток, Ω₀ = 0°.

Отметим, что в данном случае группа витков с тенью состоит из 29 витков – это 18.5%, или почти пятая часть от всей угловой дальности перелета, а по продолжительности – 32.4 суток, или 17.7%. На траектории, реализующей такой перелет (с нулевой тягой в области тени), количество витков с тенью увеличивается до 31 (127‑й–157-й витки), ΔМ_РВ = 0.45% (2.19 кг). А, например, на траектории с Ω₀ = 50° группа витков с тенью состоит из 34 витков (124-й–157-й витки), ΔМ_РВ = –0.24% (–1.15 кг). По угловой дальности это 21.7%, по продолжительности перелета 19.4%.

Из рис. 2 также видно, что почти на трети диапазона изменения Ω₀ других витков с тенью, кроме группы конечных, нет. В этом диапазоне ∆М_РВ < 0.85%. За пределами этого диапазона количество витков с тенью в начале траектории возрастает до 38 при Ω₀ = 275°–295°, причем ∆М_РВ достигает максимума 2% при Ω₀ ≈ 290°. Были проведены расчеты, в которых влияние тени из первой группы витков было полностью “отключено”, т.е. на этих участках тяга не обнулялась. Они показали, что в таком случае ∆М_РВ изменяется в диапазоне –0.146…0.842%. Т.е. вклад в ∆М_РВ участков тени на витках в конце траектории действительно невелик, и он не сильно зависит от Ω₀.

ОСОБОЕ ПОЛОЖЕНИЕ НАЧАЛЬНОЙ ОРБИТЫ ОТНОСИТЕЛЬНО ОРБИТЫ СОЛНЦА

Из табл. 3 видно, что наряду с “хорошими” столбцами есть несколько строк, в которых также все или почти все варианты “хорошие”. Это строки с Ω₀, равным 160°–190°, 350° и 0°. Среди них особо выделяется строка с Ω₀ = 180°, в которой значения ΔМ_РВ не превышают 0.56% (2.7 кг). Ясно, что дело, скорей всего, именно в начальной орбите, представляющей собой сильно вытянутый эллипс с эксцентриситетом e = 0.685 (радиус перигея 15.6, апогея 83.2 тыс. км). К тому же наибольшее наклонение на многовитковой траектории перелета на ГСО именно в начальный момент времени, в данном случае i = 13°.

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

В чем особенность положения начальной орбиты с Ω₀ = 180°, хорошо видно на рис. 4. На нем, в проекции на плоскость 〈Y, Z〉 небесной системы координат, показаны: плоскость эклиптики, в которой находится орбита Солнца, плоскость экватора, в которой находится ГСО, плоскость начальной орбиты для Ω₀ = 180° и плоскость начальной орбиты для Ω₀ = 0°. Все эти плоскости проходят через ось X (на рисунке в проекции дает точку O на пересечении осей Y, Z, направлена в сторону от зрителя).

Плоскости начальных орбит с Ω₀ = 0° и 180° расположены симметрично относительно плоскости экватора. Большие оси орбит – линии апсид – совпадают с осью X, но направлены противоположно друг другу. При движении Солнца по орбите цилиндр тени заметает окрестность плоскости эклиптики, пересекая начальные орбиты вблизи оси X, т.е. вблизи их линий апсид, в окрестности перигея и апогея. Линии апсид являются осями симметрии эллипса, и именно это является главной особенностью орбит с Ω₀ = 0° и 180°, т.к. возмущения от участков тени, расположенных вблизи оси симметрии, по-видимому, легче нивелировать.

Но между двумя этими орбитами есть различия. Во-первых, это противоположные направления больших осей. Во-вторых – разное положение плоскостей орбит по отношению к плоскости эклиптики. Для орбиты с Ω₀ = 180° угол между плоскостью орбиты и плоскостью эклиптики равен 23.5° + 13° = 36.5°, для орбиты с Ω₀ = 0° он равен 23.5° – 13° = 10.5°. А чем меньше этот угол, тем больший отрезок орбиты заметается цилиндром тени. Например, при наклонении орбиты i = 23.5° угол будет равен 0°, т.е. начальная орбита с Ω₀ = 0° будет лежать в плоскости эклиптики, и цилиндр тени будет заметать всю орбиту. В результате по затратам ΔМ_РВ строка с Ω₀ = 0° в табл. 3 оказывается хуже, чем с Ω₀ = 180°.

ВОЗМОЖНОСТЬ СУЩЕСТВОВАНИЯ НЕСКОЛЬКИХ РЕШЕНИЙ

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

На рис. 5, 6 представлены графики зависимости ∆М_РВ и ∆N от Ω₀ для двух ветвей решений с t₀ = 183. Основная ветвь, решения которой существуют во всем диапазоне изменения параметра Ω₀, помечена квадратиками: белыми – на большей части ветви, черными – где функционал больше, чем на второй ветви. Решения второй ветви существуют в небольшом диапазоне углов и дают лучшие решения для Ω₀ в диапазоне 70°–130° (помечены крестиками).

Особенностью первой ветви является то, что она не замкнута, хотя ее решения существуют для всех значений Ω₀: при приближении Ω₀ к 90° слева мы получаем одно решение, при приближении Ω₀ к 90° справа – другое решение. В результате для Ω₀ = 90° имеем три решения. Они приведены в табл. 4. Причем в данном случае все три решения существуют не только для оптимальных N (помечены звездочкой), но и для N = 157 (как на номинальной траектории).

ПРОБЛЕМА НАЧАЛЬНОГО ПРИБЛИЖЕНИЯ

При решении краевых задач возникает проблема выбора хорошего начального приближения. В рассматриваемой задаче это могут быть решения с близкими значениями даты старта t₀ или долготы восходящего узла Ω₀. Однако, как следует из рис. 5, 6, близость параметров не всегда означает близость решений: имея решение с Ω₀ < 90°, принадлежащее левой части первой ветви, получить решение с Ω₀ > 90°, принадлежащее правой части этой же ветви, не удастся. Тем более, не удастся получить решение, принадлежащее другой ветви (разве что случайно).

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

Если бы изменение тяги в области тени было в несколько раз меньше, с номинального (максимального) не до нуля, а, например до 0.9 номинального значения, то и невязки уменьшились бы на порядок, и краевую задачу, возможно, удалось бы решить. Полученное решение затем можно было бы использовать, чтобы получить решение с тягой в 0.8 от номинального значения, и т.д. Т.е. решить серию краевых задач, “выключая” тягу на теневых участках постепенно, пока не получим искомое решение.

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

На рисунке представлены данные по траекториям с оптимальным количеством витков, последние витки отмечены белыми кружками. Черными кружками отмечены витки с участками тени (при данном масштабе – каждый второй и крайние в группе). В этом примере “эволюция” решения краевой задачи довольно существенна. На траектории с максимальной тягой в области тени количество витков с участками тени – 54 (это витки с 76-го по 129-й), при общем количестве витков 157. На траектории с нулевой тягой в области тени количество витков с тенью – 24 (это витки с 50-го по 73-й), при общем количестве витков 146. То есть общее количество витков уменьшилось на 11 (с 157 до 146), количество витков с тенью – на 30 (с 54 до 24).

Но надо отметить, что это очень “плохой” вариант – различие в затратах рабочего вещества 30.3%, или 147 кг (см. табл. 3 для t₀ = 137 и Ω₀ = 50°). Если решение краевой задачи очень сильно отличается от номинального, не всегда удается получить решение таким образом. В “хороших” случаях различие невелико и по функционалу – затратам рабочего вещества, и по другим характеристикам траекторий.

ТЕНЬ КАК БАРЬЕР

В отличие от других возмущающих факторов – гравитационного влияния Луны, Солнца, вариаций геопотенциала, которые представляют собой поле сил, рассматриваемое влияние тени таковым не является. Тень ведет себя, скорее, как барьер. А проявляется это, в частности, в незамкнутости ветви решений и сильном изменении количества витков, как на рис. 5, 6. Каков механизм такого воздействия?

В [8], где рассматривались перелеты без тени, было показано, что существует много решений краевой задачи, различающихся количеством витков. Причем в большом диапазоне изменения количества витков (он тем больше, чем больше наилучшая угловая дальность перелета на ГСО, или, в терминах данной статьи – чем больше угловая дальность номинальной траектории) функционал меняется не сильно. В случае постоянно работающей и постоянной по величине тяги это эквивалентно тому, что в этом диапазоне витков почти постоянна продолжительность перелета. Например, для траекторий перелета на ГСО с рассматриваемой здесь начальной орбиты (табл. 1), с количеством витков N в диапазоне от 151 до 166 (16 траекторий), отличие функционала от номинального (табл. 2) не превышает 1%. Для N из диапазона [149–171] (23 траектории) – менее 2%, из диапазона [144–183] (40 траекторий) – менее 5%, из диапазона [141–201] (61 траектория) – менее 10%.

И здесь напрашивается аналогия со спирально закрученной пружиной будильника. Только в случае с пружиной постоянна ее длина. При скручивании пружины увеличивается количество витков, а так как общая длина не изменяется – уменьшаются длина каждого витка и его размеры. Аналогично и с многовитковой траекторией: с увеличением количества витков при почти постоянной общей продолжительности перелета продолжительность каждого витка и его размеры уменьшаются. С уменьшением количества витков их размеры, наоборот, увеличиваются. Частично за счет этого можно “обогнуть” препятствие в виде цилиндра тени. А в зависимости от того, с наружной или внутренней стороны витков находится большая часть цилиндра тени, количество витков будет уменьшаться или увеличиваться. Именно этим, по-видимому, можно объяснить эффекты, проявляющиеся на рис. 5, 6. Хотя в реальности все сложнее, поскольку витков с тенью много, цилиндр тени от витка к витку смещается в пространстве, и продолжительность перелета может изменяться очень сильно.

МИНИМИЗАЦИЯ ЗАТРАТ ПО ПАРАМЕТРУ Ω₀

Определим для каждого t₀ наилучшее значение Ω₀ и соответствующее ΔМ_РВ.

Изменяя время старта t_С ракетоносителя в интервале одних суток, можно получить весь спектр значений параметра Ω₀. Подобрать t_С можно так, чтобы момент t(Ω₀) выхода КА на нужную начальную орбиту лежал в диапазоне t₀ ± 0.5. После чего может потребоваться пауза Δt продолжительностью от 0 до 1.264 суток, чтобы старт с начальной орбиты был в перигее. В зависимости от значений t(Ω₀) и Δt можно скорректировать время t_С так, чтобы время старта с орбиты отличалось от t₀ не более чем на 0.5 суток. Например, если t(Ω₀) = 0.5, а Δt = 0.25, 0.5, 0.75, 1.0 или 1.25 – можно уменьшить t_С на сутки или двое, и время старта с орбиты будет соответственно t₀ – 0.25, t₀, t₀ + 0.25, t₀ ± 0.5, t₀ – 0.25. Из-за таких изменений времени старта значения ΔМ_РВ будут отличаться от табличных не сильно, и совсем немного – для хороших вариантов.

На рис. 8 приведен график минимальных по параметру Ω₀ значений ΔМ_РВ. Расчеты проводились для t₀ с шагом 5 суток и Ω₀ с шагом 10°. Как видно, во всем диапазоне изменения даты старта t₀ затраты ΔM_РВ < 0, т.е. затраты на перелет с обнулением тяги в области тени меньше, чем на номинальной траектории. Однако выигрыш в затратах небольшой – менее 0.5%. Для сравнения приведен также график ΔМ_РВ для Ω₀ = 180°: для большинства значений t₀ затраты ΔМ_РВ > 0, но не превышают 0.62%.

Т.о., для перелетов на ГСО с рассмотренных здесь начальных орбит “неполная” краевая задача позволила получить траектории с нулевой тягой в области тени с затратами рабочего вещества меньшими, чем на траекториях номинальных – с постоянно работающей малой тягой. При этом перелеты с начальной орбиты с фиксированным значением параметра Ω₀ = 180° лишь ненамного хуже номинальных.