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

  • Так как его можно рассматривать как симметричный объект (как например cubesat) и его САУ управления угловым положением как правило полностью твердотельная и построена на индуктивных катушках которые позволяют управлять угловым положением МКА относительно магнитного поля земли.

  • Индуктивный метод управления позволяет управлять МКА не изменяя его массу (по сравнению например с крупными спутниками которые для управления используют реактивное движение и использованием сжатого азота) что позволяет сильно упростить управления движения и как следствие упросить САУ управления космическим аппаратом.

Входные данные:

  • Начальный кватернион ориентации космического аппарата в инерциальной системе координат ( qн)

  • Начальная угловая скорость — любая в пределах 0,3 град/с;

  • Конечный кватернион в инерциальное системе координат — любой;

  • Конечная угловая скорость — 0 град/с.

Максимально допустимые значения:

  • Максимально допустимое угловое ускорение космического аппарата \epsilon_{max} — 0,01 град/с^2;

  • Максимально допустимая угловая скорость космического аппарата (\omega_{max}) — 0,5 град/с;

  • Такт формирования программного кватерниона ‑ 0,1 с.

Рисунок 1 – примерный вид малого космического аппарата про модель углового движения которого говорится в статье.
Рисунок 1 – примерный вид малого космического аппарата про модель углового движения которого говорится в статье.

Назначение системы программного управления:

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

  • Расчет целеуказаний для солнечной батареи и остр��направленной антенны;

  • Режим ориентации малого космического аппарата (МКА) определяется параметрами ориентации (скорость и кватернион), которые могут быть как постоянные, так и переменные, в заданной системе координат (инерциальной, орбитальной) в заданный временной промежуток.

Структура тракта управления ориентацией малого космического аппарата представлена на рисунке 2.

Рисунок 2 – Структура тракта управления ориентацией малого космического аппарата
Рисунок 2 – Структура тракта управления ориентацией малого космического аппарата

Задачи системы программного управления (СПУ):

  • Система программного управления должна обеспечить малые рассогласования по углам и по скорости на входе системы спутниковой ориентации (ССО);

  • Все развороты выполняются как пространственные одноосные, если нет других требований.

Изначально данный алгоритм разрабатывался в программе MATLAB, но в последствии, в целях импортозамещения было решено использовать отечественную среду разработки Engee. Данная среда математического моделирования работает онлайн и использует Julia-подобный язык программирования.

Рисунок 3 – внешний вид среды разработки Engee
Рисунок 3 – внешний вид среды разработки Engee

Алгоритм с помощью которого будем вести расчет параметров углового положения МКА построен с использованием кватернионов и представлен в виде пунктов 1-8

Рисунок 4 – геометрический смысл кватерниона
Рисунок 4 – геометрический смысл кватерниона

Кватернионы — это расширение комплексных чисел, используемое для представления пространственных поворотов при моделировании углового движения подвижных объектов. Они позволяют описывать вращения без математических проблем таких как «складывание рамок карданного подвеса». При примене��ии трёх поворотов поочерёдно возможно, что первая или вторая ось будет направлена в том же направлении, что и одна из предыдущих осей. Это приводит к потере степени свободы. Например, при пространственном повороте вокруг осей ZYX, когда угол поворота вокруг оси Y становится равным ±90 градусов, оси X и Z выравниваются, что приводит к потере одной степени свободы. 

Кватернион — это четырехмерное гиперкомплексное число вида: q = w +xi+yi+zk

Где i,j,k - гиперкомплексные единицы, удовлетворяющие следующему соотношению:

i^2 = j^2=k^2 = i*j*k = -1

Определение начального и конечного кватерниона;

Вычисление кватерниона поворота по следующим соотношениям:

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

Если \varphi_пбольше 180 градусов, то необходимо поменять знак у конечного кватерниона и будет осуществляться обратный разворот на угол равный:

\varphi_п = -(360-\varphi_{п1})

q_{k1} =- q_{k}

Если \varphi_пменьше -180 градусов, то необходимо поменять знак у конечного кватерниона и будет осуществляться прямой разворот на угол равный:

\varphi_п = (360-\varphi_{п1})

q_{k1} =- q_{k}

  • Определение моментов времени конца разгона (t1), начала (t2) и конца торможения (t3), при этом начало разгона осуществляется в момент времени t0=0:

    Определение программного кватерниона на каждом такте с шагом равным t_{шага} = 0.1 с. Для этого на каждом такте формируются обновленные значения угловой скорости и угла, угловое ускорение на каждом такте до момента времени t1 равно \epsilon_{max}, от t1 до t2 равно 0, от t2 до t3 равно \epsilon_{max}. с помощью нового значения угла и компонент оси поворота на каждом такте рассчитывается новый кватер��ион, который необходимо перемножать с начальным кватернионом для формирования программного кватерниона:

На последнем такте значение программного кватерниона присваивается заданному: q_{прог} = q_k

Построение графиков изменения компонент программного кватерниона, угла, угловой скорости и углового ускорения (рисунки 2-4).

Скрипт, реализующий данный алгоритм, был написан в Engee и представлена в листинге 1.

Листинг 1 – Код на языке Engee, реализующий алгоритм программного управления пространственной ориентацией КА.

# подключение библиотеки построения графиков

using Plots

#компоненты начального кватерниона

q1st = 1.0;

q2st = 0.0;

q3st = 0.0;

q4st = 0.0;

#компоненты конечного кватерниона

q1end = 0.0;

q2end = 1.0;

q3end = 0.0;

q4end = 0.0;



#Вычисление кватерниона поворота по формуле деления кватернионов

# представленных в виде отдельных компонентов

# начальный делится на конечный

qp1 = q1st*q1end+q2st*q2end+q3st*q3end+q4st*q4end;

qp2 = q1st*q2end-q2st*q1end-q3st*q4end+q4st*q3end;

qp3 = q1st*q3end-q3st*q1end+q2st*q4end-q4st*q2end;

qp4 = q1st*q4end-q4st*q1end-q2st*q3end+q3st*q2end;

# норма кватерниона

Normp = qp1^2+qp2^2+qp3^2+qp4^2;

#пе��есчет компонент кватернионов

fip1 = 2*acosd(qp1)          

if fip1>180 # если значение угла поворота больше 180 градусов то необходим пересчет

    fip = -(360-fip1)

q1end=-q1end; # смена знака при угле поворота более 180 градусов

q2end=-q2end;

q3end=-q3end;

q4end=-q4end;

elseif fip1<-180 # для поворота на угол менее -180 аналогичная ситуация

    fip = (360+fip1)

q1end=-q1end;  # смена знака при угле поворота более 180 градусов

q2end=-q2end;

q3end=-q3end;

q4end=-q4end;

else

    fip = fip1;

end

# расчет длины ортов

e1 = qp2/sind(fip1/2);

e2 = qp3/sind(fip1/2);

e3 = qp4/sind(fip1/2);

#Вычисление моментов времени

if fip>0

    eps_max = 0.01; # макcисмальная погрешность расчета

    w_max = 0.5;  # максимальная угловая скорость

    h = 0.1; # длина шага по времени

    t1 = w_max/eps_max;# время конца ускорения

    delta_fi = (eps_max*(t1)^2)/2;# погрешность угла поворота

    t2 = t1+((fip-2*delta_fi)/w_max); # время начала торможения

    t3 = t1+t2;# время окончания торможения

elseif fip<0

    eps_max = -0.005; # макисмальная погрешность расчета

    w_max = -0.2;  # максимальная угловая скорость

    h = 0.1;  # длина шага по времени

    t1 = w_max/eps_max; # время конца ускорения

    delta_fi = (eps_max*(t1)^2)/2; # погрешность угла поворота

    t2 = t1+((fip-2*delta_fi)/w_max);  # время начала торможения

    t3 = t1+t2; # время окончания торможения

end



# максимальное значение счетчика

N=floor(Int,round(t3/h)+1);

# инициализация выходных массивов нулями

t = zeros(N); # массив времени

Eps = zeros(N,1); # массив содержащий изменение погрешности

w = zeros(N,1); # массив содержащий угловую скорость

q = zeros(4,N);  # массив содержаний компоненты кватерниона

fi = zeros(N,1);  # массив содержаний углы ориентации космического аппарата

q_prog1 = zeros(4,N);  # компоненты программного кватерниона

for i = 2 : N

    #Увеличение времени на шаг

    t[i] = t[i-1] + h;

    #изменение углового ускорения

    if t[i]<=t1  ## диапазон времени 1

        Eps[i]=eps_max;

    elseif t[i]<=t2  # диапазон времени 2

        Eps[i]=0;

    elseif t[i]<=t3 # диапазон времени 3

        Eps[i]=-eps_max;

    else

        Eps[i]=0;

    end

    #Новые значения угловой скорости, угла и кватерниона

    w[i] = w[i-1] + h*Eps[i]; # накопление массива содержащего угловую скорость

    fi[i] = fi[i-1] + h*w[i]; # накопление массива содержащего угол поворота МКА

    q[1,i] = cosd(fi[i]/2);  # накопление массива содержашего первую компоненту кватерниона

    q[2,i] = e1*sind(fi[i]/2); # накопление массива содержашего вторую компоненту кватерниона

    q[3,i] = e2*sind(fi[i]/2); # накопление массива содержашего третью компоненту кватерниона

    q[4,i] = e3*sind(fi[i]/2); # накопление массива содержашего четвертую компоненту кватерниона

    #Вычисление программного кватерниона по формуле умножения кватернионов (начальный умножается на пересчитанный каждый такт)

    #        в бортовой программе нет никаких массивов, там только i-й элемент

    q_prog1[1,i] = q1st*q[1,i] - q2st*q[2,i] - q3st*q[3,i] - q4st*q[4,i]; # первая компонента программного кватерниона

    q_prog1[2,i] = q1st*q[2,i] + q2st*q[1,i] + q3st*q[4,i] - q4st*q[3,i]; # вторая компонен��а программного кватерниона

    q_prog1[3,i] = q1st*q[3,i] - q2st*q[4,i] + q3st*q[1,i] + q4st*q[2,i]; # третья компонента программного кватерниона

    q_prog1[4,i] = q1st*q[4,i] + q2st*q[3,i] - q3st*q[2,i] + q4st*q[1,i]; # четвертая компонента программного кватернионаи

end

# расчет погрешности программного кватерниона в конечной точке поворота МКА

delta_prog1 = abs(q1end-q_prog1[1,N]);

delta_prog2 = abs(q2end-q_prog1[2,N]);

delta_prog3 = abs(q3end-q_prog1[3,N]);

delta_prog4 = abs(q4end-q_prog1[4,N]);

# вывод на экран значений погрешности программного кватерниона

print("delta_prog1 = ",delta_prog1,"\n");

print("delta_prog2 = ",delta_prog2,"\n");

print("delta_prog3 = ",delta_prog3,"\n");

print("delta_prog4 = ",delta_prog4,"\n");

#приравнивание к конечному кватерниону

q_prog1[1,N] = q1end;

q_prog1[2,N] = q2end;

q_prog1[3,N] = q3end;

q_prog1[4,N] = q4end;

# построение парамептров поворота МКА (угловое ускорение, угловая скорость, угол)

p1 = plot(t,Eps, title="Изменение углового ускорения", xlabel="Время, с", ylabel = "Угловое ускорение, град/с^2", linewidth=3);

p2 = plot(t,w, title="Изменение угловой скорости", xlabel="Время, с", ylabel = "Угловая скорость, град/с", linewidth=3);

p3 = plot(t,fi, title="Изменение угла", xlabel="Время, с", ylabel = "Угол, град", linewidth=3);

# графики компонент программного кватерниона

p4 = plot(t,[q_prog1[1,:]], title="Первая компонента программного кватерниона", xlabel="Время, с", ylabel = "Компонента 1", linewidth=3);

p5 = plot(t,[q_prog1[2,:]], title="Вторая компонента программного кватерниона", xlabel="Время, с", ylabel = "Компонента 2", linewidth=3);

p6 = plot(t,[q_prog1[3,:],q_prog1[4,:]], title="3,4 компоненты программного кватерниона", xlabel="Время, с", ylabel = "Компоненты 3-4", linewidth=3);

# вывод графиков

display(p1)

display(p2)

display(p3)

display(p4)

display(p5)

display(p6)

## время конца участка углового ускорения

print("t1 = ",t1,"\n");

## время начала участка угловго торможения

print("t2 = ",t2,"\n");

## время конца участка угловго торможения

print("t3 = ",t3,"\n");

Результаты работы скрипта для расчета параметров ориентации МКА

Начальный кватернион q_н

Конечный кватернион q_к

Тогда угол поворота: \varphi = 180

Орты: e1 = 1, e2 = 0, e3 = 0

Время конца разгона: t1 = 50 c

Время начала торможения: t2 = 360 c

Время окончания торможения: t3 = 410 c

Рисунок 5 – циклограмма процесса изменения ориентации МКА
Рисунок 5 – циклограмма процесса изменения ориентации МКА

На рисунке 6 представлен график изменения угла со временем при прямом развороте на заданный угол. На данном рисунке можно увидеть что кривя делится на три части в соответствии с циклограммой (рисунок 5):

С момента времени 0-50 секунд происходит угловое ускорение

В момент 50-360 секунд угловое ускорение равно нулю и угол поворота увеличивается линейно, угловая скорость имеет постоянное значение.

В момент 360-410 секунд угловое ускорение имеет отрицательное значение, угловая скорость линейно уменьшается до нуля, угол стремится к заданному значению (с учетом погрешности)

Рисунок 6 – Изменение угла со временем при прямом развороте на заданный угол
Рисунок 6 – Изменение угла со временем при прямом развороте на заданный угол

На рисунке 7 представлен график изменения угловой скорости со временем при прямом развороте на заданный угол.

Рисунок 7 – Изменение угловой скорости со временем при прямом развороте на заданный угол
Рисунок 7 – Изменение угловой скорости со временем при прямом развороте на заданный угол

На рисунке 8 представлен график изменения углового ускорения со временем при прямом развороте на заданный угол. На графике видно что в момент начала поворота ускорение представляет собой положительный импульс (для создания положительной угловой скорости). Далее когда угловая скорость достигает заданного значения ускорение равно нулю. В мом��нт торможения график ускорения представляет собой отрицательный импульс равный по модулю тому импульсу который был при разгоне.

Рисунок 8 – Изменение углового ускорения со временем при прямом развороте на заданный угол
Рисунок 8 – Изменение углового ускорения со временем при прямом развороте на заданный угол

На рисунке 9-10 представлены графики изменения первой и второй компонент программного кватерниона со временем при прямом развороте на заданный угол.

Рисунок 9 – Изменение первой компоненты программного кватерниона со временем при прямом развороте на заданный угол
Рисунок 9 – Изменение первой компоненты программного кватерниона со временем при прямом развороте на заданный угол
Рисунок 10 – Изменение второй компоненты программного кватерниона со временем при прямом развороте на заданный угол
Рисунок 10 – Изменение второй компоненты программного кват��рниона со временем при прямом развороте на заданный угол

Так как МКА поворачивается в одной плоскости то 3 и 4 компоненты кватерниона будут равны нулю, что и видно по графику (рисунок 11)

Рисунок 11 – Изменение второй компоненты программного кватерниона со временем при прямом развороте на заданный угол
Рисунок 11 – Изменение второй компоненты программного кватерниона со временем при прямом развороте на заданный угол

На последнем такте значение программного кватерниона программного приравнивается к значению конечного кватерниона и тогда можно определить погрешность вычисления вычитая по модулю фактическое значение кватерниона и требуемое:

И тогда погрешность составляет:

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