Известия РАН. Механика твердого тела, 2019, № 4, стр. 48-62

1. Введение. Одно из направлений полезного применения космических тросовых систем это возможность изменения орбитальных параметров с помощью проводящих ток тросов [1]. Взаимодействие проводящих тросов с магнитным полем Земли (в перспективе с магнитными полями других планет) приводит к возникновению силы Ампера, которая, действуя на систему в целом, изменяет орбиту ее центра масс. Электродинамические тросовые системы (ЭДТС) могут функционировать в пассивном и активном режимах [1]. При движении в магнитном поле в проводнике (тросе) возникает электродвижущая сила (ЭДС) индукции, которая является причиной возникновения электрического тока определенного направления. Сила Ампера, действующая на элементарный участок троса, определяется известной формулой [1]

где I – величина тока, $B$ – вектор магнитной индукции, $\tau $ – единичный вектор, направленный по касательной к тросу.

В пассивном режиме использование ЭДТС приводит к возникновению силы (1.1), тормозящей движение системы по орбите [1–4]. Активный режим функционирования ЭДТС соответствует случаю, когда наряду с ЭДС индукции имеется еще управляющая разность потенциалов [1], которая возникает за счет бортовых источников энергии концевых тел, и с помощью которой можно изменять величину и направление тока в тросе, а значит, и направление действия силы (1.1). Тем самым можно не только тормозить орбитальное движение, но и осуществлять стабилизацию орбитальных параметров или переводить систему на более высокую орбиту [5, 6]. Перспективность применения ЭДТС в космосе подтверждает тот факт, что до настоящего времени проведено уже более десятка реальных тросовых экспериментов с проводящими ток тросами на околоземных орбитах [3].

2. Постановка задачи. Известно, что положения равновесия ЭДТС, расположенные вблизи местной вертикали, так называемые “рабочие” положения равновесия системы, неустойчивы [1, 2, 7, 8]. Неустойчивость движения системы проявляется в возрастании ее маятниковых колебаний, что в конечном итоге (в зависимости от величины тока в тросе) ведет или к переходу системы во вращение, или к большой деформации троса, когда трос сминается и становится не натянутым. И в том, и в другом случае ЭДТС прекращает выполнять свои функции по целенаправленному изменению орбитальных параметров системы. Положение еще более усложняется для наклонных орбит с относительно большим эксцентриситетом, так как в этом случае всегда реализуется пространственный случай движения, при этом имеют место угловые маятниковые колебания системы как в плоскости, так вне плоскости орбиты [9]. Для наклонных эллиптических орбит на систему начинают действовать периодические возмущения, связанные с изменением аргумента широты и с угловыми маятниковыми колебаниями троса перпендикулярно плоскости орбиты. Это существенно усложняет решение задачи стабилизации движения ЭДТС в плоскости орбиты. Поэтому вопросам стабилизации движения ЭДТС в настоящее время уделяется большое внимание [10–13], при этом используются как пассивные [10], так и активные методы стабилизации [11–13]. При реализации активных методов чаще всего используются алгоритмы стабилизации, основанные на определении поправок к номинальным значениям тока [11–13], то есть $I = {{I}_{n}} + \Delta I$, где I_n – номинальная величина тока, ΔI – управляющие стабилизирующие воздействия, определяемые на основании принципа обратной связи. Поправки ΔI в известных работах строятся в зависимости от значений производных углов, определяющих положение троса относительно местной вертикали. Причем при построении алгоритмов стабилизации предполагается, что форма троса близка к прямой линии [11–13]. Однако, как известно, при действии распределенной нагрузки при движении ЭДТС по околоземной орбите имеют место изгибные колебания троса, причем, если масса троса много меньше массы концевых тел (то есть трос можно считать невесомым), то форма троса в любой момент времени представляют собой часть окружности некоторого радиуса [1], [14]. Почти периодические изменения радиуса этой окружности как раз и характеризуют его изгибные колебания. Взаимодействие изгибных колебаний троса с маятниковыми колебаниями системы имеет сложный характер, в частности, возможны случаи синхронизации изгибных и маятниковых колебаний в системе [14]. Неучет этих особенностей в движении ЭДТС, как показано ниже в представленной работе, может привести к потере устойчивости “рабочих” положений ЭДТС при использовании регуляторов, учитывающих только угловые скорости вращения системы относительно ее положений равновесия.

В настоящей работе предлагаются и исследуются алгоритмы стабилизации движения ЭДТС, в которых регулирование осуществляется не только по угловым скоростям вращения системы относительно вертикали, но и по скорости изменения расстояния между концевыми точками системы, то есть производной по времени хорды окружности, характеризующей изгибные колебания троса. На основании анализа предлагаемого подхода даются рекомендации по наиболее рациональному выбору параметров ЭДТС и ее режимов движения, обеспечивающих устойчивое движение системы на эллиптических и наклонных орбитах на больших интервалах времени. Сформулированные выводы подтверждаются результатами численного моделирования движения ЭДТС, состоящей из двух малых космических аппаратов и проводящего ток троса. Для оценки эффективности рассматриваемого подхода используются две математические модели, учитывающие изгибные колебания троса. Первая модель получена методом Лагранжа для невесомого троса в работе [14] с использованием методики, изложенной в [1]. В этой модели в отличие от [1] рассматривается случай произвольного соотношения масс концевых тел и наклонные эллиптические орбиты. Вторую модель движения ЭДТС можно назвать моделью с распределенными параметрами, так как в ней трос представляется совокупностью материальных точек. Эта модель учитывает распределение массы троса в гравитационном поле и изменение вектора магнитной индукции в пределах троса. Подобные модели использовались во многих работах для проведения поверочных расчетов при обосновании использования более простых моделей движения космических тросовых систем [15–18].

3. Уравнения движения системы относительно центра масс с невесомым тросом. Уравнения, описывающие колебания ЭДТС относительно центра масс, записываются в правой орбитальной прямоугольной системе координат $C{{x}_{o}}{{y}_{o}}{{z}_{o}}$ (рис. 1), начало которой совпадает с центром масс системы, ось $C{{x}_{o}}$ направлена по его радиус-вектору ${{R}_{c}}$, ось $C{{y}_{o}}$ – перпендикулярно ${{R}_{c}}$ и направлена в сторону орбитального движения. При записи уравнений используется предположение, что на трос действует постоянная распределенная нагрузка от сил Ампера (1.1), то есть изменением вектора магнитной индукции B вдоль троса пренебрегается. В этом случае форма невесомого троса в любой момент времени есть часть окружности некоторого радиуса [1]. Для этого случая в монографии [1] методом Лагранжа для экваториальной круговой орбиты в обобщенных координатах $r,{\theta ,}\,{\varphi }$, где $r$ – расстояние между концевыми точками (хорда окружности), ${\theta ,}\,{\varphi }$ – углы, определяющие положение хорды (или системы координат $C{{x}_{t}}{{y}_{t}}{{z}_{t}}$ с ней связанной, рис. 1) относительно системы координат $C{{x}_{o}}{{y}_{o}}{{z}_{o}}$ (рис. 1), получены соответствующие уравнения. Эти уравнения построены для случая, когда масса одного концевого тела много больше массы другого. В работе [14] проведено обобщение уравнений [1] на случай произвольного распределения масс концевых точек и на эллиптические наклонные орбиты. В этом случае имеем [14]

где ${{m}_{e}} = {{m}_{a}}{{m}_{b}}{\text{/}}m$, ${{m}_{{a,b}}}$ – массы концевых точек (для определенности далее будем считать, что m_b – верхнее тело), $m = {{m}_{a}} + {{m}_{b}}$, ${\omega } = \dot {\vartheta } = {{({\mu /}{{p}^{3}})}^{{0.5}}}{{{\nu }}^{2}}$, ${\dot {\omega }} = \ddot {\vartheta } = - 2{\mu }{{{\nu }}^{3}}e\sin \vartheta {\text{/}}{{p}^{3}}$, ${\mu }$ – гравитационный параметр Земли, ${\nu } = 1 + e\cos \vartheta $, e – эксцентриситет орбиты, $\vartheta $ – истинная аномалия, p – параметр орбиты, ${{Q}_{{{\theta },r,{\varphi }}}}$ – обобщенные силы Ампера.

Наряду с предположением о форме троса при выводе уравнений использовались следующие допущения: 1) гравитационное поле ньютоновское; 2) длина хорды много меньше расстояния до притягивающего центра $r \ll {{R}_{c}}$; 3) концевые тела – материальные точки; 4) орбита центра масс системы – невозмущенная. Если φ = 0, то движение системы происходит в плоскости орбиты.

Используется модель магнитного поля Земли – прямой диполь [1], поэтому

где ${{B}_{0}} = {{{\mu }}_{m}}{\text{/}}R_{c}^{3}$, μ_m = 8 × 10⁶ Тл км³ – магнитный момент земного диполя, ${{e}_{R}} = {{R}_{c}}{\text{/}}{{R}_{c}}$, ${{e}_{z}}$ – единичный вектор, направленный по оси вращения Земли; $\left( {{{e}_{z}} \cdot {{e}_{R}}} \right)$ – скалярное произведение.

Обобщенные силы вычисляются из следующих выражений [14]

где $i$ – наклонение орбиты, u – аргумент широты, $\Delta = 0.5r\left( {{{m}_{b}} - {{m}_{a}}} \right){\text{/}}m$, ${\psi }$ – угол между касательной к окружности, проведенной в концевых точках, и хордой их соединяющей; $I$ – величина тока, , $\Delta {{Q}_{{\theta }}}$ = 3B₀IrΔsinicosθsinφsinu.

В отличие от [14], где использовалась более простая модель для магнитного поля $B = {{B}_{o}}{{e}_{z}}$, здесь приводятся обобщенные силы, которые соответствуют модели (3.4). Величина тока I считается положительной, если он направлен от точки m_a к точке m_b.

Угол ${\psi }$ определяется из геометрического соотношения [1]

где L – недеформированная длина троса, ${{{\gamma }}_{t}} = 2{{E}_{t}}{\psi /}\left( {2{{E}_{t}}{\psi } - B\left| I \right|L} \right)$ – относительное удлинение троса, ${{E}_{t}}$ – жесткость троса. Предполагается, что при растяжении трос подчиняется закону Гука.

Система (3.1–3.3) совместно с выражениями (3.4–3.6) описывает движение ЭДТС относительно своего центра масс с учетом изгибных колебаний троса. Изгиб троса характеризуется расстоянием между концевыми точками r и соответственно углом ${\psi }$ (3.6), значения которых определяют радиус части окружности ${{r}_{t}} = r{\text{/}}(2\sin {\psi )}$.

4. Положения относительного равновесия системы и анализ линеаризованных уравнений. Положения относительного равновесия можно определить только для экваториальной и круговой орбиты центра масс системы. Дело в том, что для наклонных и эллиптических орбит в системе существуют периодические возмущения, связанные с изменением аргумента широты $u$ и угла ${\varphi }$, что исключает возможность определения неизменных равновесных положений системы.

Для экваториальной круговой орбиты положения равновесия, близкие к вертикали, определяются из уравнений (3.1)–(3.2) и имеют вид $\left( {{\varphi } = {\dot {\varphi }} = 0} \right)$

где ${\sigma } = {{{\mu }}_{m}}I\left( {{{m}_{b}} - {{m}_{a}}} \right){\text{/}}3{\mu }{{m}_{a}}{{m}_{b}}$

Остальные переменные системы (3.1–3.3) определяются из выражений

Положения равновесия (4.1–4.2) иногда называют “рабочими” положениями равновесия [1], так как с точки зрения использования ЭДТС для изменения орбитальных параметров системы они являются наиболее предпочтительными. Если ЭДТС совершает ограниченные колебания относительно положений равновесия (4.1–4.2), которые располагаются вблизи вертикали, то равнодействующая сил Ампера (1.1) почти параллельна вектору скорости центра масс системы и изменение высоты орбиты будут близки к максимальным. Очевидно, что условие существования положений равновесия (4.1–4.2) (условие статической устойчивости) будет иметь вид $\left| {\sigma } \right| < 1$.

Линеаризованная система относительно положений равновесия (4.1–4.2) имеет стандартную форму и приводится здесь в безразмерном виде

где $y = \left( {\Delta {\theta ,\theta }',\Delta r,r',{\varphi ,\varphi '}} \right)$ – вектор отклонений обобщенных координат и скоростей от положений равновесия, ${\tau } = \Omega \,t$ – безразмерное время, $\Omega $ – угловая скорость орбитального движения. При этом матрица $A$ частных производных имеет следующий вид

В безразмерных выражениях (4.3–4.4) сохранены те же обозначения, что и в исходных уравнениях (3.1–3.3), то есть, сделаны формальные замены $r \Rightarrow r{\text{/}}L$, ${{r}_{{1,2}}} \Rightarrow {{r}_{{1,2}}}{\text{/}}L$, $\Delta r \Rightarrow \Delta r{\text{/}}L$. В линеаризованной системе (4.3–4.4) уравнения для переменных ${\varphi ,}\Delta {\dot {\varphi }}$ отделяются и имеют всегда чисто мнимые корни. Остальные собственные числа матрицы (4.4) определяются из следующего характеристического уравнения

где ${{p}_{2}} = - \left( {{{A}_{{21}}} + {{A}_{{43}}} + {{A}_{{24}}}{{A}_{{42}}}} \right)$, ${{p}_{1}} = - {{A}_{{24}}}{{A}_{{41}}}$, ${{p}_{0}} = {{A}_{{21}}}{{A}_{{43}}}$

Если ${{p}_{1}} = 2I{{B}_{o}}\left( {{{m}_{a}} - {{m}_{b}}} \right){\text{/}}({{m}_{a}}{{m}_{b}}{{{\omega }}^{2}}) = 0$, то есть когда массы концевых точек одинаковы и момент от сил Ампера отсутствует, то остальные уравнения также имеют чисто мнимые корни. Если ${{p}_{1}} \ne 0$, то решения уравнений (4.3) всегда неустойчивы, что также было отмечено в монографии [1] для частного случая ${{m}_{b}} \gg {{m}_{a}}$.

Для анализа эволюционных движений ЭДТС при постоянном токе первые четыре уравнения системы (4.3) приводятся к двум уравнениям второго порядка

В нормальных координатах система (4.6) принимает

где $V$– матрица собственных векторов матрицы $c$.

Здесь необходимо отметить, что при малых изгибах троса, когда ${{{\psi }}_{{1,2}}} \to 0$ (${{{\psi }}_{{1,2}}} > 0$), имеем ${{{\omega }}_{2}} \gg {{{\omega }}_{1}}$, так как ${{A}_{{43}}} \to \infty $ (см. (4.4)). Кроме того, первая нормальная координата совпадает с исходной переменной ${{{\xi }}_{1}} = \Delta {\theta }$, а вторая координата ${{{\xi }}_{2}} = {\delta }\Delta {\theta + }\Delta r$.

Представление системы (4.6) в виде (4.7) позволяет выделить в явном виде слагаемые, непосредственно влияющие на неустойчивость движения системы (4.4). Это компоненты $ - {\delta }{{A}_{{24}}}$ и ${\delta }{{A}_{{24}}}$ матрицы ${{b}_{n}}$ (4.7), причем эти слагаемые пропорциональны малому параметру δ, так как ${\delta } \to 0$ при ${{{\psi }}_{{1,2}}} \to 0$.

В связи с этим применим к системе (4.7) стандартный метод медленно изменяющихся амплитуд [19], применяя замену переменных

где ${{K}_{{1,2}}}$ и ${{{\gamma }}_{{1,2}}}$ – амплитуды и фазы колебаний, и рассматривая правую часть системы (4.7) как возмущение.

После усреднения по фазам нетрудно получить

где $k = 1,2$ – индекс положения равновесия (4.1).

В усредненных уравнениях для простоты сохранены те же обозначения для амплитуд колебаний ${{K}_{{1,2}}}$, что и для исходных не усредненных переменных (4.8). Величина коэффициента ${\eta }$, определяющая скорость возрастания амплитуд колебаний, есть приближенная оценка модулей вещественных частей собственных значений матрицы A. Так, например, для исходных данных, приведенных в таблице 1, для круговой (e = 0) и экваториальной (i = 0) орбиты имеем ${\eta } = 3.66 \times {{10}^{{ - 3}}}$ и $\left| {\operatorname{Re} ({{\lambda }_{{1,2,3,4}}})} \right| = 3.58 \times {{10}^{{ - 3}}}$, при этом ω₁ = 1.71, ω₂ = 13.62, δ = 3.63 × 10^–3. Частота колебаний по углу ${\varphi }$ равна 1.415.

Таким образом, если концевое тело (m_b) смещается противоположно направлению орбитального движения $\sin 2{{{\theta }}_{k}} > 0$, то неустойчивость движения системы обусловлена возрастанием амплитуд маятниковых колебаний в плоскости орбиты (по углу $\Delta {\theta }$). Если же $\sin 2{{{\theta }}_{k}} < 0$, то неустойчивость движения проявляется в возрастании амплитуд продольных (или изгибных) колебаний троса $\Delta r = {{{\xi }}_{2}} - {\delta }\Delta {\theta }$, так как ${\delta }$ – малый параметр. Причем, как следует из (4.9), знак и величина $\sin 2{{{\theta }}_{k}}$ определяется величиной тока и соотношением масс концевых тел, то есть от знака и величины момента равнодействующей сил Ампера.

5. Синтез приближенно оптимального регулятора для экваториальной орбиты. Для синтеза приближенно оптимального регулятора воспользуемся принципом динамического программирования Беллмана [20] в сочетании с методом усреднения. Метод усреднения использовался во многих работах для решения различных задач построения приближенно оптимальных управлений, например, в [21–24]. Здесь применяется процедура синтеза линейного регулятора, аналогичная алгоритму синтеза при решении задачи стабилизации движения твердого тела вокруг неподвижной точки [25] для случая, близкого к движению волчка Лагранжа. Рассматриваемый алгоритм, как показано ниже, позволит для частного случая круговой орбиты получить уравнения для регулятора в аналитической форме.

Для решения задачи синтеза регулятора рассматривается система (4.7) с управлением и малым параметром ${\varepsilon }$

где $\Delta u = \Delta I{\text{/}}{{I}_{n}}$ – управление в переменных системы (4.6), M – матрица-столбец коэффициентов при управлении, ${{M}_{1}} = {{B}_{o}}{{I}_{n}}\left( {{{m}_{b}} - {{m}_{a}}} \right){\text{/}}2{{m}_{a}}{{m}_{b}}{{{\omega }}^{2}}$, ${{M}_{2}} = {{B}_{o}}{{I}_{n}}{\text{/}}2{{m}_{e}}{{{\omega }}^{2}}$. Здесь необходимо отметить, что при выбранном способе управления (поправки к величине тока) линеаризованная система (4.3) неуправляема по переменным ${\varphi ,\varphi '}$.

Введение малого параметра ${\varepsilon }$ в систему (5.1) является формальным (масштабный коэффициент), то есть предполагается, что движение системы (5.1) слабо отличается от гармонических колебаний в нормальных координатах. Причем управление также ищется в классе малых управляющих воздействий, сопоставимых по порядку с возмущениями.

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

где $P\left( {{{K}_{1}},{{K}_{2}},\Delta u} \right) = {{s}_{1}}K_{1}^{2} + {{s}_{2}}K_{2}^{2} + {{s}_{3}}\Delta {{u}^{2}}$, ${{s}_{{1,2,3}}}$ – положительные безразмерные весовые коэффициенты, причем полагается, что $\sum {{{s}_{k}} = 1} $.

Движение системы (5.1) рассматривается на асимптотически большом интервале времени , где ${{L}_{t}} < \infty $ – некоторая константа. В этом случае функционал изменяется на величину порядка O(1).

В соответствии с принципом динамического программирования, оптимальное управление определяется из условия [20], [25]

где $W\left( {K,{\gamma },{\text{u}}} \right)$ – функция Беллмана [25], ${\text{K,}}\,{\gamma }$ – вектора амплитуд и фаз (4.8), а точка $\left( \cdot \right)$ означает скалярное произведение векторов.

Дифференциальные уравнения для переменных ${\text{K,}}\,{\gamma }$ в соответствии с методом медленно изменяющихся амплитуд [19] имеют вид

С учетом вида правых частей уравнения (5.4–5.5) и функции $P$ выражение, стоящее под знаком минимума (5.3), представляет собой квадратичный степенной полином по управлению $u$. Поэтому из условия минимума этого полинома оптимальное управление определяется в виде

Подставив управление (5.6) в выражение (5.3), получаем уравнение в частных производных для определения функции $W$

Уравнение (5.7) решается приближенно методом усреднения. В соответствии с этим решения для уравнений (5.4–5.5), (5.7) ищутся в виде асимптотических рядов

где ${{K}_{{a{v}}}},{{{\gamma }}_{{a{v}}}}$ – векторы переменных метода усреднения (Averaging method), ${{W}_{0}}({{a}_{{a{v}}}})$ – функция Беллмана, определенная в силу усредненной системы первого приближения; ${{\nu }_{1}}\left( {{{K}_{{a{v}}}},{{{\gamma }}_{{a{v}}}}} \right)$, ${{{\mu }}_{1}}\left( {{{K}_{{a{v}}}},{{{\gamma }}_{{a{v}}}}} \right)$, ${{W}_{1}}\left( {{{K}_{{a{v}}}},{{{\gamma }}_{{a{v}}}}} \right)$ – поправки к первому приближению метода усреднения, имеющие нулевое среднее.

Подставляя (5.8) в уравнение (5.7) и усредняя по фазам, получим

где $\left\langle \cdot \right\rangle $ – оператор усреднения по фазам.

При записи усредненного уравнения (5.9) для упрощения обозначений индекс $(a{v})$ для новых переменных опущен. При усреднении также учитывалось, что $\left\langle {\sum\limits_{k = 1}^2 {{{\partial {{W}_{1}}} \mathord{\left/ {\vphantom {{\partial {{W}_{1}}} {\partial {{{\gamma }}_{k}}{{{\omega }}_{k}}}}} \right. \kern-0em} {\partial {{{\gamma }}_{k}}{{{\omega }}_{k}}}}} } \right\rangle $ = 0.

Функция ${{W}_{0}}\left( {{{K}_{1}},{{K}_{2}}} \right)$, соответствующая первому приближению метода усреднения (5.9), ищется в виде положительно определенной квадратичной формы

Подставляя (5.10) в уравнение первого приближения (5.9), нетрудно найти

при этом ${{a}_{{12}}} = 0$.

Полная производная функции ${{W}_{0}}$ по времени, определенная в силу уравнений первого приближения, будет иметь вид

где $\langle P({{K}_{1}},{{K}_{2}},\Delta {{u}^{o}})\rangle $ – знакоопределенная по ${{K}_{{1,2}}}$ и положительная функция, определенная на решениях усредненной системы первого приближения.

Положительно определенная функция ${{W}_{0}}$ в соответствии с уравнением (5.12) монотонное убывает. Отсюда следует монотонное убыванию амплитуд колебаний ${{K}_{{1,2}}}$. Соответствие свойств устойчивости усредненной и исходной систем уравнений рассмотрено, например, в монографии [26] и обеспечивается общими теоремами метода усреднения.

Подставляя функцию (5.10) в выражение для оптимального управления (5.6), переходя к исходным переменным (учитывая (4.8) и ${\xi } = {{V}^{{ - 1}}}x$), получим уравнение регулятора в виде

После подстановки (5.13) в уравнения (5.4–5.5) и проведения процедуры усреднения, усредненная система с управлением будет иметь вид

где ${{{\eta }}_{k}}$ – коэффициенты, определяющие скорость уменьшения амплитуд колебаний; ${{N}_{1}} = {{M}_{1}}$, ${{N}_{2}} = {{M}_{2}} + {\delta }{{M}_{1}}$

При постоянных коэффициентах ${{{\eta }}_{k}}$ уравнения (5.14) нетрудно проинтегрировать и получить аналитическое решение

Для исходных данных, приведенных в таблице (при $e = i = 0$, ${{I}_{n}} = 0.1\,A$), и при равных весовых коэффициентах ${{s}_{{1,2,3}}} = 1{\text{/}}3$, имеем ${{k}_{{\theta }}} = - 8.30$, ${{k}_{r}} = - 1.03$, ${{{\eta }}_{1}} = - 1.38$, ${{{\eta }}_{2}}$ = –0.35.

Применение метода усреднения позволяет определить аналитические соотношения между параметрами ЭДТС, весовыми коэффициентами критерия оптимальности (5.2), коэффициентами регулятора и получающимися характеристиками переходного процесса в системе регулирования.

Если рассматривать движение усредненных систем на фазовой плоскости (K₁, K₂), то очевидно особая точка “седло” при введении управления превращается в особую точку “устойчивый узел”. Изменяя соотношения между весовыми коэффициентами ${{s}_{{1,2,3}}}$ всегда можно обеспечить требуемые характеристики переходных процессов в системе.

6. Математическая модель движения системы с распределенными параметрами. При построении математической модели движения ЭДТС (3.1–3.3) было принято ряд допущений. Здесь можно отметить: 1) невесомость троса; 2) форма троса – часть окружности; 3) неизменность вектора магнитной индукции B в пределах длины троса; 4) постоянство орбитальных параметров центра масс системы. Поэтому для сравнительного анализа здесь также используется более сложная модель движения, в которой трос представляется дискретной совокупностью материальных точек, связанных между собой упругими односторонними механическими связями. Подобные модели использовались во многих работах, например, в [15–18]. Особенностью данной модели является учет распределенных нагрузок от сил Ампера.

После дискретизации троса уравнения ЭДТС записываются в виде

где ${{r}_{k}}$, ${{V}_{k}}$ и ${{m}_{k}}$ – радиус-вектор, скорость и масса k-й материальной точки; ${{F}_{k}}$ – сила Ампера, действующая на k-ю материальную точку; T_k – силы натяжения троса, действующие между k-й и k + 1-й точками и приложенные к k-й точке; ${{G}_{k}}$ – гравитационные силы. При записи системы (6.1) полагается, что ${{m}_{1}} = {{m}_{a}}$, ${{m}_{n}} = {{m}_{b}}$.

Для определения ${{F}_{k}}$ сначала вычисляются силы Ампера, действующие на k-й участок троса

где I – величина тока, ${{B}_{k}}$ – вектор магнитной индукции, определенный в середине отрезка троса; ${{\tau }_{k}} = \left( {{{r}_{{k + 1}}} - {{r}_{k}}} \right){\text{/}}\left| {{{r}_{{k + 1}}} - {{r}_{k}}} \right|$.

Распределение сил ${{S}_{k}}$ по материальным точкам ЭДТС осуществляется следующим образом

При использовании системы (6.1) гравитационные силы G_k (центральное ньютоновское поле) и силы натяжения троса T_k определяются так же, как в работе [18], где при растяжении применялся закон Гука, и рассматривалась односторонняя механическая связь между точками.

7. Численные результаты. Сравнение результатов моделирования управляемого движения ЭДТС по различным моделям проводилось для экваториальных и наклонных, круговых и эллиптических орбит. Для уменьшения уровня управляющего воздействия $\Delta I$ и для обеспечения приблизительно равной длительности переходных процессов по различным координатам были выбраны не равные весовые коэффициенты в критерии оптимальности (5.2), удовлетворяющие следующим соотношениям: ${{s}_{1}}{\text{/}}{{s}_{3}} = 0.01$, ${{s}_{2}}{\text{/}}{{s}_{3}} = 0.1$. При этом имеем ${{k}_{{\theta }}} = - 0.85$, ${{k}_{r}} = - 0.32$, ${{{\eta }}_{1}} = - 0.14$, ${{{\eta }}_{2}} = - 0.11$. Исходные данные соответствуют таблице ($i = 0,\,\,e = 0,\,\,I = 0.1\,{\text{A}}$). Изменение направления тока ($I = - 0.1\,{\text{A}}$) практически не влияет на полученный результат.

Если рассматривать линеаризованные не усредненные уравнения (4.3), то совпадение результатов с аналитическим решением (5.15) почти идеальное по обеим координатам ($\Delta {\theta },\,\,\Delta r$). Сравнение аналитического решения (5.15) с нелинейными моделями (3.1–3.3) и (6.1) показывает, что малая погрешность определения амплитуд колебаний сохраняется по переменной $\Delta {\theta }$, а по переменной $\Delta r$ имеет место совпадение результатов по времени переходного процесса. Данный вывод иллюстрируется рис. 2, 3, где рис. 2 соответствует модели (3.1–3.3), а рис. 3 – модели (6.1). Результаты, приведенные на рис. 2, 3, соответствуют случаю, когда в начальном положении ЭДТС расположена на вертикали и по тросу начинают пропускать ток. При реализации управляемого движения начальные малые возмущения по углу ${\varphi }$, по которому система не управляема, например, ${{{\varphi }}_{0}} < 1$, в соответствии с нелинейными моделями (3.1–3.3) и (6.1) остаются ограниченными, причем их амплитуды колебаний немного уменьшаются. При использовании модели (6.1) количество материальных точек, моделирующих распределенную массу троса, была равной шести. Дальнейшее увеличение дискретизации троса, даже в два раза, практически не изменяло полученных результатов. Учет массы троса в модели (6.1) немного изменяло положение равновесие системы, однако из-за сравнительно ее малой величины (для таблицы масса троса $0.2\,{\text{кг}}$) это изменение было не значительным (0.001 по углу ${\theta }$ и 0.005 по переменной r).

Использование полученного управления при тех же исходных данных для эллиптических и наклонных орбит также возможно, однако приводит к периодическим возмущениям, связанным с изменением аргумента широты (3.5) и колебаниями по углу ${\varphi }$. Причем в соответствии с уравнениями (3.1–3.3) возмущения по углу ${\varphi }$ имеют место, если даже в начальном состоянии ЭДТС находится в плоскости орбиты. В этом случае в окрестности положений равновесия (4.1–4.2) системы возникают ограниченные колебания, близкие к периодическим. На рис. 4 приводятся зависимости, полученные по модели (6.1), и иллюстрирующие действие указанных периодических возмущений, когда $e = 0.01$ и $i = {\pi /}3$.

Изменение параметров орбитального движения (центра масс системы) за время реализации переходного процесса практически не влияет на результаты приведенных расчетов. В данных примерах почти круговых орбит изменение большой полуоси орбиты на отрезке ${\tau } \in \left[ {0,50} \right]$ была равна соответственно –24.16 км и –11.80 км, что подтверждает известный факт снижения эффективности использования ЭДТС на наклонных орбитах (приблизительно в $1{\text{/cos}}i$ раза). Орбиты, близкие к полярным, требуют отдельного рассмотрения.

Дискретность управления (обновление значения управляющего воздействия) в соответствии с (5.13) для приведенных результатов была равна $20\,{\text{с}}$, увеличении времени до 100 c не сказывалось на качестве переходных процессов в системе. Другие особенности реальных систем регулирования, например, случайный характер измерений, запаздывание, не учитывались.

Здесь необходимо отметить, что все приведенные результаты справедливы при достаточно малой величине тока I_n. Во-первых, должно выполняться условие статической устойчивости вертикальных положений равновесия $\sigma < 1$ (4.1), критическая величина тока в приведенном примере для экваториальной орбиты – $0.45\,{\text{A}}$; во-вторых, в соответствии с величиной тока I_n должна быть выбрана величина весового коэффициента ${{s}_{3}}$ (5.2), определяющая уровень управляющего воздействия.

При использовании управления без обратной связи по отклонению $r'$ (5.13), то есть при ${{k}_{r}} = 0$, система теряет устойчивость, как на экваториальных, так и на наклонных орбитах. На рис. 5 приводятся результаты моделирования, показывающие поведение системы при отсутствии регулирования по отклонению $r'$ для случая, когда $e = 0.01$ и $i = {\pi /}3$. Используются те же исходные данные, что и для рис. 4. Здесь имеет место постепенное увеличение деформации троса, что проявляется в уменьшении расстояния между концевыми точками рис. 5,a. Увеличение деформации троса сопровождается увеличением амплитуд колебаний троса по углу ${\varphi }$ вне плоскости орбиты рис. 5,b. В конечном итоге трос сминается, так как не воспринимает сжимающих усилий, и имеет сложную форму рис. 5,с. Результаты, приведенные на рис. 5, получены по математической модели (6.1). Естественно математическая модель (3.1–3.3), в которой используется предположение о форме троса в виде дуги окружности, такие эволюционные движения системы полностью описать не может.

7. Заключение. В работе рассмотрены некоторые особенности стабилизации движения ЭДТС на околоземной орбите. Предложена аналитическая процедура расчета приближенно оптимального регулятора для данной задачи, основанная на совместном применении принципа динамического программирования Беллмана и метода усреднения. Показано, что при наличии малых начальных возмущений с помощью предлагаемого подхода удается обеспечить устойчивость переходных процессов для экваториальных и наклонных орбит. Данный вывод подтверждается результатами математического моделирования по модели с распределенными параметрами, учитывающей растяжимость и массу троса. Предлагаемый приближенно оптимальный регулятор обеспечивает уменьшение амплитуд маятниковых и продольных (изгибных) колебаний ЭДТС. Приводится пример потери устойчивости движения ЭДТС при отсутствии контроля над изгибными колебаниями троса.

Работа поддержана грантом РФФИ РФ 16-41-630637.