MATLAB 利用YALMIP+Gurobi 求解线性规划 -多无人机扫描覆盖

使用要点

1. 创建决策变量
2. 设置目标函数
3. 添加约束条件
4. 参数配置
5. 求解问题

问题描述

假设M个无人机的任务是尽快覆盖一组由 P 顶点表示的多边形凸区域，假设每架无人机的最大飞行时间是有限的，并且是预先知道的。每架无人机的都配备了一个指向下方的机载摄像头，无人机需要使用它来感知 P 指定的整个区域。具体问题如下：

1. 最大限度减少覆盖区域 P 所需时间的无人机数量 m <= M
2. 以最短完成任务时间为目标规划每个 UAV 路径

变量定义

$C_{ij}$ represents the traversing cost of the edge $(i,j)$ between nodes $i$ and $j$.

for  i=1:numberOfVertices
	for j=1:numberOfVertices
		if i~=j
			C(i,j)=norm(V(i,:)-V(j,:))/(uavSpeed*1000/60);
		end
	end	
end

X i j k ∈ { 0 , 1 } {X_{ij}^{k}}\in\{0,1\} Xijk​∈{0,1} represents whether or not the $k$-th UAV is going to fly from vertex $i$ to vertex $j$.

$V_{ij}^k\in \mathbb{R}$ represent the flight speed of the $k$-th UAV while flying from vertex $i$ to vertex $j$.

$t_s\in \mathbb{R}$ be the individual setup time.

$L_k$ be the battery duration of UAV number $k$.

$m\in \mathbb{N}$ be the variable that represents the number of UAVs designed for a mission

$M\in \mathbb{N}$ be the total number of UAVs available.

$O\in \mathbb{N}$ be the constant number of UAV operators.

$N\in \mathbb{N}$ be the number of nodes of the graph.

$d_k$ represents the extra time necessary to launch UAV $k$.

每个无人机的最短时间可表示为：
$$T_k=\sum _{i=1}^N\sum _{j=1}^N\frac{C_{ij}}{V_{ij}^k}{X}_{ij}^k+d_k$$
线性规划问题定义：
$$min(\mathcal{V})$$
$$\sum _{i=1}^N\sum _{j=1}^N\frac{C_{ij}}{V_{ij}^k}{X}_{ij}^k+d_k\le \mathcal{V},k=1,\dots ,M$$
$$t_s\lceil \frac{k}{O}\rceil \sum _{j=1}^N{X}_{1j}^k=d_k,k=1,\dots ,M$$

01 | 创建决策变量

变量设置常用的三种形式：

• sdqvar() 设置实型变量
• intvar() 设置整型变量
• binvar() 设置0-1变量

举例：P = sdqvar(n,m) 表示用 n 行 m 列定义一个矩阵（或标量）P

v = sdpvar(1,1);
X = binvar(numberOfVertices,numberOfVertices,uavNumber,'full');
u = sdpvar(numberOfVertices,1);
m = intvar(1,1);
for k = 1:uavNumber
	for i = 1:numberOfVertices
		X(i,i,k) = 0;
	end
	if uav(k) == 0
		X(:,:,k) = zeros(numberOfVertices);
	end
end

02 | 设置目标函数

YALMIP 默认求解最小值问题，如果问题的目标函数是求解最大值，则需要在目标函数前加负号
目标函数即求最小值 v

03 | 添加约束条件

先初始化约束条件为空，再逐个将约束条件添加

constraints = [];
% 限制1 - 总成本
for k = 1:uavNumber
    sum1 = 0;
    for i = 1:numberOfVertices
        for j = 1:numberOfVertices
            if i ~= j
                sum1 = sum1 + C(i,j)*sum(X(i,j,k));
            end
        end
    end
    constraints = [constraints, v >= sum1 + sum(X(1,:,k))*ceil(k/O)*uavSetupTime];
    % 限制9 - 飞行时间
    constraints = [constraints, sum1 <= uavFlightTime];
end
% 限制2 - 每个顶点必须被访问一次
for j = 2:numberOfVertices
    constraints = [constraints, sum(sum(X(:,j,:))) == 1];
end
% 限制3 - 到达给定节点的节点与离开该节点的节点相同。
for k = 1:uavNumber
    for p = 1:numberOfVertices
        constraints = [constraints, sum(X(:,p,k)) - sum(X(p,:,k)) == 0];
    end
end
% 限制4&5 - 无人机数量
constraints = [constraints, sum(sum(X(1,:,:))) == m];
constraints = [constraints, m <= sum(uav)];
% 限制7 - 循环
for i = 2:numberOfVertices
    for j = 2:numberOfVertices
        if i ~= j
            constraints = [constraints, u(i) - u(j) + numberOfVertices*sum(X(i,j,:)) <= numberOfVertices - 1];
        end
    end
end
% 限制 8 - 强制每条线路仅由一架无人机沿一个方向穿越
for i = 2:2:numberOfVertices
   constraints = [constraints, 1 == sum(X(i,i+1,:)) + sum(X(i+1,i,:))];
end
% 限制 10 - 避免 UA Vs 沿着不平行于覆盖范围的边缘穿过覆盖区域行
for i = 2:2:numberOfVertices
    constraints = [constraints, sum(X(i,i+1,:)) == sum(sum(X(i,3:2:numberOfVertices,:)))];
    constraints = [constraints, sum(X(i+1,i,:)) == sum(sum(X(i+1,2:2:numberOfVertices,:)))];
end

04 | 参数设置

使用 sdpsetting() 函数进行参数设置

% 使用 gurobi 最小化 v 受限于约束变量中包含的约束
options = sdpsettings('verbose',1,'gurobi.Threads',4);

05 | 求解问题

首先使用 solvesdp() 函数求解
再使用 double() 函数将求解的变量转换为实数

solvesdp(constraints,v,options);
X = round(double(X));