2023 年 数维杯（A题）国际大学生数学建模挑战赛 |数学建模完整代码+建模过程全解全析

当大家面临着复杂的数学建模问题时，你是否曾经感到茫然无措？作为2021年美国大学生数学建模比赛的O奖得主，我为大家提供了一套优秀的解题思路，让你轻松应对各种难题。
让我们来看看数维杯A题！

问题重述

1、俯仰力矩和俯仰角变化：
推导俯仰力矩的表达式。
基于给定参数建立俯仰角变化模型。
计算 5 秒、10 秒和 20 秒时的姿态角。

2、滚转、俯仰和偏航力矩：
建立滚转、俯仰和偏航力矩的表达式。
建立姿态角变化模型。
计算 5 秒、10 秒和 20 秒时的姿态角。

3.机动特性：
设计低速和高速飞行的机动以实现平飞任务。

4、加速机动任务：
设计控制输入以实现前进加速和平飞。
考虑低速和高速飞行特性。

问题 1: 俯仰力矩和俯仰角变化

建模思路：

1. 俯仰力矩表达式：

   • 俯仰力矩主要受到共轴刚性转子、螺旋桨推进器、水平尾翼的影响。

   • 共轴刚性转子产生的气动力矩可以表示为：
     $M_{\text{rotor}}=K_{\text{rotor}}\cdot \rho \cdot A\cdot V_{\text{tip}}$

   • 螺旋桨推进器产生的推力和旋转力矩：
     $T_{\text{propeller}}=C_{\text{propeller}}$

   • 水平尾翼产生的气动力矩：
     $M_{\text{horizontal tail}}=C_{\text{horizontal tail}}\cdot q\cdot S_{\text{horizontal}}\cdot {\bar{y}}_{\text{horizontal}}$

2. 俯仰角变化模型：

   • 利用力矩和力的平衡，可以得到俯仰角变化的微分方程：
     $I_{yy}\cdot \dot{\theta }=M_{\text{rotor}}+T_{\text{propeller}}+M_{\text{horizontal tail}}$

   • 其中$I_{yy}$ 是飞行器绕 y 轴的惯性矩。

3. 数值求解：

   • 使用数值求解器（例如欧拉法）对微分方程进行离散求解，得到不同时刻的俯仰角。

4. 初始条件：

   • 使用提供的初始条件（flight altitude, flight speed, control inputs）进行模拟。

import numpy as np
from scipy.integrate import odeint

# 定义常数和初始条件
K_rotor = 0.1
rho = 1.225
A = 10.0
V_tip = 5.0
C_propeller = 0.05
C_horizontal_tail = 0.02
q = 100.0
S_horizontal = 8.0
y_horizontal = 2.0
I_yy = 100.0
initial_conditions = [0.0]  # 初始俯仰角

# 定义微分方程
def pitch_dynamics(theta, t):
    M_rotor = K_rotor * rho * A * V_tip
    T_propeller = C_propeller
    M_horizontal_tail = C_horizontal_tail * q * S_horizontal * y_horizontal

    # 俯仰角变化微分方程
    dtheta_dt = (M_rotor + T_propeller + M_horizontal_tail) / I_yy

    return dtheta_dt

# 定义时间范围
time_points = np.linspace(0, 20, 1000)
#见完整版代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

问题 2: 滚转、俯仰和偏航力矩

建模思路：

1. 滚转力矩表达式：

   • 滚转力矩主要由共轴刚性转子和水平尾翼贡献。

   • 共轴刚性转子产生的滚转力矩与问题一中的俯仰力矩相似，可以表示为：
     $M_{\text{roll}}=K_{\text{rotor}}\cdot \rho \cdot A\cdot V_{\text{tip}}$

   • 水平尾翼产生的滚转力矩：
     $M_{\text{horizontal tail}}=C_{\text{horizontal tail}}\cdot p\cdot S_{\text{horizontal}}\cdot {\bar{y}}_{\text{horizontal}}$

     • $C_{\text{horizontal tail}}$ 是水平尾翼力矩系数，$p$ 是滚转角速度，$S_{\text{horizontal}}$ 是水平尾翼面积，${\bar{y}}_{\text{horizontal}}$ 是水平尾翼相对飞行器中心的距离。

2. 俯仰力矩表达式：

   • 俯仰力矩在问题一中已经建模过，主要由共轴刚性转子和水平尾翼贡献。

3. 偏航力矩表达式：

   • 偏航力矩主要由共轴刚性转子和垂直尾翼贡献。

   • 共轴刚性转子产生的偏航力矩与滚转和俯仰情形相似：
     $M_{\text{yaw}}=K_{\text{rotor}}\cdot \rho \cdot A\cdot V_{\text{tip}}$

   • 垂直尾翼产生的偏航力矩：
     $M_{\text{vertical tail}}=C_{\text{vertical tail}}\cdot r\cdot S_{\text{vertical}}\cdot {\bar{y}}_{\text{vertical}}$

     • $C_{\text{vertical tail}}$ 是垂直尾翼力矩系数，(r) 是偏航角速度，$S_{\text{vertical}}$ 是垂直尾翼面积，${\bar{y}}_{\text{vertical}}$ 是垂直尾翼相对飞行器中心的距离。

4. 滚转、俯仰和偏航角变化模型：

   • 利用力矩和力的平衡，可以得到滚转、俯仰和偏航角变化的微分方程：
     $I_{\text{roll}}\cdot \dot{p}=M_{\text{roll}}+M_{\text{horizontal tail}}$

     $I_{\text{pitch}}\cdot \dot{q}=M_{\text{rotor}}+T_{\text{propeller}}+M_{\text{horizontal tail}}$

   $I_{\text{yaw}}\cdot \dot{r}=M_{\text{yaw}}+M_{\text{vertical tail}}$

   • 其中 $I_{\text{roll}}$、$I_{\text{pitch}}$ 和 $I_{\text{yaw}}$分别是飞行器绕 x、y 和 z 轴的惯性矩。

5. 数值求解：

   • 使用数值求解器（例如欧拉法）对微分方程进行离散求解，得到不同时刻的滚转、俯仰和偏航角。

6. 初始条件：

   • 使用提供的初始条件（flight altitude, flight speed, control inputs）进行模拟。

import numpy as np
from scipy.integrate import odeint

# 定义常数和初始条件
K_rotor = 0.1
rho = 1.225
A = 10.0
V_tip = 5.0
C_horizontal_tail = 0.02
p = 0.1
S_horizontal = 8.0
y_horizontal = 2.0
C_vertical_tail = 0.03
r = 0.05
S_vertical = 5.0
y_vertical = 1.5
I_roll = 150.0
I_pitch = 100.0
I_yaw = 80.0
initial_conditions = [0.0, 0.0, 0.0]  # 初始滚转、俯仰和偏航角

# 定义微分方程
def dynamics(variables, t):
    p, q, r = variables  # 滚转、俯仰和偏航角速度

    M_roll = K_rotor * rho * A * V_tip + C_horizontal_tail * p * S_horizontal * y_horizontal
    M_pitch = K_rotor * rho * A * V_tip + C_horizontal_tail * q * S_horizontal * y_horizontal
    M_yaw = K_rotor * rho * A * V_tip + C_vertical_tail * r * S_vertical * y_vertical

    # 滚转、俯仰和偏航角速度变化微分方程
    dp_dt = M_roll / I_roll
    dq_dt = M_pitch / I_pitch
    dr_dt = M_yaw / I_yaw

    return [dp_dt, dq_dt, dr_dt]

# 定义时间范围
time_points = np.linspace(0, 20, 1000)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

问题 3: 低速和高速飞行的机动特性设计

建模思路：

1. 低速飞行的机动特性设计：

   • 低速飞行时，主要由共轴刚性转子和螺旋桨推进器提供控制。
   • 设计共轴刚性转子的位置和螺旋桨推进器的工作能力，使飞行器实现水平飞行（零姿态角）。

2. 高速飞行的机动特性设计：

   • 高速飞行时，主要通过螺旋桨推进器、升降舵和方向舵实现控制。
   • 设计螺旋桨推进器、升降舵和方向舵的控制输入，使飞行器在高速时能够保持水平飞行（零姿态角）。

3. 设计方法：

   • 使用数值模拟和优化方法，调整共轴刚性转子、螺旋桨推进器、升降舵和方向舵的参数，以满足低速和高速飞行时的控制需求。
   • 可以采用基于物理模型的仿真工具，对不同飞行条件下的机动特性进行模拟。

4. 优化目标：

   • 低速飞行：使共轴刚性转子和螺旋桨推进器的输出满足水平飞行的需求，即姿态角为零。
   • 高速飞行：通过调整螺旋桨推进器、升降舵和方向舵的输出，使飞行器在高速时能够保持水平飞行。

5. 约束条件：

   • 考虑飞行器的物理限制，如最大速度、最大推力、最大姿态角等。

6. 数学表达式：

   • 基于飞行器的物理模型，建立低速和高速飞行时的动力学方程，并在此基础上进行优化。

7. 数值模拟：

   • 使用数值模拟工具，对设计的参数进行测试，验证低速和高速飞行时的机动特性。

具体步骤：

1. 低速飞行设计：

   • 设计共轴刚性转子的位置和螺旋桨推进器的工作能力。
   • 利用数值模拟验证在低速条件下，飞行器的机动特性是否满足水平飞行的要求。

2. 高速飞行设计：

   • 设计螺旋桨推进器、升降舵和方向舵的控制输入。
   • 利用数值模拟验证在高速条件下，飞行器的机动特性是否满足水平飞行的要求。

3. 参数调整和优化：

   • 根据数值模拟的结果，调整和优化飞行器的参数，以达到设计的低速和高速飞行机动特性。

4. 验证和测试：

   • 对优化后的飞行器进行验证和测试，确保其在实际飞行中能够实现设计的机动特性。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize

# 定义初始条件和参数
initial_conditions_low_speed = [0.0, 0.0, 0.0]  # 滚转、俯仰和偏航角
initial_conditions_high_speed = [0.0, 0.0, 0.0]  # 滚转、俯仰和偏航角
params_low_speed = [0.1, 1.0, 5.0]  # 共轴刚性转子和螺旋桨推进器参数
params_high_speed = [1.0, 2.0, 3.0, 1.0, 2.0, 3.0]  # 螺旋桨推进器、升降舵和方向舵参数

# 定义低速飞行动力学方程
def dynamics_low_speed(variables, t, params):
    # 省略动力学方程，根据具体问题补充
    return [0, 0, 0]

# 定义高速飞行动力学方程
def dynamics_high_speed(variables, t, params):
    # 省略动力学方程，根据具体问题补充
    return [0, 0, 0]

# 定义目标函数
def objective(params):
    # 低速飞行
    result_low_speed = odeint(dynamics_low_speed, initial_conditions_low_speed, time_points, args=(params[:3],))

    # 高速飞行
    result_high_speed = odeint(dynamics_high_speed, initial_conditions_high_speed, time_points, args=(params[3:],))

    # 计算目标函数，例如使得姿态角尽量接近零
    error_low_speed = np.sum(np.abs(result_low_speed[:, :3]))
    error_high_speed = np.sum(np.abs(result_high_speed[:, :3]))

    # 返回总体目标函数
    return error_low_speed + error_high_speed

# 定义约束条件
def constraint(params):
    # 可以根据需要添加约束条件
    return [0]

# 定义时间范围
time_points = np.linspace(0, 20, 1000)

# 优化参数
initial_guess = params_low_speed + params_high_speed
result = minimize(objective, initial_guess, constraints={'type': 'eq', 'fun': constraint})
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

问题四：直升机加速度机动任务

建模思路：

1. 飞行动力学建模：

   • 基于直升机的飞行动力学，建立飞行器的速度、加速度和力矩之间的关系。
   • 考虑在低速和高速飞行模式下，直升机的动力学特性可能不同，因此需要建立不同飞行模式下的动力学模型。

2. 控制输入设计：

   • 设计控制输入，如螺旋桨推进器、升降舵和方向舵的工作状态，以实现加速度机动任务。

3. 目标函数和约束条件：

   • 确定加速度机动任务的目标，例如在给定时间内达到目标速度，并制定相应的目标函数。
   • 添加物理约束条件，如最大速度、最大推力、最大姿态角等。

4. 数学表达式：

   • 建立描述直升机加速度机动任务的数学模型，考虑飞行器的动力学方程、控制输入和约束条件。

5. 数值模拟和优化：

   • 使用数值模拟工具，对设计的参数进行测试，并通过优化算法调整参数以满足加速度机动任务的要求。

具体步骤：

1. 低速和高速动力学模型建立：

   • 建立低速和高速飞行模式下直升机的飞行动力学模型，包括滚转、俯仰、偏航和加速度的方程。

2. 控制输入设计：

   • 设计螺旋桨推进器、升降舵和方向舵的控制输入，以实现加速度机动任务。
   • 在低速和高速模式下，可能需要不同的控制策略和参数。

3. 目标函数和约束条件制定：

   • 定义加速度机动任务的目标函数，可能包括最小时间达到目标速度等。
   • 添加物理约束条件，如最大速度、最大推力、最大姿态角等。

4. 数学模型建立：

   • 将飞行动力学模型、控制输入和约束条件整合为数学模型，以描述直升机在加速度机动任务中的行为。

5. 数值模拟和优化：

   • 使用数值模拟工具，对设计的参数进行测试，并通过优化算法调整参数以满足加速度机动任务的要求。
   • 可以采用常见的优化算法，如遗传算法、粒子群优化等。

6. 验证和测试：

   • 对优化后的直升机进行验证和测试，确保其在实际飞行中能够实现设计的加速度机动任务。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize

# 定义初始条件和参数
initial_conditions = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  # 滚转、俯仰、偏航角和速度、加速度
params_low_speed = [0.1, 1.0, 5.0]  # 低速模式下的控制输入参数
params_high_speed = [1.0, 2.0, 3.0, 1.0, 2.0, 3.0]  # 高速模式下的控制输入参数

# 定义低速和高速飞行动力学方程
def dynamics(variables, t, params):
    # 省略动力学方程，根据具体问题补充
    return [0, 0, 0, 0, 0, 0]

# 定义目标函数
def objective(params):
    # 低速飞行
    result_low_speed = odeint(dynamics, initial_conditions, time_points, args=(params[:3],))

    # 高速飞行
    result_high_speed = odeint(dynamics, initial_conditions, time_points, args=(params[3:],))

    # 计算目标函数，例如使得加速度尽量大
    error_low_speed = -result_low_speed[-1, 4]  # 取负号表示最大化加速度
    error_high_speed = -result_high_speed[-1, 4]

    # 返回总体目标函数
    return error_low_speed + error_high_speed

# 定义约束条件
def constraint(params):
    # 可以根据需要添加约束条件
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

更多完整的代码和思路看这里：
2023 年 数维杯（A题）国际大学生数学建模挑战赛 |数学建模完整代码+建模过程全解全析