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

Аннотация

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

Полный текст

ВВЕДЕНИЕ Возможности применения двигателей малой тяги, в первую очередь, электроракетных двигательных установок (ЭРДУ), для обеспечения космических транспортных операций начали рассматриваться еще пионерами космонавтики в начале 20-го века. Интерес к этим двигателям и их использованию сохраняется до сих пор и продиктован их основным преимуществом - высоким удельным импульсом тяги, недоступным химическим двигателям. В связи с этим применение ЭРДУ в качестве основной двигательной установки, обеспечивающей наибольший вклад в транспортные операции космического аппарата (КА), гарантирует уменьшение расхода топлива. Однако из-за низкого уровня тяги, свойственного ЭРДУ, наиболее эффективно использовать их оказывается возможным только на достаточно больших удалениях от притягивающих объектов (планет и массивных спутников), т. е. в первую очередь на гелиоцентрических участках межпланетных перелетов. Так, в случае использования маршевого ЭРДУ в сильном гравитационном поле в окрестности массивного небесного тела, располагаемое реактивное ускорение может оказаться крайне малым по отношению к гравитационному ускорению притягивающего центра и быть на уровне 10-5-10-4. В этих условиях траектории перелета в окрестности массивных тел содержат большое количество витков и часто называются задачами перелета с малой тягой. На межпланетных траекториях уровень реактивного ускорения ЭРДУ не сильно уступает притяжению Солнца, и их отношение может иметь порядок 10-2-10-1. В этих случаях можно говорить о перелете не с малой, а с конечной тягой. Именно в задачах исследования Солнечной системы на рубеже веков ЭРДУ стали широко применяться в качестве маршевых. Первыми такими аппаратами стали Deep Space 1 (пролет астероида и двух комет), Smart-1 (выход на окололунную орбиту), Hayabusa (доставка образцов грунта с астероида Итокава), Dawn (последовательный перелет к астероидам Веста и Церера). При этом подходы к решению задач оптимального управления КА с двигателем малой тяги существенно отличаются от методов, применяемых к задачам с большой тягой (соответствующим КА с химическим маршевым двигателем). Для задач оптимизации траекторий с двигателями большой тяги, ввиду малой продолжительности активных участков, общепринятым является использование в расчетах импульсного приближения [11, 15], допускающего замену участка работы двигательной установки (ДУ) мгновенным изменением скорости КА. В этом случае задача оптимизации сводится к минимизации характеристической скорости и заключается в определении последовательности импульсов - величины, направления и времени их приложения, что является существенным упрощением. Кроме того, такая задача в случае движения в центральном гравитационном поле может быть сведена к системе нелинейных уравнений и вообще не требовать использования интегрирования уравнений движения для определения траектории КА. Задача оптимизации перелетов с двигателем малой тяги не позволяет использовать импульсное приближение, так как для нее продолжительности активных и пассивных участков оказываются сопоставимыми, что приводит к необходимости управления вектором тяги не в конечном числе точек, а в каждой точке траектории. Все методы оптимизации траекторий КА с малой тягой так или иначе сталкиваются при их практической реализации с рядом свойственных этой задаче проблем. Одна из них заключается в том, что обычно существует множество траекторий, удовлетворяющих необходимым и достаточным условиям оптимальности. Кроме того, задачи оптимального управления КА с малой тягой, как правило, не имеют аналитических решений и обычно могут быть решены только с помощью численных методов, использование которых осложняется разрывностью правых частей системы дифференциальных уравнений движения КА (при включении и выключении ДУ), наличием особой точки в притягивающем центре, их знакопеременностью и осцилляцией в случае использования декартовой системы координат. Все это может приводить к накоплению значительных ошибок при длительных периодах интегрирования и невозможности достоверного определения производных от краевых условий по начальным параметрам (зачастую необходимых для решения краевой задачи). Однако необходимо заметить, что на сегодняшний день эти проблемы достаточно легко устраняются благодаря сглаживанию релейной функции тяги и использованию методов интегрирования с переменным шагом и контролем точности. Куда более сложной и важной представляется проблема существования искомого решения, так как она играет важную практическую роль при решении конкретных задач. Теория необходимых условий оптимальности может оказаться бесполезной и не приводить к желаемому результату, если не выяснен вопрос о существовании решения в заданном классе допустимых функций, поскольку все необходимые условия априори предполагают существование экстремума. В случае отсутствия решения численные алгоритмы, основанные на необходимых условиях оптимальности, очевидно не приведут к построению оптимизирующей последовательности допустимых элементов даже тогда, когда необходимые условия позволяют получить все соотношения, необходимые для работы численного метода. На сегодняшний день получены важные результаты в виде теорем существования оптимального управления в классе измеримых функций управления и непрерывных траекторий, в том числе в работах: Филиппова 1959 года [21], Ли и Маркуса [14], Cesari [26], а также см. [2, 4]. Эти теоремы применимы и к рассматриваемым в данной работе задачам. Однако, равно как в условиях оптимальности обязательным элементом является предположение о существовании оптимального управления, в теоремах существования предполагается наличие пробного управления. В данной работе проанализированы две основные постановки задачи оптимального управления КА с малой тягой (разделы 2 и 3) и показано, что наличие пробного управления в них не всегда очевидно. Для устранения этой проблемы предложен метод построения границы области существования (разделы 5 и 6), а в качестве иллюстрации приведены результаты его использования (раздел 8). 102 А. В. ИВАНЮХИН 2. МАТЕМАТИЧЕСКИЕ МОДЕЛИ КОСМИЧЕСКОГО АППАРАТА С ОГРАНИЧЕННОЙ ТЯГОЙ Обычно математическая модель космического аппарата с ограниченной тягой включает в себя две основные части: динамические уравнения движения и массовую модель. В рамках данной работы космический аппарат рассматривается как материальная точка переменной массы, движение вокруг центра масс аппарата не рассматривается. Математическая модель движения центра масс в инерциальной декартовой системе координат, связанной с барицентром системы, может быть представлена в виде следующих дифференциальных уравнений ⎧ d2x T ⎨⎪ dt2 = Ωx + m, (2.1) - ⎩ = ⎪ dm dt |T | . c где x - вектор положения КА, t - время, m - масса КА, Ω - силовая функция гравитационного поля (нижний индекс обозначает производную), T - вектор тяги, c - скорость истечения. В случае ньютоновского, центрального гравитационного поля, силовая функция имеет вид: Ω = μ/r, где μ - гравитационный параметр центрального тела, r = |x| - удаление КА от центрального тела. В общем случае управлением является вектор тяги (его величина и направление) и скорость истечения: (T, c) ∈ R3 × R, c > 0. Часто для описания ДУ вместо скорости истечения используется удельный импульс тяги, равный отношению скорости истечения к стандартному ускорению свободного падения (g0 = 9,80665 м/с2): I = c . (2.2) g0 Реактивная мощность ЭРДУ N связана с тягой и скоростью истечения следующим соотношением: N = |T |c . (2.3) 2 Различия в математических моделях движения КА связаны с режимами функционирования ЭРД и ограничениями на управление. В общем случае параметры двигательной установки, такие как мощность, тяга и скорость истечения, могут быть функциями координат и времени. Традиционно рассматриваются следующие модели функционирования ЭРД (см. [9, 11, 15]): o идеально-регулируемый двигатель ограниченной мощности - управлением является величина и направление тяги, а также величина скорости истечения, при этом они ограничены только величиной располагаемой реактивной мощности: |T |c = N χ(x, t); (2.4) 2 0 o двигатель ограниченной тяги - управлением является направление тяги и ее значение, при этом направление тяги не ограничено, а величина тяги ограничена максимально допустимым значением, скорость истечения является заданной функцией: 0 :( |T | :( T0η(x, t), c = c0w(x, t), (2.5) где N0, T0 и c0 - значения мощности, тяги и скорости истечения в начальный момент времени, а χ(x, t), η(x, t), w(x, t) - заданные функции, при этом: χ(x, t)= η(x, t)w(x, t), ∀ t, x χ(x, t) > 0, η(x, t) > 0, w(x, t) > 0, χ(x(t0), t0)= 1, η(x(t0), t0)= 1, w(x(t0), t0)=1 (конкретный вид этих функций зависит от ДУ). Задача нахождения оптимальной траектории КА с ЭРД прямого перелета, сформулированная в наиболее общем виде, может содержать следующие краевые условия в начальный и конечный момент времени: g0(x(t0), v(t0), t0)= 0, m(t0)= m0, k tk = t0 + D.t, gk (x(tk ), v(tk ), tk )= 0, m(tk ) ?: m∗ , (2.6) ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 103 где t0 и tk - начальный и конечный моменты времени, (x, v, m) - положение, скорость и масса КА, ∗ mk - минимальная конечная масса КА. Эти моменты времени могут быть заданы или выбираться оптимальным образом. Функции g0 и gk , ограничивающие положение и скорость аппарата в начальный и конечный момент, могут включать параметры орбит и положение на этих орбитах. Основным критерием оптимизации в задачах механики космического полета является конечная масса КА m(tk ) → max; (2.7) помимо этого, рассматривается и классическая для теории оптимального управления задача оптимального быстродействия D.t → min . (2.8) 3. ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ КОСМИЧЕСКИМ АППАРАТОМ Идеально-регулируемый (ИР) двигатель ограниченной мощности (ОМ-задача) является математической моделью электроракетного двигателя, в рамках которой предполагается, что на скорость истечения и тягу ЭРД накладывается единственное ограничение - величина реактивной мощности ЭРД (2.4). В рамках указанного ограничения скорость истечения и тяга ИР-двигателя могут произвольно изменяться. Оптимизация траекторий КА с ИР-двигателем рассматривалась во многих работах (например, [9, 17, 18, 28]), привлекая к себе внимание относительной простотой оптимального управления. Самостоятельно эта задача может иметь мало приложений на практике ввиду сложности реализации регулирования двигателя в большом диапазоне по тяге и скорости истечения. Тем не менее, ОМ-задача позволяет определить максимально возможную величину конечной массы КА при заданной начальной мощности ДУ, а также дает достаточно хорошее начальное приближение для постановок задач с более реальными моделями работы ЭРД. Учитывая зависимость (2.4), дифференциальное уравнение для массы КА из системы (2.1) можно переписать в следующем виде: dm T 2 dt = - 2N m2a2 - = 2N0χ(x, t) , (3.1) где a = T /m - вектор реактивного ускорения. Это уравнение имеет следующее решение (см. [9, 17, 28]): m0 где m(t)= m0 1+ N0 t J (t) , (3.2) 1 r J (t)= 2 t0 a2(τ ) χ(x(τ ),τ ) dτ, (3.3) и минимум J (t) соответствует максимальной конечной массе. Возможность записать отдельно выражение для массы (3.2) и отсутствие зависимости правых частей дифференциальных уравнений для положения и скорости от мощности приводит к разделению задачи на две независимые части (см. [9, 17, 28]): · динамическую - нахождение оптимальной программы ускорения от реактивной тяги, минимизирующей интегральный критерий (3.3) для заданного динамического маневра; · параметрическую - определение оптимальных массовых соотношений и уровня мощности, проверка удовлетворения условия на конечную массу (2.6). Такое разделение задачи существенно упрощает ее решение и оправдывает идеализацию работы ЭРД. Далее будет сформулирована динамическая часть ОМ-задачи как задача оптимального управления. Именно ее решение будет исследоваться в дальнейшем, и далее в тексте она будет именоваться просто ОМ-задачей, а весовая часть как таковая рассматриваться не будет. 104 А. В. ИВАНЮХИН 1. ОМ-задача (динамическая часть). Среди всех управлений (a ∈ R3), переводящих КА с ИР-двигателем с начального многообразия на конечное за заданное время, требуется найти управление, позволяющее минимизировать функционал (3.3) и соответствующую этому управлению траекторию перелета. Диапазон выбора дат старта ограничен. ⎧ tk r ⎪ 1 ⎪ J (tk )= ⎪ ⎪⎪ 2 t ⎨ 0 a2 χ dt → min, dx ⎪ dt ⎪ = v, dv dt = Ωx + a, (3.4) ⎪ g (x(t ), v(t ),t )= 0, g (x(t ), v(t ),t )= 0, 0 0 0 0 ⎪ k k k k ⎩⎪ tk - t0 = D.t, tmin :( t0 :( tmax. 0 0 Прежде всего необходимо отметить, что в этой задаче всегда существует оптимальное управление в классе измеримых управлений a ∈ L2[t0, tk ] и непрерывных траекторий (x, v). Такой случай рассматривался еще в вариационном исчислении, где была установлена теорема Тонелли и ее различные варианты (см. [4, 14]). Так как система линейна по управлению, а функционал квадратичен, следовательно, из слабой сходимости по управлению следует равномерная сходимость (x, v). Относительно функции χ без нарушения общности справедливо следующее: χ(x, t) ?: Cχ > 0. Тогда коэффициент 1 χ ограничен сверху, и функционал ОМ-задачи (3.4) коэрцитивен (т. е. стремится к +∞ при ⊕a⊕→ ∞), поэтому всякая минимизирующая последовательность ограничена. Кроме того, в динамической части ОМ-задачи всегда существует и пробное управление (удовлетворяющее краевым условиям), так как всегда можно провести кривую, соединяющую начальное и конечное положение, и по этой кривой восстановить управление с помощью уравнений (3.4). А учитывая вид зависимости конечной массы от функционала динамической части ОМ-задачи (3.2), из которой видно, что она не может быть меньше нуля, можно сделать заключение, что и оптимальное решение параметрической части ОМ-задачи всегда существует при отсутствии ограничений на конечную массу. Рассмотрим задачу в рамках принципа максимума (см. [19]). Функция Понтрягина имеет вид a2 H = -λJ 2χ + (pv , Ωx) + (pv , a) + (px, v) , (3.5) а терминант (краевая функция Лагранжа) gk l = (g0(x(t0), v(t0), t0), λ0) + / (x(tk ), v(tk ), tk ), λk \\ + (3.6) + (tk - t0 - D.t)λt + (tmin - t0)λ1 + (t0 - tmax)λ2 , 0 t0 0 t0 t0 где px, pv - липшицевы функции (сопряженные переменные к x и v соответственно), λJ , λt, λ1 , λ2 0 k t0 - скалярные множители, λ , λ - векторы. Условия принципа максимума принимают следующий вид: § условия оптимальности по управлению a = argmaxH(t, x(t), v(t), a!): a∗∈R3 § уравнения сопряженной системы: 1 λ a = χpv , (3.7) J 2 v dpx = -Ω p ∂χ p - , dt dpv dt xx v = -px, 2λJ ∂x (3.8) ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 105 § условия трансверсальности: ⎧ ∂ 0) ∂ k ) ⎪px(t0)= ⎪ ⎪ ∂x(t0) ∂ (g0(x(t0), v(t0), t0),λ , px(tk )= - ∂x(tk ) ∂ (gk (x(tk ), v(tk ), tk ),λ , ⎨pv (t0)= ⎪ ∂v(t0) (g0(x(t0), v(t0), t0), λ0) , pv (tk )= - ∂v(tk ) (g0(x(tk ), v(tk ), tk ), λk ) , ⎪ ⎪⎩H(t0)= - ∂ ∂t0 t0 (g0(x(t0), v(t0), t0), λ0) + λt + λ1 t0 - λ2 , H(tk )= ∂ ∂tk (gk (x(tk ), v(tk ), tk ), λk ) + λt, (3.9) § условия дополняющей нежесткости: (tmin 1 max 2 0 - t0)λt0 = 0, (t0 - t0 )λt0 = 0, (3.10) § условия нетривиальности и неотрицательности: t0 + λ t0 |λ0| + |λk | + |λt| + λ1 2 t0 + λJ > 0, λJ ?: 0, λ1 t0 ?: 0, λ2 ?: 0. (3.11) В связи с однородностью по сопряженным переменным функции Понтрягина можно ввести одно условие нормировки. В данном случае удобно ввести ее следующим образом, не зависящим от вида краевых условий g0, gk и остальных множителей: λJ = 1. (3.12) Случай λJ =0 противоречит приведенным условиям оптимальности. Так как максимум функции Понтрягина достигается при условии ∂H ∂a = pv - λJ a = 0, χ что при λJ = 0 приводит к pv (t) ≡ 0, тогда из (3.8) следует, что px(t) ≡ 0, а из (3.5) вытекает, что H(t) ≡ 0. То есть, pv (t0) = pv (tk ) = 0, px(t0) = px(tk ) = 0, H(t0) = H(tk ) = 0, и условия трансверсальности (3.9) для невырожденных краевых условий могут быть выполнены только при t0 тривиальном наборе множителей λJ = 0, λ0 = 0, λk = 0, λt = 0, λ1 t0 = 0, λ2 = 0, что противоречит условию (3.11). Значит, λJ > 0 и справедлива нормировка (3.12). 2. ОТ-задача. Как отмечалось ранее, ОМ-задача использует идеализированную математическую модель ЭРДУ. Задача о перелете КА с ЭРДУ ограниченной тяги (ОТ-задача) существенно реалистичнее, поэтому оптимизация ОТ-траекторий представляет значительный практический интерес. Однако с методической точки зрения, оптимизация ОТ-траекторий осложняется невозможностью разделения задачи на динамическую и весовую части и более жесткими ограничениями на управление (2.5). Сформулируем ОТ-задачу следующим образом. Среди всех управлений (e, δ) ∈ R3 × R, переводящих КА с начального многообразия на конечное за заданное время, требуется найти управление, позволяющее максимизировать массу КА в конечный момент времени, и соответствующую ему траекторию перелета. Диапазон выбора дат старта ограничен. ⎧ m(tk ⎪ dx dv η(x, t)T0 dm η(x, t)T 0 δ, = v, = Ωx + eδ dt dt m , = dt - w(x, t)c 0 ⎪ ⎪ ⎪ ⎨ ) → max, g0(x(t0), v(t0), t0)= 0, m(t0)= m0, gk (x(tk ), v(tk ), tk )= 0, ⎪ ⎪ tk - t0 = D.t, tmin :( t0 :( tmax, (3.13) ⎪ 0 0 ⎩⎪ |e| = 1, δ ∈ [0; 1]. где δ - функция тяги, e - единичный вектор в направлении вектора тяги. Рассмотрим эту задачу в рамках принципа максимума (см. [19]). Функция Понтрягина имеет вид (pv ,e) pm \\ H = (px, v) + (pv , Ωx) + δη(x, t)T0 m - w(x, t)c0 , (3.14) 106 А. В. ИВАНЮХИН а терминант (краевая функция Лагранжа) l = (g0(x(t0), v(t0), t0), λ0) + /gk (x(tk ), v(tk ), tk ), λk \\ + (3.15) + (m(t0) - m0)λ0 - m(tk )λk + (tk - t0 - D.t)λt + (tmin - t0)λ1 + (t0 - tmax)λ2 , m m 0 t0 0 t0 где px, pv , pm - липшицевы функции (сопряженные переменные), λ0 , λk , λt, λ1 , λ2 · скалярные m m множители, λ0, λk - векторы. Условия принципа максимума принимают следующий вид: t0 t0 o условия оптимальности по управлению (e, δ) = argmax (e∗,δ∗)∈R3×R ⎧ H(t, x(t), v(t), m(t), e!, δ!): 1, Ψ > 0, e = pv |pv | , Ψ= |pv | m pm · wc , δ = ⎨⎪ ∀δ! ∈ [0; 1], Ψ= 0, ⎩⎪ 0, Ψ < 0, (3.16) o уравнения сопряженной системы: ⎧ dpx ⎪ = -Ωxxpv - Г (pv ,e) - pm \\ ∂η - η pm ∂w δT0, ⎪ dt ⎪⎨ dpv dt ⎪ ⎪ = -px, m wc0 ∂x w2c0 ∂x (3.17) ⎪ dpm (pv ,e) ⎩ dt = m2 ηT0δ, o условия трансверсальности: ⎧ ∂ 0) ∂ k ) ⎪px(t0)= ∂x(t ) (g0(x(t0), v(t0), t0),λ , px(tk )= - ∂x(t (gk (x(tk ), v(tk ), tk ),λ , ) ⎪ 0 ⎪ ⎪ ∂ 0) k ∂ k ) ⎨⎪pv (t0)= ∂v(t ) (g0(x(t0), v(t0), t0),λ , pv (tk )= - ∂v(t (gk (x(tk ), v(tk ), tk ),λ , ) 0 k ⎪pm(t0)= λ0 , pm(tk )= λk , m m ⎪ ⎪ ∂ ⎪ 0) 1 2 ∂ ( k ) ⎪⎩H(t0)= - ∂t (g0(x(t0), v(t0), t0),λ + λt + λt0 - λt0, H(tk )= ∂t gk (x(tk ), v(tk ), tk ),λ + λt, 0 o условия дополняющей нежесткости: (tmin 1 k max 2 (3.18) 0 - t0)λt0 = 0, (t0 - t0 )λt0 = 0, (3.19) o условия нетривиальности и неотрицательности: |λ0 | + |λ0| + |λk | + |λt| + λ1 + λ2 + λk > 0, λk ?: 0, λ1 ?: 0, λ2 ?: 0. (3.20) m t0 t0 m m t0 t0 В связи с однородностью по сопряженным переменным функции Понтрягина, так же как и ранее, можно ввести одно условие нормировки, связывающее сопряженные переменные. Как видно, оптимальное управление (3.16) может иметь особенности следующего вида: 1. pv (t) ≡ 0 на некотором конечном интервале [τ0, τk ] (t0 :( τ0 < τk :( tk ), тогда из (3.17) следует, что и px(t) ≡ 0, pm(t) ≡ 0 ∀t ∈ [τ0, τk ], и значит, на всем интервале t ∈ [t0, tk ] имеет место px(t) ≡ 0, pv (t) ≡ 0, pm(t) ≡ 0, а тогда условия трансверсальности (3.18) при невырожденных m краевых условиях допускают только тривиальное решение λ0 = 0, λ0 = 0, λk = 0, λt = 0, λ1 2 k t0 = 0, λt0 = 0, λm = 0, что противоречит (3.20), а значит, оптимальная траектория не может содержать рассматриваемых участков; 2. Ψ(t) ≡ 0 на некотором конечном интервале [τ0, τk ] (t0 :( τ0 < τk :( tk ), тогда принцип максимума допускает на этом интервале особое управление соответствующее промежуточной тяге δ ∈ [0, 1], но не позволяет однозначно определить его; для определения оптимального промежуточного управления необходимо воспользоваться условием Коппа-Мойера (см. [29]), которое для центрального ньютоновского гравитационного поля имеет вид равенств (см. [7, 8]) Ψ(t)= 0, Ψ˙ (t)= 0, Ψ¨ (t)= 0, ... t)= 0, ....(t)=0 t [τ ,τ ], Ψ( Ψ ∀ ∈ 0 k ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 107 последнее из которых позволяет выразить управление с помощью условия ∂ d4 ∂H ∂δ dt4 ∂δ :( 0 (однако в данной работе особые управления не рассматриваются). m Рассмотрим подробнее случай λk = 0, так как он приводит к тому, что pm(tk )= 0, а с учетом уравнения для pm из (3.17) это значит, что pm(t) :( 0 ∀t ∈ [t0, tk ]. Следовательно, ∀t ∈ [t0, tk ] Ψ ?: 0 и при отсутствии участков особого управления δ ≡ 1, что соответствует задаче оптимального быстродействия, и заданное время перелета D.t является минимальным для рассматриваемой задачи. Кроме того, в ОТ-задаче (3.13) вопрос существования решений не так тривиален, как в случае ОМ-задачи (3.4). Однако в работах Oberle и Taubert [32], Caillau, Gergaud, Noailles, Haberkorn [25, 27] на основе теоремы Филиппова [21] делается вывод о существовании оптимального решения в ОТ-задаче. Однако во всех перечисленных работах остается недостаточно исследованным вопрос существования пробного управления и четкого выделения задач, в которых оптимальное управление заведомо существует на множестве основных параметров, таких как время перелета, тяга и скорость истечения (D.t, T0, c0). Для подтверждения существования допустимых управлений, вообще говоря, можно указать только два пути: построение области достижимости и определение области существования решения на множестве параметров системы. Оба эти подхода применимы к задачам механики полета с малой тягой. Однако построить область достижимости - значит определить все возможные положения, в которые можно попасть из начальной точки за заданное время, управляя системой допустимым образом. Очевидно, что определение такой области фазового пространства для задач механики полета с малой тягой является очень трудоемким и сложным процессом, в основном из-за «проклятия размерности» и нелинейности уравнений (2.1). Из-за этих недостатков построение множеств достижимости в рассматриваемых задачах на практике может носить только весьма ограниченный характер. Для оценки проектных параметров КА с двигателем малой тяги в ОТ-задаче и существования допустимого управления для заданного перелета часто используются решения задачи с идеально регулируемым двигателем (ОМ-задача, см. [18]), оптимального быстродействия (см. [25, 32]) и импульсные траектории (см. [5, 6]). При этом часто (см. [5, 6, 18]) эти решения используются, в первую очередь, в качестве начального приближения, и вопрос существования допустимого управления решается в них только косвенно. При этом решения с импульсной тягой не являются для ОТ-задачи допустимыми управлениями, а решения ОМ-задачи будут таковыми, только если изменение реактивного ускорения на траектории ОМ-задачи не превышает допустимых значений ОТ-задачи, т. е., если траекторию ОМ-задачи можно реализовать с моделью двигателя ограниченной тяги, что выполняется далеко не всегда. В этих условиях нельзя гарантировать успешность перехода от задачи перелета с идеально-регулируемой или импульсной тягой к задаче с ограниченной тягой. Такой проблемы не возникает при использовании задачи оптимального быстродействия, в работах [27, 32] она рассматривается в качестве гарантии существования допустимого управления - если есть решение оптимального быстродействия, то существует и пробное управление с любым временем перелета, большим минимального [25]. Однако использование задачи оптимального быстродействия в качестве проверки осуществимости перелета имеет один существенный недостаток - в процессе ее решения неизбежно будет происходить изменение угловой дальности перелета. А как отмечалось ранее, существует деление семейств экстремалей, отличающихся числом целых витков вокруг притягивающего центра. И изменение угловой дальности в процессе решения может приводить к неконтролируемому переходу от одного семейства экстремалей к другому, что осложняет анализ задачи. Чтобы гарантировать в процессе решения нахождение в рамках выбранного семейства экстремалей, необходимо зафиксировать угловую дальность или, по крайней мере, минимизировать и отслеживать ее изменения. Более конструктивным и оправданным представляется метод, основанный на том, что решение ОТ-задачи может существовать не всегда, а точнее, не для всех значений тяги и скорости истечения. Область существования решения ограничена минимальными значениями тяги и скорости 108 А. В. ИВАНЮХИН истечения, так как, во-первых, для осуществления заданного перелета требуется некоторое конечное приращение характеристической скорости, для реализации которого за определенное время D.t требуется достаточная величина реактивного ускорения, а следовательно - тяги. Во-вторых, из условий (2.6) следует, что допустимый расход рабочего тела ограничен сверху, следовательно, скорость истечения также имеет предельное минимальное значение. Более того, для каждого значения скорости истечения c ∈ [cmin; ∞] существует минимальное значение тяги Tmin и реактивной мощности Nmin, при котором существует решение: Tmin(c), Nmin(Tmin(c), c), c ∈ [cmin; ∞]. (3.21) Минимальное значение скорости истечения cmin определяется выходом на предельное значение конечной массы КА при бесконечной величине тяги (импульсное решение): v min - c = D. ln(mk /m0) , (3.22) где D.v - суммарные затраты приращения скорости для импульсного решения. Тогда можно привести следующие рассуждения относительно существования решений ОТmax задачи. Предположим, есть решение ОМ-задачи, и известно значение максимальной тяги на этом решении T OM , тогда всегда есть допустимое управление в задаче с минимальной тягой без массового расхода: max T ∈ r0; T OM l , c = ∞, m(tk )= m0. (3.23) Решив задачу на минимум тяги без массового расхода, можно продолжить это решение по скорости истечения, получив семейство решений, являющееся нижней границей области существования ОТ-задачи (3.21). После этого можно сделать заключение, что для каждого значения тяги больше минимального при заданной скорости истечения в ОТ-задаче максимизации конечной массы решение с минимальной тягой является допустимым (но не оптимальным) управлением, и из теоремы существования следует наличие оптимального управления этой задачи с тягой больше минимальной T ∗ ?: Tmin(c), T ∈ [0; T ∗]. (3.24) Отсутствие решений с тягой меньше минимальной следует из постановки ОТ-задачи с минимальной тягой. Увеличивать тягу можно до еще одного предельного случая - импульсного решения. Всякое решение с конечной тягой может рассматриваться как промежуточное в последовательности решений краевых задач с неограниченно возрастающей тягой, сходящихся к импульсному решению, как это показано в работе [31]. Таким образом, совокупность значений тяги и скорости истечения при заданном времени перелета ограничивает снизу кривая минимальной тяги Tmin(c). Метод построения этой кривой и формулировка задачи о минимуме тяги приведены далее в разделах 5 и 6. 1. МЕТОД ПРОДОЛЖЕНИЯ ПО ПАРАМЕТРУ И КРАЕВАЯ ЗАДАЧА ПРИНЦИПА МАКСИМУМА Приведенные выше соотношения (3.4)-(3.12) и (3.13)-(3.20) сводят ОМи ОТ-задачи к краевым, в которых необходимо определить недостающие значения переменных, удовлетворяющие краевым условиям и соответствующим условиям оптимальности. При этом, если начальное и конечное многообразие не имеют особых точек, то число неизвестных краевой задачи всегда совпадает с числом условий (см. [19]). В данной работе все краевые задачи принципа максимума будут решаться методом продолжения по параметру (гомотопии). Кроме того, с его помощью будет строиться граница области существования решений ОТ-задачи. В связи с этим приведем краткое описание метода. Идея метода продолжения прослеживается в литературе еще с работ Шварца (1869 г.), Леверье (1886 г.) и Пуанкаре (1892 г.). Тем не менее, для численного решения задач этим методом стали пользоваться только во второй половине прошлого века в связи с расширением возможностей ЭВМ. Для случая решения систем нелинейных уравнений в непрерывном виде этот метод был предложен Д. Ф. Давиденко [10] и М. К. Гавуриным [3], а в дальнейшем развит во многих работах (см. [13, 16, 23, 24]), в том числе и для решения краевых задач для систем обыкновенных ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 109 дифференциальных уравнений. В частности, последние 20 лет метод продолжения активно применяется в задачах механики космического полета с малой тягой (см., например, [17, 18, 25, 27]). Применение метода продолжения по параметру к краевой задаче для обыкновенных дифференциальных уравнений заключается в следующем. Будем решать нелинейную краевую задачу ⎧ dx dt ⎪ = f (x, τ, t), t ∈ [t0; tk ], ⎨⎪ g(τ, x(t0), x(tk )) = 0, ⎪ (x, τ, t) ∈ Rn × R × R, ⎪ (4.1) ⎩ g : R × Rn × Rn → Rn × Rn, с вектором неизвестных начальных значений p (p ∈ Rm, m :( n), т. е. p - часть вектора x(t0), подлежащая определению для удовлетворения краевых условий. Эту краевую задачу нужно представить в виде некоторой нелинейной системы Φ(p, τ )= 0. (4.2) При этом, если известно решение краевой задачи (4.1) при каком-либо значении параметра τ (т. е. i xi(t0)= p0, i = 1,..., m), то функция (4.2) совпадает с краевыми условиями задачи Φ(p, τ )= g(τ, x(t0), x(tk )) = 0. (4.3) i Если такое решение не известно, и начальное приближение xi(t0) = p0, i = 1,..., m, не удовлетворяет краевым условиям, функция (4.2) может быть представлена в виде Φ(p, τ )= g(τ, x(t0), x(tk )) - τb = 0, b = g(τ, x(t0), x(tk )), τ ∈ [1; 0], (4.4) где p0 - вектор начального приближения неизвестных значений, b - вектор начальных невязок, параметр τ и уравнения краевой задачи (4.1) нормированы и преобразованы соответствующим образом, чтобы обеспечить выполнение принадлежности параметра заданному отрезку (например, τ ∈ [1; 0]). Продифференцировав соотношение (4.2), легко получить производную решения от параметра ∂Φ dp ∂Φ + ∂p dτ ∂τ = 0. (4.5) Это выражение в первом случае функции Φ(p, τ ) (4.3) можно переписать как dp Г ∂g -1 ∂g = - , (4.6) а в случае (4.4) dτ ∂p ∂τ dp Г ∂g -1 ∂g = - \\ - b . (4.7) dτ ∂p ∂τ Таким образом, процесс построения неявной функции, заданной уравнением Φ(p, τ ) = 0, на интересующем отрезке изменения параметра сводится к интегрированию системы (4.5). При этом должны быть выполнены локальные условия существования неявной функции на всем отрезке продолжения (см. [13, 16, 22, 23]) rank Г ∂g , ∂p ∂g ∂τ = m, (4.8) и задача (4.1) должна иметь непрерывные производные по начальным условиям и параметру, что можно обеспечить наложением ограничения на гладкость правых частей и краевых условий (см. [13, 24]): f (x, τ, t) ∈ C1, g(τ, x(t0), x(tk )) ∈ C1. (4.9) Эти условия необходимы для локальной продолжаемости по параметру. Что касается нелокальной продолжаемости на весь участок изменения параметра, этот вопрос более сложный и никакими условиями не может быть гарантирован. Как отмечается в ряде работ [10, 23, 24], во многих случаях рационально переходить к продолжению по длине дуги или менять параметр в процессе продолжения (чередуя по необходимости τ с компонентами вектора p) для организации немонотонного изменения параметра τ вдоль линии 110 А. В. ИВАНЮХИН продолжения. Такой метод способен проходить точки разворота (предельные точки), характеризующиеся вырождением матрицы чувствительности по начальным значениям на единицу, при сохранении полного ранга расширенной матрицы rank Г ∂g ∂p = m - 1, rank Г ∂g , ∂p ∂g ∂τ = m. (4.10) Для построения кривой продолжения совершенно неважно, какая именно переменная выступает в роли параметра продолжения: им может быть и их некоторая комбинация. Производная от (4.2) по абстрактному параметру продолжения s будет иметь вид (предполагается, что Φ(p(s),τ (s)). dΦ ∂Φ dp = + ds ∂p ds ∂Φ ∂τ ∂τ ∂s (4.11) Определить s, вообще говоря, можно любым удобным способом, присоединив к системе (4.11) соответствующее условие. Часто в качестве такого параметра используется длина дуги (см. [23, 24]), тогда определяющее его уравнение имеет вид (ds)2 = (dp)2 + (dτ )2. (4.12) Недостатком такого определения параметра продолжения является нелинейность уравнения (4.12). Можно указать еще один достаточно оригинальный способ выбора параметра продолжения (в некоторых источниках именуемый «методом Драгилева»). Решение уравнений (4.5) по методу Крамера имеет вид dpi = - Δi dτ , i = 1,..., m, (4.13) ds Δ0 ds ∂Φ где Δ0 - определитель матрицы , Δi - определитель матрицы, в которой i-ый столбец матрицы ∂p ∂Φ заменен на ∂Φ . Теперь можно ввести параметр s, приняв производную от τ по s равной ∂p ∂τ определителю -Δ0. Тогда система уравнений примет следующий вид: ⎧ dpi ⎪⎨ ds dτ ⎪ = Δi, (4.14) ⎩ ds = -Δ0. Параметры, заданные таким образом (длина дуги, «метод Драгилева»), как уже отмечалось ранее, допускают немонотонное изменение τ и позволяют проходить точки разворота кривой продолжения, что оказывается полезным в некоторых случаях. Для повышения точности процесс продолжения можно организовать по схеме предиктор- корректор, такой вариант рассматривается, например, в книге [23]. Поскольку известен вид неявной функции (4.2), это выражение должно выполняться при каждом значении параметра τ. Тогда в качестве предиктора можно использовать выполнение одного или нескольких шагов выбранного для решения системы (4.6) или (4.7) метода численного интегрирования, а в качестве корректора, с целью уточнения текущего значения p, - численный метод для решения системы нелинейных уравнений (4.2). Эту систему можно решать, например, методом Ньютона, тогда итерации корректора имеют следующий вид: пока |Φ(pj ,τ )| ?: h, pj = pj-1 - ∂Φ ∂p p=p j-1 \\-1 Φ(pj-1,τ ), (4.15) где h - точность процесса продолжения, j - номер итерации метода Ньютона. Вычисление матрицы чувствительности ∂Φ/∂p и вектора ∂Φ/∂τ может производиться разными способами, например, с использованием конечных разностей, совместного интегрирования дифференциальных уравнений для элементов этой матрицы с уравнениями краевой задачи, либо с помощью метода комплексного шага для высокоточного численного дифференцирования (см. [30, 33]). ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 111 Как видно, в ОТ-задаче правые части дифференциальных уравнений содержат ступенчатую (релейную) функцию тяги δ, которая приводит к их разрывам, что противоречит условиям гладкости метода продолжения (4.9). И чтобы иметь возможность воспользоваться им, необходимо использовать сглаженные правые части, вернее, сглаженное управление, так как разрывность обусловлена только включением-выключением двигателя. Можно использовать следующую сглаженную функцию тяги: 1 δ = 1+ 2 Ψ \\ |Ψ| + ε , (4.16) где ε - параметр сглаживания, Ψ - функция переключения. Используя этот подход, всегда можно получить функцию сколь угодно близкую к релейной, устремляя параметр сглаживания к нулю на столько, на сколько это позволит устойчивость метода продолжения. Как отмечалось ранее, одной из особенностей задач механики перелета является возможность существования множества локальных экстремумов рассматриваемой задачи оптимального управления при заданном времени перелета, соответствующих различному количеству целых витков вокруг гравитирующего центра. Эти локальные экстремумы образуют изолированные семейства решений (семейства экстремалей). Для регуляризации процесса решения задачи оптимального управления и более полного анализа семейств экстремалей необходим алгоритм, обеспечивающий при заданном начальном приближении сходимость к решению с заданной угловой дальностью перелета (числом витков). Из полученного ранее оптимального управления видно, что для рассматриваемых задач нулевое начальное приближение для начальных значений сопряженных переменных соответствует нулевому реактивному ускорению на всей траектории перелета, и в конечный момент времени КА останется на начальной орбите. Для того, чтобы использовать пассивное движение аппарата на всем участке траектории, кроме начальной и конечной точки (т. е. нулевое начальное приближение для вектора сопряженных переменных), можно воспользоваться двумя видами параметризации: § модификацией гравитационного параметра - изменением параметра притягивающего центра таким образом, чтобы угловая дальность пассивного движения по начальной орбите совпадала с заданной величиной с точностью до 2π (см. [17, 18]); § изменением начальной и конечной скорости КА - заменой их на скорости, полученные из решения задачи Ламберта. Первый прием легко проиллюстрировать следующими рассуждениями. Предположим, известны начальные приближения положения и скорости КА в начальный и конечный момент времени (x0, v0) и (xk , vk ), а также время перелета. Тогда можно найти угол между этими двумя положениями ∠(x0, xk ) и, зная параметры начальной орбиты, перейти от него к средней аномалии, которая связана с гравитационным параметром следующим соотношением: μ D.M + 2πn = a 3 D. t, (4.17) 0 где D.M - угол средней аномалии между начальным и конечным положением, D.M ∈ [0; 2π], n - число целых витков, a0 - большая полуось начальной орбиты, μ - гравитационный параметр, D.t - время перелета. Из этой формулы видно, что, задавая время перелета и число целых витков, можно всегда подобрать подходящее значение гравитационного параметра. Другой вариант априорного выбора решения, соответствующего определенному числу витков вокруг центрального тела - это использование решения задачи Ламберта. Задача Ламберта относится к классическим задачам небесной механики и заключается в определении кеплеровской траектории движения по двум известным векторам положения и заданному интервалу времени, за которое КА перемещается между этими положениями (после определения орбиты находятся и векторы скорости в граничных положениях). Эта задача хорошо описана в литературе, например, в [1, 20] подробно приведены эффективные методы ее решения. Следует отметить, что при достаточно большом времени перелета задача Ламберта может быть решена для разного числа витков вокруг центрального тела, однако, в отличие от модификации гравитационного параметра, она имеет по два решения для каждого числа витков, большего нуля. Кроме того, задача Ламберта может не иметь решения при достаточно большом числе витков и 112 А. В. ИВАНЮХИН достаточно малом времени перелета, что связано с геометрическими ограничениями кеплеровского движения. Метод модификации гравитационной постоянной избавлен от этого недостатка, так как использует меньше соотношений кеплеровского движения, но по этой причине он не определяет точное значение векторов скорости и не удовлетворяет условиям по положению на правом конце траектории. Преимуществом использования задачи Ламберта является то, что с ее помощью можно удовлетворить все краевые условия по положению и скорости, используя в качестве начального приближения двухимпульсный перелет, полученный на ее основе. 2. ЗАДАЧА МИНИМИЗАЦИИ ТЯГИ Можно сформулировать задачу минимизации тяги следующим образом, аналогично работе [12], с добавлением ограничения на величину конечной массы. Задача на минимум тяги. В рамках модели ограниченной тяги (ОТ-задачи) найти управление (e, δ), переводящее КА с начального многообразия на конечное за заданное время с минимально возможной величиной начальной реактивной тяги: ⎪ ⎧ T0 → min, ⎪ dx dv ⎪ η(x, t)T0 dm η(x, t)T0 δ, - = v, ⎨⎪ dt = Ωx + dt eδ, m = dt w(x, t)c0 k g0(x(t0), v(t0), t0)= 0, m(t0)= m0, gk (x(tk ), v(tk ), tk )= 0, m(tk ) ?: m∗ , ⎪ (5.1) ⎪ t - t = D.t, tmin :( t0 :( tmax, ⎪ k 0 0 0 ⎪ ⎩ (e, δ) ∈ R3 × R, |e| = 1, δ ∈ [0; 1], где δ - функция тяги, e - единичный вектор в направлении вектора тяги. В соответствии с идеями принципа максимума введем дополнительное дифференциальное уравнение запишем функцию Понтрягина dT0 = 0, (5.2) dt (pv ,e) pm \\ H = (px, v) + (pv , Ωx) + δη(x, t)T0 и терминант (краевую функцию Лагранжа) m - w(x, t)c0 , (5.3) l = T0λT + (g0(x(t0), v(t0), t0), λ0) + /gk (x(tk ), v(tk ), tk ), λk \\ + (5.4) + (m(t0) - m0)λ0 + (m∗ - m(tk ))λk + (tk - t0 - D.t)λt + (tmin - t0)λ1 + (t0 - tmax)λ2 , m k m 0 t0 0 t0 где px, pv , pm - липшицевы функции (сопряженные переменные), λ0 , λk , λt, λ1 , λ2 · скалярные множители и λ0, λk - векторы. Условия принципа максимума принимают следующий вид: o условия оптимальности по управлению (e, δ) = argmax (e∗,δ∗)∈R3×R ⎧ m m t0 t0 H(t, x(t), v(t), m(t), e!, δ!): 1, Ψ > 0, e = pv |pv | , Ψ= |pv | m pm · wc , δ = ⎨⎪ ∀δ! ∈ [0; 1], Ψ= 0, ⎩⎪ 0, Ψ < 0, (5.5) o уравнения сопряженной системы: ⎧ dpx T0 ∂η T0 ∂η η ∂w\\ ⎪ ⎪ dt ⎪ ⎪ dp = -Ωxxpv · δ m |pv | ∂x + δ pm wc0 ∂x - , w ∂x ⎪ v ⎪⎨ dt = -px, m dp ηT0 ⎪ = δ p |, (5.6) ⎪ dt ⎪ ⎪ m2 | v ⎪ dpT - ⎩ ⎪ = δηΨ, dt ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 113 o условия трансверсальности: ⎧ ∂ ∂ ⎪px(t0)= ⎪ ⎪ ∂x(t0) (g0(x(t0), v(t0), t0), λ0) , px(tk )= - ∂x(tk ) (gk (x(tk ), v(tk ), tk ), λk ) , ( 0) ⎪ ∂ ⎪ ∂ ( k ) pv (t0)= ⎪ ∂v(t ) g0(x(t0), v(t0), t0),λ , pv (tk )= - ∂v(t ) g0(x(tk ), v(tk ), tk ),λ , ⎪ 0 k ⎨pm(t0)= λ0 , pm(tk )= λk , ⎪ p T ⎪ ⎪ (t0)= λ ⎪ ⎪ ⎪ m T , pT ∂ (tk m )= 0, ∂ ⎪⎩H(t0)= - ∂t0 t0 (g0(x(t0), v(t0), t0), λ0) + λt + λ1 t0 - λ2 , H(tk )= ∂tk (gk (x(tk ), v(tk ), tk ), λk ) + λt, (5.7) o условия дополняющей нежесткости (tmin 1 max 2 ∗ k 0 - t0)λt0 = 0, (t0 - t0 )λt0 = 0, (mk - m(tk ))λm = 0, (5.8) o условия нетривиальности и неотрицательности |λ0 | + |λ0| + |λk | + |λt| + λ1 + λ2 + λk + λT > 0, λT ?: 0, λk ?: 0. (5.9) m t0 t0 m m Как и ранее, данная задача допускает условие нормировки, связывающее сопряженные переменные, и может иметь участки особого управления, аналогичные рассмотренным ранее. Из условий трансверсальности (5.7) и дополняющей нежесткости (5.8) видно, что задача распадается на два случая: 1. конечная масса больше минимально допустимой: m(tk ) > m∗ , λk = 0, (5.10) k m o в этом случае из-за неотрицательности правой части дифференциального уравнения для сопряженной переменной к массе (5.6) и ее равенства нулю в конечной точке следует, что на всей траектории она неположительна, а функция переключения всегда неотрицательна, значит, двигатель на траектории всегда включен: ∀t ∈ [t0; tk ] pm(t) :( 0, Ψ > 0 ⇒ δ ≡ 1; (5.11) 2. конечная масса равна минимально допустимой: m(tk )= m∗ , λk ?: 0, (5.12) k m o при этом оптимальное управление определяется по общему правилу (5.5), и траектория может содержать участки нулевой тяги. У этих двух вариантов есть одна общая точка m(tk )= m∗ , λk = 0. k m Можно проинтегрировать уравнение для сопряженной переменной к тяге и записать условие для минимума тяги: tk r δηΨdt = λT ?: 0. (5.13) t0 Поскольку под интегралом стоит неотрицательная функция, то λT > 0, если на траектории есть хотя бы один активный участок. Можно показать, что условия оптимальности задач для минимума времени перелета и минимума тяги связаны через условия нормировки, т. е. они приводят к одинаковым решениям, у которых по-разному нормированы значения сопряженных переменных. Тогда справедливо равенство pT (t0)= Gpt0 (t0), (5.14) где pT (t0) и pt0 (t0) - векторы начальных значений сопряженных параметров для задачи на минимум тяги и оптимального быстродействия, G - некоторая неотрицательная константа. Справедливость (5.14) можно легко проверить простым сравнением соответствующих решений. 114 А. В. ИВАНЮХИН 1. ВЫЧИСЛЕНИЕ ГРАНИЦЫ ОБЛАСТИ СУЩЕСТВОВАНИЯ Построение решений задачи на минимум тяги на интервале c ∈ [cmin; ∞] удобно начинать со случая бесконечной скорости истечения (нулевого массового расхода). Этому случаю соответствует перелет с постоянной массой КА m(t) = m0, и, следовательно, соответствующее уравнение может быть исключено, а задача минимизации начальной тяги переписана относительно начального реактивного ускорения в следующем виде: ⎧ ⎪ ⎪ dx ⎪ ⎨⎪ dt = v, a0 → min, dv dt = Ωx + η(x, t)a0eδ, g0(x(t0), v(t0), t0)= 0, gk (x(tk ), v(tk ), tk )= 0, (6.1) ⎪ ⎪ tk ⎪ - t0 = D.t, t min 0 :( t0 :( t , max 0 ⎪⎩ с функцией Понтрягина (e, δ) ∈ R3 × R, |e| = 1, δ ∈ [0; 1], H = (px, v) + (pv , Ωx) + δη(x, t)a0 (pv , e) (6.2) и терминантом (краевой функцией Лагранжа) l = a0λa + (g0(x(t0), v(t0), t0), λ0) + /gk (x(tk ), v(tk ), tk ), λk \\ + (6.3) +(tk - t0 - D.t)λt + (tmin - t0)λ1 + (t0 - tmax)λ2 , 0 t0 0 t0 где px, pv , - липшицевы функции (сопряженные переменные), λa, λt, λ1 , λ2 o скалярные множители, λ0, λk - векторы. t0 t0 Понятно, что условия принципа максимума этой задачи являются частным случаем полученных ранее: § условия оптимальности по управлению (e, δ) = argmax (e∗,δ∗)∈R3×R H(t, x(t), v(t), e!, δ!): e = pv |pv | , δ = ( 1, |pv | > 0, ∀δ! ∈ [0; 1], |pv | = 0, (6.4) § уравнения сопряженной системы: ⎧ dpx ⎪ ⎪ dt ⎪ ∂η = -Ωxxpv - δa0|pv | ∂x, ⎨ dpv dt ⎪ ⎪ ⎪ dpa = -px, (6.5) § условия трансверсальности ⎧ ∂ ⎩ dt = -δη|pv |, 0) ∂ k ) ⎪px(t0)= ∂x(t ) (g0(x(t0), v(t0), t0),λ , px(tk )= - ∂x(t (gk (x(tk ), v(tk ), tk ),λ , ) ⎪ 0 ⎪ ⎪ ∂ 0) k ∂ k ) ⎨⎪pv (t0)= ∂v(t ) (g0(x(t0), v(t0), t0),λ , pv (tk )= - ∂v(t (g0(x(tk ), v(tk ), tk ),λ , ) 0 k ⎪pa(t0)= λa, pa(tk )= 0, ⎪ ⎪ ⎪ ∂ 0) 1 2 ∂ ( k ) ⎪⎩H(t0)= - ∂t (g0(x(t0), v(t0), t0),λ + λt + λt0 - λt0, H(tk )= ∂t gk (x(tk ), v(tk ), tk ),λ + λt, 0 § условия дополняющей нежесткости (tmin 1 k max 2 (6.6) 0 - t0)λt0 = 0, (t0 - t0 )λt0 = 0, (6.7) § условия нетривиальности и неотрицательности t0 + λ t0 |λ0| + |λk | + |λt| + λ1 2 t0 + λa > 0, λa ?: 0, λ1 t0 ?: 0, λ2 ?: 0. (6.8) ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 115 Очевидно, что оптимальная траектория для данной задачи не имеет участков особого управления. Действительно, если pv (t) ≡ 0 на некотором конечном интервале [τ0, τk ] (t0 :( τ0 < τk :( tk ), тогда из (6.5) следует, что и px(t) ≡ 0, pa(t) ≡ 0 ∀t ∈ [τ0, τk ], а значит, на всем интервале t ∈ [t0, tk ] px(t) ≡ 0, pv (t) ≡ 0, pa(t) ≡ 0, и тогда условия трансверсальности (6.6) при невырожt0 λ2 денных краевых условиях допускают только тривиальное решение λ0 = 0, λk = 0, λt = 0, λ1 = 0, t0 = 0, λa = 0, что противоречит (6.8) и, следовательно, оптимальной траектории δ ≡ 1. Условие для минимума начального ускорения можно записать отдельно как tk r η|pv |dt = λa ?: 0, (6.9) t0 причем λa > 0 для любой активной траектории. Для решения краевой задачи (6.1), (6.5)-(6.8) можно использовать метод продолжения по параметру с начальным приближением в виде ОМ-задачи (3.4). Для этого необходимо параметризовать правые части системы следующим образом: ⎧ dx ⎪ ⎪ dt ⎪ ⎪ dv ⎪ = v, \\ a0 = Ωx + ⎪⎨ dt η(x, t)τ + χ(x, t) |pv | (1 - τ ) pv , x dp ∂η p2 ∂χ τ ∈ [1; 0], (6.10) ⎪ = - (Ωxx, pv )- ⎪ v τ - a0|pv |(1 - τ ), dt ⎪ ⎪ ⎪ dpv ⎪⎩ dt = -px, ∂x 2 ∂x где τ - параметр продолжения. Тогда при τ = 1 уравнения (6.10) совпадают с уравнениями оптимального движения КА в ОМ-задаче (3.4), (3.8), а при τ = 0 - с уравнениями оптимального движения для задачи на минимум реактивного ускорения (6.1), (6.5). Краевые условия для этих задач совпадают, за исключением условия на минимум ускорения. Для его удовлетворения и обеспечения непрерывности продолжения необходимо воспользоваться возможностью произвольной нормировки сопряженных переменных в задаче (6.1) и выбрать значение константы λa в условии (6.9) равным соответствующему интегралу на траектории ОМ-задачи: λa = tk r v η|pOM |dt, (6.11) t0 v где |pOM | - модуль сопряженной переменной к скорости на траектории ОМ-задачи. Это никак не противоречат условиям оптимальности, так как этот интеграл всегда неотрицателен, а на активной траектории положителен. Вектор начального приближения в этой задаче, удовлетворяющий краевым условиям, может быть взят из решения ОМ-задачи и иметь вид (px(t0), pv (t0), a0)= (pOM (t0), pOM (t0), 0) . (6.12) x v Дифференциальные уравнения метода продолжения в данном случае имеют форму (4.7). Их интегрирование и дает решение задачи на минимум ускорения. Для продолжения по скорости истечения c ∈ [cmin; ∞] и представления задачи минимизации начальной тяги как однопараметрического семейства используется аналогичная методика со скоk ростью истечения в качестве параметра продолжения, где вектор Φ из (4.2) составлен из краевых условий (5.1), (5.7), (5.13), а вектор неизвестных параметров имеет вид (px(t0), pv (t0), pm(t0), T0). При этом понятно, что никакая дополнительная параметризация задачи не требуется, и для получения всех решений с конечной массой, принадлежащих отрезку m(tk ) ∈ [m∗; m0], достаточно проинтегрировать соответствующую систему (4.5) для задачи минимизации тяги при начальных 116 А. В. ИВАНЮХИН pa значениях вектора p в точке c-1 = 0, соответствующих задачи минимизации начального ускорения: ⎛pT ⎛ x (t0)⎞ ⎜ ⎜pT ⎟ ⎜ ⎜ v (t0)⎟ ⎜ x(t0) ⎞ 0 pa(t ) ⎟ v ⎟ ⎟ ⎟ ⎜ ⎟ = ⎜ tk ⎟ , λT = λa, (6.13) ⎜pT ⎜ ⎜ m(t0)⎟ a0 { ⎜- v ⎟ η|pa|dt⎟ T ⎝ ⎠ ⎝ 0 m0 t0 ⎠ m0a0 где индексы T и a соответствуют задачам на минимум тяги и ускорения. Как уже отмечалось ранее, на пути этого продолжения есть два характерных участка: без выхода на ограничение по конечной массе и с выходом на ограничение по конечной массе. В первом случае все решения являются гладкими и не содержат пассивных участков, во втором - могут содержать пассивные участки, и правые части внутренней задачи Коши метода продолжения (5.1), (5.6) будут претерпевать разрыв в моменты включения и выключения двигателя, что приводит к необходимости использования сглаженной формы функции тяги (4.16); в остальном продолжение на этих участках ничем не отличается. 2. ПЕРЕХОД К ЗАДАЧЕ С ЗАДАННЫМИ ЗНАЧЕНИЯМИ ПАРАМЕТРОВ Задача с минимальной тягой, помимо самостоятельного интереса и определения границы области существования, дает хорошее начальное приближение для задачи с заданным значением величины тяги. После получения зависимости минимальной тяги от скорости истечения Tmin(c), c ∈ [cmin, ∞], можно организовать переход к случаю с большей тягой, используя метод продолжения - продолжая по ее величине до заданного значения T ∈ [Tmin(c),T ∗]. Такой переход может быть реализован только для сглаженного управления (4.16), и каждая точка такой кривой продолжения будет решением ОТ-задачи для текущей величины тяги, с учетом погрешностей, вносимых сглаживанием. Уравнения, описывающие процесс продолжения по величине тяги, состоят из: o внутренней задачи Коши - сглаженной дифференциальной системы ОТ-задачи (3.13), (4.16) и присоединенных к ней уравнений сопряженной системы (3.17); o внешней задачи Коши - соответствующей случаю (4.6), где система краевых условий составлена из исходных условий ОТ-задачи (3.13) и условий трансверсальности принципа максимума (3.18), а сама система дифференциальных уравнений внешней задачи Коши получается ее дифференцированием; в начальной точке вектор неизвестных соответствует случаю минимальной тяги. Среди всех условий, входящих в краевую задачу, отдельно необходимо остановиться на услоm вии pm(tk ) = λk ?: 0 из (3.18), так как на решении ОТ-задачи оно может быть выполнено в виде равенства только в случае работы двигателя на всей траектории перелета, что соответствует задачи с минимальной тягой. А при продолжении по величине тяги необходимо обеспечить возможность появления пассивных участков, что может быть, как следует из вида оптимального m управления (3.16), только при λk > 0. Следовательно, для организации непрерывного продолжеm ния по тяге необходимо обеспечить увеличение λk вместе с ростом тяги. Поскольку в ОТ-задаче еще не было выбрано условие нормировки, это можно сделать достаточно произвольным образом, m например, положив его равным текущему изменению тяги λk = T0 - Tmin(c), или в нормированном λk к 1 случае m = (T0 - Tmin(c))/(T ∗ o Tmin(c)). Это обеспечит непрерывный сход с решения на минимум тяги для ОТ-задачи. 3. ЧИСЛЕННЫЕ ПРИМЕРЫ Элементарной единицей всех миссий в рамках Солнечной системы является прямой межпланетный перелет. Рассмотрим несколько таких перелетов: Земля-Марс и Земля-Меркурий. В качестве инерциальной декартовой системы координат в расчетах использовалась гелиоцентрическая эклиптическая система J2000. Движение КА рассматривалось в ньютоновском гравитационном поле, силовая функция имела вид: Ω= μ/r, где μ - гравитационный параметр центрального тела, ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 117 Начальное реактивное ускорение, мм/с2 0.4 2.0 0.355 0.31 0.265 0.5≤ 0.7≤ 0.9≤ 0.22 1.2 0.175 0.13 0.3≤ 1.0 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Удельный импульс, с Рис. 1. Перелет Земля-Марс. Область существования решений. r = |x| - удаление КА от центрального тела. Для определения движения планет использовалось готовое эфемеридное обеспечение JPL DE405/LE405. Краевые условия, записанные ранее в общем виде, соответствуют случаю равенства координат и скоростей КА и планет в начальный и конечный моменты. Они имеют следующий вид: ⎪⎪ 0 0 0 k k k ⎧x(t ) - x (t )= 0, x(t ) - x (t )= 0, ⎨ v(t0) - v0(t0)= 0, v(tk ) - vk (tk )= 0, ⎪ k ⎩⎪ m(t0) - m0 = 0, m∗ - m(tk ) :( 0, (8.1) где (xi, vi) - положение и скорость начальной (i = 0) и конечной (i = k) точки (планеты). В этом случае условия трансверсальности по дате старта (для tmin < t0 < tmax) при фиксиро- 0 0 ванном времени перелета имеют следующий вид: o для ОМ-задачи χ(x0, t0) (pv (t0), pv (t0))- χ(xk , tk ) (pv (tk ), pv (tk )) = 0; (8.2) o для ОТ-задачи (на минимум тяги и максимум конечной массы) δ(t0)η(x0, t0)T0Ψ(t0) - δ(tk )η(xk , tk )T0Ψ(tk )= 0. (8.3) При этом в приведенных условиях не учитывается тот факт, что для КА и планет используются, вообще говоря, разные модели движения. Так, движение КА рассматривается в рамках ограниченной задачи двух тел, в то время как движение планет, определяемое по эфемеридам JPL DE405/LE405, учитывает значительно больше факторов, таких как наличие других планет и тел солнечной системы, релятивистские эффекты и т. д. Однако в большинстве случаев такое упрощение является несущественным. Для простоты расчетов реактивная мощность, тяга и скорость истечения считаются постоянными (χ(x, t)= 1, η(x, t)= 1, w(x, t)= 1) . Пример 1. Перелет Земля-Марс с оптимальной датой отлета от Земли, выбираемой из синодического периода 23.04.2019-11.06.2021 и временем перелета 380 суток. На рисунке 1 показана зависимость минимального начального реактивного ускорения от скорости истечения для этой задачи. Как отмечалось ранее, область существования решений для ограниченной конечной массы имеет предельные случаи, соответствующие бесконечному удельному импульсу и бесконечной тяге, что можно видеть на рис. 1. Участок монотонного убывания начального реактивного ускорения при уменьшении удельного импульса соответствует отсутствию ограничения на конечную массу, он построен вплоть до значения относительной конечной массы 0,3, на всей этой линии оптимальные решения не имеют пассивных участков на траектории перелета. От нее отходят линии, соответствующие случаям выхода решения с минимальной тягой на 118 А. В. ИВАНЮХИН Рис. 2. Перелет Земля-Марс. T0 = 1,2Tmin Рис. 3. Перелет Земля-Марс. T0 = 2Tmin ограничение по конечной массе: на каждой из них она остается постоянной и равна на приведенном графике 0,9, 0,7, 0,5 и 0,3 от начальной массы. При этом на траектории есть пассивные участки, для получения этих решений в расчетах использовался параметр сглаживания, равный 10-5. Эти линии с минимальной тягой ограничивают области соответствующих решений (на графике они выделены заливкой). Таким образом, все возможные решения ОТ-задачи данного перелета с продолжительностью меньше или равной рассматриваемой и угловой дальностью из отрезка [0; 2π] находятся над этими линиями, и область над графиком является областью существования решений. Как видно на рисунке 1, она ограничена только снизу, сверху нет ограничения ни по удельному импульсу, ни по тяге. В частности, для рассматриваемой задачи минимальное начальное реактивное ускорение при бесконечном удельном импульсе (постоянной массе КА) равняется 0,21255 мм/с2. Функционал ОМ-задачи равен 0,31912 м2/с3. Чтобы проследить эволюцию управления при переходе с границы области существования в ее внутреннюю часть, на рис. 2 и 3 приведены проекции траекторий на плоскость XY решений с начальным реактивным ускорением больше минимального в 1,2 и 2 раза соответственно и удельным импульсом 3100 сек. (что соответствует следующим значениям начального реактивного ускорения: 0,1951 мм/с2, 0,2342 мм/с2, 0,3903 мм/с2). Эти решения отмечены на рис. 1 точками. Конечная масса КА для случая минимального реактивного ускорения равна 0,7892 от начальной (в этом случае двигатель работает постоянно); при увеличении ускорения в 1,2 раза она увеличивается до 0,8264, а в 2 раза - до 0,829. Видно, что приведенные решения имеют в начале и конце активные участки, что говорит о наличии оптимальной даты старта для такого времени перелета при заданных параметрах, что в общем случае не обязательно, так как при увеличении тяги оптимальнее может оказаться не сдвигать дату старта, а уменьшать общее время перелета, тогда в начале или конце будет образовываться пассивный участок (такое решение будет оптимально как по дате старта, так и по времени перелета). Оптимальные даты старта для случаев 1, 1,2, 2 соответственно равны: 14.04.2020, 27.04.2020, 13.05.2020. Пример 2. Перелет Земля-Меркурий с датой отлета от Земли 19.03.2020 и временем перелета 2140 суток, 13 полных витков вокруг Солнца. Такие продолжительные перелеты особенно ярко показывают многоэкстремальный характер решений подобных задач. Очевидно, что могут существовать решения, отличающиеся числом целых витков, и заранее неизвестно, какое из них лучше. ОБЛАСТЬ СУЩЕСТВОВАНИЯ РЕШЕНИЙ В ЗАДАЧЕ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ КОСМИЧЕСКИМ АППАРАТОМ 119 Начальное реактивное ускорение, мм/с2 0.18 0.16 0.14 0.3≤ 2.0 0.5≤ 0.7≤ 0.12 0.1 0.08 1.2 1.0 0.06 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Удельный импульс, с Рис. 4. Перелет Земля-Меркурий. Область существования решений. Рис. 5. Перелет Земля-Меркурий. T0 = 1,2Tmin Рис. 6. Перелет Земля-Меркурий. T0 = 2Tmin Решения с тягой больше минимальной были получены с помощью продолжения по величине тяги с параметром сглаживания 10-5. На рис. 4 показана область существования решений, аналогичная приведенной ранее. Как видно, ее структура схожа с предыдущей, однако сход «ветвей» с постоянной конечной массой более резкий. В этом случае минимальное начальное реактивное ускорение при бесконечном удельном импульсе (постоянной массе КА) равняется 0,10711 мм/с2, а функционал ОМ-задачи равен 0,51608 м2/с3. Так же, как и в случае Марса, для этого перелета интересно проследить эволюцию управления, вида траектории и изменение относительной конечной массы при увеличении тяги. На рис. 5 и 6 приведены вид проекции траектории на плоскость XY при увеличении тяги до 1,2 и 2 от минимальной с удельным импульсом 3100 сек. На начальной траектории (с минимальной тягой) начальное реактивное ускорение равно 0,082 мм/с2, а относительная конечная масса КА - 0,5014, при увеличении тяги в 1,2 раза она составила 0,5311, а при росте тяги в 2 раза - 0,5465. По рис. 5 и 6 можно видеть, что с увеличением тяги перестройка решения происходит таким образом, что в конце образуется пассивный участок, в связи с чем сокращается время перелета и определяется ее оптимальное значение. Также на рис. 6 прослеживается очевидная локализация 120 А. В. ИВАНЮХИН активных участков в районах апоцентра и перицентра траектории перелета, что иллюстрирует эволюцию в сторону импульсного решения с их приложением в этих точках при стремлении тяги к бесконечности.