✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。
🍎完整代码获取 定制创新 论文复现点击：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：做科研，博学之、审问之、慎思之、明辨之、笃行之，是为：博学慎思，明辨笃行。
🔥 内容介绍
一、引言
在现代无人机应用中，多无人机协同作业已成为提升任务效率和拓展应用范围的关键手段。然而，多无人机在共享空间中飞行时，路径规划与防撞问题至关重要。基于多智能体 Q - Learning 强化学习的方法为解决这些问题提供了有效途径。该方法借助 Q - Learning 算法的优势，让每架无人机通过维护独立的 Q 值表进行学习，并通过奖励函数的巧妙设计实现协同防撞，从而在复杂环境中实现高效的协同路径规划。
二、Q - Learning 算法核心原理
（一）Q 值表与状态 - 动作价值函数
Q - Learning 作为值迭代强化学习算法，其核心在于维护 Q 值表，即状态 - 动作价值函数 Q(s,a)。这里的 s 代表智能体所处的状态，a 表示智能体可执行的动作。Q 值表示智能体在某状态 s 下执行某动作 a 能获得的长期累积奖励期望。例如，在无人机的飞行场景中，状态 s 可以是无人机当前的位置、速度、方向等信息的组合，动作 a 可以是向前飞行、向左转向、向右转向等操作。Q 值表的初始值可以随机设定，随着智能体与环境的不断交互，Q 值将逐步更新。
（二）Q 值更新与最优策略收敛
智能体通过不断与环境交互来更新 Q 值。其更新公式为：
三、多智能体场景下的协同实现
（一）独立 Q 表维护
在多智能体场景下，每架无人机作为一个独立的智能体，维护自己的 Q 表。这是因为每架无人机的飞行任务、位置和状态可能不同，需要根据自身情况进行决策。例如，在一个搜索救援任务中，有的无人机负责搜索特定区域，有的无人机负责运输救援物资，它们各自的起始点、目标点和飞行限制不同，因此需要独立的 Q 表来进行学习和决策。
（二）考虑其他无人机位置的决策
每架无人机在做出决策时，不仅要考虑自身的状态，还要考虑其他无人机的位置。这是实现协同防撞的关键。例如，通过传感器或通信设备，无人机可以获取周围其他无人机的位置信息，并将其纳入到自身的状态表示中。这样，当某架无人机在某状态下选择动作时，会综合考虑自身飞行方向、速度以及周围其他无人机的位置，避免与其他无人机发生碰撞。
（三）奖励函数约束实现协同防撞
为了实现多无人机的协同防撞，设计特殊的奖励函数至关重要。奖励函数应包含对碰撞行为的惩罚和对安全协同飞行的奖励。例如：
- 碰撞惩罚
：当某架无人机与其他无人机发生碰撞或过于接近（超过安全距离）时，给予一个较大的负奖励，如 -100。这将促使无人机在学习过程中尽量避免靠近其他无人机，从而实现防撞。
- 安全协同奖励
：当无人机能够在不碰撞的前提下，按照协同任务要求飞行，如保持一定的编队形式或共同完成目标搜索，给予正奖励。例如，当无人机在协同搜索任务中，按照预定的搜索模式飞行且未与其他无人机发生冲突时，给予奖励值 10。通过这样的奖励函数，引导无人机在追求自身最优路径的同时，实现与其他无人机的协同飞行，避免碰撞。
四、多无人机协同路径规划与防撞实现步骤
（一）环境建模
- 空间与障碍物建模
：将多无人机飞行的空间进行建模，划分成网格或其他合适的空间单元。标记出空间中的障碍物位置，如建筑物、山脉等，这些障碍物区域是无人机不能进入的。
（四）路径规划与防撞决策
- 路径规划
：在训练完成后，每架无人机根据收敛的 Q 表，从起始状态开始，每次选择 Q 值最大的动作，逐步规划出一条从起始点到目标点的路径。
- 实时防撞调整
：在实际飞行过程中，无人机实时获取自身状态和其他无人机的位置信息，根据 Q 表和当前状态做出决策。如果检测到与其他无人机可能发生碰撞，通过调整动作，选择能避免碰撞且尽量靠近规划路径的动作，实现实时防撞。
⛳️ 运行结果
📣 部分代码
function [t, x, y, phase, theta_hist] = slip_simulation(m, g, l0, K, ...
theta_neutral, vx_desired, F_thrust, x0, y0, vx0, vy0, t_max, debug)
if nargin < 13
debug = false;
end
% Initialize
x_pos = x0;
y_pos = y0;
vx = vx0;
vy = vy0;
current_phase = 'swing';
xc = NaN; % contact point
theta_current = theta_neutral; % current foot placement angle
K_raibert = 0.2; % Raibert controller gain (controlling horizontal speed)
stance_start_time = NaN; % track when stance began
% Storage
t = [];
x = [];
y = [];
phase = [];
theta_hist = [];
time = 0;
dt = 0.001;
while time < t_max
if strcmp(current_phase, 'swing')
% Swing dynamics: free fall
state = [x_pos; y_pos; vx; vy];
[t_seg, state_seg] = ode45(@swing_dynamics, [time, time+dt], state);
% Extract final state
x_pos = state_seg(end, 1);
y_pos = state_seg(end, 2);
vx = state_seg(end, 3);
vy = state_seg(end, 4);
% Compute desired theta using Raibert controller
theta_desired = theta_neutral + K_raibert * (vx - vx_desired);
% Smoothly adjust theta toward target
theta_rate = 2.0; % rad/s - how fast theta can change
delta_theta = theta_desired - theta_current;
theta_change = theta_rate * dt * sign(delta_theta) * min(1, abs(delta_theta));
theta_current = theta_current + theta_change;
% Debug print theta and vx during swing
if debug && mod(length(t), 50) == 0
fprintf(' [SWING t=%.3f] theta_current=%.2f deg, theta_desired=%.2f deg, vx=%.4f m/s, vx_desired=%.4f m/s\n', ...
time, rad2deg(theta_current), rad2deg(theta_desired), vx, vx_desired);
end
% Check for touchdown
touchdown_height = l0 * cos(theta_current);
if y_pos <= touchdown_height && vy < 0
if debug
fprintf('TOUCHDOWN at t=%.3f: x=%.3f, y=%.3f, vx=%.3f, vy=%.3f, theta=%.1f deg\n', ...
time, x_pos, y_pos, vx, vy, rad2deg(theta_current));
end
% Calculate contact point (foot placement)
xc = x_pos + l0 * sin(theta_current);
if debug
l_check = sqrt((x_pos - xc)^2 + y_pos^2);
fprintf(' Contact at xc=%.3f, spring length at touchdown: l=%.3f (l0=%.3f)\n', ...
xc, l_check, l0);
fprintf(' Expected touchdown height: y = %.3f, actual y = %.3f\n', touchdown_height, y_pos);
end
current_phase = 'stance';
stance_start_time = time; % record when stance began
end
% Store swing data
t = [t; t_seg];
x = [x; state_seg(:, 1)];
y = [y; state_seg(:, 2)];
phase = [phase; ones(length(t_seg), 1)]; % 1 = swing
theta_hist = [theta_hist; theta_current * ones(length(t_seg), 1)];
time = t_seg(end);
else % stance
% Stance dynamics with thrust force
state = [x_pos; y_pos; vx; vy];
[t_seg, state_seg] = ode45(@(t, s) stance_dynamics(t, s, m, g, K, l0, xc, F_thrust, stance_start_time, debug, time), ...
[time, time+dt], state);
% Extract final state
x_pos = state_seg(end, 1);
y_pos = state_seg(end, 2);
vx = state_seg(end, 3);
vy = state_seg(end, 4);
% Detailed stance debugging
if debug
l = sqrt((x_pos - xc)^2 + y_pos^2);
spring_compression = l0 - l;
% Force components
F_spring_x = K * (l0 - l) * (x_pos - xc) / l;
F_spring_y = K * (l0 - l) * y_pos / l;
F_gravity_y = -m * g;
% Accelerations
xddot = F_spring_x / m;
yddot = F_spring_y / m + F_gravity_y / m;
% Print every 10 steps to avoid spam
if mod(length(t), 10) == 0
fprintf(' [STANCE t=%.3f] theta=%.2f deg, vx=%.4f m/s, x=%.3f, y=%.3f, vy=%.3f\n', ...
time, rad2deg(theta_current), vx, x_pos, y_pos, vy);
fprintf(' l=%.3f, compression=%.3f, xc=%.3f\n', l, spring_compression, xc);
fprintf(' F_spring_x=%.2f, F_spring_y=%.2f, F_grav=%.2f\n', ...
F_spring_x, F_spring_y, F_gravity_y);
fprintf(' xddot=%.2f, yddot=%.2f\n', xddot, yddot);
end
end
% Check for ground penetration during stance
l = sqrt((x_pos - xc)^2 + y_pos^2);
if y_pos < 0
if debug
fprintf(' WARNING: Ground penetration at t=%.3f: y=%.3f, l=%.3f\n', ...
time, y_pos, l);
end
end
% Check for liftoff
l = sqrt((x_pos - xc)^2 + y_pos^2);
if l >= l0 && vy > 0
if debug
fprintf('LIFTOFF at t=%.3f: x=%.3f, y=%.3f, vx=%.3f, vy=%.3f, l=%.3f\n', ...
time, x_pos, y_pos, vx, vy, l);
end
current_phase = 'swing';
xc = NaN;
stance_start_time = NaN;
end
% Store stance data
t = [t; t_seg];
x = [x; state_seg(:, 1)];
y = [y; state_seg(:, 2)];
phase = [phase; 2*ones(length(t_seg), 1)]; % 2 = stance
theta_hist = [theta_hist; theta_current * ones(length(t_seg), 1)];
% Debug: Check for min y during this stance segment
if debug && min(state_seg(:, 2)) < -0.01
fprintf(' Ground penetration detected in stance segment! min_y=%.3f\n', ...
min(state_seg(:, 2)));
end
time = t_seg(end);
end
end
% Nested functions for dynamics
function dstate = swing_dynamics(~, s)
dstate = zeros(4, 1);
dstate(1) = s(3); % xdot
dstate(2) = s(4); % ydot
dstate(3) = 0; % xddot = 0
dstate(4) = -g; % yddot = -g
end
function dstate = stance_dynamics(t_current, s, m, g, K, l0, xc, F_thrust, stance_start_time, ~, ~)
dstate = zeros(4, 1);
l = sqrt((s(1) - xc)^2 + s(2)^2);
% Determine if we're in second half of stance
% Estimate stance duration based on spring oscillation period
t_in_stance = t_current - stance_start_time;
stance_half_period = pi * sqrt(m / K); % half period of spring oscillation
% Apply thrust only in second half of stance
if t_in_stance > stance_half_period
thrust_active = F_thrust;
else
thrust_active = 0;
end
% Check if mass is at ground level
epsilon = 0.001; % ground contact threshold
if s(2) <= epsilon
% Mass in contact with ground - constrain y = 0, compute ground force
dstate(1) = s(3); % xdot
dstate(2) = 0; % ydot = 0
dstate(3) = K/m * (l0 - l) * (s(1) - xc) / l;
dstate(4) = 0;
else
% Mass above ground - free 2D motion under spring + gravity
dstate(1) = s(3); % xdot
dstate(2) = s(4); % ydot
dstate(3) = K/m * (l0 - l) * (s(1) - xc) / l;
dstate(4) = K/m * (l0 - l) * s(2) / l - g + thrust_active / m;
end
end
end
🔗 参考文献
[1]胡细兵.基于强化学习算法的最优潮流研究[D].华南理工大学,2011.
