1. 固定翼无人机非线性动力学建模基础
固定翼无人机的动力学建模是飞行控制系统的基石,其核心在于准确描述飞行器在三维空间中的运动特性与受力关系。这个建模过程需要融合刚体运动学与空气动力学两大理论体系,通过严谨的数学推导建立非线性微分方程组。
1.1 坐标系定义与转换
在飞行器建模中,我们主要使用两种坐标系:
-
机体坐标系(Body Frame):
- 原点:飞行器质心
- X轴:指向机头方向
- Y轴:指向右机翼
- Z轴:根据右手定则确定(向下)
-
地面坐标系(Earth Frame):
- 通常采用北-东-地(NED)坐标系
- 原点:地面某固定点
- X轴:指向正北
- Y轴:指向正东
- Z轴:指向地心
两个坐标系之间的转换通过欧拉角实现,具体转换矩阵为:
code复制R = [cosθcosψ sinφsinθcosψ-cosφsinψ cosφsinθcosψ+sinφsinψ;
cosθsinψ sinφsinθsinψ+cosφcosψ cosφsinθsinψ-sinφcosψ;
-sinθ sinφcosθ cosφcosθ]
其中φ、θ、ψ分别为滚转角、俯仰角和偏航角。
1.2 刚体运动学方程
基于牛顿-欧拉方程,我们可以建立飞行器的六自由度运动方程:
平移动力学方程:
code复制m(dV/dt + ω×V) = F
其中:
- m为飞行器质量
- V = [u v w]^T为机体坐标系下的速度矢量
- ω = [p q r]^T为角速度矢量
- F为外力矢量
旋转动力学方程:
code复制I(dω/dt) + ω×(Iω) = M
其中:
- I为惯性矩阵
- M为外力矩矢量
1.3 气动力与力矩建模
气动力和力矩是飞行器运动的主要驱动力,可以表示为:
气动力:
code复制F_aero = [X Y Z]^T = qS[-C_D C_Y -C_L]^T
其中:
- q = 0.5ρV^2为动压
- S为参考面积
- C_D, C_Y, C_L分别为阻力、侧力和升力系数
气动力矩:
code复制M_aero = [L M N]^T = qS[bC_l cC_m bC_n]^T
其中:
- b为翼展
- c为平均气动弦长
- C_l, C_m, C_n分别为滚转、俯仰和偏航力矩系数
这些气动系数通常是攻角α、侧滑角β、控制面偏转δ以及无量纲角速度p̂、q̂、r̂的复杂函数,需要通过风洞试验或计算流体力学(CFD)获得。
2. 非线性动力学模型构建
2.1 状态变量与输入变量定义
完整的非线性模型需要明确定义状态变量和输入变量:
状态变量:
code复制x = [u v w p q r φ θ ψ x y z]^T
包含:
- 线速度(u,v,w)
- 角速度(p,q,r)
- 欧拉角(φ,θ,ψ)
- 位置(x,y,z)
控制输入:
code复制u = [δ_e δ_a δ_r δ_t]^T
包含:
- 升降舵偏角δ_e
- 副翼偏角δ_a
- 方向舵偏角δ_r
- 油门开度δ_t
2.2 完整非线性方程推导
结合运动学方程和气动力模型,可以得到完整的非线性动力学方程:
线加速度方程:
code复制u̇ = rv - qw - gsinθ + (X_aero + X_thrust)/m
v̇ = pw - ru + gcosθsinφ + Y_aero/m
ẇ = qu - pv + gcosθcosφ + Z_aero/m
角加速度方程:
code复制ṗ = (I_zL + I_xzN - (I_z(I_z-I_y)+I_xz²)qr + I_xz(I_x-I_y+I_z)pq)/(I_xI_z-I_xz²)
q̇ = (M + (I_z-I_x)pr - I_xz(p²-r²))/I_y
ṙ = (I_xzL + I_xN + (I_x(I_x-I_y)+I_xz²)pq - I_xz(I_x-I_y+I_z)qr)/(I_xI_z-I_xz²)
姿态角速率方程:
code复制φ̇ = p + qsinφtanθ + rcosφtanθ
θ̇ = qcosφ - rsinφ
ψ̇ = (qsinφ + rcosφ)/cosθ
位置方程:
code复制ẋ = ucosθcosψ + v(sinφsinθcosψ-cosφsinψ) + w(cosφsinθcosψ+sinφsinψ)
ẏ = ucosθsinψ + v(sinφsinθsinψ+cosφcosψ) + w(cosφsinθsinψ-sinφcosψ)
ż = -usinθ + vsinφcosθ + wcosφcosθ
2.3 气动系数建模
气动系数的准确建模对非线性模型的精度至关重要。通常采用如下形式:
code复制C_L = C_L0 + C_Lαα + C_Lq(q̂) + C_Lδeδe
C_D = C_D0 + kC_L² + C_Dδeδe
C_Y = C_Yββ + C_Yp(p̂) + C_Yr(r̂) + C_Yδaδa + C_Yδrδr
C_l = C_lββ + C_lp(p̂) + C_lr(r̂) + C_lδaδa + C_lδrδr
C_m = C_m0 + C_mαα + C_mq(q̂) + C_mδeδe
C_n = C_nββ + C_np(p̂) + C_nr(r̂) + C_nδaδa + C_nδrδr
其中:
- α = arctan(w/u)为攻角
- β = arcsin(v/V)为侧滑角
- V = √(u²+v²+w²)为空速
- p̂ = pb/(2V), q̂ = qc/(2V), r̂ = rb/(2V)为无量纲角速度
3. 小扰动线性化方法
3.1 线性化理论基础
小扰动线性化是基于泰勒展开的一种方法,假设飞行器在标称状态附近做小幅度运动。对于非线性系统:
code复制ẋ = f(x,u)
在平衡点(x0,u0)附近进行一阶泰勒展开:
code复制Δẋ ≈ AΔx + BΔu
其中:
- Δx = x - x0
- Δu = u - u0
- A = ∂f/∂x|(x0,u0)为状态矩阵
- B = ∂f/∂u|(x0,u0)为控制矩阵
3.2 巡航状态平衡点分析
对于固定翼无人机,典型的巡航状态平衡条件为:
- 直线水平飞行(φ=0)
- 恒定空速(V=V0)
- 恒定高度(h=h0)
- 零侧滑(β=0)
- 零滚转和偏航角速度(p=r=0)
- 恒定俯仰角速度(q=0)
在这些条件下,可以解出平衡状态下的控制输入u0。
3.3 状态空间模型推导
通过计算雅可比矩阵,可以得到线性化的状态空间模型。对于固定翼无人机,通常可以分解为纵向和横侧向两个子系统。
纵向状态空间模型:
状态变量:Δu, Δw, Δq, Δθ, Δh
输入变量:Δδe, Δδt
横侧向状态空间模型:
状态变量:Δv, Δp, Δr, Δφ, Δψ
输入变量:Δδa, Δδr
完整的6阶状态空间模型可以表示为:
code复制Δẋ = AΔx + BΔu
Δy = CΔx + DΔu
其中Δx = [Δu Δv Δw Δp Δq Δr]^T为状态扰动向量。
4. MATLAB实现与仿真
4.1 非线性模型实现
在MATLAB中,我们可以通过编写函数文件实现非线性模型:
matlab复制function dxdt = nonlinear_model(t,x,u,params)
% 提取参数
m = params.m;
Ix = params.Ix;
Iy = params.Iy;
Iz = params.Iz;
Ixz = params.Ixz;
S = params.S;
b = params.b;
c = params.c;
rho = params.rho;
g = params.g;
% 提取状态变量
u = x(1); v = x(2); w = x(3);
p = x(4); q = x(5); r = x(6);
phi = x(7); theta = x(8); psi = x(9);
% 计算气动角度
V = sqrt(u^2 + v^2 + w^2);
alpha = atan2(w,u);
beta = asin(v/V);
% 计算气动系数
[C_D,C_Y,C_L,C_l,C_m,C_n] = compute_aero_coeffs(alpha,beta,p,q,r,u,params);
% 计算气动力和力矩
q_bar = 0.5*rho*V^2;
X_aero = q_bar*S*(-C_D);
Y_aero = q_bar*S*C_Y;
Z_aero = q_bar*S*(-C_L);
L_aero = q_bar*S*b*C_l;
M_aero = q_bar*S*c*C_m;
N_aero = q_bar*S*b*C_n;
% 计算推力
T = compute_thrust(u(4),V,params);
% 动力学方程
dudt = r*v - q*w - g*sin(theta) + (X_aero + T)/m;
dvdt = p*w - r*u + g*cos(theta)*sin(phi) + Y_aero/m;
dwdt = q*u - p*v + g*cos(theta)*cos(phi) + Z_aero/m;
% 角动力学方程
dpdt = (Iz*L_aero + Ixz*N_aero - (Iz*(Iz-Iy)+Ixz^2)*q*r + Ixz*(Ix-Iy+Iz)*p*q)/(Ix*Iz-Ixz^2);
dqdt = (M_aero + (Iz-Ix)*p*r - Ixz*(p^2-r^2))/Iy;
drdt = (Ixz*L_aero + Ix*N_aero + (Ix*(Ix-Iy)+Ixz^2)*p*q - Ixz*(Ix-Iy+Iz)*q*r)/(Ix*Iz-Ixz^2);
% 运动学方程
dphidt = p + q*sin(phi)*tan(theta) + r*cos(phi)*tan(theta);
dthetadt = q*cos(phi) - r*sin(phi);
dpsidt = (q*sin(phi) + r*cos(phi))/cos(theta);
% 位置方程
dxdt = u*cos(theta)*cos(psi) + v*(sin(phi)*sin(theta)*cos(psi)-cos(phi)*sin(psi)) + w*(cos(phi)*sin(theta)*cos(psi)+sin(phi)*sin(psi));
dydt = u*cos(theta)*sin(psi) + v*(sin(phi)*sin(theta)*sin(psi)+cos(phi)*cos(psi)) + w*(cos(phi)*sin(theta)*sin(psi)-sin(phi)*cos(psi));
dhdt = -u*sin(theta) + v*sin(phi)*cos(theta) + w*cos(phi)*cos(theta);
dxdt = [dudt; dvdt; dwdt; dpdt; dqdt; drdt; dphidt; dthetadt; dpsidt; dxdt; dydt; dhdt];
end
4.2 线性化实现
在MATLAB中,可以使用符号计算工具箱进行线性化:
matlab复制function [A,B] = linearize_model(x0,u0,params)
% 创建符号变量
syms u v w p q r phi theta psi real
syms de da dr dt real
x_sym = [u; v; w; p; q; r; phi; theta; psi];
u_sym = [de; da; dr; dt];
% 计算符号形式的导数
dxdt_sym = nonlinear_model_sym(0,x_sym,u_sym,params);
% 计算雅可比矩阵
A_sym = jacobian(dxdt_sym,x_sym);
B_sym = jacobian(dxdt_sym,u_sym);
% 在平衡点求值
A = double(subs(A_sym,[x_sym;u_sym],[x0;u0]));
B = double(subs(B_sym,[x_sym;u_sym],[x0;u0]));
% 提取6阶状态空间模型
A = A(1:6,1:6);
B = B(1:6,:);
end
4.3 LQR控制器设计
基于线性化模型,可以设计LQR控制器:
matlab复制function [K,S] = design_lqr(A,B,Q,R)
% 解决连续时间代数Riccati方程
[K,S] = lqr(A,B,Q,R);
end
% 使用示例
Q = diag([10 10 10 1 1 1]); % 状态权重矩阵
R = diag([1 1 1 1]); % 控制权重矩阵
[K,~] = design_lqr(A,B,Q,R);
5. 仿真结果与分析
5.1 非线性模型验证
通过对比非线性模型与已知飞行数据的响应,可以验证模型的准确性。典型的验证场景包括:
- 阶跃升降舵响应
- 阶跃副翼响应
- 双脉冲方向舵输入
- 风扰动响应
5.2 线性化模型精度分析
通过比较线性化模型和非线性模型对小扰动的响应,评估线性化模型的精度。通常考察:
- 模型在平衡点附近的预测能力
- 不同扰动幅度下的误差变化
- 主要模态(短周期、长周期等)的频率和阻尼比
5.3 LQR控制性能评估
对设计的LQR控制器进行闭环仿真,评估以下性能指标:
- 阶跃响应的超调量和调节时间
- 抗风扰能力
- 控制能量消耗
- 鲁棒性分析
6. 实际应用中的注意事项
6.1 模型不确定性处理
实际应用中需要考虑:
- 气动参数的不确定性
- 质量特性变化(如燃油消耗)
- 传感器噪声和延迟
- 执行器动态特性
可以采用鲁棒控制或自适应控制方法增强系统鲁棒性。
6.2 实时实现考虑
在嵌入式系统实现时需注意:
- 离散化采样频率选择
- 计算复杂度优化
- 数值稳定性问题
- 故障检测与处理
6.3 飞行测试验证
建议的验证流程:
- 硬件在环仿真(HIL)
- 受限飞行测试(限制扰动幅度)
- 逐步扩大飞行包线
- 全状态飞行测试
在实际飞行测试中,应准备完善的安全措施,包括应急切换机制和飞行终止系统。