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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Кватернион — это четырехмерное гиперкомплексное число вида:

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

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

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

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

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

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

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

  

  Определение программного кватерниона на каждом такте с шагом равным Для этого на каждом такте формируются обновленные значения угловой скорости и угла, угловое ускорение на каждом такте до момента времени t1 равно , от t1 до t2 равно 0, от t2 до t3 равно с помощью нового значения угла и компонент оси поворота на каждом такте рассчитывается новый кватернион, который необходимо перемножать с начальным кватернионом для формирования программного кватерниона:

  

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

Построение графиков изменения компонент программного кватерниона, угла, угловой скорости и углового ускорения (рисунки 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");

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

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

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

Тогда угол поворота:

Орты:

Время конца разгона:

Время начала торможения:

Время окончания торможения:

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

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

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

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

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

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

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

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

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

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

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