Математическая модель движения БПЛА
В работе рассматривается движение в вертикальной плоскости. Системы координат, которые применяются в дальнейших выводах представлены на рисунке 1.
Рис.1. Системы координат, используемые при разработке математической модели БПЛА
На летательный аппарат действуют сила тяжести, тяга двигателя, аэродинамические силы и силы, необходимые для управления. Эти силы определены как: mg, P, F_a. Силы, которые созданы для управления БПЛА более подробно рассмотрены далее.
Математическая модель БПЛА представляет собой системой линейных однородных дифференциальных уравнений 2 – ого порядка:
, (1)
где
- перемещение ЛА вдоль оси Х (по горизонтали) системы координат ;
– перемещение ЛА вдоль оси Y (по вертикали) системы координат ;
– угол разворота ЛА относительно оси Z (угол тангажа) системы координат ;
- масса БПЛА;
- момент инерции БПЛА относительно оси Z системы координат ;
- сумма сил, которые приложены к БПЛА, действующих вдоль оси X системы координат ;
- сумма сил, которые действуют вдоль оси Y системы координат ;
– сумма моментов, которые приложены к БПЛА относительно оси Z системы координат .
Детально распишем правую часть системы уравнений (1):
, (2)
где
- сила тяги БПЛА;
- сила аэродинамического сопротивления и подъёмная сила;
- сила на рулевой поверхности;
L1 – плечо силы
– сумма моментов, которые приложены к БПЛА относительно оси Z системы координат ;
момент, созданный аэродинамическими силами, которые приложены к фюзеляжу БПЛА;
Проекция аэродинамической силы на ось X выглядит следующим образом:
,
где
аэродинамический коэффициент лобового сопротивления;
плотность воздуха;
скорость полёта;
площадь аэродинамических поверхностей ЛА;
- коэффициент.
Рис.2. Схема действия аэродинамических сил на крыло
Проекция аэродинамической силы на ось y имеет следующий вид:
где
подъёмная сила крыла;
аэродинамическая сила, которая создана элеронами;
аэродинамический коэффициент крыла;
площадь крыльев;
? - угол поворота элерона;
аэродинамический коэффициент элерона;
площадь элеронов;
- коэффициент;
- коэффициент.
Рис.3. Схема действия моментов от аэродинамических сил
Сумма моментов корпуса рассчитывается согласно формуле:
где
момент аэродинамических сил, действующих на 1 часть корпуса ЛА;
момент аэродинамических сил, действующих на 2 часть корпуса ЛА;
угол тангажа;
плечо действия сил ;
– аэродинамические коэффициенты лобового сопротивления 1 и 2 – ой части фюзеляжа;
– площади 1 – ой и 2 – ой части корпуса соответственно;
– аэродинамический коэффициент, стабилизирующий БПЛА по углу .
Примем угол довольно малым ( <30°), тогда . Так же считаем, что в математической модели
После преобразования системы (2) выглядит следующим образом:
(3)
Основные формулы и расчеты описываются в приложении 1.
1.1 Приведение системы дифференциальных уравнений (1) к форме Коши
Преобразуем исходную систему (1). Для этого введем дополнительные переменные:
Тогда система (1) примет вид:
(4)
Система (4) является системой дифференциальных уравнений первого порядка.
1.2 Векторно-матричная форма модели БПЛА
Запишем систему уравнений (4) в векторно-матричном форме:
(5)
Учитывая выполнение условия: , формула (5) приобретет вид:
(6)
Таким образом, движение БПЛА можно рассматривать в виде уравнения в векторно-матричной форме, представленной ниже:
(7)
где,
- вектор состояния динамической системы;
- вектор управления;
Коэффициенты А1 (динамическая матрица системы) и B1 (матрица при управлении):
A1=
B1=
Размерность вектора X1 (6*1), вектора управления U1 (3*1).
Расчетная модель пространственного движения БПЛА
Математическая модель БПЛА представляет собой систему линейных однородных дифференциальных уравнений 2 – ого порядка:
, (8)
где
- перемещение ЛА вдоль оси Z системы координат ;
– угол разворота ЛА относительно оси X (угол крена) системы координат ;
– угол разворота ЛА относительно оси У (угол рыскания) системы координат ;
- масса БПЛА;
- момент инерции БПЛА относительно оси Х системы координат ;
- момент инерции БПЛА относительно оси У системы координат ;
- сумма сил, приложенных к БПЛА, действующих вдоль оси Z системы координат ;
– сумма моментов, приложенных к БПЛА относительно оси X системы координат ;
– сумма моментов, которые приложены к БПЛА относительно оси У системы координат ;
Рис.4. Схема БПЛА, вид сверху
Рис.5. Схема БПЛА, вид по направлению полета
Детально распишем правую часть системы уравнений (8):
, (9)
где
- управляющая сила по оси z;
- управляющая сила по углу ? (крен);
- управляющая сила по углу ? (курс);
L2 – плечо силы
L3 – плечо силы
момент, который возникает из-за действия аэродинамических сил, создаваемыми крыльями ЛА;
момент, который возникает из-за действия аэродинамических сил, создаваемыми корпусом БПЛА;
Рис.6. Схема БПЛА, вид по направлению полета
Сумма моментов? M?_xkр рассчитывается по формуле:
M_xkр= M_(xkр 1)-M_(xkр 2)=1/2 ??V_x?^2 ?(L_кр1 C_(x 1) S_кр1-L_кр2 C_(x 2) S_кр2)??(K_? ) V_x ??K_? ?;
где,
M_хkр1- момент аэродинамических сил, которые действуют на 1 часть крыла БПЛА;
M_xkр2 - момент аэродинамических сил, которые действуют на 2 часть крыла БПЛА;
? - угол крена;
L_кр1,L_кр2 - плечо действия сил F_хkр1 и F_xkр2;
C_x1,C_x2 - аэродинамическая характеристика первой и второй части крыла;
S_кр1,S_кр2 - площадь частей 1 и 2;
K_? – коэффициент.
Полагаем, что угол ? достаточно мал (? <30°) тогда, cos ? ? 1. Так же предполагаем, что в математической модели V_x=V_(x const).
Рис.7. Схема БПЛА, вид сверху
Сумма моментов? M?_yk вычисляется по формуле:
M_yk= M_(yk 1)-M_(yk 2)=1/2 ??V_x?^2 ?(L_корп1 C_(x 1) S_к1-L_корп2 C_(x 2) S_к2)??(K_? ) V_x ??K_? ?;
где,
M_yk1- момент аэродинамических сил, воздействующих на 1 часть корпуса БПЛА
M_yk2 - момент аэродинамических сил, воздействующих на 2 часть корпуса БПЛА;
F_Zкорп1=1/2 ??V_x?^2 C_(x 1) S_к1 - аэродинамическая сила, которая возникает на первой части корпуса;
F_Zкорп2=1/2 ??V_x?^2 C_(x 2) S_к2 - аэродинамическая сила, возникающая на второй части корпуса;
? - угол курса (рыскания);
L_корп1,L_корп2 - плечо действия сил F_zk1 и F_zk2;
C_x1,C_x2 - аэродинамическая характеристика частей 1и 2 корпуса;
S_1,S_2 - площадь частей 1 и 2;
K_? – коэффициент.
Подразумеваем, что угол ? довольно мал (?<30°). Отсюда следует, что cos? ? 1. Так же предполагаем, что в математической модели V_x=V_(x const).
Поучаем преобразованную систему :
, (10)
Основные расчетные формулы и вычисления приведены в приложении 1.
2.1 Приведение системы дифференциальных уравнений (8) к форме Коши
Преобразуем систему (8). Для этого введем дополнительные переменные:
Тогда исходная система (8) примет следующий вид:
(11)
Система (11) представляет собой систему дифференциальных уравнений первого порядка.
2.2 Векторно-матричная форма модели БПЛА
Преобразуем систему уравнений (13) в векторно-матричную форму:
(12)
Тогда, движение БПЛА запишется в виде уравнения в векторно-матричной форме:
(13)
где
- вектор состояния динамической системы;
- вектор управления;
Коэффициенты А2 (динамическая матрица системы) и B2 (матрица при управлении):
A2=
B2=
Размерность вектора X2 (6*1), вектора управления U2 (3*1).
2.3. Уравнение пространственного движения
Таким образом, в работе найдено уравнение движения в вертикальной плоскости (6). Присвоим элементам данного уравнения индекс 1:
(14)
Затем построена математическая модель поперечного движения (12). Присвоим элементам данного уравнения индекс 2.
(15)
Соединим эти два вида движений в одно уравнение. Тогда векторно-матричное уравнение пространственного движения примет вид:
(16)
Матрицы уравнения (16) состоят из матриц уравнений (14) и (15), полученных выше:
X=[¦(X_1@X_2 )], U=[¦(U_1@U_2 )], A=[¦(A_1&0@0&A_2 )], B=[¦(B_1&0@0&B_2 )].
3. Формирование оптимального управления
Итак, векторно-матричное уравнение движения БПЛА представлено в виде:
Неоходимо подобрать такое управление, которое будет обеспечивать , где
Для синтеза оптимального регуляторов линейных стационарных систем в Matlab в блоке ControlSystemToolbox есть функция , вычисляющая оптимальную матрицу коэффициентов регулирования c квадратичным функционалом качества:
J=?_0^(t_k)-{X^T QX+U^T RU}dt,
здесь - квадратичный критерий качества, и - квадратные весовые матрицы, Q - неотрицательно определенная симметрическая матрица размера (n*n), R - положительно определенная симметрическая матрица (q*q).
4. Схема системы управления
Рис.8. Структурная схема САУ БПЛА
На рис.8 в блоке «Формирование Xтреб» рассчитываются параметры требуемой траектории Xтреб. После этого в схеме есть суммирующий элемент, на выходе которого стоит вектор рассогласований (отклонения от требуемой траектории ). По сигналу рассогласования в регуляторе рассчитывается вектор сигнала управления U (U=k•?). Затем в схеме стоит объект управления «ОУ». На сумматор поступает выходной сигнал Х через обратную связь.
5. Требуемая траектория
На рис.9 показан график требуемой траектории в вертикальной плоскости:
Рис.9. Требуемая траектория движения БПЛА в плоскости ОХУаэр
Численные значения точек заданной траектории движения показаны в таблице 1.
При x>=17000 м груз сбрасывается.
Таблица 1
№ , м
, м
0 0 5000
1 3000 5000
2 4000 2500
3 17000 2500
4 21000 2500
5 22000 5000
6 23000 5000
Рис.10. Требуемая траектория движения БПЛА в плоскости ОХZаэр
Таблица 2
№ , м
z, м
7 28000 0
8 33900 2500
9 35850 0
10 59000 0
Рис.11. Траектория движения БПЛА в трехмерном пространстве
Масса БПЛА рассчитывается согласно формуле: ,
где
- начальная масса БПЛА при взлете;
- промежуток времени от начала полёта до настоящего момента;
- расход топлива в единицу времени;
- масса отделяемого груза;
6. Этап отделения груза
На данном этапе, учитывая отделение груза векторно-матричная форма модели БПЛА приобретет вид:
, (17)
где
- вектор внешней силы, который действует на БПЛА во время отделении груза.
Распишем элементы данного вектора по отдельности:
1) - проекция аэродинамической силы на ось X, которая создана открытыми створками грузового люка БПЛА,
где
C_(x створок) – аэродинамический коэффициент открытых створок грузового люка;
?- плотность воздуха;
V_x- скорость полёта;
S_створок- площадь аэродинамических поверхностей открытых створок грузового люка БПЛА;
2) - сила, которая возникает в точках подвеса груза при отделении груза от БПЛА,
где
,
где
- сила, которая приложена к точке подвеса 1
- сила, которая приложена к точке подвеса 2
3) - момент, который возникает в точках подвеса груза по причине действия сил при отделении груза от БПЛА.
- момент, который возникает при отделении груза
Рис.12. Схема сил и моментов, действующих на БПЛА при отделении груза
Рис.13. Циклограмма действия сил и момента при отделении груза
После анализа сил и моментов, которые действуют на БПЛА во время отделении груза, была смоделирована математическая модель в системе Simulink.
7. Моделирование математической модели БПЛА в Simulink
7.1 Программная реализация в системе Мatlab
В системе Matlab были расчитыны процесы управляемого движения БПЛА.
Для этого в системе Matlab был создан m-file. В данном файле были прописаны все численные параметры БПЛА, введены численные значения требуемой траектории. Затем разработана векторно – матричная модель БПЛА, учитывая этап отделения груза. В итоге, проведен синтез оптимального параметрического управления с использованием функции lqr в системе Matlab. В результате получена оптимальная матрица коэффициентов K в регуляторе .
Далее представлена программа, реализованная в Matlab:
Mv=5000;
L1=6.09;
L2=11.875;
L3=6.09;
Kx=1.8;
Ky=7.63;
Ko=-534.24;
m1=2500;
mdv=363;
l1=5.8;
l3=2.7;
Jz=2*m1*l1*l1;
Jx=2*m1*l1*l1+2*mdv*l3*l3;
Jy=2*m1*l1*l1;
V=60;
Mgr=0;
T=10;
dMgr=(Mgr/1.0)*T;
Mgr_sbr=0.0;
dMt=85/(60*60);
%параметры требуемой траектории
xnach=0;
Vpost=V;
h=5000;
Dy=0;
alfa=0;
Dalfa=0;
z=0;
Dz=0;
beta=0;
Dbeta=0;
gamma=0;
Dgamma=0;
Dx=V;
y=h;
y1=2500;
x1=3000;
x2=4000;
x3=17000;
x4=21000;
x5=22000;
x6=23000;
%завершение параметров требуемой траектории
A1=[0,1,0,0,0,0; 0,(-1/Mv)*Kx,0,0,0,0; 0,0,0,1,0,0; 0,0,0,0,0,0; 0,0,0,0,0,1; 0,0,0,0,(1/Jz)*Ko,0];
A2=[0,1,0,0,0,0; 0,0,0,0,0,0; 0,0,0,1,0,0; 0,0,0,0,0,0; 0,0,0,0,0,1; 0,0,0,0,0,0];
B1=[0,0,0; 1/Mv,0,0; 0,0,0; 0,(1/Mv)*Ky*V,0; 0,0,0; 0,0,(-1/Jz)*L1];
B2=[0,0,0; -1/Mv,0,0; 0,0,0; 0,L2/Jx,0; 0,0,0; 0,0,L3/Jy];
A=[A1,zeros(6); zeros(6),A2]
B=[B1,zeros(6,3); zeros(6,3),B2]
C=eye(12);
D=zeros(12,6);
H=ss(A,B,C,D);
Q=1*eye(12);
Q(1,1)=1000000000000;
Q(3,3)=10000000000000;
Q(5,5)=10000000000000;
Q(7,7)=1000000000000;
Q(9,9)=10000000000000;
Q(11,11)=10000000000000;
R=1*eye(6);
N=0;
[K,S,E]=lqr(A,B,Q,R,N)
7.1.1. Анализ системы на управляемость, наблюдаемость и устойчивость
Чтобы система была управляемой, матрица управляемости должна иметь полный ранг.
rank(ctrb(A,B))
ans =
12
Ранг матрицы равен 12, значит, система управляема.
Согласно критерию наблюдаемости, матрица наблюдаемости должна иметь полный ранг.
rank(obsv(A,C))
ans =
12
Ранг матрицы равен 12, отсюда следует, что система наблюдаема.
Для проверки системы на устойчивость нужно найти корни характеристического уравнения. Причем действительные части корней уравнения должны быть отрицательными.
