2023年第二十届五一杯数学建模竞赛A题保姆级思路代码（gzh数模孵化园）

2023年第二十届五一数学建模竞赛题目

A题：无人机定点投放问题

by:公从号：数模孵化园

拉个个交流裙：715266881，里面会第一时间分享和继续更新各种

问题一：

问题1： 假设无人机以平行于水平面的方式飞行，在空中投放物资（物资为球形，半径20cm，重量50kg）到达地面指定位置。

（1）建立数学模型，给出无人机投放距离（投放物资时无人机与地面物资指定落地点之间的直线距离）与无人机飞行高度、飞行速度、空气阻力等之间的关系。

（2）假设无人机的飞行高度为300m，飞行速度为300km/h，风速为5m/s，风向与水平面平行。建立数学模型，分别给出无人机飞行方向与风向相同（夹角为0度）、相反（夹角为180度）、垂直（夹角为90度）情况下无人机的投放距离。

这里完全就是物理题了，但是这里的物理题是有很大的发挥空间的。

下面给出具体做法：

这个是重点！我们可以将无人机投放物资看作是一个自由落体运动。考虑到无人机的飞行高度、飞行速度和空气阻力，我们可以建立如下的数学模型：

假设：

h 为无人机的飞行高度；

v0 为无人机的水平飞行速度；

m 为物资的质量（50kg）；

R 为物资的半径（0.2m）；

ρ 为空气密度（约为1.225 kg/m³）；

Cd 为物资的阻力系数（假设为0.47，球形物体）；

g 为重力加速度（约为9.81 m/s²）；

t 为物资从投放到落地所用的时间；

D 为投放距离；

vw 为风速；

θ 为无人机飞行方向与风向之间的夹角。

根据自由落体运动，我们可以得到：

t = sqrt(2h/g)

考虑空气阻力，我们可以得到阻力公式：

Fd = 0.5 * ρ * v² * Cd * A

其中 A 为物资的横截面积：A = πR²

物资受到的净水平力为：

Fnet = m * a = m * (v0 + vw * cos(θ)) - Fd

根据牛顿第二定律：

a = Fnet / m

物资的水平位移可以通过以下公式计算：

D = (v0 + vw * cos(θ)) * t

（2）在给定的情况下，我们可以分别计算不同夹角下的投放距离：

无人机飞行高度：h = 300m

无人机飞行速度：v0 = 300km/h = 83.33m/s

风速：vw = 5m/s

情况1：夹角为0度（飞行方向与风向相同）

θ = 0°

t = sqrt(2h/g) = sqrt(2*300/9.81) ≈ 7.8s

D = (v0 + vw * cos(θ)) * t = (83.33 + 5) * 7.8 ≈ 688.2m

情况2：夹角为180度（飞行方向与风向相反）

θ = 180°

t = sqrt(2h/g) = sqrt(2*300/9.81) ≈ 7.8s

D = (v0 + vw * cos(θ)) * t = (83.33 - 5) * 7.8 ≈ 610.6m

情况3：夹角为90度（飞行方向与风向垂直）

θ = 90°

t = sqrt(2h/g) = sqrt(2*300/9.81) ≈ 7.8 s

D = (v0 + vw * cos(θ)) * t = (83.33 + 5 * cos(90°)) * 7.8 ≈ 649.4m

综上所述，在三种不同风向情况下，无人机的投放距离分别为：

飞行方向与风向相同（夹角为0度）：D ≈ 688.2m

飞行方向与风向相反（夹角为180度）：D ≈ 610.6m

飞行方向与风向垂直（夹角为90度）：D ≈ 649.4m

这里随便算的，比较简单，想弄复杂点可以根据实际情况调整无人机的飞行高度、速度等参数。

问题二：

无人机不仅能定点投放物资，而且还可以通过安装在无人机前端的发射筒发射爆炸物疏通河道。其大致过程是：无人机首先水平飞行接近障碍物所处区域，然后俯冲找准时机发射爆炸物，发射结束后无人机随即拉升飞走。现有一处河流被冰块堆积阻断，需要用无人机发射爆炸物（爆炸物为球形，半径8cm，重量5kg）对目标进行爆破。假设无人机的初始点到目标的水平距离为10000m。受环境影响，无人机必须俯冲发射，并且发射方向与无人机的飞行方向一致。

1. 建立数学模型，给出无人机发射距离（发射点与目标的直线距离）与无人机的飞行高度、飞行速度、俯冲角度及发射速度等因素之间的关系。

问题一说，重点是投放物资，可以看成自由落体，这里就不行了。这里是爆破，可以看成斜抛运动。

这一小题可以给出两种思路：

思路一：考虑到无人机的飞行高度、飞行速度、俯冲角度和发射速度，我们可以建立如下的数学模型：

假设：

h0 为无人机的初始飞行高度；

v0 为无人机的水平飞行速度；

m 为爆炸物的质量（5kg）；

R 为爆炸物的半径（0.08m）；

g 为重力加速度（约为9.81 m/s²）；

α 为无人机的俯冲角度；

vf 为发射速度；

D 为发射距离。

根据斜抛运动，我们可以得到爆炸物的水平速度和垂直速度分量：

水平速度分量：vx = v0 * cos(α) + vf * cos(α)

垂直速度分量：vy = v0 * sin(α) + vf * sin(α) - g * t

爆炸物的水平位移为：

Dx = vx * t

爆炸物的垂直位移为：

Dy = vy * t - 0.5 * g * t²

由于发射点与目标的直线距离为发射距离 D，我们有：

D = sqrt(Dx² + Dy²)

将水平位移和垂直位移的公式代入 D 的公式中，我们可以得到：

D = sqrt((v0 * cos(α) + vf * cos(α))² * t² + (v0 * sin(α) + vf * sin(α) - g * t)² * t² - 0.5 * g * t^4)

思路2：

在这个情况下，我们可以将无人机发射爆炸物的过程分为两个阶段：第一阶段是无人机的俯冲过程，第二阶段是爆炸物的飞行过程。

第一阶段：无人机俯冲过程

假设：

h0 为无人机的初始飞行高度；

v0 为无人机的水平飞行速度；

α 为无人机的俯冲角度；

t1 为无人机俯冲所需时间。

在无人机俯冲过程中，我们可以得到：

无人机的水平速度：v0x = v0 * cos(α)

无人机的垂直速度：v0y = v0 * sin(α) - g * t1

在 t1 时间内，无人机的水平位移和垂直位移分别为：

Dx1 = v0x * t1

Dy1 = v0y * t1 - 0.5 * g * t1²

第二阶段：爆炸物飞行过程

假设：

vf 为爆炸物的发射速度；

t2 为爆炸物飞行所需时间。

在爆炸物飞行过程中，我们可以得到：

爆炸物的水平速度：vfx = (v0 * cos(α) + vf * cos(α))

爆炸物的垂直速度：vfy = (v0 * sin(α) + vf * sin(α) - g * t2)

在 t2 时间内，爆炸物的水平位移和垂直位移分别为：

Dx2 = vfx * t2

Dy2 = vfy * t2 - 0.5 * g * t2²

综合两个阶段，我们可以得到总的水平位移和垂直位移：

Dx = Dx1 + Dx2

Dy = Dy1 + Dy2

由于发射点与目标的直线距离为发射距离 D，我们有：

D = sqrt(Dx² + Dy²)

这个模型同样需要知道无人机的飞行高度、飞行速度、俯冲角度和发射速度

（2）假设风速为6m/s，无人机接近目标时的飞行高度为800m、飞行速度为300km/h，爆炸物的发射速度为600km/h（相对于无人机的速度）。要求发射爆炸物时无人机与目标的距离在1000 m-3000 m之间，且无人机的高度不低于300m，请给出无人机击中目标的发射策略。

首先，我们将给定的速度单位统一为米/秒：

无人机飞行速度：v0 = 300km/h = 83.33m/s

爆炸物发射速度：vf = 600km/h = 166.67m/s

风速：vw = 6m/s

思路1：斜抛运动

我们要求无人机的发射距离 D 在1000m到3000m之间，以及无人机的高度不低于300m。为了满足这些条件，我们可以通过调整无人机的俯冲角度α来实现。

在这个问题中，我们需要考虑风速对爆炸物的影响。我们可以将风速分解为与无人机飞行方向平行和垂直的分量：

水平分量：vw_x = vw * cos(α)

垂直分量：vw_y = vw * sin(α)

将风速的分量考虑在内，爆炸物的水平速度和垂直速度分量为：

水平速度分量：vx = v0 * cos(α) + vf * cos(α) + vw_x

垂直速度分量：vy = v0 * sin(α) + vf * sin(α) - g * t + vw_y

此时，我们需要遍历可能的俯冲角度α，找到一个合适的角度使得发射距离 D 在1000m到3000m之间，并且无人机的高度不低于300m。

思路2：分阶段分析

与思路1类似，我们需要考虑风速对无人机和爆炸物的影响。首先，我们需要考虑风速对无人机俯冲过程的影响。将风速分解为与无人机飞行方向平行和垂直的分量：

水平分量：vw_x1 = vw * cos(α)

垂直分量：vw_y1 = vw * sin(α)

在无人机俯冲过程中，考虑风速的影响，我们可以得到：

Dx1 = (v0 * cos(α) + vw_x1) * t1

Dy1 = (v0 * sin(α) + vw_y1 - g * t1) * t1 - 0.5 * g * t1²

接下来，我们需要考虑风速对爆炸物飞行过程的影响。将风速分解为与爆炸物飞行方向平行和垂直的分量：

水平分量：vw_x2 = vw * cos(α)

垂直分量：vw_y2 = vw * sin(α)

在爆炸物飞行过程中，考虑风速的影响，我们可以得到：

Dx2 = (v0 * cos(α) + vf * cos(α) + vw_x2) * t2

Dy2 = (v0 * sin(α) + vf * sin(α) - g * t2 + vw_y2) * t2 - 0.5 * g * t2²

综合两个阶段，我们可以得到总的水平位移和垂直位移：

Dx = Dx1 + Dx2

Dy = Dy1 + Dy2

由于发射点与目标的直线距离为发射距离 D，我们有：

D = sqrt(Dx² + Dy²)

同样，我们需要遍历可能的俯冲角度α，找到一个合适的角度使得发射距离 D 在1000m到3000m之间，并且无人机的高度不低于300m。

划重点：：：：：：：：：

由于涉及到多个变量，包括俯冲角度α，以及俯冲时间t1和爆炸物飞行时间t2，我们需要使用数值方法（例如遍历、优化算法等）来寻找一个合适的策略使得发射距离 D 在1000m到3000m之间，同时满足无人机的高度不低于300m。这里要使用算法了！！

这里我们使用遗传算法。我们的目标是优化无人机的俯冲角度α、俯冲时间t1和爆炸物飞行时间t2。

下面给出具体代码：

import numpy as np

import random

from scipy.optimize import minimize



# 参数

v0 = 83.33

vf = 166.67

vw = 6

g = 9.81

h0 = 800



# 目标函数：发射距离

def target_function(params):

    alpha, t1, t2 = params



    vw_x = vw * np.cos(alpha)

    vw_y = vw * np.sin(alpha)



    Dx1 = (v0 * np.cos(alpha) + vw_x) * t1

    Dy1 = (v0 * np.sin(alpha) + vw_y - g * t1) * t1 - 0.5 * g * t1**2



    Dx2 = (v0 * np.cos(alpha) + vf * np.cos(alpha) + vw_x) * t2

    Dy2 = (v0 * np.sin(alpha) + vf * np.sin(alpha) - g * t2 + vw_y) * t2 - 0.5 * g * t2**2



    D = np.sqrt((Dx1 + Dx2)**2 + (Dy1 + Dy2)**2)

    return D



# 约束条件

def constraint(params):

    alpha, t1, t2 = params

    h = h0 - (v0 * np.sin(alpha) - g * t1) * t1 - 0.5 * g * t1**2

    return h - 300



# 遗传算法

from pyeasyga import pyeasyga



# 初始化遗传算法

ga = pyeasyga.GeneticAlgorithm(seed_data=None,

                               population_size=200,

                               generations=100,

                               crossover_probability=0.8,

                               mutation_probability=0.2,

                               elitism=True,

                               maximise_fitness=False)



# 适应度函数

def fitness_function(params):

    D = target_function(params)

    h = constraint(params)

    return abs(D) if 1000 <= D <= 3000 and h >= 0 else 10000



ga.fitness_function = fitness_function



# 创建初始种群

def create_individual(data):

    alpha = random.uniform(0, np.pi / 2)

    t1 = random.uniform(0, 30)

    t2 = random.uniform(0, 30)

    return (alpha, t1, t2)



ga.create_individual = create_individual



# 交叉

def crossover(parent_1, parent_2):

    idx = random.randrange(1, len(parent_1))

    child_1 = parent_1[:idx] + parent_2[idx:]

    child_2 = parent_2[:idx] + parent_1[idx:]

    return child_1, child_2



ga.crossover_function = crossover



# 突变

def mutate(individual):

    mutate_idx = random.randrange(len(individual))

    if mutate_idx == 0:

        individual[ mutate_idx] = random.uniform(0, np.pi / 2)

else:

individual[mutate_idx] = random.uniform(0, 30)

ga.mutate_function = mutate



选择

def selection(population):

return random.choice(population)



ga.selection_function = selection



运行遗传算法

ga.run()



输出最佳结果

best_individual = ga.best_individual()

best_params = best_individual[1]

best_D = target_function(best_params)

best_h = constraint(best_params)



print("Best Parameters: ", best_params)

print("Best Distance: ", best_D)

print("Best Height: ", best_h)

问题三：问题3：无人机发射爆炸物命中目标的精度与无人机飞行的稳定性有很大关系。相同条件下，无人机发射爆炸物时越稳定，命中目标的精度越高。开始俯冲后，无人机操控员需要不断调整无人机的飞行姿态以修正风向、风速对无人机的影响。

1. 在飞行速度、发射速度一定的情况下，综合考虑各种因素，建立数学模型，量化无人机飞行的稳定性，给出稳定性与命中精度之间的关系，并利用数值仿真等方法对无人机的稳定性进行分析验证。

为了量化无人机飞行的稳定性，我们可以考虑引入一个稳定性指标S，这个指标将综合无人机在飞行过程中所受到的气动力、气动矩以及外部环境因素（如风速、风向）的影响。一个高稳定性的无人机将具有较小的S值。假设无人机的稳定性指标S与命中精度P之间存在线性关系，即：

P = k * S + b

其中，k 和 b 是常数，可以通过实验数据进行拟合。

为了建立稳定性指标S，我们可以考虑以下因素：

气动力：当无人机受到较大的气动力时，其稳定性可能较差。气动力F可以表示为：

F = 0.5 * ρ * v² * A * Cd

其中，ρ 是空气密度，v 是无人机的飞行速度，A 是无人机受到气流影响的表面积，Cd 是阻力系数。

气动矩：气动矩M反映了气动力对无人机姿态的影响。气动矩可以表示为：

M = 0.5 * ρ * v² * A * Cl * l

其中，Cl 是力矩系数，l 是力臂。

外部环境因素：我们可以考虑风速w和风向θ对无人机稳定性的影响。风速和风向可以与无人机的飞行速度和飞行方向结合，得到相对风速和相对风向。

综合以上因素，我们可以定义稳定性指标S为：

S = k1 * F + k2 * M + k3 * w * cos(θ)

其中，k1、k2 和 k3 是权重系数，可以根据实际情况进行调整。

蒙特卡洛模拟稳定性检验步骤：

定义相关参数和变量：

N：模拟次数

ρ：空气密度

A：无人机受到气流影响的表面积

Cd：阻力系数

Cl：力矩系数

l：力臂

k1, k2, k3：权重系数

v：无人机的飞行速度

θ：风向

w：风速

S：稳定性指标

P：命中精度

初始化参数值，例如：

N = 10000

ρ = 1.225 kg/m³

A = 1.0 m²

Cd = 0.1

Cl = 0.1

l = 1.0 m

k1 = 1.0

k2 = 1.0

k3 = 1.0

进行蒙特卡洛模拟：

对于 i = 1 到 N：

随机生成飞行速度v、风向θ和风速w的值

计算气动力F、气动矩M和相对风速w相对风向θ

计算稳定性指标S：S = k1 * F + k2 * M + k3 * w * cos(θ)

计算命中精度P：P = k * S + b

将计算得到的S和P添加到对应的结果列表中

分析结果：

计算S和P的平均值、方差等统计量

绘制S和P的散点图，观察它们之间的关系

如果需要，调整权重系数k1、k2和k3，然后重复步骤3和4，以研究各因素对稳定性的贡献

代码：

import numpy as np

import matplotlib.pyplot as plt



# 初始化参数

N = 10000

rho = 1.225

A = 1.0

Cd = 0.1

Cl = 0.1

l = 1.0

k1 = 1.0

k2 = 1.0

k3 = 1.0

k = -1.0

b = 1.0



# 结果列表

S_list = []

P_list = []



# 蒙特卡洛模拟

for i in range(N):

    v = np.random.uniform(50, 150)  # 飞行速度 (m/s)

    theta = np.random.uniform(0, 2 * np.pi)  # 风向 (rad)

    w = np.random.uniform(0, 20)  # 风速 (m/s)



    F = 0.5 * rho * v**2 * A * Cd

    M = 0.5 * rho * v**2 * A * Cl * l

    S = k1 * F + k2 * M + k3 * w * np.cos(theta)



    P = k * S + b



    S_list.append(S)

    P_list.append(P)



# 分析结果

S_mean = np.mean(S_list)

P_mean = np.mean(P_list)

S_std = np.std(S_list)

P_std = np.std(P_list)



print("S Mean:", S_mean)

print("P Mean:", P_mean)

print("S Std:", S_std)

print("P Std:", P_std)



plt.scatter(S_list, P_list, alpha=0.5)

plt.xlabel("Stability Index (S)")

plt.ylabel("Hit Accuracy (P)")

plt.show()

1. 假设风速为6m/s，无人机的飞行速度范围为300 km/h -400 km/h，爆炸物的发射速度为500km/h（相对于无人机的速度）。无人机在800m高度开始俯冲，初始俯冲角度为45°，发射爆炸物时的飞行高度不低于300m，请给出为尽量保持无人机稳定而采取的飞行姿态最优调整策略。

这里是优化策略，理论上可以从非常多的角度去切入，下面给出几个角度，优化一下就行。

在俯冲过程中，无人机的飞行速度应保持在300-400 km/h范围内。这可以通过在俯冲过程中适当调整发动机推力来实现。

在俯冲过程中，操控员应密切关注风速和风向的变化，以便对无人机的姿态进行实时调整。当风速增大或风向发生改变时，操控员应调整无人机的姿态以减小风对其稳定性的影响。

在发射爆炸物时，操控员应保持无人机的飞行高度不低于300m。这可以通过在接近目标时控制俯冲角度并适时调整发动机推力来实现。

在发射爆炸物前，操控员应尽量将无人机调整到相对稳定的姿态，以减小发射过程中受到的干扰。这包括在接近目标时适时减小俯冲角度，以减小气动力对无人机的影响。

在爆炸物发射后，操控员应迅速调整无人机的姿态，以减小发射过程中对无人机稳定性的影响。这可以通过在发射后立即拉升无人机，并适当调整其姿态以应对风速和风向变化来实现。

在整个过程中，操控员应密切关注无人机的各项指标，包括飞行速度、飞行高度、俯冲角度等，以便及时作出调整。

下面给出我自己的实操：

确定影响无人机稳定性的关键因素。假设影响无人机稳定性的主要因素为：飞行速度v、俯冲角度θ、风速w以及风向ϕ。

我们可以采用基于梯度下降的优化算法来寻找最优调整策略。以下是基于这种方法的数学建模过程、算法和代码：

定义一个评估无人机稳定性的目标函数f(v, θ, w, ϕ)，该函数应反映v、θ、w和ϕ对无人机稳定性的影响。这里我们假设目标函数为：

f(v, θ, w, ϕ) = a * v^2 + b * θ^2 + c * w^2 + d * ϕ^2，其中a、b、c和d为权重系数。

初始化飞行速度v、俯冲角度θ、风速w和风向ϕ的值。

计算目标函数关于各参数的梯度，并根据梯度值更新参数：

dv = 2 * a * v

dθ = 2 * b * θ

dw = 2 * c * w

dϕ = 2 * d * ϕ

v = v - lr * dv

θ = θ - lr * dθ

w = w - lr * dw

ϕ = ϕ - lr * dϕ

其中lr为学习率。

重复步骤3，直到达到预定的迭代次数或目标函数收敛

代码：

import numpy as np



def objective_function(v, theta, w, phi, a, b, c, d):

    return a * v**2 + b * theta**2 + c * w**2 + d * phi**2



def gradient_descent(v_init, theta_init, w_init, phi_init, a, b, c, d, lr, max_iter):

    v = v_init

    theta = theta_init

    w = w_init

    phi = phi_init

   

    for _ in range(max_iter):

        dv = 2 * a * v

        dtheta = 2 * b * theta

        dw = 2 * c * w

        dphi = 2 * d * phi

       

        v -= lr * dv

        theta -= lr * dtheta

        w -= lr * dw

        phi -= lr * dphi

       

        if objective_function(v, theta, w, phi, a, b, c, d) < 1e-6:

            break

   

    return v, theta, w, phi



# 初始化参数

v_init = 300 / 3.6

theta_init = np.radians(45)

w_init = 6

phi_init = 0

a, b, c, d = 1, 1, 1, 1

lr = 0.001

max_iter = 1000



# 执行梯度下降

v_opt, theta_opt, w_opt, phi_opt = gradient_descent(v_init, theta_init, w_init, phi_init, a, b, c, d, lr, max_iter)



# 将速度和角度转换回对应的单位

v_opt_kmh = v_opt * 3.6

theta_opt_deg = np.degrees(theta_opt)



print(f"Optimal flight speed: {v_opt_kmh:.2f} km/h")

print(f"Optimal dive angle: {theta_opt_deg:.2f} degrees")

print(f"Optimal wind speed: {w_opt:.2f} m/s")

print(f"Optimal wind direction: {phi_opt:.2f} radians")

作者Nus计算机在读硕士，热爱数模，欢迎大家关注，我会持续在裙里第一时间更新各类数模比赛思路。问题的结尾我再说两句：我的思路都是一天内完成的，难免会有错误，甚至很多错误，很多小伙伴很贴心地支出了，我这边不一一回复了，答疑和反馈请戳这里：

数模孵化园思路反馈答疑问题汇总